FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2026-10-10 14:24:39
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 21204687 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 3171191 times.
✓ Branch 1 taken 18033496 times.
 21204687 if (a == 0)
41 3171191 return b;
42
2/2
✓ Branch 0 taken 4097 times.
✓ Branch 1 taken 18029399 times.
 18033496 if (b == 0)
43 4097 return a;
44 18029399 za = ff_ctzll(a);
45 18029399 zb = ff_ctzll(b);
46 18029399 k = FFMIN(za, zb);
47 18029399 u = llabs(a >> za);
48 18029399 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 90042585 times.
✓ Branch 1 taken 18029399 times.
 108071984 while (u != v) {
50
2/2
✓ Branch 0 taken 55633619 times.
✓ Branch 1 taken 34408966 times.
 90042585 if (u > v)
51 55633619 FFSWAP(int64_t, v, u);
52 90042585 v -= u;
53 90042585 v >>= ff_ctzll(v);
54 }
55 18029399 return (uint64_t)u << k;
56 }
57
58 25192328 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 25192328 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 25192328 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 25192328 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 25192328 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 25192328 times.
 25192328 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 739911 times.
✓ Branch 1 taken 24452417 times.
 25192328 if (rnd & AV_ROUND_PASS_MINMAX) {
69
4/4
✓ Branch 0 taken 738115 times.
✓ Branch 1 taken 1796 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 738114 times.
 739911 if (a == INT64_MIN || a == INT64_MAX)
70 1797 return a;
71 738114 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 138785 times.
✓ Branch 1 taken 25051746 times.
 25190531 if (a < 0)
75 138785 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 23710821 times.
✓ Branch 1 taken 1340925 times.
 25051746 if (rnd == AV_ROUND_NEAR_INF)
78 23710821 r = c / 2;
79
2/2
✓ Branch 0 taken 364054 times.
✓ Branch 1 taken 976871 times.
 1340925 else if (rnd & 1)
80 364054 r = c - 1;
81
82
4/4
✓ Branch 0 taken 24880347 times.
✓ Branch 1 taken 171399 times.
✓ Branch 2 taken 24767283 times.
✓ Branch 3 taken 113064 times.
 25051746 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 24681133 times.
✓ Branch 1 taken 86150 times.
 24767283 if (a <= INT_MAX)
84 24681133 return (a * b + r) / c;
85 else {
86 86150 int64_t ad = a / c;
87 86150 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 1327 times.
✓ Branch 1 taken 84823 times.
✓ Branch 2 taken 1327 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 19 times.
✓ Branch 5 taken 1308 times.
 86150 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 19 return INT64_MIN;
90 86131 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 284463 uint64_t a0 = a & 0xFFFFFFFF;
95 284463 uint64_t a1 = a >> 32;
96 284463 uint64_t b0 = b & 0xFFFFFFFF;
97 284463 uint64_t b1 = b >> 32;
98 284463 uint64_t t1 = a0 * b1 + a1 * b0;
99 284463 uint64_t t1a = t1 << 32;
100 int i;
101
102 284463 a0 = a0 * b0 + t1a;
103 284463 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 284463 a0 += r;
105 284463 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 18205632 times.
✓ Branch 1 taken 284463 times.
 18490095 for (i = 63; i >= 0; i--) {
108 18205632 a1 += a1 + ((a0 >> i) & 1);
109 18205632 t1 += t1;
110
2/2
✓ Branch 0 taken 2364633 times.
✓ Branch 1 taken 15840999 times.
 18205632 if (c <= a1) {
111 2364633 a1 -= c;
112 2364633 t1++;
113 }
114 }
115
2/2
✓ Branch 0 taken 134 times.
✓ Branch 1 taken 284329 times.
 284463 if (t1 > INT64_MAX)
116 134 return INT64_MIN;
117 284329 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 1473056 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 1473056 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 22751030 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 22751030 int64_t b = bq.num * (int64_t)cq.den;
138 22751030 int64_t c = cq.num * (int64_t)bq.den;
139 22751030 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 21499724 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 21499724 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 1569282 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 1569282 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 1569282 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 476266 times.
✓ Branch 1 taken 1093016 times.
✓ Branch 2 taken 9326 times.
✓ Branch 3 taken 1559956 times.
✓ Branch 4 taken 1561732 times.
✓ Branch 5 taken 7550 times.
 1569282 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 1561732 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 299 int64_t av_compare_mod(uint64_t a, uint64_t b, uint64_t mod)
161 {
162 299 int64_t c = (a - b) & (mod - 1);
163
2/2
✓ Branch 0 taken 212 times.
✓ Branch 1 taken 87 times.
 299 if (c > (mod >> 1))
164 212 c -= mod;
165 299 return c;
166 }
167
168 367988 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 367988 times.
 367988 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 367988 times.
 367988 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1593 times.
✓ Branch 1 taken 366395 times.
✓ Branch 2 taken 17018 times.
✓ Branch 3 taken 349377 times.
✓ Branch 4 taken 295312 times.
✓ Branch 5 taken 54065 times.
 367988 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 313923 simple_round:
176 314574 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 314574 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 54065 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 54065 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 54065 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 651 times.
✓ Branch 3 taken 53414 times.
 54065 if (*last < 2*a - b || *last > 2*b - a)
183 651 goto simple_round;
184
185 53414 this = av_clip64(*last, a, b);
186 53414 *last = this + duration;
187
188 53414 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 2616746 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 2606742 times.
 2616746 if (inc != 1)
196 10004 inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 2616746 m = inc_tb.num * (int64_t)ts_tb.den;
199 2616746 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 998965 times.
✓ Branch 1 taken 1617781 times.
✓ Branch 2 taken 998965 times.
✗ Branch 3 not taken.
 2616746 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 998965 return ts + m / d;
203
2/2
✓ Branch 0 taken 2823 times.
✓ Branch 1 taken 1614958 times.
 1617781 if (m < d)
204 2823 return ts;
205
206 {
207 1614958 int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208 1614958 int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209
210
3/6
✓ Branch 0 taken 1614958 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1614958 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1614958 times.
 1614958 if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211 ✗ return ts;
212
213 1614958 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }
215 }
216
217 87824530 static inline double eval_poly(const double *coeff, int size, double x) {
218 87824530 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 834120065 times.
✓ Branch 1 taken 87824530 times.
 921944595 for (i = size-2; i >= 0; --i) {
221 834120065 sum *= x;
222 834120065 sum += coeff[i];
223 }
224 87824530 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 43914814 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 2549 times.
✓ Branch 1 taken 43912265 times.
 43914814 if (x == 0)
307 2549 return 1.0;
308 43912265 x = fabs(x);
309
2/2
✓ Branch 0 taken 43876770 times.
✓ Branch 1 taken 35495 times.
 43912265 if (x <= 15) {
310 43876770 y = x * x;
311 43876770 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 35495 y = 1 / x - 1.0 / 15;
315 35495 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 35495 factor = exp(x) / sqrt(x);
317 35495 return factor * r;
318 }
319 }
320