rational.c (5332B)
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 << (62 - 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 122 return a; 123 } 124 125 int av_nearer_q(AVRational q, AVRational q1, AVRational q2) 126 { 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; 130 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); 133 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); 136 137 return ((x_up > q.num) - (x_down < q.num)) * av_cmp_q(q2, q1); 138 } 139 140 int av_find_nearest_q_idx(AVRational q, const AVRational* q_list) 141 { 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) 145 nearest_q_idx = i; 146 147 return nearest_q_idx; 148 } 149 150 uint32_t av_q2intfloat(AVRational q) { 151 int64_t n; 152 int shift; 153 int sign = 0; 154 155 if (q.den < 0) { 156 q.den *= -1; 157 q.num *= -1; 158 } 159 if (q.num < 0) { 160 q.num *= -1; 161 sign = 1; 162 } 163 164 if (!q.num && !q.den) return 0xFFC00000; 165 if (!q.num) return 0; 166 if (!q.den) return 0x7F800000 | (q.num & 0x80000000); 167 168 shift = 23 + av_log2(q.den) - av_log2(q.num); 169 if (shift >= 0) n = av_rescale(q.num, 1LL<<shift, q.den); 170 else n = av_rescale(q.num, 1, ((int64_t)q.den) << -shift); 171 172 shift -= n >= (1<<24); 173 shift += n < (1<<23); 174 175 if (shift >= 0) n = av_rescale(q.num, 1LL<<shift, q.den); 176 else n = av_rescale(q.num, 1, ((int64_t)q.den) << -shift); 177 178 av_assert1(n < (1<<24)); 179 av_assert1(n >= (1<<23)); 180 181 return sign<<31 | (150-shift)<<23 | (n - (1<<23)); 182 } 183 184 AVRational av_gcd_q(AVRational a, AVRational b, int max_den, AVRational def) 185 { 186 int64_t gcd, lcm; 187 188 gcd = av_gcd(a.den, b.den); 189 lcm = (a.den / gcd) * b.den; 190 return lcm < max_den ? av_make_q(av_gcd(a.num, b.num), lcm) : def; 191 }