1 /* 2 * Copyright (c) 2018 Iván Santa María <ghevan@gmail.com> 3 * 4 * This program is free software; you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation; either version 2 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program; if not, write to the Free Software 16 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 17 */ 18 19 #ifndef VC_ADDITIONAL_MATH_H 20 #define VC_ADDITIONAL_MATH_H 21 22 #include <config-vc.h> 23 24 #if defined HAVE_VC 25 26 #include <Vc/Vc> 27 #include <Vc/IO> 28 29 class VcExtraMath 30 { 31 public: 32 // vectorized erf function, precision 1e-5 erf(Vc::float_v x)33 static inline Vc::float_v erf(Vc::float_v x) { 34 Vc::float_v xa = abs(x); 35 Vc::float_m precisionLimit(xa >= 9.3f); // wrong result for any number beyond this 36 xa(precisionLimit) = 0; 37 Vc::float_v sign(Vc::One); 38 Vc::float_m invertMask = x < 0.f; 39 sign(invertMask) = -1.f; 40 41 // CONSTANTS 42 float a1 = 0.254829592; 43 float a2 = -0.284496736; 44 float a3 = 1.421413741; 45 float a4 = -1.453152027; 46 float a5 = 1.061405429; 47 float p = 0.3275911; 48 49 Vc::float_v t = 1.0f / (1.0f + p * xa); 50 Vc::float_v y = 1.0f - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * exp(-xa * xa); 51 y(precisionLimit) = 1.0f; 52 return sign * y; 53 } 54 }; 55 #endif /* defined HAVE_VC */ 56 57 58 #endif // VC_ADDITIONAL_MATH_H 59