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