00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00028
00029 #include <limits.h>
00030
00031 #include "common.h"
00032 #include "mathematics.h"
00033 #include "rational.h"
00034
00035 int av_reduce(int *dst_nom, int *dst_den, int64_t nom, int64_t den, int64_t max){
00036 AVRational a0={0,1}, a1={1,0};
00037 int sign= (nom<0) ^ (den<0);
00038 int64_t gcd= ff_gcd(FFABS(nom), FFABS(den));
00039
00040 if(gcd){
00041 nom = FFABS(nom)/gcd;
00042 den = FFABS(den)/gcd;
00043 }
00044 if(nom<=max && den<=max){
00045 a1= (AVRational){nom, den};
00046 den=0;
00047 }
00048
00049 while(den){
00050 uint64_t x = nom / den;
00051 int64_t next_den= nom - den*x;
00052 int64_t a2n= x*a1.num + a0.num;
00053 int64_t a2d= x*a1.den + a0.den;
00054
00055 if(a2n > max || a2d > max){
00056 if(a1.num) x= (max - a0.num) / a1.num;
00057 if(a1.den) x= FFMIN(x, (max - a0.den) / a1.den);
00058
00059 if (den*(2*x*a1.den + a0.den) > nom*a1.den)
00060 a1 = (AVRational){x*a1.num + a0.num, x*a1.den + a0.den};
00061 break;
00062 }
00063
00064 a0= a1;
00065 a1= (AVRational){a2n, a2d};
00066 nom= den;
00067 den= next_den;
00068 }
00069 assert(ff_gcd(a1.num, a1.den) <= 1U);
00070
00071 *dst_nom = sign ? -a1.num : a1.num;
00072 *dst_den = a1.den;
00073
00074 return den==0;
00075 }
00076
00077 AVRational av_mul_q(AVRational b, AVRational c){
00078 av_reduce(&b.num, &b.den, b.num * (int64_t)c.num, b.den * (int64_t)c.den, INT_MAX);
00079 return b;
00080 }
00081
00082 AVRational av_div_q(AVRational b, AVRational c){
00083 return av_mul_q(b, (AVRational){c.den, c.num});
00084 }
00085
00086 AVRational av_add_q(AVRational b, AVRational c){
00087 av_reduce(&b.num, &b.den, b.num * (int64_t)c.den + c.num * (int64_t)b.den, b.den * (int64_t)c.den, INT_MAX);
00088 return b;
00089 }
00090
00091 AVRational av_sub_q(AVRational b, AVRational c){
00092 return av_add_q(b, (AVRational){-c.num, c.den});
00093 }
00094
00095 AVRational av_d2q(double d, int max){
00096 AVRational a;
00097 #define LOG2 0.69314718055994530941723212145817656807550013436025
00098 int exponent= FFMAX( (int)(log(fabs(d) + 1e-20)/LOG2), 0);
00099 int64_t den= 1LL << (61 - exponent);
00100 av_reduce(&a.num, &a.den, (int64_t)(d * den + 0.5), den, max);
00101
00102 return a;
00103 }