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 }
185 
187 {
188  int64_t gcd, lcm;
189 
190  gcd = av_gcd(a.den, b.den);
191  lcm = (a.den / gcd) * b.den;
192  return lcm < max_den ? av_make_q(av_gcd(a.num, b.num), lcm) : def;
193 }
q1
static const uint8_t q1[256]
Definition: twofish.c:96
av_div_q
AVRational av_div_q(AVRational b, AVRational c)
Divide one rational by another.
Definition: rational.c:88
rational.h
b
#define b
Definition: input.c:41
max
#define max(a, b)
Definition: cuda_runtime.h:33
mathematics.h
av_sub_q
AVRational av_sub_q(AVRational b, AVRational c)
Subtract one rational from another.
Definition: rational.c:101
av_gcd
int64_t av_gcd(int64_t a, int64_t b)
Compute the greatest common divisor of two integer operands.
Definition: mathematics.c:37
av_q2intfloat
uint32_t av_q2intfloat(AVRational q)
Convert an AVRational to a IEEE 32-bit float expressed in fixed-point format.
Definition: rational.c:152
U
#define U(x)
Definition: vp56_arith.h:37
AV_ROUND_UP
@ AV_ROUND_UP
Round toward +infinity.
Definition: mathematics.h:83
av_reduce
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
AVRational::num
int num
Numerator.
Definition: rational.h:59
a1
#define a1
Definition: regdef.h:47
avassert.h
limits.h
FFABS
#define FFABS(a)
Absolute value, Note, INT_MIN / INT64_MIN result in undefined behavior as they are not representable ...
Definition: common.h:72
if
if(ret)
Definition: filter_design.txt:179
AVRational
Rational number (pair of numerator and denominator).
Definition: rational.h:58
isnan
#define isnan(x)
Definition: libm.h:340
c
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
av_rescale_rnd
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
FFMAX
#define FFMAX(a, b)
Definition: common.h:94
av_make_q
static AVRational av_make_q(int num, int den)
Create an AVRational.
Definition: rational.h:71
a
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:41
FFMIN
#define FFMIN(a, b)
Definition: common.h:96
a0
#define a0
Definition: regdef.h:46
av_find_nearest_q_idx
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
av_assert2
#define av_assert2(cond)
assert() equivalent, that does lie in speed critical code.
Definition: avassert.h:64
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:269
av_gcd_q
AVRational av_gcd_q(AVRational a, AVRational b, int max_den, AVRational def)
Return the best rational so that a and b are multiple of it.
Definition: rational.c:186
common.h
av_assert1
#define av_assert1(cond)
assert() equivalent, that does not lie in speed critical code.
Definition: avassert.h:53
AV_ROUND_DOWN
@ AV_ROUND_DOWN
Round toward -infinity.
Definition: mathematics.h:82
av_d2q
AVRational av_d2q(double d, int max)
Convert a double precision floating point number to a rational.
Definition: rational.c:106
av_rescale
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
av_cmp_q
static int av_cmp_q(AVRational a, AVRational b)
Compare two rationals.
Definition: rational.h:89
AVRational::den
int den
Denominator.
Definition: rational.h:60
av_mul_q
AVRational av_mul_q(AVRational b, AVRational c)
Multiply two rationals.
Definition: rational.c:80
shift
static int shift(int a, int b)
Definition: sonic.c:82
av_nearer_q
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
av_add_q
AVRational av_add_q(AVRational b, AVRational c)
Add two rationals.
Definition: rational.c:93
av_log2
int av_log2(unsigned v)
Definition: intmath.c:26