FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2026-09-14 03:06:03
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 13414267 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 1615108 times.
✓ Branch 1 taken 11799159 times.
 13414267 if (a == 0)
41 1615108 return b;
42
2/2
✓ Branch 0 taken 3823 times.
✓ Branch 1 taken 11795336 times.
 11799159 if (b == 0)
43 3823 return a;
44 11795336 za = ff_ctzll(a);
45 11795336 zb = ff_ctzll(b);
46 11795336 k = FFMIN(za, zb);
47 11795336 u = llabs(a >> za);
48 11795336 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 64044393 times.
✓ Branch 1 taken 11795336 times.
 75839729 while (u != v) {
50
2/2
✓ Branch 0 taken 37803488 times.
✓ Branch 1 taken 26240905 times.
 64044393 if (u > v)
51 37803488 FFSWAP(int64_t, v, u);
52 64044393 v -= u;
53 64044393 v >>= ff_ctzll(v);
54 }
55 11795336 return (uint64_t)u << k;
56 }
57
58 17731171 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 17731171 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 17731171 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 17731171 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 17731171 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 17731171 times.
 17731171 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 737682 times.
✓ Branch 1 taken 16993489 times.
 17731171 if (rnd & AV_ROUND_PASS_MINMAX) {
69
4/4
✓ Branch 0 taken 736107 times.
✓ Branch 1 taken 1575 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 736106 times.
 737682 if (a == INT64_MIN || a == INT64_MAX)
70 1576 return a;
71 736106 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 128087 times.
✓ Branch 1 taken 17601508 times.
 17729595 if (a < 0)
75 128087 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 16273594 times.
✓ Branch 1 taken 1327914 times.
 17601508 if (rnd == AV_ROUND_NEAR_INF)
78 16273594 r = c / 2;
79
2/2
✓ Branch 0 taken 359387 times.
✓ Branch 1 taken 968527 times.
 1327914 else if (rnd & 1)
80 359387 r = c - 1;
81
82
4/4
✓ Branch 0 taken 17444732 times.
✓ Branch 1 taken 156776 times.
✓ Branch 2 taken 17332276 times.
✓ Branch 3 taken 112456 times.
 17601508 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 17246431 times.
✓ Branch 1 taken 85845 times.
 17332276 if (a <= INT_MAX)
84 17246431 return (a * b + r) / c;
85 else {
86 85845 int64_t ad = a / c;
87 85845 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 1260 times.
✓ Branch 1 taken 84585 times.
✓ Branch 2 taken 1260 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1259 times.
 85845 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 1 return INT64_MIN;
90 85844 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 269232 uint64_t a0 = a & 0xFFFFFFFF;
95 269232 uint64_t a1 = a >> 32;
96 269232 uint64_t b0 = b & 0xFFFFFFFF;
97 269232 uint64_t b1 = b >> 32;
98 269232 uint64_t t1 = a0 * b1 + a1 * b0;
99 269232 uint64_t t1a = t1 << 32;
100 int i;
101
102 269232 a0 = a0 * b0 + t1a;
103 269232 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 269232 a0 += r;
105 269232 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 17230848 times.
✓ Branch 1 taken 269232 times.
 17500080 for (i = 63; i >= 0; i--) {
108 17230848 a1 += a1 + ((a0 >> i) & 1);
109 17230848 t1 += t1;
110
2/2
✓ Branch 0 taken 2171666 times.
✓ Branch 1 taken 15059182 times.
 17230848 if (c <= a1) {
111 2171666 a1 -= c;
112 2171666 t1++;
113 }
114 }
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 269232 times.
 269232 if (t1 > INT64_MAX)
116 ✗ return INT64_MIN;
117 269232 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 1258598 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 1258598 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 15523830 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 15523830 int64_t b = bq.num * (int64_t)cq.den;
138 15523830 int64_t c = cq.num * (int64_t)bq.den;
139 15523830 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 14278963 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 14278963 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 1079301 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 1079301 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 1079301 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 481315 times.
✓ Branch 1 taken 597986 times.
✓ Branch 2 taken 5818 times.
✓ Branch 3 taken 1073483 times.
✓ Branch 4 taken 1071683 times.
✓ Branch 5 taken 7618 times.
 1079301 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 1071683 return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
153
2/2
✓ Branch 0 taken 6344 times.
✓ Branch 1 taken 1274 times.
 7618 if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
154 6344 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 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 361138 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 361138 times.
 361138 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 361138 times.
 361138 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1577 times.
✓ Branch 1 taken 359561 times.
✓ Branch 2 taken 16998 times.
✓ Branch 3 taken 342563 times.
✓ Branch 4 taken 288917 times.
✓ Branch 5 taken 53646 times.
 361138 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 307492 simple_round:
176 308142 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 308142 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 53646 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 53646 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 53646 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 650 times.
✓ Branch 3 taken 52996 times.
 53646 if (*last < 2*a - b || *last > 2*b - a)
183 650 goto simple_round;
184
185 52996 this = av_clip64(*last, a, b);
186 52996 *last = this + duration;
187
188 52996 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 589423 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 579419 times.
 589423 if (inc != 1)
196 10004 inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 589423 m = inc_tb.num * (int64_t)ts_tb.den;
199 589423 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 575576 times.
✓ Branch 1 taken 13847 times.
✓ Branch 2 taken 575576 times.
✗ Branch 3 not taken.
 589423 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 575576 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 85271996 static inline double eval_poly(const double *coeff, int size, double x) {
218 85271996 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 809877448 times.
✓ Branch 1 taken 85271996 times.
 895149444 for (i = size-2; i >= 0; --i) {
221 809877448 sum *= x;
222 809877448 sum += coeff[i];
223 }
224 85271996 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 42638284 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 2286 times.
✓ Branch 1 taken 42635998 times.
 42638284 if (x == 0)
307 2286 return 1.0;
308 42635998 x = fabs(x);
309
2/2
✓ Branch 0 taken 42601579 times.
✓ Branch 1 taken 34419 times.
 42635998 if (x <= 15) {
310 42601579 y = x * x;
311 42601579 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 34419 y = 1 / x - 1.0 / 15;
315 34419 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 34419 factor = exp(x) / sqrt(x);
317 34419 return factor * r;
318 }
319 }
320