blob: 35ee08877f3497379d383b06a7c424a1c940d1cd
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 | |
80 | AVRational av_mul_q(AVRational b, AVRational c) |
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 | |
88 | AVRational av_div_q(AVRational b, AVRational c) |
89 | { |
90 | return av_mul_q(b, (AVRational) { c.den, c.num }); |
91 | } |
92 | |
93 | AVRational av_add_q(AVRational b, AVRational c) { |
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 | |
101 | AVRational av_sub_q(AVRational b, AVRational c) |
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 | |
127 | int av_nearer_q(AVRational q, AVRational q1, AVRational q2) |
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 | |
142 | int av_find_nearest_q_idx(AVRational q, const AVRational* q_list) |
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 | |
152 | uint32_t av_q2intfloat(AVRational q) { |
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 |