FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2024-07-26 21:54:09
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 1457501 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 534504 times.
✓ Branch 1 taken 922997 times.
 1457501 if (a == 0)
41 534504 return b;
42
2/2
✓ Branch 0 taken 3260 times.
✓ Branch 1 taken 919737 times.
 922997 if (b == 0)
43 3260 return a;
44 919737 za = ff_ctzll(a);
45 919737 zb = ff_ctzll(b);
46 919737 k = FFMIN(za, zb);
47 919737 u = llabs(a >> za);
48 919737 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 3880313 times.
✓ Branch 1 taken 919737 times.
 4800050 while (u != v) {
50
2/2
✓ Branch 0 taken 1467598 times.
✓ Branch 1 taken 2412715 times.
 3880313 if (u > v)
51 1467598 FFSWAP(int64_t, v, u);
52 3880313 v -= u;
53 3880313 v >>= ff_ctzll(v);
54 }
55 919737 return (uint64_t)u << k;
56 }
57
58 14179214 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 14179214 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 14179214 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 14179214 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 14179214 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 14179214 times.
 14179214 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 437318 times.
✓ Branch 1 taken 13741896 times.
 14179214 if (rnd & AV_ROUND_PASS_MINMAX) {
69
3/4
✓ Branch 0 taken 435998 times.
✓ Branch 1 taken 1320 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 435998 times.
 437318 if (a == INT64_MIN || a == INT64_MAX)
70 1320 return a;
71 435998 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 105188 times.
✓ Branch 1 taken 14072706 times.
 14177894 if (a < 0)
75 105188 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 12909747 times.
✓ Branch 1 taken 1162959 times.
 14072706 if (rnd == AV_ROUND_NEAR_INF)
78 12909747 r = c / 2;
79
2/2
✓ Branch 0 taken 283590 times.
✓ Branch 1 taken 879369 times.
 1162959 else if (rnd & 1)
80 283590 r = c - 1;
81
82
4/4
✓ Branch 0 taken 13939684 times.
✓ Branch 1 taken 133022 times.
✓ Branch 2 taken 13849724 times.
✓ Branch 3 taken 89960 times.
 14072706 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 13742890 times.
✓ Branch 1 taken 106834 times.
 13849724 if (a <= INT_MAX)
84 13742890 return (a * b + r) / c;
85 else {
86 106834 int64_t ad = a / c;
87 106834 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 1364 times.
✓ Branch 1 taken 105470 times.
✓ Branch 2 taken 1364 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1363 times.
 106834 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 1 return INT64_MIN;
90 106833 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 222982 uint64_t a0 = a & 0xFFFFFFFF;
95 222982 uint64_t a1 = a >> 32;
96 222982 uint64_t b0 = b & 0xFFFFFFFF;
97 222982 uint64_t b1 = b >> 32;
98 222982 uint64_t t1 = a0 * b1 + a1 * b0;
99 222982 uint64_t t1a = t1 << 32;
100 int i;
101
102 222982 a0 = a0 * b0 + t1a;
103 222982 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 222982 a0 += r;
105 222982 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 14270848 times.
✓ Branch 1 taken 222982 times.
 14493830 for (i = 63; i >= 0; i--) {
108 14270848 a1 += a1 + ((a0 >> i) & 1);
109 14270848 t1 += t1;
110
2/2
✓ Branch 0 taken 1712359 times.
✓ Branch 1 taken 12558489 times.
 14270848 if (c <= a1) {
111 1712359 a1 -= c;
112 1712359 t1++;
113 }
114 }
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 222982 times.
 222982 if (t1 > INT64_MAX)
116 ✗ return INT64_MIN;
117 222982 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 1071039 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 1071039 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 12241580 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 12241580 int64_t b = bq.num * (int64_t)cq.den;
138 12241580 int64_t c = cq.num * (int64_t)bq.den;
139 12241580 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 11402785 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 11402785 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 896082 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 896082 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 896082 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 405146 times.
✓ Branch 1 taken 490936 times.
✓ Branch 2 taken 3185 times.
✓ Branch 3 taken 892897 times.
✓ Branch 4 taken 888132 times.
✓ Branch 5 taken 7950 times.
 896082 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 888132 return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
153
2/2
✓ Branch 0 taken 6696 times.
✓ Branch 1 taken 1254 times.
 7950 if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
154 6696 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 303040 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 303040 times.
 303040 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 303040 times.
 303040 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1813 times.
✓ Branch 1 taken 301227 times.
✓ Branch 2 taken 16456 times.
✓ Branch 3 taken 284771 times.
✓ Branch 4 taken 239783 times.
✓ Branch 5 taken 44988 times.
 303040 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 258052 simple_round:
176 258685 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 258685 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 44988 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 44988 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 44988 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 633 times.
✓ Branch 3 taken 44355 times.
 44988 if (*last < 2*a - b || *last > 2*b - a)
183 633 goto simple_round;
184
185 44355 this = av_clip64(*last, a, b);
186 44355 *last = this + duration;
187
188 44355 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 493726 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 493726 times.
 493726 if (inc != 1)
196 ✗ inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 493726 m = inc_tb.num * (int64_t)ts_tb.den;
199 493726 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 480785 times.
✓ Branch 1 taken 12941 times.
✓ Branch 2 taken 480785 times.
✗ Branch 3 not taken.
 493726 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 480785 return ts + m / d;
203
2/2
✓ Branch 0 taken 2822 times.
✓ Branch 1 taken 10119 times.
 12941 if (m < d)
204 2822 return ts;
205
206 {
207 10119 int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208 10119 int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209
210
3/6
✓ Branch 0 taken 10119 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10119 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 10119 times.
 10119 if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211 ✗ return ts;
212
213 10119 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }
215 }
216
217 84983556 static inline double eval_poly(const double *coeff, int size, double x) {
218 84983556 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 807198192 times.
✓ Branch 1 taken 84983556 times.
 892181748 for (i = size-2; i >= 0; --i) {
221 807198192 sum *= x;
222 807198192 sum += coeff[i];
223 }
224 84983556 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 42493569 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 1791 times.
✓ Branch 1 taken 42491778 times.
 42493569 if (x == 0)
307 1791 return 1.0;
308 42491778 x = fabs(x);
309
2/2
✓ Branch 0 taken 42467513 times.
✓ Branch 1 taken 24265 times.
 42491778 if (x <= 15) {
310 42467513 y = x * x;
311 42467513 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 24265 y = 1 / x - 1.0 / 15;
315 24265 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 24265 factor = exp(x) / sqrt(x);
317 24265 return factor * r;
318 }
319 }
320