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