libavutil/mathematics.c

   1 /*
   2  * Copyright (c) 2005-2012 Michael Niedermayer <michaelni@gmx.at>
   3  *
   4  * This file is part of FFmpeg.
   5  *
   6  * FFmpeg is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Lesser General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2.1 of the License, or (at your option) any later version.
  10  *
  11  * FFmpeg is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public
  17  * License along with FFmpeg; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  19  */
  20
  21 /**
  22  * @file
  23  * miscellaneous math routines and tables
  24  */
  25
  26 #include <stdint.h>
  27 #include <limits.h>
  28
  29 #include "avutil.h"
  30 #include "mathematics.h"
  31 #include "libavutil/intmath.h"
  32 #include "libavutil/common.h"
  33 #include "avassert.h"
  34
  35 /* Stein's binary GCD algorithm:
  36  * https://en.wikipedia.org/wiki/Binary_GCD_algorithm */
  37 int64_t av_gcd(int64_t a, int64_t b) {
  38     int za, zb, k;
  39     int64_t u, v;
  40     if (a == 0)
  41         return b;
  42     if (b == 0)
  43         return a;
  44     za = ff_ctzll(a);
  45     zb = ff_ctzll(b);
  46     k  = FFMIN(za, zb);
  47     u = llabs(a >> za);
  48     v = llabs(b >> zb);
  49     while (u != v) {
  50         if (u > v)
  51             FFSWAP(int64_t, v, u);
  52         v -= u;
  53         v >>= ff_ctzll(v);
  54     }
  55     return (uint64_t)u << k;
  56 }
  57
  58 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
  59 {
  60     int64_t r = 0;
  61     av_assert2(c > 0);
  62     av_assert2(b >=0);
  63     av_assert2((unsigned)(rnd&~AV_ROUND_PASS_MINMAX)<=5 && (rnd&~AV_ROUND_PASS_MINMAX)!=4);
  64
  65     if (c <= 0 || b < 0 || !((unsigned)(rnd&~AV_ROUND_PASS_MINMAX)<=5 && (rnd&~AV_ROUND_PASS_MINMAX)!=4))
  66         return INT64_MIN;
  67
  68     if (rnd & AV_ROUND_PASS_MINMAX) {
  69         if (a == INT64_MIN || a == INT64_MAX)
  70             return a;
  71         rnd -= AV_ROUND_PASS_MINMAX;
  72     }
  73
  74     if (a < 0)
  75         return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
  76
  77     if (rnd == AV_ROUND_NEAR_INF)
  78         r = c / 2;
  79     else if (rnd & 1)
  80         r = c - 1;
  81
  82     if (b <= INT_MAX && c <= INT_MAX) {
  83         if (a <= INT_MAX)
  84             return (a * b + r) / c;
  85         else {
  86             int64_t ad = a / c;
  87             int64_t a2 = (a % c * b + r) / c;
  88             if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
  89                 return INT64_MIN;
  90             return ad * b + a2;
  91         }
  92     } else {
  93 #if 1
  94         uint64_t a0  = a & 0xFFFFFFFF;
  95         uint64_t a1  = a >> 32;
  96         uint64_t b0  = b & 0xFFFFFFFF;
  97         uint64_t b1  = b >> 32;
  98         uint64_t t1  = a0 * b1 + a1 * b0;
  99         uint64_t t1a = t1 << 32;
 100         int i;
 101
 102         a0  = a0 * b0 + t1a;
 103         a1  = a1 * b1 + (t1 >> 32) + (a0 < t1a);
 104         a0 += r;
 105         a1 += a0 < r;
 106
 107         for (i = 63; i >= 0; i--) {
 108             a1 += a1 + ((a0 >> i) & 1);
 109             t1 += t1;
 110             if (c <= a1) {
 111                 a1 -= c;
 112                 t1++;
 113             }
 114         }
 115         if (t1 > INT64_MAX)
 116             return INT64_MIN;
 117         return t1;
 118 #else
 119         /* reference code doing (a*b + r) / c, requires libavutil/integer.h */
 120         AVInteger ai;
 121         ai = av_mul_i(av_int2i(a), av_int2i(b));
 122         ai = av_add_i(ai, av_int2i(r));
 123
 124         return av_i2int(av_div_i(ai, av_int2i(c)));
 125 #endif
 126     }
 127 }
 128
 129 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
 130 {
 131     return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
 132 }
 133
 134 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
 135                          enum AVRounding rnd)
 136 {
 137     int64_t b = bq.num * (int64_t)cq.den;
 138     int64_t c = cq.num * (int64_t)bq.den;
 139     return av_rescale_rnd(a, b, c, rnd);
 140 }
 141
 142 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
 143 {
 144     return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
 145 }
 146
 147 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
 148 {
 149     int64_t a = tb_a.num * (int64_t)tb_b.den;
 150     int64_t b = tb_b.num * (int64_t)tb_a.den;
 151     if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
 152         return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
 153     if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
 154         return -1;
 155     if (av_rescale_rnd(ts_b, b, a, AV_ROUND_DOWN) < ts_a)
 156         return 1;
 157     return 0;
 158 }
 159
 160 int64_t av_compare_mod(uint64_t a, uint64_t b, uint64_t mod)
 161 {
 162     int64_t c = (a - b) & (mod - 1);
 163     if (c > (mod >> 1))
 164         c -= mod;
 165     return c;
 166 }
 167
 168 int64_t av_rescale_delta(AVRational in_tb, int64_t in_ts,  AVRational fs_tb, int duration, int64_t *last, AVRational out_tb){
 169     int64_t a, b, this;
 170
 171     av_assert0(in_ts != AV_NOPTS_VALUE);
 172     av_assert0(duration >= 0);
 173
 174     if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
 175 simple_round:
 176         *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
 177         return av_rescale_q(in_ts, in_tb, out_tb);
 178     }
 179
 180     a =  av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN)   >>1;
 181     b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP  )+1)>>1;
 182     if (*last < 2*a - b || *last > 2*b - a)
 183         goto simple_round;
 184
 185     this = av_clip64(*last, a, b);
 186     *last = this + duration;
 187
 188     return av_rescale_q(this, fs_tb, out_tb);
 189 }
 190
 191 int64_t av_add_stable(AVRational ts_tb, int64_t ts, AVRational inc_tb, int64_t inc)
 192 {
 193     int64_t m, d;
 194
 195     if (inc != 1)
 196         inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
 197
 198     m = inc_tb.num * (int64_t)ts_tb.den;
 199     d = inc_tb.den * (int64_t)ts_tb.num;
 200
 201     if (m % d == 0 && ts <= INT64_MAX - m / d)
 202         return ts + m / d;
 203     if (m < d)
 204         return ts;
 205
 206     {
 207         int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
 208         int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
 209
 210         if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
 211             return ts;
 212
 213         return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
 214     }
 215 }
 216
 217 static inline double eval_poly(const double *coeff, int size, double x) {
 218     double sum = coeff[size-1];
 219     int i;
 220     for (i = size-2; i >= 0; --i) {
 221         sum *= x;
 222         sum += coeff[i];
 223     }
 224     return sum;
 225 }
 226
 227 /**
 228  * 0th order modified bessel function of the first kind.
 229  * Algorithm taken from the Boost project, source:
 230  * https://searchcode.com/codesearch/view/14918379/
 231  * Use, modification and distribution are subject to the
 232  * Boost Software License, Version 1.0 (see notice below).
 233  * Boost Software License - Version 1.0 - August 17th, 2003
 234 Permission is hereby granted, free of charge, to any person or organization
 235 obtaining a copy of the software and accompanying documentation covered by
 236 this license (the "Software") to use, reproduce, display, distribute,
 237 execute, and transmit the Software, and to prepare derivative works of the
 238 Software, and to permit third-parties to whom the Software is furnished to
 239 do so, all subject to the following:
 240
 241 The copyright notices in the Software and this entire statement, including
 242 the above license grant, this restriction and the following disclaimer,
 243 must be included in all copies of the Software, in whole or in part, and
 244 all derivative works of the Software, unless such copies or derivative
 245 works are solely in the form of machine-executable object code generated by
 246 a source language processor.
 247
 248 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 249 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 250 FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
 251 SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
 252 FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
 253 ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 254 DEALINGS IN THE SOFTWARE.
 255  */
 256
 257 double av_bessel_i0(double x) {
 258 // Modified Bessel function of the first kind of order zero
 259 // minimax rational approximations on intervals, see
 260 // Blair and Edwards, Chalk River Report AECL-4928, 1974
 261     static const double p1[] = {
 262         -2.2335582639474375249e+15,
 263         -5.5050369673018427753e+14,
 264         -3.2940087627407749166e+13,
 265         -8.4925101247114157499e+11,
 266         -1.1912746104985237192e+10,
 267         -1.0313066708737980747e+08,
 268         -5.9545626019847898221e+05,
 269         -2.4125195876041896775e+03,
 270         -7.0935347449210549190e+00,
 271         -1.5453977791786851041e-02,
 272         -2.5172644670688975051e-05,
 273         -3.0517226450451067446e-08,
 274         -2.6843448573468483278e-11,
 275         -1.5982226675653184646e-14,
 276         -5.2487866627945699800e-18,
 277     };
 278     static const double q1[] = {
 279         -2.2335582639474375245e+15,
 280          7.8858692566751002988e+12,
 281         -1.2207067397808979846e+10,
 282          1.0377081058062166144e+07,
 283         -4.8527560179962773045e+03,
 284          1.0,
 285     };
 286     static const double p2[] = {
 287         -2.2210262233306573296e-04,
 288          1.3067392038106924055e-02,
 289         -4.4700805721174453923e-01,
 290          5.5674518371240761397e+00,
 291         -2.3517945679239481621e+01,
 292          3.1611322818701131207e+01,
 293         -9.6090021968656180000e+00,
 294     };
 295     static const double q2[] = {
 296         -5.5194330231005480228e-04,
 297          3.2547697594819615062e-02,
 298         -1.1151759188741312645e+00,
 299          1.3982595353892851542e+01,
 300         -6.0228002066743340583e+01,
 301          8.5539563258012929600e+01,
 302         -3.1446690275135491500e+01,
 303         1.0,
 304     };
 305     double y, r, factor;
 306     if (x == 0)
 307         return 1.0;
 308     x = fabs(x);
 309     if (x <= 15) {
 310         y = x * x;
 311         return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
 312     }
 313     else {
 314         y = 1 / x - 1.0 / 15;
 315         r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
 316         factor = exp(x) / sqrt(x);
 317         return factor * r;
 318     }
 319 }