FFmpeg: libavutil/mathematics.c Source File

FFmpeg
mathematics.c
Go to the documentation of this file.
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  int64_t av_gcd(int64_t a, int64_t b) {
38  int za, zb, k;
39  int64_t u, v;
40  if (a == 0)
41  return b;
42  if (b == 0)
43  return a;
44  za = ff_ctzll(a);
45  zb = ff_ctzll(b);
46  k = FFMIN(za, zb);
47  u = llabs(a >> za);
48  v = llabs(b >> zb);
49  while (u != v) {
50  if (u > v)
51  FFSWAP(int64_t, v, u);
52  v -= u;
53  v >>= ff_ctzll(v);
54  }
55  return (uint64_t)u << k;
56 }
57 
58  int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
59 {
60  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  if (c <= 0 || b < 0 || !((unsigned)(rnd&~AV_ROUND_PASS_MINMAX)<=5 && (rnd&~AV_ROUND_PASS_MINMAX)!=4))
66  return INT64_MIN;
67 
68  if (rnd & AV_ROUND_PASS_MINMAX) {
69  if (a == INT64_MIN || a == INT64_MAX)
70  return a;
71  rnd -= AV_ROUND_PASS_MINMAX;
72  }
73 
74  if (a < 0)
75  return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
76 
77  if (rnd == AV_ROUND_NEAR_INF)
78  r = c / 2;
79  else if (rnd & 1)
80  r = c - 1;
81 
82  if (b <= INT_MAX && c <= INT_MAX) {
83  if (a <= INT_MAX)
84  return (a * b + r) / c;
85  else {
86  int64_t ad = a / c;
87  int64_t a2 = (a % c * b + r) / c;
88  if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
89  return INT64_MIN;
90  return ad * b + a2;
91  }
92  } else {
93 #if 1
94  uint64_t a0 = a & 0xFFFFFFFF;
95  uint64_t a1 = a >> 32;
96  uint64_t b0 = b & 0xFFFFFFFF;
97  uint64_t b1 = b >> 32;
98  uint64_t t1 = a0 * b1 + a1 * b0;
99  uint64_t t1a = t1 << 32;
100  int i;
101 
102  a0 = a0 * b0 + t1a;
103  a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
104  a0 += r;
105  a1 += a0 < r;
106 
107  for (i = 63; i >= 0; i--) {
108  a1 += a1 + ((a0 >> i) & 1);
109  t1 += t1;
110  if (c <= a1) {
111  a1 -= c;
112  t1++;
113  }
114  }
115  if (t1 > INT64_MAX)
116  return INT64_MIN;
117  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  int64_t av_rescale(int64_t a, int64_t b, int64_t c)
130 {
131  return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
132 }
133 
134  int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
135  enum AVRounding rnd)
136 {
137  int64_t b = bq.num * (int64_t)cq.den;
138  int64_t c = cq.num * (int64_t)bq.den;
139  return av_rescale_rnd(a, b, c, rnd);
140 }
141 
142  int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
143 {
144  return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
145 }
146 
147  int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
148 {
149  int64_t a = tb_a.num * (int64_t)tb_b.den;
150  int64_t b = tb_b.num * (int64_t)tb_a.den;
151  if ((FFABS64U(ts_a)|a|FFABS64U(ts_b)|b) <= INT_MAX)
152  return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
153  if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
154  return -1;
155  if (av_rescale_rnd(ts_b, b, a, AV_ROUND_DOWN) < ts_a)
156  return 1;
157  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  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  av_assert0(in_ts != AV_NOPTS_VALUE);
172  av_assert0(duration >= 0);
173 
174  if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
175 simple_round:
176  *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
177  return av_rescale_q(in_ts, in_tb, out_tb);
178  }
179 
180  a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
181  b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
182  if (*last < 2*a - b || *last > 2*b - a)
183  goto simple_round;
184 
185  this = av_clip64(*last, a, b);
186  *last = this + duration;
187 
188  return av_rescale_q(this, fs_tb, out_tb);
189 }
190 
191  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  if (inc != 1)
196  inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
197 
198  m = inc_tb.num * (int64_t)ts_tb.den;
199  d = inc_tb.den * (int64_t)ts_tb.num;
200 
201  if (m % d == 0 && ts <= INT64_MAX - m / d)
202  return ts + m / d;
203  if (m < d)
204  return ts;
205 
206  {
207  int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
208  int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
209 
210  if (old == INT64_MAX || old == AV_NOPTS_VALUE || old_ts == AV_NOPTS_VALUE)
211  return ts;
212 
213  return av_sat_add64(av_rescale_q(old + 1, inc_tb, ts_tb), ts - old_ts);
214  }
215 }
216 
217  static inline double eval_poly(const double *coeff, int size, double x) {
218  double sum = coeff[size-1];
219  int i;
220  for (i = size-2; i >= 0; --i) {
221  sum *= x;
222  sum += coeff[i];
223  }
224  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  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  if (x == 0)
307  return 1.0;
308  x = fabs(x);
309  if (x <= 15) {
310  y = x * x;
311  return eval_poly(p1, FF_ARRAY_ELEMS(p1), y) / eval_poly(q1, FF_ARRAY_ELEMS(q1), y);
312  }
313  else {
314  y = 1 / x - 1.0 / 15;
315  r = eval_poly(p2, FF_ARRAY_ELEMS(p2), y) / eval_poly(q2, FF_ARRAY_ELEMS(q2), y);
316  factor = exp(x) / sqrt(x);
317  return factor * r;
318  }
319 }
AV_ROUND_UP
@ AV_ROUND_UP
Round toward +infinity.
Definition: mathematics.h:134
q1
static const uint8_t q1[256]
Definition: twofish.c:100
r
const char * r
Definition: vf_curves.c:127
av_compare_ts
int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
Compare two timestamps each in its own time base.
Definition: mathematics.c:147
av_add_stable
int64_t av_add_stable(AVRational ts_tb, int64_t ts, AVRational inc_tb, int64_t inc)
Add a value to a timestamp.
Definition: mathematics.c:191
av_bessel_i0
double av_bessel_i0(double x)
0th order modified bessel function of the first kind.
Definition: mathematics.c:257
int64_t
long long int64_t
Definition: coverity.c:34
u
#define u(width, name, range_min, range_max)
Definition: cbs_apv.c:68
b
#define b
Definition: input.c:42
AVRounding
AVRounding
Rounding methods.
Definition: mathematics.h:130
mathematics.h
FFMAX
#define FFMAX(a, b)
Definition: macros.h:47
av_i2int
int64_t av_i2int(AVInteger a)
Convert the given AVInteger to an int64_t.
Definition: integer.c:160
av_gcd
int64_t av_gcd(int64_t a, int64_t b)
Compute the greatest common divisor of two integer operands.
Definition: mathematics.c:37
b1
static double b1(void *priv, double x, double y)
Definition: vf_xfade.c:2034
a2
static double a2(void *priv, double x, double y)
Definition: vf_xfade.c:2030
AVRational::num
int num
Numerator.
Definition: rational.h:59
av_clip64
#define av_clip64
Definition: common.h:103
avassert.h
rnd
#define rnd()
Definition: checkasm.h:190
FF_ARRAY_ELEMS
#define FF_ARRAY_ELEMS(a)
Definition: sinewin_tablegen.c:29
duration
int64_t duration
Definition: movenc.c:65
av_int2i
AVInteger av_int2i(int64_t a)
Convert the given int64_t to an AVInteger.
Definition: integer.c:149
av_assert0
#define av_assert0(cond)
assert() equivalent, that is always enabled.
Definition: avassert.h:41
limits.h
eval_poly
static double eval_poly(const double *coeff, int size, double x)
Definition: mathematics.c:217
av_rescale_q
int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
Rescale a 64-bit integer by 2 rational numbers.
Definition: mathematics.c:142
AV_ROUND_NEAR_INF
@ AV_ROUND_NEAR_INF
Round to nearest and halfway cases away from zero.
Definition: mathematics.h:135
AV_ROUND_PASS_MINMAX
@ AV_ROUND_PASS_MINMAX
Flag telling rescaling functions to pass INT64_MIN/MAX through unchanged, avoiding special cases for ...
Definition: mathematics.h:159
av_rescale_delta
int64_t av_rescale_delta(AVRational in_tb, int64_t in_ts, AVRational fs_tb, int duration, int64_t *last, AVRational out_tb)
Rescale a timestamp while preserving known durations.
Definition: mathematics.c:168
av_add_i
AVInteger av_add_i(AVInteger a, AVInteger b)
Definition: integer.c:36
fabs
static __device__ float fabs(float a)
Definition: cuda_runtime.h:182
AVRational
Rational number (pair of numerator and denominator).
Definition: rational.h:58
ff_ctzll
#define ff_ctzll
Definition: intmath.h:125
exp
int8_t exp
Definition: eval.c:73
c
Undefined Behavior In the C some operations are like signed integer dereferencing freed accessing outside allocated Undefined Behavior must not occur in a C it is not safe even if the output of undefined operations is unused The unsafety may seem nit picking but Optimizing compilers have in fact optimized code on the assumption that no undefined Behavior occurs Optimizing code based on wrong assumptions can and has in some cases lead to effects beyond the output of computations The signed integer overflow problem in speed critical code Code which is highly optimized and works with signed integers sometimes has the problem that often the output of the computation does not c
Definition: undefined.txt:32
inc
static int inc(int num, int period)
Definition: perlin.c:34
AV_ROUND_DOWN
@ AV_ROUND_DOWN
Round toward -infinity.
Definition: mathematics.h:133
av_rescale_rnd
int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
Rescale a 64-bit integer with specified rounding.
Definition: mathematics.c:58
size
int size
Definition: twinvq_data.h:10344
AV_NOPTS_VALUE
#define AV_NOPTS_VALUE
Undefined timestamp value.
Definition: avutil.h:247
av_mul_i
AVInteger av_mul_i(AVInteger a, AVInteger b)
Definition: integer.c:66
a
The reader does not expect b to be semantically here and if the code is changed by maybe adding a a division or other the signedness will almost certainly be mistaken To avoid this confusion a new type was SUINT is the C unsigned type but it holds a signed int to use the same example SUINT a
Definition: undefined.txt:41
a0
static double a0(void *priv, double x, double y)
Definition: vf_xfade.c:2028
av_assert2
#define av_assert2(cond)
assert() equivalent, that does lie in speed critical code.
Definition: avassert.h:68
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:256
AVInteger
Definition: integer.h:36
common.h
FFMIN
#define FFMIN(a, b)
Definition: macros.h:49
av_rescale
int64_t av_rescale(int64_t a, int64_t b, int64_t c)
Rescale a 64-bit integer with rounding to nearest.
Definition: mathematics.c:129
mod
static int mod(int a, int b)
Modulo operation with only positive remainders.
Definition: vf_v360.c:755
av_compare_mod
int64_t av_compare_mod(uint64_t a, uint64_t b, uint64_t mod)
Compare the remainders of two integer operands divided by a common divisor.
Definition: mathematics.c:160
FFSWAP
#define FFSWAP(type, a, b)
Definition: macros.h:52
av_div_i
AVInteger av_div_i(AVInteger a, AVInteger b)
Return a/b.
Definition: integer.c:143
av_sat_add64
#define av_sat_add64
Definition: common.h:139
AVRational::den
int den
Denominator.
Definition: rational.h:60
av_mul_q
AVRational av_mul_q(AVRational b, AVRational c)
Multiply two rationals.
Definition: rational.c:80
factor
static const int factor[16]
Definition: vf_pp7.c:80
avutil.h
coeff
static const double coeff[2][5]
Definition: vf_owdenoise.c:80
b0
static double b0(void *priv, double x, double y)
Definition: vf_xfade.c:2033
a1
static double a1(void *priv, double x, double y)
Definition: vf_xfade.c:2029
FFABS64U
#define FFABS64U(a)
Definition: common.h:92
av_rescale_q_rnd
int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq, enum AVRounding rnd)
Rescale a 64-bit integer by 2 rational numbers with specified rounding.
Definition: mathematics.c:134
intmath.h
Generated on Fri Oct 17 2025 19:23:37 for FFmpeg by   doxygen 1.8.17

AltStyle によって変換されたページ (->オリジナル) /