rational.c 6.37 KB
Newer Older
Michael Niedermayer's avatar
Michael Niedermayer committed
1
/*
2
 * rational numbers
Michael Niedermayer's avatar
Michael Niedermayer committed
3 4
 * Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at>
 *
5 6 7
 * This file is part of FFmpeg.
 *
 * FFmpeg is free software; you can redistribute it and/or
Michael Niedermayer's avatar
Michael Niedermayer committed
8 9
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
10
 * version 2.1 of the License, or (at your option) any later version.
Michael Niedermayer's avatar
Michael Niedermayer committed
11
 *
12
 * FFmpeg is distributed in the hope that it will be useful,
Michael Niedermayer's avatar
Michael Niedermayer committed
13 14 15 16 17
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * 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
Michael Niedermayer's avatar
Michael Niedermayer committed
20
 */
21

Michael Niedermayer's avatar
Michael Niedermayer committed
22
/**
23
 * @file
24
 * rational numbers
Michael Niedermayer's avatar
Michael Niedermayer committed
25 26 27
 * @author Michael Niedermayer <michaelni@gmx.at>
 */

28
#include "avassert.h"
Michael Niedermayer's avatar
Michael Niedermayer committed
29
#include <limits.h>
30

Michael Niedermayer's avatar
Michael Niedermayer committed
31
#include "common.h"
32
#include "mathematics.h"
Michael Niedermayer's avatar
Michael Niedermayer committed
33 34
#include "rational.h"

35 36 37 38 39 40 41 42 43 44
int av_reduce(int *dst_num, int *dst_den,
              int64_t num, int64_t den, int64_t max)
{
    AVRational a0 = { 0, 1 }, a1 = { 1, 0 };
    int sign = (num < 0) ^ (den < 0);
    int64_t gcd = av_gcd(FFABS(num), FFABS(den));

    if (gcd) {
        num = FFABS(num) / gcd;
        den = FFABS(den) / gcd;
45
    }
46 47 48
    if (num <= max && den <= max) {
        a1 = (AVRational) { num, den };
        den = 0;
49
    }
50

51 52 53 54 55
    while (den) {
        uint64_t x        = num / den;
        int64_t next_den  = num - den * x;
        int64_t a2n       = x * a1.num + a0.num;
        int64_t a2d       = x * a1.den + a0.den;
56

57 58 59
        if (a2n > max || a2d > max) {
            if (a1.num) x =          (max - a0.num) / a1.num;
            if (a1.den) x = FFMIN(x, (max - a0.den) / a1.den);
60

61 62
            if (den * (2 * x * a1.den + a0.den) > num * a1.den)
                a1 = (AVRational) { x * a1.num + a0.num, x * a1.den + a0.den };
63 64
            break;
        }
65

66 67 68 69
        a0  = a1;
        a1  = (AVRational) { a2n, a2d };
        num = den;
        den = next_den;
70
    }
71
    av_assert2(av_gcd(a1.num, a1.den) <= 1U);
72

73
    *dst_num = sign ? -a1.num : a1.num;
74
    *dst_den = a1.den;
75

76
    return den == 0;
77 78
}

79 80 81 82 83
AVRational av_mul_q(AVRational b, AVRational c)
{
    av_reduce(&b.num, &b.den,
               b.num * (int64_t) c.num,
               b.den * (int64_t) c.den, INT_MAX);
Michael Niedermayer's avatar
Michael Niedermayer committed
84 85 86
    return b;
}

87 88 89
AVRational av_div_q(AVRational b, AVRational c)
{
    return av_mul_q(b, (AVRational) { c.den, c.num });
Michael Niedermayer's avatar
Michael Niedermayer committed
90 91
}

92 93 94 95 96
AVRational av_add_q(AVRational b, AVRational c) {
    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);
Michael Niedermayer's avatar
Michael Niedermayer committed
97 98 99
    return b;
}

100 101 102
AVRational av_sub_q(AVRational b, AVRational c)
{
    return av_add_q(b, (AVRational) { -c.num, c.den });
Michael Niedermayer's avatar
Michael Niedermayer committed
103 104
}

105 106
AVRational av_d2q(double d, int max)
{
Michael Niedermayer's avatar
Michael Niedermayer committed
107
    AVRational a;
Michael Niedermayer's avatar
Michael Niedermayer committed
108
#define LOG2  0.69314718055994530941723212145817656807550013436025
109 110
    int exponent;
    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 116
    exponent = FFMAX( (int)(log(fabs(d) + 1e-20)/LOG2), 0);
    den = 1LL << (61 - exponent);
117 118
    // (int64_t)rint() and llrint() do not work with gcc on ia64 and sparc64
    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);
Michael Niedermayer's avatar
Michael Niedermayer committed
121 122 123

    return a;
}
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142

int av_nearer_q(AVRational q, AVRational q1, AVRational q2)
{
    /* n/d is q, a/b is the median between q1 and q2 */
    int64_t a = q1.num * (int64_t)q2.den + q2.num * (int64_t)q1.den;
    int64_t b = 2 * (int64_t)q1.den * q2.den;

    /* rnd_up(a*d/b) > n => a*d/b > n */
    int64_t x_up = av_rescale_rnd(a, q.den, b, AV_ROUND_UP);

    /* rnd_down(a*d/b) < n => a*d/b < n */
    int64_t x_down = av_rescale_rnd(a, q.den, b, AV_ROUND_DOWN);

    return ((x_up > q.num) - (x_down < q.num)) * av_cmp_q(q2, q1);
}

int av_find_nearest_q_idx(AVRational q, const AVRational* q_list)
{
    int i, nearest_q_idx = 0;
143
    for (i = 0; q_list[i].den; i++)
144 145 146 147 148
        if (av_nearer_q(q, q_list[i], q_list[nearest_q_idx]) > 0)
            nearest_q_idx = i;

    return nearest_q_idx;
}
149 150

#ifdef TEST
151 152
int main(void)
{
153
    AVRational a,b,r;
154 155 156 157 158 159 160 161 162 163 164
    for (a.num = -2; a.num <= 2; a.num++) {
        for (a.den = -2; a.den <= 2; a.den++) {
            for (b.num = -2; b.num <= 2; b.num++) {
                for (b.den = -2; b.den <= 2; b.den++) {
                    int c = av_cmp_q(a,b);
                    double d = av_q2d(a) == av_q2d(b) ?
                               0 : (av_q2d(a) - av_q2d(b));
                    if (d > 0)       d = 1;
                    else if (d < 0)  d = -1;
                    else if (d != d) d = INT_MIN;
                    if (c != d)
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 168
                    r = av_sub_q(av_add_q(b,a), b);
                    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);
170 171 172 173
                }
            }
        }
    }
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203

    for (a.num = 1; a.num <= 10; a.num++) {
        for (a.den = 1; a.den <= 10; a.den++) {
            if (av_gcd(a.num, a.den) > 1)
                continue;
            for (b.num = 1; b.num <= 10; b.num++) {
                for (b.den = 1; b.den <= 10; b.den++) {
                    int start;
                    if (av_gcd(b.num, b.den) > 1)
                        continue;
                    if (av_cmp_q(b, a) < 0)
                        continue;
                    for (start = 0; start < 10 ; start++) {
                        int acc= start;
                        int i;

                        for (i = 0; i<100; i++) {
                            int exact = start + av_rescale_q(i+1, b, a);
                            acc = av_add_stable(a, acc, b, 1);
                            if (FFABS(acc - exact) > 2) {
                                av_log(NULL, AV_LOG_ERROR, "%d/%d %d/%d, %d %d\n", a.num,
                                       a.den, b.num, b.den, acc, exact);
                                return 1;
                            }
                        }
                    }
                }
            }
        }
    }
204
    return 0;
205 206
}
#endif