FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libavutil/mathematics.c
Date: 2024-11-20 23:03:26
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 1544522 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 562491 times.
✓ Branch 1 taken 982031 times.
 1544522 if (a == 0)
41 562491 return b;
42
2/2
✓ Branch 0 taken 3309 times.
✓ Branch 1 taken 978722 times.
 982031 if (b == 0)
43 3309 return a;
44 978722 za = ff_ctzll(a);
45 978722 zb = ff_ctzll(b);
46 978722 k = FFMIN(za, zb);
47 978722 u = llabs(a >> za);
48 978722 v = llabs(b >> zb);
49
2/2
✓ Branch 0 taken 3898710 times.
✓ Branch 1 taken 978722 times.
 4877432 while (u != v) {
50
2/2
✓ Branch 0 taken 1472827 times.
✓ Branch 1 taken 2425883 times.
 3898710 if (u > v)
51 1472827 FFSWAP(int64_t, v, u);
52 3898710 v -= u;
53 3898710 v >>= ff_ctzll(v);
54 }
55 978722 return (uint64_t)u << k;
56 }
57
58 14769095 int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60 14769095 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 14769095 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 14769095 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 14769095 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 14769095 times.
 14769095 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 448316 times.
✓ Branch 1 taken 14320779 times.
 14769095 if (rnd & AV_ROUND_PASS_MINMAX) {
69
3/4
✓ Branch 0 taken 446988 times.
✓ Branch 1 taken 1328 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 446988 times.
 448316 if (a == INT64_MIN || a == INT64_MAX)
70 1328 return a;
71 446988 rnd -= AV_ROUND_PASS_MINMAX;
72 }
73
74
2/2
✓ Branch 0 taken 106090 times.
✓ Branch 1 taken 14661677 times.
 14767767 if (a < 0)
75 106090 return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76
77
2/2
✓ Branch 0 taken 13484621 times.
✓ Branch 1 taken 1177056 times.
 14661677 if (rnd == AV_ROUND_NEAR_INF)
78 13484621 r = c / 2;
79
2/2
✓ Branch 0 taken 283649 times.
✓ Branch 1 taken 893407 times.
 1177056 else if (rnd & 1)
80 283649 r = c - 1;
81
82
4/4
✓ Branch 0 taken 14528485 times.
✓ Branch 1 taken 133192 times.
✓ Branch 2 taken 14417069 times.
✓ Branch 3 taken 111416 times.
 14661677 if (b <= INT_MAX && c <= INT_MAX) {
83
2/2
✓ Branch 0 taken 14331681 times.
✓ Branch 1 taken 85388 times.
 14417069 if (a <= INT_MAX)
84 14331681 return (a * b + r) / c;
85 else {
86 85388 int64_t ad = a / c;
87 85388 int64_t a2 = (a % c * b + r) / c;
88
5/6
✓ Branch 0 taken 1364 times.
✓ Branch 1 taken 84024 times.
✓ Branch 2 taken 1364 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1363 times.
 85388 if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89 1 return INT64_MIN;
90 85387 return ad * b + a2;
91 }
92 } else {
93 #if 1
94 244608 uint64_t a0 = a & 0xFFFFFFFF;
95 244608 uint64_t a1 = a >> 32;
96 244608 uint64_t b0 = b & 0xFFFFFFFF;
97 244608 uint64_t b1 = b >> 32;
98 244608 uint64_t t1 = a0 * b1 + a1 * b0;
99 244608 uint64_t t1a = t1 << 32;
100 int i;
101
102 244608 a0 = a0 * b0 + t1a;
103 244608 a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104 244608 a0 += r;
105 244608 a1 += a0 < r;
106
107
2/2
✓ Branch 0 taken 15654912 times.
✓ Branch 1 taken 244608 times.
 15899520 for (i = 63; i >= 0; i--) {
108 15654912 a1 += a1 + ((a0 >> i) & 1);
109 15654912 t1 += t1;
110
2/2
✓ Branch 0 taken 1888506 times.
✓ Branch 1 taken 13766406 times.
 15654912 if (c <= a1) {
111 1888506 a1 -= c;
112 1888506 t1++;
113 }
114 }
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 244608 times.
 244608 if (t1 > INT64_MAX)
116 ✗ return INT64_MIN;
117 244608 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 1078555 int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131 1078555 return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133
134 12808997 int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135 enum AVRounding rnd)
136 {
137 12808997 int64_t b = bq.num * (int64_t)cq.den;
138 12808997 int64_t c = cq.num * (int64_t)bq.den;
139 12808997 return av_rescale_rnd(a, b, c, rnd);
140 }
141
142 11959153 int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144 11959153 return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146
147 944122 int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149 944122 int64_t a = tb_a.num * (int64_t)tb_b.den;
150 944122 int64_t b = tb_b.num * (int64_t)tb_a.den;
151
6/6
✓ Branch 0 taken 418298 times.
✓ Branch 1 taken 525824 times.
✓ Branch 2 taken 5423 times.
✓ Branch 3 taken 938699 times.
✓ Branch 4 taken 936168 times.
✓ Branch 5 taken 7954 times.
 944122 if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152 936168 return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
153
2/2
✓ Branch 0 taken 6700 times.
✓ Branch 1 taken 1254 times.
 7954 if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
154 6700 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 303140 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 303140 times.
 303140 av_assert0(in_ts != AV_NOPTS_VALUE);
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 303140 times.
 303140 av_assert0(duration >= 0);
173
174
6/6
✓ Branch 0 taken 1817 times.
✓ Branch 1 taken 301323 times.
✓ Branch 2 taken 16456 times.
✓ Branch 3 taken 284867 times.
✓ Branch 4 taken 239872 times.
✓ Branch 5 taken 44995 times.
 303140 if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 258145 simple_round:
176 258778 *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177 258778 return av_rescale_q(in_ts, in_tb, out_tb);
178 }
179
180 44995 a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181 44995 b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182
3/4
✓ Branch 0 taken 44995 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 633 times.
✓ Branch 3 taken 44362 times.
 44995 if (*last < 2*a - b || *last > 2*b - a)
183 633 goto simple_round;
184
185 44362 this = av_clip64(*last, a, b);
186 44362 *last = this + duration;
187
188 44362 return av_rescale_q(this, fs_tb, out_tb);
189 }
190
191 505489 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 505489 times.
 505489 if (inc != 1)
196 ✗ inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197
198 505489 m = inc_tb.num * (int64_t)ts_tb.den;
199 505489 d = inc_tb.den * (int64_t)ts_tb.num;
200
201
3/4
✓ Branch 0 taken 492526 times.
✓ Branch 1 taken 12963 times.
✓ Branch 2 taken 492526 times.
✗ Branch 3 not taken.
 505489 if (m % d == 0 && ts <= INT64_MAX - m / d)
202 492526 return ts + m / d;
203
2/2
✓ Branch 0 taken 2822 times.
✓ Branch 1 taken 10141 times.
 12963 if (m < d)
204 2822 return ts;
205
206 {
207 10141 int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208 10141 int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209
210
3/6
✓ Branch 0 taken 10141 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10141 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 10141 times.
 10141 if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211 ✗ return ts;
212
213 10141 return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214 }
215 }
216
217 84991860 static inline double eval_poly(const double *coeff, int size, double x) {
218 84991860 double sum = coeff[size-1];
219 int i;
220
2/2
✓ Branch 0 taken 807274332 times.
✓ Branch 1 taken 84991860 times.
 892266192 for (i = size-2; i >= 0; --i) {
221 807274332 sum *= x;
222 807274332 sum += coeff[i];
223 }
224 84991860 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 42497739 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 1809 times.
✓ Branch 1 taken 42495930 times.
 42497739 if (x == 0)
307 1809 return 1.0;
308 42495930 x = fabs(x);
309
2/2
✓ Branch 0 taken 42471207 times.
✓ Branch 1 taken 24723 times.
 42495930 if (x <= 15) {
310 42471207 y = x * x;
311 42471207 return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312 }
313 else {
314 24723 y = 1 / x - 1.0 / 15;
315 24723 r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316 24723 factor = exp(x) / sqrt(x);
317 24723 return factor * r;
318 }
319 }
320