mathematics.c (10002B)
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 }