3 * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at>
5 * This file is part of FFmpeg.
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.
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.
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
25 * @author Michael Niedermayer <michaelni@gmx.at>
32 #include "mathematics.h"
35 int av_reduce(int *dst_num
, int *dst_den
,
36 int64_t num
, int64_t den
, int64_t max
)
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
));
43 num
= FFABS(num
) / gcd
;
44 den
= FFABS(den
) / gcd
;
46 if (num
<= max
&& den
<= max
) {
47 a1
= (AVRational
) { num
, 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
;
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
);
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
};
67 a1
= (AVRational
) { a2n
, a2d
};
71 av_assert2(av_gcd(a1
.num
, a1
.den
) <= 1U);
73 *dst_num
= sign
? -a1
.num
: a1
.num
;
79 AVRational
av_mul_q(AVRational b
, AVRational c
)
81 av_reduce(&b
.num
, &b
.den
,
82 b
.num
* (int64_t) c
.num
,
83 b
.den
* (int64_t) c
.den
, INT_MAX
);
87 AVRational
av_div_q(AVRational b
, AVRational c
)
89 return av_mul_q(b
, (AVRational
) { c
.den
, c
.num
});
92 AVRational
av_add_q(AVRational b
, AVRational c
) {
93 av_reduce(&b
.num
, &b
.den
,
94 b
.num
* (int64_t) c
.den
+
95 c
.num
* (int64_t) b
.den
,
96 b
.den
* (int64_t) c
.den
, INT_MAX
);
100 AVRational
av_sub_q(AVRational b
, AVRational c
)
102 return av_add_q(b
, (AVRational
) { -c
.num
, c
.den
});
105 AVRational
av_d2q(double d
, int max
)
108 #define LOG2 0.69314718055994530941723212145817656807550013436025
112 return (AVRational
) { 0,0 };
113 if (fabs(d
) > INT_MAX
+ 3LL)
114 return (AVRational
) { d
< 0 ? -1 : 1, 0 };
115 exponent
= FFMAX( (int)(log(fabs(d
) + 1e-20)/LOG2
), 0);
116 den
= 1LL << (61 - exponent
);
117 // (int64_t)rint() and llrint() do not work with gcc on ia64 and sparc64
118 av_reduce(&a
.num
, &a
.den
, floor(d
* den
+ 0.5), den
, max
);
119 if ((!a
.num
|| !a
.den
) && d
&& max
>0 && max
<INT_MAX
)
120 av_reduce(&a
.num
, &a
.den
, floor(d
* den
+ 0.5), den
, INT_MAX
);
125 int av_nearer_q(AVRational q
, AVRational q1
, AVRational q2
)
127 /* n/d is q, a/b is the median between q1 and q2 */
128 int64_t a
= q1
.num
* (int64_t)q2
.den
+ q2
.num
* (int64_t)q1
.den
;
129 int64_t b
= 2 * (int64_t)q1
.den
* q2
.den
;
131 /* rnd_up(a*d/b) > n => a*d/b > n */
132 int64_t x_up
= av_rescale_rnd(a
, q
.den
, b
, AV_ROUND_UP
);
134 /* rnd_down(a*d/b) < n => a*d/b < n */
135 int64_t x_down
= av_rescale_rnd(a
, q
.den
, b
, AV_ROUND_DOWN
);
137 return ((x_up
> q
.num
) - (x_down
< q
.num
)) * av_cmp_q(q2
, q1
);
140 int av_find_nearest_q_idx(AVRational q
, const AVRational
* q_list
)
142 int i
, nearest_q_idx
= 0;
143 for (i
= 0; q_list
[i
].den
; i
++)
144 if (av_nearer_q(q
, q_list
[i
], q_list
[nearest_q_idx
]) > 0)
147 return nearest_q_idx
;
154 for (a
.num
= -2; a
.num
<= 2; a
.num
++) {
155 for (a
.den
= -2; a
.den
<= 2; a
.den
++) {
156 for (b
.num
= -2; b
.num
<= 2; b
.num
++) {
157 for (b
.den
= -2; b
.den
<= 2; b
.den
++) {
158 int c
= av_cmp_q(a
,b
);
159 double d
= av_q2d(a
) == av_q2d(b
) ?
160 0 : (av_q2d(a
) - av_q2d(b
));
162 else if (d
< 0) d
= -1;
163 else if (d
!= d
) d
= INT_MIN
;
165 av_log(NULL
, AV_LOG_ERROR
, "%d/%d %d/%d, %d %f\n", a
.num
,
166 a
.den
, b
.num
, b
.den
, c
,d
);
167 r
= av_sub_q(av_add_q(b
,a
), b
);
168 if(b
.den
&& (r
.num
*a
.den
!= a
.num
*r
.den
|| !r
.num
!= !a
.num
|| !r
.den
!= !a
.den
))
169 av_log(NULL
, AV_LOG_ERROR
, "%d/%d ", r
.num
, r
.den
);
175 for (a
.num
= 1; a
.num
<= 10; a
.num
++) {
176 for (a
.den
= 1; a
.den
<= 10; a
.den
++) {
177 if (av_gcd(a
.num
, a
.den
) > 1)
179 for (b
.num
= 1; b
.num
<= 10; b
.num
++) {
180 for (b
.den
= 1; b
.den
<= 10; b
.den
++) {
182 if (av_gcd(b
.num
, b
.den
) > 1)
184 if (av_cmp_q(b
, a
) < 0)
186 for (start
= 0; start
< 10 ; start
++) {
190 for (i
= 0; i
<100; i
++) {
191 int exact
= start
+ av_rescale_q(i
+1, b
, a
);
192 acc
= av_add_stable(a
, acc
, b
, 1);
193 if (FFABS(acc
- exact
) > 2) {
194 av_log(NULL
, AV_LOG_ERROR
, "%d/%d %d/%d, %d %d\n", a
.num
,
195 a
.den
, b
.num
, b
.den
, acc
, exact
);