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