00001 /* 00002 * rational numbers 00003 * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at> 00004 * 00005 * This file is part of FFmpeg. 00006 * 00007 * FFmpeg is free software; you can redistribute it and/or 00008 * modify it under the terms of the GNU Lesser General Public 00009 * License as published by the Free Software Foundation; either 00010 * version 2.1 of the License, or (at your option) any later version. 00011 * 00012 * FFmpeg is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00015 * Lesser General Public License for more details. 00016 * 00017 * You should have received a copy of the GNU Lesser General Public 00018 * License along with FFmpeg; if not, write to the Free Software 00019 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00020 */ 00021 00028 #include "avassert.h" 00029 //#include <math.h> 00030 #include <limits.h> 00031 00032 #include "common.h" 00033 #include "mathematics.h" 00034 #include "rational.h" 00035 00036 int av_reduce(int *dst_num, int *dst_den, 00037 int64_t num, int64_t den, int64_t max) 00038 { 00039 AVRational a0 = { 0, 1 }, a1 = { 1, 0 }; 00040 int sign = (num < 0) ^ (den < 0); 00041 int64_t gcd = av_gcd(FFABS(num), FFABS(den)); 00042 00043 if (gcd) { 00044 num = FFABS(num) / gcd; 00045 den = FFABS(den) / gcd; 00046 } 00047 if (num <= max && den <= max) { 00048 a1 = (AVRational) { num, den }; 00049 den = 0; 00050 } 00051 00052 while (den) { 00053 uint64_t x = num / den; 00054 int64_t next_den = num - den * x; 00055 int64_t a2n = x * a1.num + a0.num; 00056 int64_t a2d = x * a1.den + a0.den; 00057 00058 if (a2n > max || a2d > max) { 00059 if (a1.num) x = (max - a0.num) / a1.num; 00060 if (a1.den) x = FFMIN(x, (max - a0.den) / a1.den); 00061 00062 if (den * (2 * x * a1.den + a0.den) > num * a1.den) 00063 a1 = (AVRational) { x * a1.num + a0.num, x * a1.den + a0.den }; 00064 break; 00065 } 00066 00067 a0 = a1; 00068 a1 = (AVRational) { a2n, a2d }; 00069 num = den; 00070 den = next_den; 00071 } 00072 av_assert2(av_gcd(a1.num, a1.den) <= 1U); 00073 00074 *dst_num = sign ? -a1.num : a1.num; 00075 *dst_den = a1.den; 00076 00077 return den == 0; 00078 } 00079 00080 AVRational av_mul_q(AVRational b, AVRational c) 00081 { 00082 av_reduce(&b.num, &b.den, 00083 b.num * (int64_t) c.num, 00084 b.den * (int64_t) c.den, INT_MAX); 00085 return b; 00086 } 00087 00088 AVRational av_div_q(AVRational b, AVRational c) 00089 { 00090 return av_mul_q(b, (AVRational) { c.den, c.num }); 00091 } 00092 00093 AVRational av_add_q(AVRational b, AVRational c) { 00094 av_reduce(&b.num, &b.den, 00095 b.num * (int64_t) c.den + 00096 c.num * (int64_t) b.den, 00097 b.den * (int64_t) c.den, INT_MAX); 00098 return b; 00099 } 00100 00101 AVRational av_sub_q(AVRational b, AVRational c) 00102 { 00103 return av_add_q(b, (AVRational) { -c.num, c.den }); 00104 } 00105 00106 AVRational av_d2q(double d, int max) 00107 { 00108 AVRational a; 00109 #define LOG2 0.69314718055994530941723212145817656807550013436025 00110 int exponent; 00111 int64_t den; 00112 if (isnan(d)) 00113 return (AVRational) { 0,0 }; 00114 if (isinf(d)) 00115 return (AVRational) { d < 0 ? -1 : 1, 0 }; 00116 exponent = FFMAX( (int)(log(fabs(d) + 1e-20)/LOG2), 0); 00117 den = 1LL << (61 - exponent); 00118 av_reduce(&a.num, &a.den, (int64_t)(d * den + 0.5), den, max); 00119 00120 return a; 00121 } 00122 00123 int av_nearer_q(AVRational q, AVRational q1, AVRational q2) 00124 { 00125 /* n/d is q, a/b is the median between q1 and q2 */ 00126 int64_t a = q1.num * (int64_t)q2.den + q2.num * (int64_t)q1.den; 00127 int64_t b = 2 * (int64_t)q1.den * q2.den; 00128 00129 /* rnd_up(a*d/b) > n => a*d/b > n */ 00130 int64_t x_up = av_rescale_rnd(a, q.den, b, AV_ROUND_UP); 00131 00132 /* rnd_down(a*d/b) < n => a*d/b < n */ 00133 int64_t x_down = av_rescale_rnd(a, q.den, b, AV_ROUND_DOWN); 00134 00135 return ((x_up > q.num) - (x_down < q.num)) * av_cmp_q(q2, q1); 00136 } 00137 00138 int av_find_nearest_q_idx(AVRational q, const AVRational* q_list) 00139 { 00140 int i, nearest_q_idx = 0; 00141 for (i = 0; q_list[i].den; i++) 00142 if (av_nearer_q(q, q_list[i], q_list[nearest_q_idx]) > 0) 00143 nearest_q_idx = i; 00144 00145 return nearest_q_idx; 00146 } 00147 00148 #ifdef TEST 00149 int main(void) 00150 { 00151 AVRational a,b; 00152 for (a.num = -2; a.num <= 2; a.num++) { 00153 for (a.den = -2; a.den <= 2; a.den++) { 00154 for (b.num = -2; b.num <= 2; b.num++) { 00155 for (b.den = -2; b.den <= 2; b.den++) { 00156 int c = av_cmp_q(a,b); 00157 double d = av_q2d(a) == av_q2d(b) ? 00158 0 : (av_q2d(a) - av_q2d(b)); 00159 if (d > 0) d = 1; 00160 else if (d < 0) d = -1; 00161 else if (d != d) d = INT_MIN; 00162 if (c != d) 00163 av_log(0, AV_LOG_ERROR, "%d/%d %d/%d, %d %f\n", a.num, 00164 a.den, b.num, b.den, c,d); 00165 } 00166 } 00167 } 00168 } 00169 return 0; 00170 } 00171 #endif