FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2026-09-27 07:50:39
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 13420390 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 1619512 times.
✓ Branch 1 taken 11800878 times.
 13420390 if (a == 0)
41 1619512 return b;
42
2/2
✓ Branch 0 taken 3846 times.
✓ Branch 1 taken 11797032 times.
 11800878 if (b == 0)
43 3846 return a;
44 11797032 za = ff_ctzll(a);
45 11797032 zb = ff_ctzll(b);
46 11797032 k = FFMIN(za, zb);
47 11797032 u = llabs(a >> za);
48 11797032 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 64050435 times.
✓ Branch 1 taken 11797032 times.
 75847467 while (u != v) {
50
2/2
✓ Branch 0 taken 37804436 times.
✓ Branch 1 taken 26245999 times.
 64050435 if (u > v)
51 37804436 FFSWAP(int64_t, v, u);
52 64050435 v -= u;
53 64050435 v >>= ff_ctzll(v);
54 }
55 11797032 return (uint64_t)u << k;
56 }
57
58 17671647 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 17671647 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 17671647 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 17671647 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 17671647 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 17671647 times.
 17671647 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 720471 times.
✓ Branch 1 taken 16951176 times.
 17671647 if (rnd & AV_ROUND_PASS_MINMAX) {
69
4/4
✓ Branch 0 taken 718885 times.
✓ Branch 1 taken 1586 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 718884 times.
 720471 if (a == INT64_MIN || a == INT64_MAX)
70 1587 return a;
71 718884 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 127813 times.
✓ Branch 1 taken 17542247 times.
 17670060 if (a < 0)
75 127813 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 16210344 times.
✓ Branch 1 taken 1331903 times.
 17542247 if (rnd == AV_ROUND_NEAR_INF)
78 16210344 r = c / 2;
79
2/2
✓ Branch 0 taken 359886 times.
✓ Branch 1 taken 972017 times.
 1331903 else if (rnd & 1)
80 359886 r = c - 1;
81
82
4/4
✓ Branch 0 taken 17385715 times.
✓ Branch 1 taken 156532 times.
✓ Branch 2 taken 17272929 times.
✓ Branch 3 taken 112786 times.
 17542247 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 17187147 times.
✓ Branch 1 taken 85782 times.
 17272929 if (a <= INT_MAX)
84 17187147 return (a * b + r) / c;
85 else {
86 85782 int64_t ad = a / c;
87 85782 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 1261 times.
✓ Branch 1 taken 84521 times.
✓ Branch 2 taken 1261 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1260 times.
 85782 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 1 return INT64_MIN;
90 85781 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 269318 uint64_t a0 = a & 0xFFFFFFFF;
95 269318 uint64_t a1 = a >> 32;
96 269318 uint64_t b0 = b & 0xFFFFFFFF;
97 269318 uint64_t b1 = b >> 32;
98 269318 uint64_t t1 = a0 * b1 + a1 * b0;
99 269318 uint64_t t1a = t1 << 32;
100 int i;
101
102 269318 a0 = a0 * b0 + t1a;
103 269318 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 269318 a0 += r;
105 269318 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 17236352 times.
✓ Branch 1 taken 269318 times.
 17505670 for (i = 63; i >= 0; i--) {
108 17236352 a1 += a1 + ((a0 >> i) & 1);
109 17236352 t1 += t1;
110
2/2
✓ Branch 0 taken 2171226 times.
✓ Branch 1 taken 15065126 times.
 17236352 if (c <= a1) {
111 2171226 a1 -= c;
112 2171226 t1++;
113 }
114 }
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 269318 times.
 269318 if (t1 > INT64_MAX)
116 ✗ return INT64_MIN;
117 269318 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 1259737 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 1259737 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 15460242 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 15460242 int64_t b = bq.num * (int64_t)cq.den;
138 15460242 int64_t c = cq.num * (int64_t)bq.den;
139 15460242 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 14231796 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 14231796 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 1048928 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 1048928 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 1048928 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 463997 times.
✓ Branch 1 taken 584931 times.
✓ Branch 2 taken 5815 times.
✓ Branch 3 taken 1043113 times.
✓ Branch 4 taken 1041375 times.
✓ Branch 5 taken 7553 times.
 1048928 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 1041375 return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
153
2/2
✓ Branch 0 taken 6280 times.
✓ Branch 1 taken 1273 times.
 7553 if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
154 6280 return -1;
155
2/2
✓ Branch 0 taken 202 times.
✓ Branch 1 taken 1071 times.
 1273 if (av_rescale_rnd(ts_b, b, a, AV_ROUND_DOWN) < ts_a)
156 202 return 1;
157 1071 return 0;
158 }
159
160 124 int64_t av_compare_mod(uint64_t a, uint64_t b, uint64_t mod)
161 {
162 124 int64_t c = (a - b) & (mod - 1);
163
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 84 times.
 124 if (c > (mod >> 1))
164 40 c -= mod;
165 124 return c;
166 }
167
168 365375 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 365375 times.
 365375 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 365375 times.
 365375 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1581 times.
✓ Branch 1 taken 363794 times.
✓ Branch 2 taken 16998 times.
✓ Branch 3 taken 346796 times.
✓ Branch 4 taken 293108 times.
✓ Branch 5 taken 53688 times.
 365375 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 311687 simple_round:
176 312338 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 312338 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 53688 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 53688 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 53688 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 651 times.
✓ Branch 3 taken 53037 times.
 53688 if (*last < 2*a - b || *last > 2*b - a)
183 651 goto simple_round;
184
185 53037 this = av_clip64(*last, a, b);
186 53037 *last = this + duration;
187
188 53037 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 592145 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 582141 times.
 592145 if (inc != 1)
196 10004 inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 592145 m = inc_tb.num * (int64_t)ts_tb.den;
199 592145 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 578366 times.
✓ Branch 1 taken 13779 times.
✓ Branch 2 taken 578366 times.
✗ Branch 3 not taken.
 592145 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 578366 return ts + m / d;
203
2/2
✓ Branch 0 taken 2823 times.
✓ Branch 1 taken 10956 times.
 13779 if (m < d)
204 2823 return ts;
205
206 {
207 10956 int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208 10956 int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209
210
3/6
✓ Branch 0 taken 10956 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10956 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 10956 times.
 10956 if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211 ✗ return ts;
212
213 10956 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }
215 }
216
217 85282530 static inline double eval_poly(const double *coeff, int size, double x) {
218 85282530 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 809975451 times.
✓ Branch 1 taken 85282530 times.
 895257981 for (i = size-2; i >= 0; --i) {
221 809975451 sum *= x;
222 809975451 sum += coeff[i];
223 }
224 85282530 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 42643572 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 2307 times.
✓ Branch 1 taken 42641265 times.
 42643572 if (x == 0)
307 2307 return 1.0;
308 42641265 x = fabs(x);
309
2/2
✓ Branch 0 taken 42606501 times.
✓ Branch 1 taken 34764 times.
 42641265 if (x <= 15) {
310 42606501 y = x * x;
311 42606501 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 34764 y = 1 / x - 1.0 / 15;
315 34764 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 34764 factor = exp(x) / sqrt(x);
317 34764 return factor * r;
318 }
319 }
320