FFmpeg coverage


Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2024-02-16 17:37:06
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 1496342 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 580703 times.
✓ Branch 1 taken 915639 times.
1496342 if (a == 0)
41 580703 return b;
42
2/2
✓ Branch 0 taken 3239 times.
✓ Branch 1 taken 912400 times.
915639 if (b == 0)
43 3239 return a;
44 912400 za = ff_ctzll(a);
45 912400 zb = ff_ctzll(b);
46 912400 k = FFMIN(za, zb);
47 912400 u = llabs(a >> za);
48 912400 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 3851738 times.
✓ Branch 1 taken 912400 times.
4764138 while (u != v) {
50
2/2
✓ Branch 0 taken 1463071 times.
✓ Branch 1 taken 2388667 times.
3851738 if (u > v)
51 1463071 FFSWAP(int64_t, v, u);
52 3851738 v -= u;
53 3851738 v >>= ff_ctzll(v);
54 }
55 912400 return (uint64_t)u << k;
56 }
57
58 15379383 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 15379383 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 15379383 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 15379383 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 15379383 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 15379383 times.
15379383 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 484285 times.
✓ Branch 1 taken 14895098 times.
15379383 if (rnd & AV_ROUND_PASS_MINMAX) {
69
3/4
✓ Branch 0 taken 482992 times.
✓ Branch 1 taken 1293 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 482992 times.
484285 if (a == INT64_MIN || a == INT64_MAX)
70 1293 return a;
71 482992 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 105275 times.
✓ Branch 1 taken 15272815 times.
15378090 if (a < 0)
75 105275 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 14016279 times.
✓ Branch 1 taken 1256536 times.
15272815 if (rnd == AV_ROUND_NEAR_INF)
78 14016279 r = c / 2;
79
2/2
✓ Branch 0 taken 329664 times.
✓ Branch 1 taken 926872 times.
1256536 else if (rnd & 1)
80 329664 r = c - 1;
81
82
4/4
✓ Branch 0 taken 15133613 times.
✓ Branch 1 taken 139202 times.
✓ Branch 2 taken 15109392 times.
✓ Branch 3 taken 24221 times.
15272815 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 15097518 times.
✓ Branch 1 taken 11874 times.
15109392 if (a <= INT_MAX)
84 15097518 return (a * b + r) / c;
85 else {
86 11874 int64_t ad = a / c;
87 11874 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 8439 times.
✓ Branch 1 taken 3435 times.
✓ Branch 2 taken 8439 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 8438 times.
11874 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 1 return INT64_MIN;
90 11873 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 163423 uint64_t a0 = a & 0xFFFFFFFF;
95 163423 uint64_t a1 = a >> 32;
96 163423 uint64_t b0 = b & 0xFFFFFFFF;
97 163423 uint64_t b1 = b >> 32;
98 163423 uint64_t t1 = a0 * b1 + a1 * b0;
99 163423 uint64_t t1a = t1 << 32;
100 int i;
101
102 163423 a0 = a0 * b0 + t1a;
103 163423 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 163423 a0 += r;
105 163423 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 10459072 times.
✓ Branch 1 taken 163423 times.
10622495 for (i = 63; i >= 0; i--) {
108 10459072 a1 += a1 + ((a0 >> i) & 1);
109 10459072 t1 += t1;
110
2/2
✓ Branch 0 taken 1307598 times.
✓ Branch 1 taken 9151474 times.
10459072 if (c <= a1) {
111 1307598 a1 -= c;
112 1307598 t1++;
113 }
114 }
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 163423 times.
163423 if (t1 > INT64_MAX)
116 return INT64_MIN;
117 163423 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 932673 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 932673 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 13533332 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 13533332 int64_t b = bq.num * (int64_t)cq.den;
138 13533332 int64_t c = cq.num * (int64_t)bq.den;
139 13533332 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 12600689 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 12600689 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 1018797 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 1018797 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 1018797 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 452143 times.
✓ Branch 1 taken 566654 times.
✓ Branch 2 taken 2857 times.
✓ Branch 3 taken 1015940 times.
✓ Branch 4 taken 1010844 times.
✓ Branch 5 taken 7953 times.
1018797 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 1010844 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 354275 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 354275 times.
354275 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 354275 times.
354275 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1806 times.
✓ Branch 1 taken 352469 times.
✓ Branch 2 taken 16456 times.
✓ Branch 3 taken 336013 times.
✓ Branch 4 taken 290992 times.
✓ Branch 5 taken 45021 times.
354275 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 309254 simple_round:
176 309887 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 309887 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 45021 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 45021 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 45021 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 633 times.
✓ Branch 3 taken 44388 times.
45021 if (*last < 2*a - b || *last > 2*b - a)
183 633 goto simple_round;
184
185 44388 this = av_clip64(*last, a, b);
186 44388 *last = this + duration;
187
188 44388 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 538751 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 538751 times.
538751 if (inc != 1)
196 inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 538751 m = inc_tb.num * (int64_t)ts_tb.den;
199 538751 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 523489 times.
✓ Branch 1 taken 15262 times.
✓ Branch 2 taken 523489 times.
✗ Branch 3 not taken.
538751 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 523489 return ts + m / d;
203
2/2
✓ Branch 0 taken 2822 times.
✓ Branch 1 taken 12440 times.
15262 if (m < d)
204 2822 return ts;
205
206 {
207 12440 int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208 12440 int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209
210
3/6
✓ Branch 0 taken 12440 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 12440 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 12440 times.
12440 if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211 return ts;
212
213 12440 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }
215 }
216
217 84791316 static inline double eval_poly(const double *coeff, int size, double x) {
218 84791316 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 805410966 times.
✓ Branch 1 taken 84791316 times.
890202282 for (i = size-2; i >= 0; --i) {
221 805410966 sum *= x;
222 805410966 sum += coeff[i];
223 }
224 84791316 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 42397115 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 1457 times.
✓ Branch 1 taken 42395658 times.
42397115 if (x == 0)
307 1457 return 1.0;
308 42395658 x = fabs(x);
309
2/2
✓ Branch 0 taken 42377902 times.
✓ Branch 1 taken 17756 times.
42395658 if (x <= 15) {
310 42377902 y = x * x;
311 42377902 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 17756 y = 1 / x - 1.0 / 15;
315 17756 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 17756 factor = exp(x) / sqrt(x);
317 17756 return factor * r;
318 }
319 }
320