FFmpeg
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
libm.h
Go to the documentation of this file.
1 /*
2  * erf function: Copyright (c) 2006 John Maddock
3  * This file is part of FFmpeg.
4  *
5  * FFmpeg is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * FFmpeg is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with FFmpeg; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 /**
21  * @file
22  * Replacements for frequently missing libm functions
23  */
24 
25 #ifndef AVUTIL_LIBM_H
26 #define AVUTIL_LIBM_H
27 
28 #include <math.h>
29 #include "config.h"
30 #include "attributes.h"
31 #if !(HAVE_COPYSIGN && HAVE_HYPOT && HAVE_ISFINITE && HAVE_ISINF && HAVE_ISNAN)
32 #include "intfloat.h"
33 #endif
34 #include "mathematics.h"
35 
36 #if HAVE_MIPSFPU && HAVE_INLINE_ASM
38 #endif /* HAVE_MIPSFPU && HAVE_INLINE_ASM*/
39 
40 #if !HAVE_ATANF
41 #undef atanf
42 #define atanf(x) ((float)atan(x))
43 #endif /* HAVE_ATANF */
44 
45 #if !HAVE_ATAN2F
46 #undef atan2f
47 #define atan2f(y, x) ((float)atan2(y, x))
48 #endif /* HAVE_ATAN2F */
49 
50 #if !HAVE_POWF
51 #undef powf
52 #define powf(x, y) ((float)pow(x, y))
53 #endif /* HAVE_POWF */
54 
55 #if !HAVE_CBRT
56 static av_always_inline double cbrt(double x)
57 {
58  return x < 0 ? -pow(-x, 1.0 / 3.0) : pow(x, 1.0 / 3.0);
59 }
60 #endif /* HAVE_CBRT */
61 
62 #if !HAVE_CBRTF
63 static av_always_inline float cbrtf(float x)
64 {
65  return x < 0 ? -powf(-x, 1.0 / 3.0) : powf(x, 1.0 / 3.0);
66 }
67 #endif /* HAVE_CBRTF */
68 
69 #if !HAVE_COPYSIGN
70 static av_always_inline double copysign(double x, double y)
71 {
72  uint64_t vx = av_double2int(x);
73  uint64_t vy = av_double2int(y);
74  return av_int2double((vx & UINT64_C(0x7fffffffffffffff)) | (vy & UINT64_C(0x8000000000000000)));
75 }
76 #endif /* HAVE_COPYSIGN */
77 
78 #if !HAVE_COSF
79 #undef cosf
80 #define cosf(x) ((float)cos(x))
81 #endif /* HAVE_COSF */
82 
83 #if !HAVE_ERF
84 static inline double ff_eval_poly(const double *coeff, int size, double x) {
85  double sum = coeff[size-1];
86  int i;
87  for (i = size-2; i >= 0; --i) {
88  sum *= x;
89  sum += coeff[i];
90  }
91  return sum;
92 }
93 
94 /**
95  * erf function
96  * Algorithm taken from the Boost project, source:
97  * http://www.boost.org/doc/libs/1_46_1/boost/math/special_functions/erf.hpp
98  * Use, modification and distribution are subject to the
99  * Boost Software License, Version 1.0 (see notice below).
100  * Boost Software License - Version 1.0 - August 17th, 2003
101 Permission is hereby granted, free of charge, to any person or organization
102 obtaining a copy of the software and accompanying documentation covered by
103 this license (the "Software") to use, reproduce, display, distribute,
104 execute, and transmit the Software, and to prepare derivative works of the
105 Software, and to permit third-parties to whom the Software is furnished to
106 do so, all subject to the following:
107 
108 The copyright notices in the Software and this entire statement, including
109 the above license grant, this restriction and the following disclaimer,
110 must be included in all copies of the Software, in whole or in part, and
111 all derivative works of the Software, unless such copies or derivative
112 works are solely in the form of machine-executable object code generated by
113 a source language processor.
114 
115 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
116 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
117 FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
118 SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
119 FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
120 ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
121 DEALINGS IN THE SOFTWARE.
122  */
123 static inline double erf(double z)
124 {
125 #ifndef FF_ARRAY_ELEMS
126 #define FF_ARRAY_ELEMS(a) (sizeof(a) / sizeof((a)[0]))
127 #endif
128  double result;
129 
130  /* handle the symmetry: erf(-x) = -erf(x) */
131  if (z < 0)
132  return -erf(-z);
133 
134  /* branch based on range of z, and pick appropriate approximation */
135  if (z == 0)
136  return 0;
137  else if (z < 1e-10)
138  return z * 1.125 + z * 0.003379167095512573896158903121545171688;
139  else if (z < 0.5) {
140  // Maximum Deviation Found: 1.561e-17
141  // Expected Error Term: 1.561e-17
142  // Maximum Relative Change in Control Points: 1.155e-04
143  // Max Error found at double precision = 2.961182e-17
144 
145  static const double y = 1.044948577880859375;
146  static const double p[] = {
147  0.0834305892146531832907,
148  -0.338165134459360935041,
149  -0.0509990735146777432841,
150  -0.00772758345802133288487,
151  -0.000322780120964605683831,
152  };
153  static const double q[] = {
154  1,
155  0.455004033050794024546,
156  0.0875222600142252549554,
157  0.00858571925074406212772,
158  0.000370900071787748000569,
159  };
160  double zz = z * z;
161  return z * (y + ff_eval_poly(p, FF_ARRAY_ELEMS(p), zz) / ff_eval_poly(q, FF_ARRAY_ELEMS(q), zz));
162  }
163  /* here onwards compute erfc */
164  else if (z < 1.5) {
165  // Maximum Deviation Found: 3.702e-17
166  // Expected Error Term: 3.702e-17
167  // Maximum Relative Change in Control Points: 2.845e-04
168  // Max Error found at double precision = 4.841816e-17
169  static const double y = 0.405935764312744140625;
170  static const double p[] = {
171  -0.098090592216281240205,
172  0.178114665841120341155,
173  0.191003695796775433986,
174  0.0888900368967884466578,
175  0.0195049001251218801359,
176  0.00180424538297014223957,
177  };
178  static const double q[] = {
179  1,
180  1.84759070983002217845,
181  1.42628004845511324508,
182  0.578052804889902404909,
183  0.12385097467900864233,
184  0.0113385233577001411017,
185  0.337511472483094676155e-5,
186  };
187  result = y + ff_eval_poly(p, FF_ARRAY_ELEMS(p), z - 0.5) / ff_eval_poly(q, FF_ARRAY_ELEMS(q), z - 0.5);
188  result *= exp(-z * z) / z;
189  return 1 - result;
190  }
191  else if (z < 2.5) {
192  // Max Error found at double precision = 6.599585e-18
193  // Maximum Deviation Found: 3.909e-18
194  // Expected Error Term: 3.909e-18
195  // Maximum Relative Change in Control Points: 9.886e-05
196  static const double y = 0.50672817230224609375;
197  static const double p[] = {
198  -0.0243500476207698441272,
199  0.0386540375035707201728,
200  0.04394818964209516296,
201  0.0175679436311802092299,
202  0.00323962406290842133584,
203  0.000235839115596880717416,
204  };
205  static const double q[] = {
206  1,
207  1.53991494948552447182,
208  0.982403709157920235114,
209  0.325732924782444448493,
210  0.0563921837420478160373,
211  0.00410369723978904575884,
212  };
213  result = y + ff_eval_poly(p, FF_ARRAY_ELEMS(p), z - 1.5) / ff_eval_poly(q, FF_ARRAY_ELEMS(q), z - 1.5);
214  result *= exp(-z * z) / z;
215  return 1 - result;
216  }
217  else if (z < 4.5) {
218  // Maximum Deviation Found: 1.512e-17
219  // Expected Error Term: 1.512e-17
220  // Maximum Relative Change in Control Points: 2.222e-04
221  // Max Error found at double precision = 2.062515e-17
222  static const double y = 0.5405750274658203125;
223  static const double p[] = {
224  0.00295276716530971662634,
225  0.0137384425896355332126,
226  0.00840807615555585383007,
227  0.00212825620914618649141,
228  0.000250269961544794627958,
229  0.113212406648847561139e-4,
230  };
231  static const double q[] = {
232  1,
233  1.04217814166938418171,
234  0.442597659481563127003,
235  0.0958492726301061423444,
236  0.0105982906484876531489,
237  0.000479411269521714493907,
238  };
239  result = y + ff_eval_poly(p, FF_ARRAY_ELEMS(p), z - 3.5) / ff_eval_poly(q, FF_ARRAY_ELEMS(q), z - 3.5);
240  result *= exp(-z * z) / z;
241  return 1 - result;
242  }
243  /* differ from Boost here, the claim of underflow of erfc(x) past 5.8 is
244  * slightly incorrect, change to 5.92
245  * (really somewhere between 5.9125 and 5.925 is when it saturates) */
246  else if (z < 5.92) {
247  // Max Error found at double precision = 2.997958e-17
248  // Maximum Deviation Found: 2.860e-17
249  // Expected Error Term: 2.859e-17
250  // Maximum Relative Change in Control Points: 1.357e-05
251  static const double y = 0.5579090118408203125;
252  static const double p[] = {
253  0.00628057170626964891937,
254  0.0175389834052493308818,
255  -0.212652252872804219852,
256  -0.687717681153649930619,
257  -2.5518551727311523996,
258  -3.22729451764143718517,
259  -2.8175401114513378771,
260  };
261  static const double q[] = {
262  1,
263  2.79257750980575282228,
264  11.0567237927800161565,
265  15.930646027911794143,
266  22.9367376522880577224,
267  13.5064170191802889145,
268  5.48409182238641741584,
269  };
270  result = y + ff_eval_poly(p, FF_ARRAY_ELEMS(p), 1 / z) / ff_eval_poly(q, FF_ARRAY_ELEMS(q), 1 / z);
271  result *= exp(-z * z) / z;
272  return 1 - result;
273  }
274  /* handle the nan case, but don't use isnan for max portability */
275  else if (z != z)
276  return z;
277  /* finally return saturated result */
278  else
279  return 1;
280 }
281 #endif /* HAVE_ERF */
282 
283 #if !HAVE_EXPF
284 #undef expf
285 #define expf(x) ((float)exp(x))
286 #endif /* HAVE_EXPF */
287 
288 #if !HAVE_EXP2
289 #undef exp2
290 #define exp2(x) exp((x) * M_LN2)
291 #endif /* HAVE_EXP2 */
292 
293 #if !HAVE_EXP2F
294 #undef exp2f
295 #define exp2f(x) ((float)exp2(x))
296 #endif /* HAVE_EXP2F */
297 
298 #if !HAVE_ISINF
299 #undef isinf
300 /* Note: these do not follow the BSD/Apple/GNU convention of returning -1 for
301 -Inf, +1 for Inf, 0 otherwise, but merely follow the POSIX/ISO mandated spec of
302 returning a non-zero value for +/-Inf, 0 otherwise. */
304 {
305  uint32_t v = av_float2int(x);
306  if ((v & 0x7f800000) != 0x7f800000)
307  return 0;
308  return !(v & 0x007fffff);
309 }
310 
312 {
313  uint64_t v = av_double2int(x);
314  if ((v & 0x7ff0000000000000) != 0x7ff0000000000000)
315  return 0;
316  return !(v & 0x000fffffffffffff);
317 }
318 
319 #define isinf(x) \
320  (sizeof(x) == sizeof(float) \
321  ? avpriv_isinff(x) \
322  : avpriv_isinf(x))
323 #endif /* HAVE_ISINF */
324 
325 #if !HAVE_ISNAN
327 {
328  uint32_t v = av_float2int(x);
329  if ((v & 0x7f800000) != 0x7f800000)
330  return 0;
331  return v & 0x007fffff;
332 }
333 
335 {
336  uint64_t v = av_double2int(x);
337  if ((v & 0x7ff0000000000000) != 0x7ff0000000000000)
338  return 0;
339  return (v & 0x000fffffffffffff) && 1;
340 }
341 
342 #define isnan(x) \
343  (sizeof(x) == sizeof(float) \
344  ? avpriv_isnanf(x) \
345  : avpriv_isnan(x))
346 #endif /* HAVE_ISNAN */
347 
348 #if !HAVE_ISFINITE
350 {
351  uint32_t v = av_float2int(x);
352  return (v & 0x7f800000) != 0x7f800000;
353 }
354 
356 {
357  uint64_t v = av_double2int(x);
358  return (v & 0x7ff0000000000000) != 0x7ff0000000000000;
359 }
360 
361 #define isfinite(x) \
362  (sizeof(x) == sizeof(float) \
363  ? avpriv_isfinitef(x) \
364  : avpriv_isfinite(x))
365 #endif /* HAVE_ISFINITE */
366 
367 #if !HAVE_HYPOT
368 static inline av_const double hypot(double x, double y)
369 {
370  double temp;
371  x = fabs(x);
372  y = fabs(y);
373 
374  if (isinf(x) || isinf(y))
375  return av_int2double(0x7ff0000000000000);
376  if (x == 0 || y == 0)
377  return x + y;
378  if (x < y) {
379  temp = x;
380  x = y;
381  y = temp;
382  }
383 
384  y = y/x;
385  return x*sqrt(1 + y*y);
386 }
387 #endif /* HAVE_HYPOT */
388 
389 #if !HAVE_LDEXPF
390 #undef ldexpf
391 #define ldexpf(x, exp) ((float)ldexp(x, exp))
392 #endif /* HAVE_LDEXPF */
393 
394 #if !HAVE_LLRINT
395 #undef llrint
396 #define llrint(x) ((long long)rint(x))
397 #endif /* HAVE_LLRINT */
398 
399 #if !HAVE_LLRINTF
400 #undef llrintf
401 #define llrintf(x) ((long long)rint(x))
402 #endif /* HAVE_LLRINT */
403 
404 #if !HAVE_LOG2
405 #undef log2
406 #define log2(x) (log(x) * 1.44269504088896340736)
407 #endif /* HAVE_LOG2 */
408 
409 #if !HAVE_LOG2F
410 #undef log2f
411 #define log2f(x) ((float)log2(x))
412 #endif /* HAVE_LOG2F */
413 
414 #if !HAVE_LOG10F
415 #undef log10f
416 #define log10f(x) ((float)log10(x))
417 #endif /* HAVE_LOG10F */
418 
419 #if !HAVE_SINF
420 #undef sinf
421 #define sinf(x) ((float)sin(x))
422 #endif /* HAVE_SINF */
423 
424 #if !HAVE_RINT
425 static inline double rint(double x)
426 {
427  return x >= 0 ? floor(x + 0.5) : ceil(x - 0.5);
428 }
429 #endif /* HAVE_RINT */
430 
431 #if !HAVE_LRINT
432 static av_always_inline av_const long int lrint(double x)
433 {
434  return rint(x);
435 }
436 #endif /* HAVE_LRINT */
437 
438 #if !HAVE_LRINTF
439 static av_always_inline av_const long int lrintf(float x)
440 {
441  return (int)(rint(x));
442 }
443 #endif /* HAVE_LRINTF */
444 
445 #if !HAVE_ROUND
446 static av_always_inline av_const double round(double x)
447 {
448  return (x > 0) ? floor(x + 0.5) : ceil(x - 0.5);
449 }
450 #endif /* HAVE_ROUND */
451 
452 #if !HAVE_ROUNDF
453 static av_always_inline av_const float roundf(float x)
454 {
455  return (x > 0) ? floor(x + 0.5) : ceil(x - 0.5);
456 }
457 #endif /* HAVE_ROUNDF */
458 
459 #if !HAVE_TRUNC
460 static av_always_inline av_const double trunc(double x)
461 {
462  return (x > 0) ? floor(x) : ceil(x);
463 }
464 #endif /* HAVE_TRUNC */
465 
466 #if !HAVE_TRUNCF
467 static av_always_inline av_const float truncf(float x)
468 {
469  return (x > 0) ? floor(x) : ceil(x);
470 }
471 #endif /* HAVE_TRUNCF */
472 
473 #endif /* AVUTIL_LIBM_H */
av_int2double
static av_always_inline double av_int2double(uint64_t i)
Reinterpret a 64-bit integer as a double.
Definition: intfloat.h:60
truncf
static av_always_inline av_const float truncf(float x)
Definition: libm.h:467
av_const
#define av_const
Definition: attributes.h:84
lrintf
static av_always_inline av_const long int lrintf(float x)
Definition: libm.h:439
mathematics.h
av_float2int
static av_always_inline uint32_t av_float2int(float f)
Reinterpret a float as a 32-bit integer.
Definition: intfloat.h:50
intfloat.h
roundf
static av_always_inline av_const float roundf(float x)
Definition: libm.h:453
cbrt
static av_always_inline double cbrt(double x)
Definition: libm.h:56
ceil
static __device__ float ceil(float a)
Definition: cuda_runtime.h:176
floor
static __device__ float floor(float a)
Definition: cuda_runtime.h:173
avpriv_isinff
static av_always_inline av_const int avpriv_isinff(float x)
Definition: libm.h:303
copysign
static av_always_inline double copysign(double x, double y)
Definition: libm.h:70
result
and forward the result(frame or status change) to the corresponding input. If nothing is possible
fabs
static __device__ float fabs(float a)
Definition: cuda_runtime.h:182
isinf
#define isinf(x)
Definition: libm.h:319
avpriv_isnan
static av_always_inline av_const int avpriv_isnan(double x)
Definition: libm.h:334
rint
static double rint(double x)
Definition: libm.h:425
exp
int8_t exp
Definition: eval.c:73
ff_eval_poly
static double ff_eval_poly(const double *coeff, int size, double x)
Definition: libm.h:84
trunc
static av_always_inline av_const double trunc(double x)
Definition: libm.h:460
powf
#define powf(x, y)
Definition: libm.h:52
hypot
static av_const double hypot(double x, double y)
Definition: libm.h:368
size
int size
Definition: twinvq_data.h:10344
libm_mips.h
attributes.h
avpriv_isinf
static av_always_inline av_const int avpriv_isinf(double x)
Definition: libm.h:311
av_double2int
static av_always_inline uint64_t av_double2int(double f)
Reinterpret a double as a 64-bit integer.
Definition: intfloat.h:70
FF_ARRAY_ELEMS
#define FF_ARRAY_ELEMS(a)
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:256
round
static av_always_inline av_const double round(double x)
Definition: libm.h:446
av_always_inline
#define av_always_inline
Definition: attributes.h:49
cbrtf
static av_always_inline float cbrtf(float x)
Definition: libm.h:63
erf
static double erf(double z)
erf function Algorithm taken from the Boost project, source: http://www.boost.org/doc/libs/1_46_1/boo...
Definition: libm.h:123
avpriv_isfinite
static av_always_inline av_const int avpriv_isfinite(double x)
Definition: libm.h:355
temp
else temp
Definition: vf_mcdeint.c:263
avpriv_isfinitef
static av_always_inline av_const int avpriv_isfinitef(float x)
Definition: libm.h:349
lrint
static av_always_inline av_const long int lrint(double x)
Definition: libm.h:432
avpriv_isnanf
static av_always_inline av_const int avpriv_isnanf(float x)
Definition: libm.h:326
coeff
static const double coeff[2][5]
Definition: vf_owdenoise.c:80