FFmpeg
rational.c
Go to the documentation of this file.
1 /*
2  * rational numbers
3  * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at>
4  *
5  * This file is part of FFmpeg.
6  *
7  * FFmpeg is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 2.1 of the License, or (at your option) any later version.
11  *
12  * FFmpeg is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with FFmpeg; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20  */
21 
22 /**
23  * @file
24  * rational numbers
25  * @author Michael Niedermayer <michaelni@gmx.at>
26  */
27 
28 #include "avassert.h"
29 #include <limits.h>
30 
31 #include "common.h"
32 #include "mathematics.h"
33 #include "rational.h"
34 
35 int av_reduce(int *dst_num, int *dst_den,
36  int64_t num, int64_t den, int64_t max)
37 {
38  AVRational a0 = { 0, 1 }, a1 = { 1, 0 };
39  int sign = (num < 0) ^ (den < 0);
40  int64_t gcd = av_gcd(FFABS(num), FFABS(den));
41 
42  if (gcd) {
43  num = FFABS(num) / gcd;
44  den = FFABS(den) / gcd;
45  }
46  if (num <= max && den <= max) {
47  a1 = (AVRational) { num, den };
48  den = 0;
49  }
50 
51  while (den) {
52  uint64_t x = num / den;
53  int64_t next_den = num - den * x;
54  int64_t a2n = x * a1.num + a0.num;
55  int64_t a2d = x * a1.den + a0.den;
56 
57  if (a2n > max || a2d > max) {
58  if (a1.num) x = (max - a0.num) / a1.num;
59  if (a1.den) x = FFMIN(x, (max - a0.den) / a1.den);
60 
61  if (den * (2 * x * a1.den + a0.den) > num * a1.den)
62  a1 = (AVRational) { x * a1.num + a0.num, x * a1.den + a0.den };
63  break;
64  }
65 
66  a0 = a1;
67  a1 = (AVRational) { a2n, a2d };
68  num = den;
69  den = next_den;
70  }
71  av_assert2(av_gcd(a1.num, a1.den) <= 1U);
72  av_assert2(a1.num <= max && a1.den <= max);
73 
74  *dst_num = sign ? -a1.num : a1.num;
75  *dst_den = a1.den;
76 
77  return den == 0;
78 }
79 
81 {
82  av_reduce(&b.num, &b.den,
83  b.num * (int64_t) c.num,
84  b.den * (int64_t) c.den, INT_MAX);
85  return b;
86 }
87 
89 {
90  return av_mul_q(b, (AVRational) { c.den, c.num });
91 }
92 
94  av_reduce(&b.num, &b.den,
95  b.num * (int64_t) c.den +
96  c.num * (int64_t) b.den,
97  b.den * (int64_t) c.den, INT_MAX);
98  return b;
99 }
100 
102 {
103  return av_add_q(b, (AVRational) { -c.num, c.den });
104 }
105 
106 AVRational av_d2q(double d, int max)
107 {
108  AVRational a;
109  int exponent;
110  int64_t den;
111  if (isnan(d))
112  return (AVRational) { 0,0 };
113  if (fabs(d) > INT_MAX + 3LL)
114  return (AVRational) { d < 0 ? -1 : 1, 0 };
115  frexp(d, &exponent);
116  exponent = FFMAX(exponent-1, 0);
117  den = 1LL << (61 - exponent);
118  // (int64_t)rint() and llrint() do not work with gcc on ia64 and sparc64,
119  // see Ticket2713 for affected gcc/glibc versions
120  av_reduce(&a.num, &a.den, floor(d * den + 0.5), den, max);
121  if ((!a.num || !a.den) && d && max>0 && max<INT_MAX)
122  av_reduce(&a.num, &a.den, floor(d * den + 0.5), den, INT_MAX);
123 
124  return a;
125 }
126 
128 {
129  /* n/d is q, a/b is the median between q1 and q2 */
130  int64_t a = q1.num * (int64_t)q2.den + q2.num * (int64_t)q1.den;
131  int64_t b = 2 * (int64_t)q1.den * q2.den;
132 
133  /* rnd_up(a*d/b) > n => a*d/b > n */
134  int64_t x_up = av_rescale_rnd(a, q.den, b, AV_ROUND_UP);
135 
136  /* rnd_down(a*d/b) < n => a*d/b < n */
137  int64_t x_down = av_rescale_rnd(a, q.den, b, AV_ROUND_DOWN);
138 
139  return ((x_up > q.num) - (x_down < q.num)) * av_cmp_q(q2, q1);
140 }
141 
143 {
144  int i, nearest_q_idx = 0;
145  for (i = 0; q_list[i].den; i++)
146  if (av_nearer_q(q, q_list[i], q_list[nearest_q_idx]) > 0)
147  nearest_q_idx = i;
148 
149  return nearest_q_idx;
150 }
151 
153  int64_t n;
154  int shift;
155  int sign = 0;
156 
157  if (q.den < 0) {
158  q.den *= -1;
159  q.num *= -1;
160  }
161  if (q.num < 0) {
162  q.num *= -1;
163  sign = 1;
164  }
165 
166  if (!q.num && !q.den) return 0xFFC00000;
167  if (!q.num) return 0;
168  if (!q.den) return 0x7F800000 | (q.num & 0x80000000);
169 
170  shift = 23 + av_log2(q.den) - av_log2(q.num);
171  if (shift >= 0) n = av_rescale(q.num, 1LL<<shift, q.den);
172  else n = av_rescale(q.num, 1, ((int64_t)q.den) << -shift);
173 
174  shift -= n >= (1<<24);
175  shift += n < (1<<23);
176 
177  if (shift >= 0) n = av_rescale(q.num, 1LL<<shift, q.den);
178  else n = av_rescale(q.num, 1, ((int64_t)q.den) << -shift);
179 
180  av_assert1(n < (1<<24));
181  av_assert1(n >= (1<<23));
182 
183  return sign<<31 | (150-shift)<<23 | (n - (1<<23));
184 }
static int shift(int a, int b)
Definition: sonic.c:82
AVRational av_div_q(AVRational b, AVRational c)
Divide one rational by another.
Definition: rational.c:88
#define a0
Definition: regdef.h:46
int num
Numerator.
Definition: rational.h:59
int av_log2(unsigned v)
Definition: intmath.c:26
The reader does not expect b to be semantically here and if the code is changed by maybe adding a a division or other the signedness will almost certainly be mistaken To avoid this confusion a new type was SUINT is the C unsigned type but it holds a signed int to use the same example SUINT a
Definition: undefined.txt:36
#define a1
Definition: regdef.h:47
static const uint8_t q1[256]
Definition: twofish.c:96
uint32_t av_q2intfloat(AVRational q)
Convert an AVRational to a IEEE 32-bit float expressed in fixed-point format.
Definition: rational.c:152
AVRational av_sub_q(AVRational b, AVRational c)
Subtract one rational from another.
Definition: rational.c:101
Round toward +infinity.
Definition: mathematics.h:83
#define av_assert2(cond)
assert() equivalent, that does lie in speed critical code.
Definition: avassert.h:64
Undefined Behavior In the C some operations are like signed integer dereferencing freed accessing outside allocated Undefined Behavior must not occur in a C it is not safe even if the output of undefined operations is unused The unsafety may seem nit picking but Optimizing compilers have in fact optimized code on the assumption that no undefined Behavior occurs Optimizing code based on wrong assumptions can and has in some cases lead to effects beyond the output of computations The signed integer overflow problem in speed critical code Code which is highly optimized and works with signed integers sometimes has the problem that often the output of the computation does not c
Definition: undefined.txt:32
#define max(a, b)
Definition: cuda_runtime.h:33
int av_reduce(int *dst_num, int *dst_den, int64_t num, int64_t den, int64_t max)
Reduce a fraction.
Definition: rational.c:35
#define U(x)
Definition: vp56_arith.h:37
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:259
simple assert() macros that are a bit more flexible than ISO C assert().
int64_t av_gcd(int64_t a, int64_t b)
Compute the greatest common divisor of two integer operands.
Definition: mathematics.c:37
#define FFMAX(a, b)
Definition: common.h:94
#define b
Definition: input.c:41
#define av_assert1(cond)
assert() equivalent, that does not lie in speed critical code.
Definition: avassert.h:53
int64_t av_rescale(int64_t a, int64_t b, int64_t c)
Rescale a 64-bit integer with rounding to nearest.
Definition: mathematics.c:129
#define FFMIN(a, b)
Definition: common.h:96
#define FFABS(a)
Absolute value, Note, INT_MIN / INT64_MIN result in undefined behavior as they are not representable ...
Definition: common.h:72
int n
Definition: avisynth_c.h:760
int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
Rescale a 64-bit integer with specified rounding.
Definition: mathematics.c:58
if(ret)
int av_find_nearest_q_idx(AVRational q, const AVRational *q_list)
Find the value in a list of rationals nearest a given reference rational.
Definition: rational.c:142
Rational number (pair of numerator and denominator).
Definition: rational.h:58
#define isnan(x)
Definition: libm.h:340
Round toward -infinity.
Definition: mathematics.h:82
AVRational av_d2q(double d, int max)
Convert a double precision floating point number to a rational.
Definition: rational.c:106
static int av_cmp_q(AVRational a, AVRational b)
Compare two rationals.
Definition: rational.h:89
common internal and external API header
Utilties for rational number calculation.
AVRational av_add_q(AVRational b, AVRational c)
Add two rationals.
Definition: rational.c:93
int av_nearer_q(AVRational q, AVRational q1, AVRational q2)
Find which of the two rationals is closer to another rational.
Definition: rational.c:127
int den
Denominator.
Definition: rational.h:60
AVRational av_mul_q(AVRational b, AVRational c)
Multiply two rationals.
Definition: rational.c:80