math: p_erf: Implement error function

src/math/p_erf.c

@@ -1,6 +1,44 @@
 #include <pal.h>
 /**
 *
+* Developer Certificate of Origin
+* Version 1.1
+*
+* Copyright (C) 2004, 2006 The Linux Foundation and its contributors.
+* 660 York Street, Suite 102,
+* San Francisco, CA 94110 USA
+*
+* Everyone is permitted to copy and distribute verbatim copies of this
+* license document, but changing it is not allowed.
+*
+* Developer's Certificate of Origin 1.1
+*
+* By making a contribution to this project, I certify that:
+*
+* (a) The contribution was created in whole or in part by me and I
+*     have the right to submit it under the open source license
+*     indicated in the file; or
+*
+* (b) The contribution is based upon previous work that, to the best
+*     of my knowledge, is covered under an appropriate open source
+*     license and I have the right under that license to submit that
+*     work with modifications, whether created in whole or in part
+*     by me, under the same open source license (unless I am
+*     permitted to submit under a different license), as indicated
+*     in the file; or
+*
+* (c) The contribution was provided directly to me by some other
+*     person who certified (a), (b) or (c) and I have not modified
+*     it.
+*
+* (d) I understand and agree that this project and the contribution
+*     are public and that a record of the contribution (including all
+*     personal information I submit with it, including my sign-off) is
+*     maintained indefinitely and may be redistributed consistent with
+*     this project or the open source license(s) involved.
+*
+*     Signed-off-by: Adam Szewczyk <[email protected]>
+* 
 * Calculate error function.
 *
 *
@@ -20,15 +58,17 @@ void p_erf_f32(const float *a, float *c, int n)
 { 
     int i; int sign; float x;
     for ( i = 0; i < n; i++)
-        { 
-	     const float *pa = (a+i); 
-         float *pc = (c+i);
-
-        // the computation is based on formula 7.1.28 from Abramowitz and Stegun for a >= 0
+      { 
+       const float *pa = (a+i); 
+             float *pc = (c+i);
 
-        if (*pa >= 0) sign = 1.0; else sign = -1.0;                 // store sign of input variable
-        if (*pa < 0) x = *pa * (-1.0f); else x = *pa;              // this is equivalent of x=abs(a)
+       // the computation is based on formula 7.1.28 from Abramowitz and Stegun for a >= 0
 
+       if (*pa >= 0) sign = 1.0; else sign = -1.0;       // store sign of input variable
+       if (*pa < 0) x = *pa * (-1.0f); else x = *pa;     // this is equivalent of x=abs(a)
+       if (x >= 3.6f) *pc = 1.0f * sign;                 // c = 1 * sign if a >=3.6 else calculate error function
+       else
+	{
         // constants for formula 7.1.28 from Abramowitz and Stegun
         float a1 =  0.0705230784f;
         float a2 =  0.0422820123f;
@@ -42,4 +82,5 @@ void p_erf_f32(const float *a, float *c, int n)
         float z = 1.0f - (t * t * t * t * t * t * t * t * t * t * t * t * t * t * t * t) ;
         *pc = z*sign;
         }
+      }
 }