FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2026-05-02 17:52:23
Exec Total Coverage
Lines: 122 125 97.6%
Functions: 11 11 100.0%
Branches: 79 92 85.9%
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 55182761 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 26605563 times.
✓ Branch 1 taken 28577198 times.
 55182761 if (a == 0)
41 26605563 return b;
42
2/2
✓ Branch 0 taken 3813 times.
✓ Branch 1 taken 28573385 times.
 28577198 if (b == 0)
43 3813 return a;
44 28573385 za = ff_ctzll(a);
45 28573385 zb = ff_ctzll(b);
46 28573385 k = FFMIN(za, zb);
47 28573385 u = llabs(a >> za);
48 28573385 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 168134391 times.
✓ Branch 1 taken 28573385 times.
 196707776 while (u != v) {
50
2/2
✓ Branch 0 taken 61448170 times.
✓ Branch 1 taken 106686221 times.
 168134391 if (u > v)
51 61448170 FFSWAP(int64_t, v, u);
52 168134391 v -= u;
53 168134391 v >>= ff_ctzll(v);
54 }
55 28573385 return (uint64_t)u << k;
56 }
57
58 16749390 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 16749390 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 16749390 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 16749390 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 16749390 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 16749390 times.
 16749390 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 711046 times.
✓ Branch 1 taken 16038344 times.
 16749390 if (rnd & AV_ROUND_PASS_MINMAX) {
69
4/4
✓ Branch 0 taken 709585 times.
✓ Branch 1 taken 1461 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 709584 times.
 711046 if (a == INT64_MIN || a == INT64_MAX)
70 1462 return a;
71 709584 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 115087 times.
✓ Branch 1 taken 16632841 times.
 16747928 if (a < 0)
75 115087 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 15345416 times.
✓ Branch 1 taken 1287425 times.
 16632841 if (rnd == AV_ROUND_NEAR_INF)
78 15345416 r = c / 2;
79
2/2
✓ Branch 0 taken 339430 times.
✓ Branch 1 taken 947995 times.
 1287425 else if (rnd & 1)
80 339430 r = c - 1;
81
82
4/4
✓ Branch 0 taken 16479292 times.
✓ Branch 1 taken 153549 times.
✓ Branch 2 taken 16367945 times.
✓ Branch 3 taken 111347 times.
 16632841 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 16282180 times.
✓ Branch 1 taken 85765 times.
 16367945 if (a <= INT_MAX)
84 16282180 return (a * b + r) / c;
85 else {
86 85765 int64_t ad = a / c;
87 85765 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 1254 times.
✓ Branch 1 taken 84511 times.
✓ Branch 2 taken 1254 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1253 times.
 85765 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 1 return INT64_MIN;
90 85764 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 264896 uint64_t a0 = a & 0xFFFFFFFF;
95 264896 uint64_t a1 = a >> 32;
96 264896 uint64_t b0 = b & 0xFFFFFFFF;
97 264896 uint64_t b1 = b >> 32;
98 264896 uint64_t t1 = a0 * b1 + a1 * b0;
99 264896 uint64_t t1a = t1 << 32;
100 int i;
101
102 264896 a0 = a0 * b0 + t1a;
103 264896 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 264896 a0 += r;
105 264896 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 16953344 times.
✓ Branch 1 taken 264896 times.
 17218240 for (i = 63; i >= 0; i--) {
108 16953344 a1 += a1 + ((a0 >> i) & 1);
109 16953344 t1 += t1;
110
2/2
✓ Branch 0 taken 2115419 times.
✓ Branch 1 taken 14837925 times.
 16953344 if (c <= a1) {
111 2115419 a1 -= c;
112 2115419 t1++;
113 }
114 }
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 264896 times.
 264896 if (t1 > INT64_MAX)
116 ✗ return INT64_MIN;
117 264896 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 1171609 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 1171609 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 14663562 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 14663562 int64_t b = bq.num * (int64_t)cq.den;
138 14663562 int64_t c = cq.num * (int64_t)bq.den;
139 14663562 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 13464296 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 13464296 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 1035100 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 1035100 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 1035100 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 454536 times.
✓ Branch 1 taken 580564 times.
✓ Branch 2 taken 5320 times.
✓ Branch 3 taken 1029780 times.
✓ Branch 4 taken 1027194 times.
✓ Branch 5 taken 7906 times.
 1035100 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 1027194 return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
153
2/2
✓ Branch 0 taken 6648 times.
✓ Branch 1 taken 1258 times.
 7906 if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
154 6648 return -1;
155
2/2
✓ Branch 0 taken 194 times.
✓ Branch 1 taken 1064 times.
 1258 if (av_rescale_rnd(ts_b, b, a, AV_ROUND_DOWN) < ts_a)
156 194 return 1;
157 1064 return 0;
158 }
159
160 3 int64_t av_compare_mod(uint64_t a, uint64_t b, uint64_t mod)
161 {
162 3 int64_t c = (a - b) & (mod - 1);
163
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
 3 if (c > (mod >> 1))
164 1 c -= mod;
165 3 return c;
166 }
167
168 346862 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 346862 times.
 346862 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 346862 times.
 346862 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1859 times.
✓ Branch 1 taken 345003 times.
✓ Branch 2 taken 16656 times.
✓ Branch 3 taken 328347 times.
✓ Branch 4 taken 274894 times.
✓ Branch 5 taken 53453 times.
 346862 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 293409 simple_round:
176 294040 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 294040 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 53453 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 53453 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 53453 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 631 times.
✓ Branch 3 taken 52822 times.
 53453 if (*last < 2*a - b || *last > 2*b - a)
183 631 goto simple_round;
184
185 52822 this = av_clip64(*last, a, b);
186 52822 *last = this + duration;
187
188 52822 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 569163 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 559159 times.
 569163 if (inc != 1)
196 10004 inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 569163 m = inc_tb.num * (int64_t)ts_tb.den;
199 569163 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 555316 times.
✓ Branch 1 taken 13847 times.
✓ Branch 2 taken 555316 times.
✗ Branch 3 not taken.
 569163 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 555316 return ts + m / d;
203
2/2
✓ Branch 0 taken 2823 times.
✓ Branch 1 taken 11024 times.
 13847 if (m < d)
204 2823 return ts;
205
206 {
207 11024 int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208 11024 int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209
210
3/6
✓ Branch 0 taken 11024 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 11024 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 11024 times.
 11024 if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211 ✗ return ts;
212
213 11024 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }
215 }
216
217 85012934 static inline double eval_poly(const double *coeff, int size, double x) {
218 85012934 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 807469033 times.
✓ Branch 1 taken 85012934 times.
 892481967 for (i = size-2; i >= 0; --i) {
221 807469033 sum *= x;
222 807469033 sum += coeff[i];
223 }
224 85012934 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 42508318 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 1851 times.
✓ Branch 1 taken 42506467 times.
 42508318 if (x == 0)
307 1851 return 1.0;
308 42506467 x = fabs(x);
309
2/2
✓ Branch 0 taken 42480827 times.
✓ Branch 1 taken 25640 times.
 42506467 if (x <= 15) {
310 42480827 y = x * x;
311 42480827 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 25640 y = 1 / x - 1.0 / 15;
315 25640 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 25640 factor = exp(x) / sqrt(x);
317 25640 return factor * r;
318 }
319 }
320