FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2026-08-28 05:02:41
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 13408440 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 1612135 times.
✓ Branch 1 taken 11796305 times.
 13408440 if (a == 0)
41 1612135 return b;
42
2/2
✓ Branch 0 taken 3823 times.
✓ Branch 1 taken 11792482 times.
 11796305 if (b == 0)
43 3823 return a;
44 11792482 za = ff_ctzll(a);
45 11792482 zb = ff_ctzll(b);
46 11792482 k = FFMIN(za, zb);
47 11792482 u = llabs(a >> za);
48 11792482 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 64039610 times.
✓ Branch 1 taken 11792482 times.
 75832092 while (u != v) {
50
2/2
✓ Branch 0 taken 37802637 times.
✓ Branch 1 taken 26236973 times.
 64039610 if (u > v)
51 37802637 FFSWAP(int64_t, v, u);
52 64039610 v -= u;
53 64039610 v >>= ff_ctzll(v);
54 }
55 11792482 return (uint64_t)u << k;
56 }
57
58 17706687 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 17706687 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 17706687 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 17706687 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 17706687 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 17706687 times.
 17706687 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 736422 times.
✓ Branch 1 taken 16970265 times.
 17706687 if (rnd & AV_ROUND_PASS_MINMAX) {
69
4/4
✓ Branch 0 taken 734867 times.
✓ Branch 1 taken 1555 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 734866 times.
 736422 if (a == INT64_MIN || a == INT64_MAX)
70 1556 return a;
71 734866 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 127901 times.
✓ Branch 1 taken 17577230 times.
 17705131 if (a < 0)
75 127901 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 16250566 times.
✓ Branch 1 taken 1326664 times.
 17577230 if (rnd == AV_ROUND_NEAR_INF)
78 16250566 r = c / 2;
79
2/2
✓ Branch 0 taken 359177 times.
✓ Branch 1 taken 967487 times.
 1326664 else if (rnd & 1)
80 359177 r = c - 1;
81
82
4/4
✓ Branch 0 taken 17421122 times.
✓ Branch 1 taken 156108 times.
✓ Branch 2 taken 17308660 times.
✓ Branch 3 taken 112462 times.
 17577230 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 17222832 times.
✓ Branch 1 taken 85828 times.
 17308660 if (a <= INT_MAX)
84 17222832 return (a * b + r) / c;
85 else {
86 85828 int64_t ad = a / c;
87 85828 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 1253 times.
✓ Branch 1 taken 84575 times.
✓ Branch 2 taken 1253 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1252 times.
 85828 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 1 return INT64_MIN;
90 85827 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 268570 uint64_t a0 = a & 0xFFFFFFFF;
95 268570 uint64_t a1 = a >> 32;
96 268570 uint64_t b0 = b & 0xFFFFFFFF;
97 268570 uint64_t b1 = b >> 32;
98 268570 uint64_t t1 = a0 * b1 + a1 * b0;
99 268570 uint64_t t1a = t1 << 32;
100 int i;
101
102 268570 a0 = a0 * b0 + t1a;
103 268570 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 268570 a0 += r;
105 268570 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 17188480 times.
✓ Branch 1 taken 268570 times.
 17457050 for (i = 63; i >= 0; i--) {
108 17188480 a1 += a1 + ((a0 >> i) & 1);
109 17188480 t1 += t1;
110
2/2
✓ Branch 0 taken 2166875 times.
✓ Branch 1 taken 15021605 times.
 17188480 if (c <= a1) {
111 2166875 a1 -= c;
112 2166875 t1++;
113 }
114 }
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 268570 times.
 268570 if (t1 > INT64_MAX)
116 ✗ return INT64_MIN;
117 268570 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 1258149 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 1258149 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 15501055 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 15501055 int64_t b = bq.num * (int64_t)cq.den;
138 15501055 int64_t c = cq.num * (int64_t)bq.den;
139 15501055 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 14257624 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 14257624 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 1077083 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 1077083 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 1077083 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 480059 times.
✓ Branch 1 taken 597024 times.
✓ Branch 2 taken 5636 times.
✓ Branch 3 taken 1071447 times.
✓ Branch 4 taken 1069472 times.
✓ Branch 5 taken 7611 times.
 1077083 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 1069472 return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
153
2/2
✓ Branch 0 taken 6340 times.
✓ Branch 1 taken 1271 times.
 7611 if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
154 6340 return -1;
155
2/2
✓ Branch 0 taken 202 times.
✓ Branch 1 taken 1069 times.
 1271 if (av_rescale_rnd(ts_b, b, a, AV_ROUND_DOWN) < ts_a)
156 202 return 1;
157 1069 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 360816 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 360816 times.
 360816 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 360816 times.
 360816 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1572 times.
✓ Branch 1 taken 359244 times.
✓ Branch 2 taken 16995 times.
✓ Branch 3 taken 342249 times.
✓ Branch 4 taken 288602 times.
✓ Branch 5 taken 53647 times.
 360816 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 307169 simple_round:
176 307819 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 307819 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 53647 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 53647 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 53647 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 650 times.
✓ Branch 3 taken 52997 times.
 53647 if (*last < 2*a - b || *last > 2*b - a)
183 650 goto simple_round;
184
185 52997 this = av_clip64(*last, a, b);
186 52997 *last = this + duration;
187
188 52997 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 588132 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 578128 times.
 588132 if (inc != 1)
196 10004 inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 588132 m = inc_tb.num * (int64_t)ts_tb.den;
199 588132 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 574288 times.
✓ Branch 1 taken 13844 times.
✓ Branch 2 taken 574288 times.
✗ Branch 3 not taken.
 588132 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 574288 return ts + m / d;
203
2/2
✓ Branch 0 taken 2823 times.
✓ Branch 1 taken 11021 times.
 13844 if (m < d)
204 2823 return ts;
205
206 {
207 11021 int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208 11021 int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209
210
3/6
✓ Branch 0 taken 11021 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 11021 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 11021 times.
 11021 if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211 ✗ return ts;
212
213 11021 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }
215 }
216
217 85190780 static inline double eval_poly(const double *coeff, int size, double x) {
218 85190780 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 809122456 times.
✓ Branch 1 taken 85190780 times.
 894313236 for (i = size-2; i >= 0; --i) {
221 809122456 sum *= x;
222 809122456 sum += coeff[i];
223 }
224 85190780 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 42597532 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 2142 times.
✓ Branch 1 taken 42595390 times.
 42597532 if (x == 0)
307 2142 return 1.0;
308 42595390 x = fabs(x);
309
2/2
✓ Branch 0 taken 42563731 times.
✓ Branch 1 taken 31659 times.
 42595390 if (x <= 15) {
310 42563731 y = x * x;
311 42563731 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 31659 y = 1 / x - 1.0 / 15;
315 31659 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 31659 factor = exp(x) / sqrt(x);
317 31659 return factor * r;
318 }
319 }
320