FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2026-11-01 05:14:56
Exec Total Coverage
Lines: 123 125 98.4%
Functions: 11 11 100.0%
Branches: 80 92 87.0%
Line Branch Exec Source
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 33473293 int64_t av_gcd(int64_t a, int64_t b) {
38 int za, zb, k;
39 int64_t u, v;
40
2/2
✓ Branch 0 taken 4603792 times.
✓ Branch 1 taken 28869501 times.
 33473293 if (a == 0)
41 4603792 return b;
42
2/2
✓ Branch 0 taken 4404 times.
✓ Branch 1 taken 28865097 times.
 28869501 if (b == 0)
43 4404 return a;
44 28865097 za = ff_ctzll(a);
45 28865097 zb = ff_ctzll(b);
46 28865097 k = FFMIN(za, zb);
47 28865097 u = llabs(a >> za);
48 28865097 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 146972554 times.
✓ Branch 1 taken 28865097 times.
 175837651 while (u != v) {
50
2/2
✓ Branch 0 taken 91069331 times.
✓ Branch 1 taken 55903223 times.
 146972554 if (u > v)
51 91069331 FFSWAP(int64_t, v, u);
52 146972554 v -= u;
53 146972554 v >>= ff_ctzll(v);
54 }
55 28865097 return (uint64_t)u << k;
56 }
57
58 25292645 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 25292645 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
4/8
✓ Branch 0 taken 25292645 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 25292645 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 25292645 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 25292645 times.
 25292645 if (c <= 0 || b < 0 || !((unsigned)(rnd&~AV_ROUND_PASS_MINMAX)<=5 && (rnd&~AV_ROUND_PASS_MINMAX)!=4))
66 ✗ return INT64_MIN;
67
68
2/2
✓ Branch 0 taken 744742 times.
✓ Branch 1 taken 24547903 times.
 25292645 if (rnd & AV_ROUND_PASS_MINMAX) {
69
4/4
✓ Branch 0 taken 742763 times.
✓ Branch 1 taken 1979 times.
✓ Branch 2 taken 34 times.
✓ Branch 3 taken 742729 times.
 744742 if (a == INT64_MIN || a == INT64_MAX)
70 2013 return a;
71 742729 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 143555 times.
✓ Branch 1 taken 25147077 times.
 25290632 if (a < 0)
75 143555 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 23797882 times.
✓ Branch 1 taken 1349195 times.
 25147077 if (rnd == AV_ROUND_NEAR_INF)
78 23797882 r = c / 2;
79
2/2
✓ Branch 0 taken 365096 times.
✓ Branch 1 taken 984099 times.
 1349195 else if (rnd & 1)
80 365096 r = c - 1;
81
82
4/4
✓ Branch 0 taken 24974810 times.
✓ Branch 1 taken 172267 times.
✓ Branch 2 taken 24861747 times.
✓ Branch 3 taken 113063 times.
 25147077 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 24775377 times.
✓ Branch 1 taken 86370 times.
 24861747 if (a <= INT_MAX)
84 24775377 return (a * b + r) / c;
85 else {
86 86370 int64_t ad = a / c;
87 86370 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 1327 times.
✓ Branch 1 taken 85043 times.
✓ Branch 2 taken 1327 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 19 times.
✓ Branch 5 taken 1308 times.
 86370 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 19 return INT64_MIN;
90 86351 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 285330 uint64_t a0 = a & 0xFFFFFFFF;
95 285330 uint64_t a1 = a >> 32;
96 285330 uint64_t b0 = b & 0xFFFFFFFF;
97 285330 uint64_t b1 = b >> 32;
98 285330 uint64_t t1 = a0 * b1 + a1 * b0;
99 285330 uint64_t t1a = t1 << 32;
100 int i;
101
102 285330 a0 = a0 * b0 + t1a;
103 285330 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 285330 a0 += r;
105 285330 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 18261120 times.
✓ Branch 1 taken 285330 times.
 18546450 for (i = 63; i >= 0; i--) {
108 18261120 a1 += a1 + ((a0 >> i) & 1);
109 18261120 t1 += t1;
110
2/2
✓ Branch 0 taken 2368169 times.
✓ Branch 1 taken 15892951 times.
 18261120 if (c <= a1) {
111 2368169 a1 -= c;
112 2368169 t1++;
113 }
114 }
115
2/2
✓ Branch 0 taken 134 times.
✓ Branch 1 taken 285196 times.
 285330 if (t1 > INT64_MAX)
116 134 return INT64_MIN;
117 285196 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 1481897 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 1481897 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 22832907 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 22832907 int64_t b = bq.num * (int64_t)cq.den;
138 22832907 int64_t c = cq.num * (int64_t)bq.den;
139 22832907 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 21573398 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 21573398 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 1566032 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 1566032 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 1566032 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 479140 times.
✓ Branch 1 taken 1086892 times.
✓ Branch 2 taken 11336 times.
✓ Branch 3 taken 1554696 times.
✓ Branch 4 taken 1558482 times.
✓ Branch 5 taken 7550 times.
 1566032 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 1558482 return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
153
2/2
✓ Branch 0 taken 6276 times.
✓ Branch 1 taken 1274 times.
 7550 if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
154 6276 return -1;
155
2/2
✓ Branch 0 taken 202 times.
✓ Branch 1 taken 1072 times.
 1274 if (av_rescale_rnd(ts_b, b, a, AV_ROUND_DOWN) < ts_a)
156 202 return 1;
157 1072 return 0;
158 }
159
160 321 int64_t av_compare_mod(uint64_t a, uint64_t b, uint64_t mod)
161 {
162 321 int64_t c = (a - b) & (mod - 1);
163
2/2
✓ Branch 0 taken 220 times.
✓ Branch 1 taken 101 times.
 321 if (c > (mod >> 1))
164 220 c -= mod;
165 321 return c;
166 }
167
168 369086 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 369086 times.
 369086 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 369086 times.
 369086 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1647 times.
✓ Branch 1 taken 367439 times.
✓ Branch 2 taken 17123 times.
✓ Branch 3 taken 350316 times.
✓ Branch 4 taken 296138 times.
✓ Branch 5 taken 54178 times.
 369086 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 314908 simple_round:
176 315559 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 315559 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 54178 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 54178 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 54178 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 651 times.
✓ Branch 3 taken 53527 times.
 54178 if (*last < 2*a - b || *last > 2*b - a)
183 651 goto simple_round;
184
185 53527 this = av_clip64(*last, a, b);
186 53527 *last = this + duration;
187
188 53527 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 2624725 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
2/2
✓ Branch 0 taken 10004 times.
✓ Branch 1 taken 2614721 times.
 2624725 if (inc != 1)
196 10004 inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 2624725 m = inc_tb.num * (int64_t)ts_tb.den;
199 2624725 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 1005223 times.
✓ Branch 1 taken 1619502 times.
✓ Branch 2 taken 1005223 times.
✗ Branch 3 not taken.
 2624725 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 1005223 return ts + m / d;
203
2/2
✓ Branch 0 taken 2823 times.
✓ Branch 1 taken 1616679 times.
 1619502 if (m < d)
204 2823 return ts;
205
206 {
207 1616679 int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208 1616679 int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209
210
3/6
✓ Branch 0 taken 1616679 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1616679 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1616679 times.
 1616679 if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211 ✗ return ts;
212
213 1616679 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }
215 }
216
217 88003622 static inline double eval_poly(const double *coeff, int size, double x) {
218 88003622 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 835783279 times.
✓ Branch 1 taken 88003622 times.
 923786901 for (i = size-2; i >= 0; --i) {
221 835783279 sum *= x;
222 835783279 sum += coeff[i];
223 }
224 88003622 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 44004688 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
2/2
✓ Branch 0 taken 2877 times.
✓ Branch 1 taken 44001811 times.
 44004688 if (x == 0)
307 2877 return 1.0;
308 44001811 x = fabs(x);
309
2/2
✓ Branch 0 taken 43959956 times.
✓ Branch 1 taken 41855 times.
 44001811 if (x <= 15) {
310 43959956 y = x * x;
311 43959956 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 41855 y = 1 / x - 1.0 / 15;
315 41855 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 41855 factor = exp(x) / sqrt(x);
317 41855 return factor * r;
318 }
319 }
320