FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2025-02-18 20:34:46
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 1849854 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 630027 times.
✓ Branch 1 taken 1219827 times.
 1849854 if (a == 0)
41 630027 return b;
42
2/2
✓ Branch 0 taken 3540 times.
✓ Branch 1 taken 1216287 times.
 1219827 if (b == 0)
43 3540 return a;
44 1216287 za = ff_ctzll(a);
45 1216287 zb = ff_ctzll(b);
46 1216287 k = FFMIN(za, zb);
47 1216287 u = llabs(a >> za);
48 1216287 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 4826525 times.
✓ Branch 1 taken 1216287 times.
 6042812 while (u != v) {
50
2/2
✓ Branch 0 taken 2118182 times.
✓ Branch 1 taken 2708343 times.
 4826525 if (u > v)
51 2118182 FFSWAP(int64_t, v, u);
52 4826525 v -= u;
53 4826525 v >>= ff_ctzll(v);
54 }
55 1216287 return (uint64_t)u << k;
56 }
57
58 14814953 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 14814953 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 14814953 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 14814953 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 14814953 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 14814953 times.
 14814953 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 450710 times.
✓ Branch 1 taken 14364243 times.
 14814953 if (rnd & AV_ROUND_PASS_MINMAX) {
69
3/4
✓ Branch 0 taken 449320 times.
✓ Branch 1 taken 1390 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 449320 times.
 450710 if (a == INT64_MIN || a == INT64_MAX)
70 1390 return a;
71 449320 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 107983 times.
✓ Branch 1 taken 14705580 times.
 14813563 if (a < 0)
75 107983 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 13529449 times.
✓ Branch 1 taken 1176131 times.
 14705580 if (rnd == AV_ROUND_NEAR_INF)
78 13529449 r = c / 2;
79
2/2
✓ Branch 0 taken 283713 times.
✓ Branch 1 taken 892418 times.
 1176131 else if (rnd & 1)
80 283713 r = c - 1;
81
82
4/4
✓ Branch 0 taken 14572441 times.
✓ Branch 1 taken 133139 times.
✓ Branch 2 taken 14461037 times.
✓ Branch 3 taken 111404 times.
 14705580 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 14375780 times.
✓ Branch 1 taken 85257 times.
 14461037 if (a <= INT_MAX)
84 14375780 return (a * b + r) / c;
85 else {
86 85257 int64_t ad = a / c;
87 85257 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 1247 times.
✓ Branch 1 taken 84010 times.
✓ Branch 2 taken 1247 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1246 times.
 85257 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 1 return INT64_MIN;
90 85256 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 244543 uint64_t a0 = a & 0xFFFFFFFF;
95 244543 uint64_t a1 = a >> 32;
96 244543 uint64_t b0 = b & 0xFFFFFFFF;
97 244543 uint64_t b1 = b >> 32;
98 244543 uint64_t t1 = a0 * b1 + a1 * b0;
99 244543 uint64_t t1a = t1 << 32;
100 int i;
101
102 244543 a0 = a0 * b0 + t1a;
103 244543 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 244543 a0 += r;
105 244543 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 15650752 times.
✓ Branch 1 taken 244543 times.
 15895295 for (i = 63; i >= 0; i--) {
108 15650752 a1 += a1 + ((a0 >> i) & 1);
109 15650752 t1 += t1;
110
2/2
✓ Branch 0 taken 1888022 times.
✓ Branch 1 taken 13762730 times.
 15650752 if (c <= a1) {
111 1888022 a1 -= c;
112 1888022 t1++;
113 }
114 }
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 244543 times.
 244543 if (t1 > INT64_MAX)
116 ✗ return INT64_MIN;
117 244543 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 1080292 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 1080292 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 12853606 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 12853606 int64_t b = bq.num * (int64_t)cq.den;
138 12853606 int64_t c = cq.num * (int64_t)bq.den;
139 12853606 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 11999912 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 11999912 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 951637 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 951637 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 951637 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 420703 times.
✓ Branch 1 taken 530934 times.
✓ Branch 2 taken 5476 times.
✓ Branch 3 taken 946161 times.
✓ Branch 4 taken 943684 times.
✓ Branch 5 taken 7953 times.
 951637 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 943684 return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
153
2/2
✓ Branch 0 taken 6699 times.
✓ Branch 1 taken 1254 times.
 7953 if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
154 6699 return -1;
155
2/2
✓ Branch 0 taken 193 times.
✓ Branch 1 taken 1061 times.
 1254 if (av_rescale_rnd(ts_b, b, a, AV_ROUND_DOWN) < ts_a)
156 193 return 1;
157 1061 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 303867 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 303867 times.
 303867 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 303867 times.
 303867 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1815 times.
✓ Branch 1 taken 302052 times.
✓ Branch 2 taken 16640 times.
✓ Branch 3 taken 285412 times.
✓ Branch 4 taken 240423 times.
✓ Branch 5 taken 44989 times.
 303867 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 258878 simple_round:
176 259509 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 259509 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 44989 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 44989 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 44989 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 631 times.
✓ Branch 3 taken 44358 times.
 44989 if (*last < 2*a - b || *last > 2*b - a)
183 631 goto simple_round;
184
185 44358 this = av_clip64(*last, a, b);
186 44358 *last = this + duration;
187
188 44358 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 507810 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 507810 times.
 507810 if (inc != 1)
196 ✗ inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 507810 m = inc_tb.num * (int64_t)ts_tb.den;
199 507810 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 494499 times.
✓ Branch 1 taken 13311 times.
✓ Branch 2 taken 494499 times.
✗ Branch 3 not taken.
 507810 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 494499 return ts + m / d;
203
2/2
✓ Branch 0 taken 2822 times.
✓ Branch 1 taken 10489 times.
 13311 if (m < d)
204 2822 return ts;
205
206 {
207 10489 int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208 10489 int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209
210
3/6
✓ Branch 0 taken 10489 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10489 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 10489 times.
 10489 if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211 ✗ return ts;
212
213 10489 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }
215 }
216
217 85003292 static inline double eval_poly(const double *coeff, int size, double x) {
218 85003292 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 807379726 times.
✓ Branch 1 taken 85003292 times.
 892383018 for (i = size-2; i >= 0; --i) {
221 807379726 sum *= x;
222 807379726 sum += coeff[i];
223 }
224 85003292 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 42503478 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 1832 times.
✓ Branch 1 taken 42501646 times.
 42503478 if (x == 0)
307 1832 return 1.0;
308 42501646 x = fabs(x);
309
2/2
✓ Branch 0 taken 42476388 times.
✓ Branch 1 taken 25258 times.
 42501646 if (x <= 15) {
310 42476388 y = x * x;
311 42476388 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 25258 y = 1 / x - 1.0 / 15;
315 25258 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 25258 factor = exp(x) / sqrt(x);
317 25258 return factor * r;
318 }
319 }
320