rational.c 5.1 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
    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 81 82 83 84
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
85 86 87
    return b;
}

88 89 90
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
91 92
}

93 94 95 96 97
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
98 99 100
    return b;
}

101 102 103
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
104 105
}

106 107
AVRational av_d2q(double d, int max)
{
Michael Niedermayer's avatar
Michael Niedermayer committed
108
    AVRational a;
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
    frexp(d, &exponent);
    exponent = FFMAX(exponent-1, 0);
117
    den = 1LL << (61 - exponent);
118 119
    // (int64_t)rint() and llrint() do not work with gcc on ia64 and sparc64,
    // see Ticket2713 for affected gcc/glibc versions
120
    av_reduce(&a.num, &a.den, floor(d * den + 0.5), den, max);
121
    if ((!a.num || !a.den) && d && max>0 && max<INT_MAX)
122
        av_reduce(&a.num, &a.den, floor(d * den + 0.5), den, INT_MAX);
Michael Niedermayer's avatar
Michael Niedermayer committed
123 124 125

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

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;
145
    for (i = 0; q_list[i].den; i++)
146 147 148 149 150
        if (av_nearer_q(q, q_list[i], q_list[nearest_q_idx]) > 0)
            nearest_q_idx = i;

    return nearest_q_idx;
}
151

152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
uint32_t av_q2intfloat(AVRational q) {
    int64_t n;
    int shift;
    int sign = 0;

    if (q.den < 0) {
        q.den *= -1;
        q.num *= -1;
    }
    if (q.num < 0) {
        q.num *= -1;
        sign = 1;
    }

    if (!q.num && !q.den) return 0xFFC00000;
    if (!q.num) return 0;
    if (!q.den) return 0x7F800000 | (q.num & 0x80000000);

    shift = 23 + av_log2(q.den) - av_log2(q.num);
    if (shift >= 0) n = av_rescale(q.num, 1LL<<shift, q.den);
    else            n = av_rescale(q.num, 1, ((int64_t)q.den) << -shift);

    shift -= n >= (1<<24);
    shift += n <  (1<<23);

    if (shift >= 0) n = av_rescale(q.num, 1LL<<shift, q.den);
    else            n = av_rescale(q.num, 1, ((int64_t)q.den) << -shift);

    av_assert1(n <  (1<<24));
    av_assert1(n >= (1<<23));

    return sign<<31 | (150-shift)<<23 | (n - (1<<23));
}