Print this page
5261 libm should stop using synonyms.h
5298 fabs is 0-sized, confuses dis(1) and others
Reviewed by: Josef 'Jeff' Sipek <jeffpc@josefsipek.net>
Approved by: Gordon Ross <gwr@nexenta.com>


  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #pragma weak erf = __erf
  31 #pragma weak erfc = __erfc
  32 
  33 /* INDENT OFF */
  34 /*
  35  * double erf(double x)
  36  * double erfc(double x)
  37  *                           x
  38  *                    2      |\
  39  *     erf(x)  =  ---------  | exp(-t*t)dt
  40  *                 sqrt(pi) \|
  41  *                           0
  42  *
  43  *     erfc(x) =  1-erf(x)
  44  *  Note that
  45  *              erf(-x) = -erf(x)
  46  *              erfc(-x) = 2 - erfc(x)
  47  *
  48  * Method:
  49  *      1. For |x| in [0, 0.84375]
  50  *          erf(x)  = x + x*R(x^2)
  51  *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]


 110  *              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
 111  *                        x*sqrt(pi)
 112  *         We use rational approximation to approximate
 113  *              g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
 114  *         Here is the error bound for R1/S1 and R2/S2
 115  *              |R1/S1 - f(x)|  < 2**(-62.57)
 116  *              |R2/S2 - f(x)|  < 2**(-61.52)
 117  *
 118  *      5. For inf > x >= 28
 119  *              erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
 120  *              erfc(x) = tiny*tiny (raise underflow) if x > 0
 121  *                      = 2 - tiny if x<0
 122  *
 123  *      7. Special case:
 124  *              erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
 125  *              erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
 126  *      erfc/erf(NaN) is NaN
 127  */
 128 /* INDENT ON */
 129 
 130 #include "libm_synonyms.h"      /* __erf, __erfc, __exp */
 131 #include "libm_macros.h"
 132 #include <math.h>
 133 
 134 static const double xxx[] = {
 135 /* tiny */      1e-300,
 136 /* half */      5.00000000000000000000e-01,     /* 3FE00000, 00000000 */
 137 /* one */       1.00000000000000000000e+00,     /* 3FF00000, 00000000 */
 138 /* two */       2.00000000000000000000e+00,     /* 40000000, 00000000 */
 139 /* erx */       8.45062911510467529297e-01,     /* 3FEB0AC1, 60000000 */
 140 /*
 141  * Coefficients for approximation to  erf on [0,0.84375]
 142  */
 143 /* efx */        1.28379167095512586316e-01,    /* 3FC06EBA, 8214DB69 */
 144 /* efx8 */       1.02703333676410069053e+00,    /* 3FF06EBA, 8214DB69 */
 145 /* pp0 */        1.28379167095512558561e-01,    /* 3FC06EBA, 8214DB68 */
 146 /* pp1 */       -3.25042107247001499370e-01,    /* BFD4CD7D, 691CB913 */
 147 /* pp2 */       -2.84817495755985104766e-02,    /* BF9D2A51, DBD7194F */
 148 /* pp3 */       -5.77027029648944159157e-03,    /* BF77A291, 236668E4 */
 149 /* pp4 */       -2.37630166566501626084e-05,    /* BEF8EAD6, 120016AC */
 150 /* qq1 */        3.97917223959155352819e-01,    /* 3FD97779, CDDADC09 */




  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #pragma weak __erf = erf
  31 #pragma weak __erfc = erfc
  32 
  33 /* INDENT OFF */
  34 /*
  35  * double erf(double x)
  36  * double erfc(double x)
  37  *                           x
  38  *                    2      |\
  39  *     erf(x)  =  ---------  | exp(-t*t)dt
  40  *                 sqrt(pi) \|
  41  *                           0
  42  *
  43  *     erfc(x) =  1-erf(x)
  44  *  Note that
  45  *              erf(-x) = -erf(x)
  46  *              erfc(-x) = 2 - erfc(x)
  47  *
  48  * Method:
  49  *      1. For |x| in [0, 0.84375]
  50  *          erf(x)  = x + x*R(x^2)
  51  *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]


 110  *              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
 111  *                        x*sqrt(pi)
 112  *         We use rational approximation to approximate
 113  *              g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
 114  *         Here is the error bound for R1/S1 and R2/S2
 115  *              |R1/S1 - f(x)|  < 2**(-62.57)
 116  *              |R2/S2 - f(x)|  < 2**(-61.52)
 117  *
 118  *      5. For inf > x >= 28
 119  *              erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
 120  *              erfc(x) = tiny*tiny (raise underflow) if x > 0
 121  *                      = 2 - tiny if x<0
 122  *
 123  *      7. Special case:
 124  *              erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
 125  *              erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
 126  *      erfc/erf(NaN) is NaN
 127  */
 128 /* INDENT ON */
 129 

 130 #include "libm_macros.h"
 131 #include <math.h>
 132 
 133 static const double xxx[] = {
 134 /* tiny */      1e-300,
 135 /* half */      5.00000000000000000000e-01,     /* 3FE00000, 00000000 */
 136 /* one */       1.00000000000000000000e+00,     /* 3FF00000, 00000000 */
 137 /* two */       2.00000000000000000000e+00,     /* 40000000, 00000000 */
 138 /* erx */       8.45062911510467529297e-01,     /* 3FEB0AC1, 60000000 */
 139 /*
 140  * Coefficients for approximation to  erf on [0,0.84375]
 141  */
 142 /* efx */        1.28379167095512586316e-01,    /* 3FC06EBA, 8214DB69 */
 143 /* efx8 */       1.02703333676410069053e+00,    /* 3FF06EBA, 8214DB69 */
 144 /* pp0 */        1.28379167095512558561e-01,    /* 3FC06EBA, 8214DB68 */
 145 /* pp1 */       -3.25042107247001499370e-01,    /* BFD4CD7D, 691CB913 */
 146 /* pp2 */       -2.84817495755985104766e-02,    /* BF9D2A51, DBD7194F */
 147 /* pp3 */       -5.77027029648944159157e-03,    /* BF77A291, 236668E4 */
 148 /* pp4 */       -2.37630166566501626084e-05,    /* BEF8EAD6, 120016AC */
 149 /* qq1 */        3.97917223959155352819e-01,    /* 3FD97779, CDDADC09 */