libF77/erf.c

   1 /*      $Id: erf.c,v 1.2 2008/02/26 19:54:41 ragge Exp $        */
   2 /*
   3  * Copyright(C) Caldera International Inc. 2001-2002. All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  *
   9  * Redistributions of source code and documentation must retain the above
  10  * copyright notice, this list of conditions and the following disclaimer.
  11  * Redistributions in binary form must reproduce the above copyright
  12  * notice, this list of conditionsand the following disclaimer in the
  13  * documentation and/or other materials provided with the distribution.
  14  * All advertising materials mentioning features or use of this software
  15  * must display the following acknowledgement:
  16  *      This product includes software developed or owned by Caldera
  17  *      International, Inc.
  18  * Neither the name of Caldera International, Inc. nor the names of other
  19  * contributors may be used to endorse or promote products derived from
  20  * this software without specific prior written permission.
  21  *
  22  * USE OF THE SOFTWARE PROVIDED FOR UNDER THIS LICENSE BY CALDERA
  23  * INTERNATIONAL, INC. AND CONTRIBUTORS ``AS IS'' AND ANY EXPRESS OR
  24  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  25  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  26  * DISCLAIMED.  IN NO EVENT SHALL CALDERA INTERNATIONAL, INC. BE LIABLE
  27  * FOR ANY DIRECT, INDIRECT INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  28  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  29  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  30  * HOWEVER CAUSED AND ON ANY THEORY OFLIABILITY, WHETHER IN CONTRACT,
  31  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
  32  * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  33  * POSSIBILITY OF SUCH DAMAGE.
  34  */
  35 /*
  36         C program for floating point error function
  37
  38         erf(x) returns the error function of its argument
  39         erfc(x) returns 1.0-erf(x)
  40
  41         erf(x) is defined by
  42         ${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$
  43
  44         the entry for erfc is provided because of the
  45         extreme loss of relative accuracy if erf(x) is
  46         called for large x and the result subtracted
  47         from 1. (e.g. for x= 10, 12 places are lost).
  48
  49         There are no error returns.
  50
  51         Calls exp.
  52
  53         Coefficients for large x are #5667 from Hart & Cheney (18.72D).
  54 */
  55
  56 #define M 7
  57 #define N 9
  58 int errno;
  59 static double torp = 1.1283791670955125738961589031;
  60 static double p1[] = {
  61         0.804373630960840172832162e5,
  62         0.740407142710151470082064e4,
  63         0.301782788536507577809226e4,
  64         0.380140318123903008244444e2,
  65         0.143383842191748205576712e2,
  66         -.288805137207594084924010e0,
  67         0.007547728033418631287834e0,
  68 };
  69 static double q1[]  = {
  70         0.804373630960840172826266e5,
  71         0.342165257924628539769006e5,
  72         0.637960017324428279487120e4,
  73         0.658070155459240506326937e3,
  74         0.380190713951939403753468e2,
  75         0.100000000000000000000000e1,
  76         0.0,
  77 };
  78 static double p2[]  = {
  79         0.18263348842295112592168999e4,
  80         0.28980293292167655611275846e4,
  81         0.2320439590251635247384768711e4,
  82         0.1143262070703886173606073338e4,
  83         0.3685196154710010637133875746e3,
  84         0.7708161730368428609781633646e2,
  85         0.9675807882987265400604202961e1,
  86         0.5641877825507397413087057563e0,
  87         0.0,
  88 };
  89 static double q2[]  = {
  90         0.18263348842295112595576438e4,
  91         0.495882756472114071495438422e4,
  92         0.60895424232724435504633068e4,
  93         0.4429612803883682726711528526e4,
  94         0.2094384367789539593790281779e4,
  95         0.6617361207107653469211984771e3,
  96         0.1371255960500622202878443578e3,
  97         0.1714980943627607849376131193e2,
  98         1.0,
  99 };
 100
 101 double
 102 erf(arg) double arg;{
 103         double erfc();
 104         int sign;
 105         double argsq;
 106         double d, n;
 107         int i;
 108
 109         errno = 0;
 110         sign = 1;
 111         if(arg < 0.){
 112                 arg = -arg;
 113                 sign = -1;
 114         }
 115         if(arg < 0.5){
 116                 argsq = arg*arg;
 117                 for(n=0,d=0,i=M-1; i>=0; i--){
 118                         n = n*argsq + p1[i];
 119                         d = d*argsq + q1[i];
 120                 }
 121                 return(sign*torp*arg*n/d);
 122         }
 123         if(arg >= 10.)
 124                 return(sign*1.);
 125         return(sign*(1. - erfc(arg)));
 126 }
 127
 128 double
 129 erfc(arg) double arg;{
 130         double erf();
 131         double exp();
 132         double n, d;
 133         int i;
 134
 135         errno = 0;
 136         if(arg < 0.)
 137                 return(2. - erfc(-arg));
 138 /*
 139         if(arg < 0.5)
 140                 return(1. - erf(arg));
 141 */
 142         if(arg >= 10.)
 143                 return(0.);
 144
 145         for(n=0,d=0,i=N-1; i>=0; i--){
 146                 n = n*arg + p2[i];
 147                 d = d*arg + q2[i];
 148         }
 149         return(exp(-arg*arg)*n/d);
 150 }