Open64: osprey/libm/mips/erf.c Source File

00001 /*
00002  * Copyright 2003, 2004, 2005, 2006 PathScale, Inc.  All Rights Reserved.
00003  */
00004 
00005 
00006 /*
00007 
00008   Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
00009 
00010   This program is free software; you can redistribute it and/or modify it
00011   under the terms of version 2.1 of the GNU Lesser General Public License 
00012   as published by the Free Software Foundation.
00013 
00014   This program is distributed in the hope that it would be useful, but
00015   WITHOUT ANY WARRANTY; without even the implied warranty of
00016   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
00017 
00018   Further, this software is distributed without any warranty that it is
00019   free of the rightful claim of any third person regarding infringement 
00020   or the like.  Any license provided herein, whether implied or 
00021   otherwise, applies only to this software file.  Patent licenses, if
00022   any, provided herein do not apply to combinations of this program with 
00023   other software, or any other product whatsoever.  
00024 
00025   You should have received a copy of the GNU Lesser General Public 
00026   License along with this program; if not, write the Free Software 
00027   Foundation, Inc., 59 Temple Place - Suite 330, Boston MA 02111-1307, 
00028   USA.
00029 
00030   Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
00031   Mountain View, CA 94043, or:
00032 
00033   http://www.sgi.com
00034 
00035   For further information regarding this notice, see:
00036 
00037   http://oss.sgi.com/projects/GenInfo/NoticeExplan
00038 
00039 */
00040 
00041 /* ====================================================================
00042  * ====================================================================
00043  *
00044  * Module: erf.c
00045  * $Revision: 1.5 $
00046  * $Date: 04/12/21 14:58:21-08:00 $
00047  * $Author: bos@eng-25.internal.keyresearch.com $
00048  * $Source: /home/bos/bk/kpro64-pending/libm/mips/SCCS/s.erf.c $
00049  *
00050  * Revision history:
00051  *  28-May-93 - Original Version
00052  *
00053  * Description: source code for erf and erfc functions
00054  *
00055  * ====================================================================
00056  * ====================================================================
00057  */
00058 
00059 /* $Header: /home/bos/bk/kpro64-pending/libm/mips/erf.c 1.5 04/12/21 14:58:21-08:00 bos@eng-25.internal.keyresearch.com $ */
00060 
00061 /*
00062   C program for floating point error function
00063 
00064   erf(x) returns the error function of its argument
00065   erfc(x) returns 1.0 - erf(x)
00066 
00067   erf(x) is defined by
00068   ${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$
00069 
00070   the entry for erfc is provided because of the
00071   extreme loss of relative accuracy if erf(x) is
00072   called for large x and the result subtracted
00073   from 1. (e.g. for x= 5, 12 places are lost).
00074 
00075   Calls exp.
00076 
00077   Based on the ideas in Hart & Cheney's book, Computer Approximations,
00078   John Wiley and Sons, 1968, Section 6.7.
00079 
00080   Uses rational approximations for erf on [0.0, 2.0], computing erfc(x)
00081   as 1.0 - erf(x).
00082   Beyond 2.0, computes erfc by exp(-x*x)*p(1/(x*x))/q(1/(x*x))/x
00083   and erf(x) = 1.0 - erfc(x).
00084 
00085 
00086 */
00087 
00088 #ifdef _CALL_MATHERR
00089 #include <stdio.h>
00090 #include <math.h>
00091 #endif
00092 
00093 #include "libm.h"
00094 
00095 #if defined(mips) && !defined(__GNUC__)
00096 extern  double  erf(double);
00097 extern  double  erfc(double);
00098 
00099 #pragma weak erf = __erf
00100 #pragma weak erfc = __erfc
00101 #endif
00102 
00103 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00104 extern double __erf(double);
00105 extern double __erfc(double);
00106 #pragma weak erf
00107 #pragma weak erfc
00108 double erf( double arg ) {
00109   return __erf( arg );
00110 }
00111 double erfc( double arg ) {
00112   return __erfc( arg );
00113 }
00114 #elif defined(__GNUC__)
00115 extern  double  __erf(double);
00116 
00117 double    erf() __attribute__ ((weak, alias ("__erf")));
00118 
00119 extern  double  __erfc(double);
00120 
00121 double    erfc() __attribute__ ((weak, alias ("__erfc")));
00122 
00123 #endif
00124 
00125 extern  double __exp(double);
00126 
00127 static const  du  Qnan = 
00128 {D(QNANHI, QNANLO)};
00129 
00130 static du Rsqrtpi = {D(0x3fe20dd7, 0x50429b6d)}; /* 1.0/sqrt(pi) */
00131 
00132 static du Ulimit1 = {D(0x4017afb4, 0x8dc96626)}; /* upper limit for erf */
00133 
00134 static du Llimit2 = {D(0xc017744f, 0x8f74e94a)}; /* lower limit for erfc */
00135 
00136 static du Ulimit2 = {D(0x403b39dc, 0x41e48bfc)}; /* upper limit for erfc */
00137 
00138 static double poly(double);
00139 
00140 /* coefficients for approximating erf on [-.75, .75] */
00141 
00142 static du p1[] = {
00143 D(0x3fc06eba, 0x8214db69),
00144 D(0xbfd49cb4, 0x2e1bc6cd),
00145 D(0xbfa23997, 0x1bd010b5),
00146 D(0xbf7ec3b8, 0x695750e1),
00147 D(0xbf3008ae, 0xd8b328e9),
00148 D(0xbef0d02f, 0xdfcc8cc0),
00149 D(0x3e800fef, 0x1ad6e676),
00150 };
00151 
00152 static du q1[]  = {
00153 D(0x3ff00000, 0x00000000),
00154 D(0x3fdaf37d, 0xbd166903),
00155 D(0x3fb3db3d, 0x1babb8d7),
00156 D(0x3f802445, 0xd3709b73),
00157 D(0x3f3d4818, 0xc17333c6),
00158 D(0x3ee80b2a, 0x2eab6d6a),
00159 D(0x00000000, 0x00000000),
00160 };
00161 
00162 /* coefficients for approximating erf on [.75, 1.25] */
00163 
00164 static du p2[] = {
00165 D(0x3fc42261, 0x62fbddd5),
00166 D(0x3feaf767, 0xa741088b),
00167 D(0x3fda911f, 0x096fbc26),
00168 D(0xbfc03d80, 0x8b1137e1),
00169 D(0x3fbbacfa, 0x66a0d1d5),
00170 D(0x3f9cfff4, 0x709b5c7d),
00171 D(0xbf5f504f, 0xfe135596),
00172 D(0x3f68d5ea, 0x7faa17b4),
00173 D(0x3f323c70, 0x817aefdb),
00174 };
00175 
00176 static du q2[] = {
00177 D(0x3ff00000, 0x00000000),
00178 D(0x3fe63821, 0x150312cf),
00179 D(0x3fe3e2f0, 0x81ee743e),
00180 D(0x3fd2a8e6, 0xefe9fae2),
00181 D(0x3fc0bc42, 0xfb9b210a),
00182 D(0x3fa5740f, 0x9a722ffb),
00183 D(0x3f83f746, 0x57cfbd94),
00184 D(0x3f5e1c3e, 0xb6d4d49f),
00185 };
00186 
00187 /* coefficients for approximating erf on [1.25, 1.75] */
00188 
00189 static du p3[] = {
00190 D(0x3fa15aaa, 0x8ec85205),
00191 D(0x3feeea55, 0x57137ae0),
00192 D(0x3fbe7237, 0x26b824a9),
00193 D(0xbfb9d1a8, 0xf32fd923),
00194 D(0x3fb77c31, 0x3b778138),
00195 D(0xbfa02178, 0x6e76caa3),
00196 D(0x3f822691, 0x77ac5924),
00197 D(0xbf46d6ce, 0xa8d3d9e2),
00198 };
00199 
00200 static du q3[] = {
00201 D(0x3ff00000, 0x00000000),
00202 D(0x3fe4dd02, 0xa28f292f),
00203 D(0x3fe2a536, 0x3e16bc6f),
00204 D(0x3fcc9a1d, 0x0eb16b1c),
00205 D(0x3fb9731a, 0x8b41322c),
00206 D(0x3f948ac7, 0x6117721e),
00207 D(0x3f73235a, 0x07cc2fd9),
00208 };
00209 
00210 /* coefficients for approximating erf on [1.75, 2.0] */
00211 
00212 static du p4[] = {
00213 D(0x3f806784, 0x42cc256f),
00214 D(0x3fefbe61, 0xeef4cf6a),
00215 D(0x3fa12ceb, 0x37ff9baf),
00216 D(0xbf9a06d6, 0x144eb107),
00217 D(0x3f9ab0c8, 0x0cca31bd),
00218 D(0xbf802a99, 0x5ed29ad3),
00219 D(0x3f5b1f65, 0x7b1d0822),
00220 };
00221 
00222 static du q4[] = {
00223 D(0x3ff00000, 0x00000000),
00224 D(0x3ff1e094, 0x50d81c1d),
00225 D(0x3feb91fe, 0xb4228831),
00226 D(0x3fd930cb, 0xc83fe50a),
00227 D(0x3fbf0110, 0x1cb1561b),
00228 D(0x3f91650c, 0xa6ca26a8),
00229 };
00230 
00231 /* coefficients for approximating erfc on [2.0, inf) in powers
00232  * of 1/(x*x)
00233  */
00234 
00235 static du p5[]  = {
00236 D(0xbfdfffff, 0xfffffff0),
00237 D(0xc0326677, 0xc6c50dc8),
00238 D(0xc06f48fc, 0x52b0853c),
00239 D(0xc098dbdb, 0xc3a76fc9),
00240 D(0xc0b33dee, 0x8ef175aa),
00241 D(0xc0bb7e1f, 0x01b13b59),
00242 D(0xc0aeacc7, 0xc40e8bc7),
00243 D(0xc08134cd, 0x96535f33),
00244 };
00245 
00246 static du q5[]  = {
00247 D(0x3ff00000, 0x00000000),
00248 D(0x40432677, 0xc6c50cd9),
00249 D(0x40815219, 0x63facf1f),
00250 D(0x40ae55a4, 0x447cc2b3),
00251 D(0x40cb6c9f, 0x19945427),
00252 D(0x40d94ebc, 0x318e1d47),
00253 D(0x40d618aa, 0xd7fba8cd),
00254 D(0x40be9acc, 0xf88372df),
00255 D(0x40859695, 0xf5e67d06),
00256 };
00257 
00258 /* coefficients for polynomial approximation of expm1 on +/- 1.8e-4 */
00259 
00260 static const du P[] =
00261 {
00262 D(0x3ff00000, 0x00000000),
00263 D(0x3fe00000, 0x00000000),
00264 D(0x3fc55555, 0x55555555),
00265 D(0x3fa55555, 0x55df1a6d),
00266 D(0x3f811111, 0x11fd3e5f),
00267 };
00268 
00269 
00270 double
00271 __erf( double arg )
00272 {
00273 double __erfc(double);
00274 double sign;
00275 double argsq;
00276 double d, n, z;
00277 #ifdef _CALL_MATHERR
00278 struct exception  exstruct;
00279 #endif
00280 
00281   if( arg != arg )
00282   {
00283     /* arg is a NaN */
00284 
00285 #ifdef _CALL_MATHERR
00286 
00287     exstruct.type = DOMAIN;
00288     exstruct.name = "erf";
00289     exstruct.arg1 = arg;
00290     exstruct.retval = Qnan.d;
00291 
00292     if ( matherr( &exstruct ) == 0 )
00293     {
00294       fprintf(stderr, "domain error in erf\n");
00295       SETERRNO(EDOM);
00296     }
00297 
00298     return ( exstruct.retval );
00299 #else
00300     NAN_SETERRNO(EDOM);
00301 
00302     return ( Qnan.d );
00303 #endif
00304   }
00305 
00306   sign = 1.0;
00307 
00308   if( arg < 0. )
00309   {
00310     arg = -arg;
00311     sign = -1.0;
00312   }
00313 
00314   if ( arg < 0.75 )
00315   {
00316     argsq = arg*arg;
00317 
00318 
00319     n = (((((p1[6].d*argsq + p1[5].d)*argsq + p1[4].d)*argsq +
00320         p1[3].d)*argsq + p1[2].d)*argsq + p1[1].d)*
00321         argsq + p1[0].d;
00322 
00323     d = ((((q1[5].d*argsq + q1[4].d)*argsq + q1[3].d)*argsq +
00324         q1[2].d)*argsq + q1[1].d)*argsq + q1[0].d;
00325 
00326     return ( sign*(arg + arg*n/d) );
00327 
00328   }
00329   else if ( arg < 1.25 )
00330   {
00331     z = arg - 1.0;
00332 
00333     n = ((((((p2[8].d*z + p2[7].d)*z + p2[6].d)*z + p2[5].d)*z +
00334       p2[4].d)*z + p2[3].d)*z + p2[2].d)*z;
00335     d = ((((((q2[7].d*z + q2[6].d)*z + q2[5].d)*z + q2[4].d)*z +
00336       q2[3].d)*z + q2[2].d)*z + q2[1].d)*z + q2[0].d;
00337 
00338     return ( sign*(p2[1].d + n/d) );
00339   }
00340   else if ( arg < 1.75 )
00341   {
00342     z = arg - 1.5;
00343 
00344     n = (((((p3[7].d*z + p3[6].d)*z +
00345       p3[5].d)*z + p3[4].d)*z + p3[3].d)*z + p3[2].d)*z;
00346     d = (((((q3[6].d*z + q3[5].d)*z + q3[4].d)*z +
00347       q3[3].d)*z + q3[2].d)*z + q3[1].d)*z + q3[0].d;
00348 
00349     return ( sign*(p3[1].d + n/d) );
00350   }
00351   else if ( arg < 2.0 )
00352   {
00353     z = arg - 1.875;
00354 
00355     n = ((((p4[6].d*z + p4[5].d)*z +
00356       p4[4].d)*z + p4[3].d)*z + p4[2].d)*z;
00357     d = ((((q4[5].d*z + q4[4].d)*z + q4[3].d)*z +
00358       q4[2].d)*z + q4[1].d)*z + q4[0].d;
00359 
00360     return ( sign*(p4[1].d + n/d) );
00361   }
00362 
00363   if ( fabs(arg) <= Ulimit1.d )
00364     return( sign*(1.0 - __erfc(arg)) );
00365   else
00366     return ( sign );
00367 }
00368 
00369 double
00370 __erfc( double arg )
00371 {
00372 double argsq;
00373 double z;
00374 double zsq;
00375 double s, f, f1;
00376 double n, d;
00377 double result;
00378 #ifdef _CALL_MATHERR
00379 struct exception  exstruct;
00380 #endif
00381 
00382   if( arg != arg )
00383   {
00384     /* arg is a NaN */
00385 
00386 #ifdef _CALL_MATHERR
00387 
00388     exstruct.type = DOMAIN;
00389     exstruct.name = "erfc";
00390     exstruct.arg1 = arg;
00391     exstruct.retval = Qnan.d;
00392 
00393     if ( matherr( &exstruct ) == 0 )
00394     {
00395       fprintf(stderr, "domain error in erfc\n");
00396       SETERRNO(EDOM);
00397     }
00398 
00399     return ( exstruct.retval );
00400 #else
00401     NAN_SETERRNO(EDOM);
00402 
00403     return ( Qnan.d );
00404 #endif
00405   }
00406 
00407   if ( arg < Llimit2.d )
00408     return ( 2.0 );
00409 
00410   if ( arg <= -2.0 )
00411     return ( 2.0 - __erfc(-arg) );
00412 
00413   if ( arg < 0.25 )
00414     return( 1.0 - __erf(arg) );
00415 
00416   result = 0.0; /* default */
00417 
00418   if ( arg < 0.75 )
00419   {
00420     argsq = arg*arg;
00421 
00422     n = (((((p1[6].d*argsq + p1[5].d)*argsq + p1[4].d)*argsq +
00423         p1[3].d)*argsq + p1[2].d)*argsq + p1[1].d)*
00424         argsq + p1[0].d;
00425 
00426     d = ((((q1[5].d*argsq + q1[4].d)*argsq + q1[3].d)*argsq +
00427         q1[2].d)*argsq + q1[1].d)*argsq + q1[0].d;
00428 
00429     return ( 0.5 + ((0.5 - arg) - arg*n/d) );
00430   }
00431   else if ( arg < 1.25 )
00432   {
00433     z = arg - 1.0;
00434 
00435     n = ((((((p2[8].d*z + p2[7].d)*z + p2[6].d)*z + p2[5].d)*z +
00436       p2[4].d)*z + p2[3].d)*z + p2[2].d)*z;
00437     d = ((((((q2[7].d*z + q2[6].d)*z + q2[5].d)*z + q2[4].d)*z +
00438       q2[3].d)*z + q2[2].d)*z + q2[1].d)*z + q2[0].d;
00439 
00440     return ( p2[0].d - n/d );
00441   }
00442   else if ( arg < 1.75 )
00443   {
00444     z = arg - 1.5;
00445 
00446     n = (((((p3[7].d*z + p3[6].d)*z +
00447       p3[5].d)*z + p3[4].d)*z + p3[3].d)*z + p3[2].d)*z;
00448     d = (((((q3[6].d*z + q3[5].d)*z +  q3[4].d)*z +
00449       q3[3].d)*z + q3[2].d)*z + q3[1].d)*z + q3[0].d;
00450 
00451     return ( p3[0].d - n/d );
00452   }
00453   else if ( arg < 2.0 )
00454   {
00455     z = arg - 1.875;
00456 
00457     n = ((((p4[6].d*z + p4[5].d)*z +
00458       p4[4].d)*z + p4[3].d)*z + p4[2].d)*z;
00459     d = ((((q4[5].d*z + q4[4].d)*z + q4[3].d)*z +
00460       q4[2].d)*z + q4[1].d)*z + q4[0].d;
00461 
00462     return ( p4[0].d - n/d );
00463   }
00464   else if ( arg <= Ulimit2.d )
00465   {
00466     zsq = 1.0/(arg*arg);
00467 
00468     n = (((((((p5[7].d*zsq + p5[6].d)*zsq + p5[5].d)*zsq +
00469       p5[4].d)*zsq + p5[3].d)*zsq + p5[2].d)*zsq + p5[1].d)*zsq +
00470       p5[0].d)*zsq;
00471 
00472     d = (((((((q5[8].d*zsq + q5[7].d)*zsq + q5[6].d)*zsq + q5[5].d)*zsq +
00473       q5[4].d)*zsq + q5[3].d)*zsq + q5[2].d)*zsq + q5[1].d)*zsq +
00474       q5[0].d;
00475 
00476     s = (float)arg;
00477 
00478     /* Now the result is exp(-arg*arg)*(1.0 + n/d)/(arg*sqrt(pi)) */
00479 
00480     if ( arg <= 26.0 )
00481     {
00482       /* To avoid loss of precision in computing exp(-arg*arg),
00483        * rewrite -arg*arg as -s*s + s*(s-arg) + arg*(s-arg)
00484        * Note that -s*s is exact.
00485        */
00486 
00487       /* Compute the second exp in the expression locally,
00488        * since the argument is small.
00489        */
00490 
00491       f1 = __exp(-s*s);
00492       f = f1 + f1*poly(s*(s-arg) + arg*(s-arg));
00493 
00494       return ( Rsqrtpi.d*(f + f*n/d)/arg );
00495     }
00496     else
00497     {
00498       /* Beyond 26.0, we have to worry about underflow in
00499        * computing exp(-arg*arg); instead we'll compute 
00500        * exp(-arg*arg/2.0) and multiply by it twice.
00501        * The result here may underflow, depending on 
00502        * whether the processor supports denormals or not.
00503        */
00504 
00505 
00506       s = (float)arg;
00507 
00508       f1 = __exp(-0.5*s*s);
00509       f = f1 + f1*poly(0.5*s*(s-arg) + 0.5*arg*(s-arg));
00510 
00511       result = Rsqrtpi.d*(f + f*n/d)/arg*f;
00512     }
00513   }
00514 
00515   if ( result == 0.0 )
00516 
00517 #ifdef _CALL_MATHERR
00518 
00519   {
00520     exstruct.type = UNDERFLOW;
00521     exstruct.name = "erfc";
00522     exstruct.arg1 = arg;
00523     exstruct.retval = 0.0;
00524 
00525     if ( matherr( &exstruct ) == 0 )
00526     {
00527       fprintf(stderr, "underflow error in erfc\n");
00528       SETERRNO(ERANGE);
00529     }
00530 
00531     return ( exstruct.retval );
00532   }
00533 #else
00534   {
00535     SETERRNO(ERANGE);
00536 
00537     return ( 0.0 );
00538   }
00539 #endif
00540 
00541   return ( result );
00542 }
00543 
00544 static
00545 double
00546 poly( double x )
00547 {
00548 double  result;
00549 
00550   /* returns expm1(x) for |x| <= 1.8e-4 */
00551 
00552   result = (((P[4].d*x + P[3].d)*x + P[2].d)*x + P[1].d)*x*x + x;
00553 
00554   return ( result );
00555 }
00556 
00557 #ifdef NO_LONG_DOUBLE
00558 
00559 #if defined(BUILD_OS_DARWIN) /* Mach-O doesn't support aliases */
00560 extern long double __erfl(long double);
00561 extern long double __erfcl(long double);
00562 long double erfl( long double x ) { 
00563   return ( (long double)__erf((double)x) );
00564 }
00565 
00566 long double erfcl( long double x ) {  
00567   return ( (long double)__erfc((double)x) );
00568 }
00569 #elif defined(__GNUC__)
00570 extern  long double  __erfl(long double);
00571 
00572 long double    erfl() __attribute__ ((weak, alias ("__erfl")));
00573 
00574 extern  long double  __erfcl(long double);
00575 
00576 long double    erfcl() __attribute__ ((weak, alias ("__erfcl")));
00577 
00578 #endif
00579 
00580 long double
00581 __erfl( long double x )
00582 { 
00583   return ( (long double)__erf((double)x) );
00584 }
00585 
00586 long double
00587 __erfcl( long double x )
00588 { 
00589   return ( (long double)__erfc((double)x) );
00590 }
00591 
00592 #endif
00593