FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2025-10-10 03:51:19
Exec Total Coverage
Lines: 116 125 92.8%
Functions: 10 11 90.9%
Branches: 75 92 81.5%
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 2069553 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 764312 times.
✓ Branch 1 taken 1305241 times.
 2069553 if (a == 0)
41 764312 return b;
42
2/2
✓ Branch 0 taken 3542 times.
✓ Branch 1 taken 1301699 times.
 1305241 if (b == 0)
43 3542 return a;
44 1301699 za = ff_ctzll(a);
45 1301699 zb = ff_ctzll(b);
46 1301699 k = FFMIN(za, zb);
47 1301699 u = llabs(a >> za);
48 1301699 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 5156496 times.
✓ Branch 1 taken 1301699 times.
 6458195 while (u != v) {
50
2/2
✓ Branch 0 taken 2230077 times.
✓ Branch 1 taken 2926419 times.
 5156496 if (u > v)
51 2230077 FFSWAP(int64_t, v, u);
52 5156496 v -= u;
53 5156496 v >>= ff_ctzll(v);
54 }
55 1301699 return (uint64_t)u << k;
56 }
57
58 15196555 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 15196555 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 15196555 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 15196555 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 15196555 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 15196555 times.
 15196555 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 666477 times.
✓ Branch 1 taken 14530078 times.
 15196555 if (rnd & AV_ROUND_PASS_MINMAX) {
69
3/4
✓ Branch 0 taken 665049 times.
✓ Branch 1 taken 1428 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 665049 times.
 666477 if (a == INT64_MIN || a == INT64_MAX)
70 1428 return a;
71 665049 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 108538 times.
✓ Branch 1 taken 15086589 times.
 15195127 if (a < 0)
75 108538 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 13871614 times.
✓ Branch 1 taken 1214975 times.
 15086589 if (rnd == AV_ROUND_NEAR_INF)
78 13871614 r = c / 2;
79
2/2
✓ Branch 0 taken 298968 times.
✓ Branch 1 taken 916007 times.
 1214975 else if (rnd & 1)
80 298968 r = c - 1;
81
82
4/4
✓ Branch 0 taken 14951070 times.
✓ Branch 1 taken 135519 times.
✓ Branch 2 taken 14839683 times.
✓ Branch 3 taken 111387 times.
 15086589 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 14754044 times.
✓ Branch 1 taken 85639 times.
 14839683 if (a <= INT_MAX)
84 14754044 return (a * b + r) / c;
85 else {
86 85639 int64_t ad = a / c;
87 85639 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 1251 times.
✓ Branch 1 taken 84388 times.
✓ Branch 2 taken 1251 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1250 times.
 85639 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 1 return INT64_MIN;
90 85638 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 246906 uint64_t a0 = a & 0xFFFFFFFF;
95 246906 uint64_t a1 = a >> 32;
96 246906 uint64_t b0 = b & 0xFFFFFFFF;
97 246906 uint64_t b1 = b >> 32;
98 246906 uint64_t t1 = a0 * b1 + a1 * b0;
99 246906 uint64_t t1a = t1 << 32;
100 int i;
101
102 246906 a0 = a0 * b0 + t1a;
103 246906 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 246906 a0 += r;
105 246906 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 15801984 times.
✓ Branch 1 taken 246906 times.
 16048890 for (i = 63; i >= 0; i--) {
108 15801984 a1 += a1 + ((a0 >> i) & 1);
109 15801984 t1 += t1;
110
2/2
✓ Branch 0 taken 1914514 times.
✓ Branch 1 taken 13887470 times.
 15801984 if (c <= a1) {
111 1914514 a1 -= c;
112 1914514 t1++;
113 }
114 }
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 246906 times.
 246906 if (t1 > INT64_MAX)
116 ✗ return INT64_MIN;
117 246906 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 1093612 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 1093612 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 13212249 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 13212249 int64_t b = bq.num * (int64_t)cq.den;
138 13212249 int64_t c = cq.num * (int64_t)bq.den;
139 13212249 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 12113028 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 12113028 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 953834 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 953834 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 953834 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 416149 times.
✓ Branch 1 taken 537685 times.
✓ Branch 2 taken 5062 times.
✓ Branch 3 taken 948772 times.
✓ Branch 4 taken 945921 times.
✓ Branch 5 taken 7913 times.
 953834 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 945921 return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
153
2/2
✓ Branch 0 taken 6656 times.
✓ Branch 1 taken 1257 times.
 7913 if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
154 6656 return -1;
155
2/2
✓ Branch 0 taken 193 times.
✓ Branch 1 taken 1064 times.
 1257 if (av_rescale_rnd(ts_b, b, a, AV_ROUND_DOWN) < ts_a)
156 193 return 1;
157 1064 return 0;
158 }
159
160 ✗ int64_t av_compare_mod(uint64_t a, uint64_t b, uint64_t mod)
161 {
162 ✗ int64_t c = (a - b) & (mod - 1);
163 ✗ if (c > (mod >> 1))
164 ✗ c -= mod;
165 ✗ return c;
166 }
167
168 309756 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 309756 times.
 309756 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 309756 times.
 309756 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1821 times.
✓ Branch 1 taken 307935 times.
✓ Branch 2 taken 16640 times.
✓ Branch 3 taken 291295 times.
✓ Branch 4 taken 238112 times.
✓ Branch 5 taken 53183 times.
 309756 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 256573 simple_round:
176 257204 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 257204 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 53183 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 53183 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 53183 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 631 times.
✓ Branch 3 taken 52552 times.
 53183 if (*last < 2*a - b || *last > 2*b - a)
183 631 goto simple_round;
184
185 52552 this = av_clip64(*last, a, b);
186 52552 *last = this + duration;
187
188 52552 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 521608 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 521608 times.
 521608 if (inc != 1)
196 ✗ inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 521608 m = inc_tb.num * (int64_t)ts_tb.den;
199 521608 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 508298 times.
✓ Branch 1 taken 13310 times.
✓ Branch 2 taken 508298 times.
✗ Branch 3 not taken.
 521608 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 508298 return ts + m / d;
203
2/2
✓ Branch 0 taken 2822 times.
✓ Branch 1 taken 10488 times.
 13310 if (m < d)
204 2822 return ts;
205
206 {
207 10488 int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208 10488 int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209
210
3/6
✓ Branch 0 taken 10488 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10488 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 10488 times.
 10488 if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211 ✗ return ts;
212
213 10488 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }
215 }
216
217 85006422 static inline double eval_poly(const double *coeff, int size, double x) {
218 85006422 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 807408321 times.
✓ Branch 1 taken 85006422 times.
 892414743 for (i = size-2; i >= 0; --i) {
221 807408321 sum *= x;
222 807408321 sum += coeff[i];
223 }
224 85006422 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 42505051 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 1840 times.
✓ Branch 1 taken 42503211 times.
 42505051 if (x == 0)
307 1840 return 1.0;
308 42503211 x = fabs(x);
309
2/2
✓ Branch 0 taken 42477763 times.
✓ Branch 1 taken 25448 times.
 42503211 if (x <= 15) {
310 42477763 y = x * x;
311 42477763 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 25448 y = 1 / x - 1.0 / 15;
315 25448 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 25448 factor = exp(x) / sqrt(x);
317 25448 return factor * r;
318 }
319 }
320