Print this page
5261 libm should stop using synonyms.h

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/LD/erfl.c
          +++ new/usr/src/lib/libm/common/LD/erfl.c
↓ open down ↓ 24 lines elided ↑ open up ↑
  25   25  /*
  26   26   * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27   27   * Use is subject to license terms.
  28   28   */
  29   29  
  30   30  /* long double function erf,erfc (long double x)
  31   31   * K.C. Ng, September, 1989.
  32   32   *                           x
  33   33   *                    2      |\
  34   34   *     erf(x)  =  ---------  | exp(-t*t)dt
  35      - *                 sqrt(pi) \| 
       35 + *                 sqrt(pi) \|
  36   36   *                           0
  37   37   *
  38   38   *     erfc(x) =  1-erf(x)
  39   39   *
  40   40   * method:
  41   41   *      Since erf(-x) = -erf(x), we assume x>=0.
  42   42   *      For x near 0, we have the expansion
  43      - * 
       43 + *
  44   44   *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....).
  45   45   *
  46   46   *      Since 2/sqrt(pi) = 1.128379167095512573896158903121545171688,
  47   47   *      we use x + x*P(x^2) to approximate erf(x). This formula will
  48   48   *      guarantee the error less than one ulp where x is not too far
  49   49   *      away from 0. We note that erf(x)=x at x = 0.6174...... After
  50      - *      some experiment, we choose the following approximation on 
       50 + *      some experiment, we choose the following approximation on
  51   51   *      interval [0,0.84375].
  52   52   *
  53   53   *      For x in [0,0.84375]
  54   54   *                 2                2        4               40
  55   55   *         P =  P(x ) = (p0 + p1 * x + p2 * x + ... + p20 * x  )
  56   56   *
  57   57   *         erf(x)  = x + x*P
  58   58   *         erfc(x) = 1 - erf(x)           if x<=0.25
  59   59   *                 = 0.5 + ((0.5-x)-x*P)  if x in [0.25,0.84375]
  60   60   *      precision: |P(x^2)-(erf(x)-x)/x| <= 2**-122.50
  61   61   *
  62      - *      For x in [0.84375,1.25], let s = x - 1, and 
       62 + *      For x in [0.84375,1.25], let s = x - 1, and
  63   63   *      c = 0.84506291151 rounded to single (24 bits)
  64   64   *         erf(x)  = c  + P1(s)/Q1(s)
  65      - *         erfc(x) = (1-c)  - P1(s)/Q1(s) 
       65 + *         erfc(x) = (1-c)  - P1(s)/Q1(s)
  66   66   *      precision: |P1/Q1 - (erf(x)-c)| <= 2**-118.41
  67      - *      
  68   67   *
  69      - *      For x in [1.25,1.75], let s = x - 1.5, and 
       68 + *
       69 + *      For x in [1.25,1.75], let s = x - 1.5, and
  70   70   *      c = 0.95478588343 rounded to single (24 bits)
  71   71   *         erf(x)  = c  + P2(s)/Q2(s)
  72      - *         erfc(x) = (1-c)  - P2(s)/Q2(s) 
       72 + *         erfc(x) = (1-c)  - P2(s)/Q2(s)
  73   73   *      precision: |P1/Q1 - (erf(x)-c)| <= 2**-123.83
  74   74   *
  75   75   *
  76   76   *      For x in [1.75,16/3]
  77   77   *         erfc(x) = exp(-x*x)*(1/x)*R1(1/x)/S1(1/x)
  78   78   *         erf(x)  = 1 - erfc(x)
  79   79   *      precision: absolute error of R1/S1 is bounded by 2**-124.03
  80   80   *
  81   81   *      For x in [16/3,107]
  82   82   *         erfc(x) = exp(-x*x)*(1/x)*R2(1/x)/S2(1/x)
  83   83   *         erf(x)  = 1 - erfc(x) (if x>=9 simple return erf(x)=1 with inexact)
  84   84   *      precision: absolute error of R2/S2 is bounded by 2**-120.07
  85   85   *
  86   86   *      Else if inf > x >= 107
  87      - *         erf(x)  = 1 with inexact 
       87 + *         erf(x)  = 1 with inexact
  88   88   *         erfc(x) = 0 with underflow
  89      - *      
       89 + *
  90   90   *      Special case:
  91   91   *         erf(inf)  = 1
  92   92   *         erfc(inf) = 0
  93   93   */
  94   94  
  95      -#pragma weak erfl = __erfl
  96      -#pragma weak erfcl = __erfcl
       95 +#pragma weak __erfl = erfl
       96 +#pragma weak __erfcl = erfcl
  97   97  
  98   98  #include "libm.h"
  99   99  #include "longdouble.h"
 100  100  
 101  101  static long double
 102  102  tiny        = 1e-40L,
 103  103  nearunfl    = 1e-4000L,
 104  104  half        = 0.5L,
 105  105  one         = 1.0L,
 106  106  onehalf     = 1.5L,
↓ open down ↓ 18 lines elided ↑ open up ↑
 125  125    -6.711366846653439036162105104991433380926e-0012L,
 126  126     4.463224090341893165100275380693843116240e-0013L,
 127  127    -2.783513452582658245422635662559779162312e-0014L,
 128  128     1.634227412586960195251346878863754661546e-0015L,
 129  129    -9.060782672889577722765711455623117802795e-0017L,
 130  130     4.741341801266246873412159213893613602354e-0018L,
 131  131    -2.272417596497826188374846636534317381203e-0019L,
 132  132     8.069088733716068462496835658928566920933e-0021L,
 133  133  };
 134  134  
 135      -/* 
      135 +/*
 136  136   * Rational erf(x) = ((float)0.84506291151) + P1(x-1)/Q1(x-1) on [0.84375,1.25]
 137  137   */
 138  138  static long double C1   = (long double)((float)0.84506291151);
 139  139  static long double P1[] = {     /*  12 top coeffs */
 140  140    -2.362118560752659955654364917390741930316e-0003L,
 141  141     4.129623379624420034078926610650759979146e-0001L,
 142  142    -3.973857505403547283109417923182669976904e-0002L,
 143  143     4.357503184084022439763567513078036755183e-0002L,
 144  144     8.015593623388421371247676683754171456950e-0002L,
 145  145    -1.034459310403352486685467221776778474602e-0002L,
↓ open down ↓ 11 lines elided ↑ open up ↑
 157  157     1.702332233546316765818144723063881095577e-0001L,
 158  158     7.498098377690553934266423088708614219356e-0002L,
 159  159     2.050154396918178697056927234366372760310e-0002L,
 160  160     7.012988534031999899054782333851905939379e-0003L,
 161  161     1.149904787014400354649843451234570731076e-0003L,
 162  162     3.185620255011299476196039491205159718620e-0004L,
 163  163     1.273405072153008775426376193374105840517e-0005L,
 164  164     4.753866999959432971956781228148402971454e-0006L,
 165  165    -1.002287602111660026053981728549540200683e-0006L,
 166  166  };
 167      -/* 
      167 +/*
 168  168   * Rational erf(x) = ((float)0.95478588343) + P2(x-1.5)/Q2(x-1.5)
 169  169   * on [1.25,1.75]
 170  170   */
 171  171  static long double C2   = (long double)((float)0.95478588343);
 172  172  static long double P2[] = {     /*  12 top coeffs */
 173  173     1.131926304864446730135126164594785863512e-0002L,
 174  174     1.273617996967754151544330055186210322832e-0001L,
 175  175    -8.169980734667512519897816907190281143423e-0002L,
 176  176     9.512267486090321197833634271787944271746e-0002L,
 177  177    -2.394251569804872160005274999735914368170e-0002L,
↓ open down ↓ 89 lines elided ↑ open up ↑
 267  267     3.792278350536686111444869752624492443659e+0005L,
 268  268     1.257750606950115799965366001773094058720e+0005L,
 269  269     3.410830600242369370645608634643620355058e+0004L,
 270  270     7.513984469742343134851326863175067271240e+0003L,
 271  271     1.313296320593190002554779998138695507840e+0003L,
 272  272     1.773972700887629157006326333696896516769e+0002L,
 273  273     1.670876451822586800422009013880457094162e+0001L,
 274  274     1.000L,
 275  275  };
 276  276  
 277      -long double erfl(x) 
      277 +long double erfl(x)
 278  278  long double x;
 279  279  {
 280  280          long double erfcl(long double),s,y,t;
 281  281  
 282  282          if (!finitel(x)) {
 283  283              if (x != x) return x+x;     /* NaN */
 284  284              return copysignl(one,x);    /* return +-1.0 is x=Inf */
 285  285          }
 286  286  
 287  287          y = fabsl(x);
↓ open down ↓ 4 lines elided ↑ open up ↑
 292  292              return  x+x*t;
 293  293          }
 294  294          if (y<=1.25L) {
 295  295              s = y-one;
 296  296              t = C1+__poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1));
 297  297              return (signbitl(x))? -t: t;
 298  298          } else if (y<=1.75L) {
 299  299              s = y-onehalf;
 300  300              t = C2+__poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2));
 301  301              return (signbitl(x))? -t: t;
 302      -        } 
      302 +        }
 303  303          if (y<=9.0L) t = erfcl(y); else t = tiny;
 304  304          return (signbitl(x))? t-one: one-t;
 305  305  }
 306  306  
 307  307  long double erfcl(x)
 308  308  long double x;
 309  309  {
 310  310          long double erfl(long double),s,y,t;
 311  311  
 312  312          if (!finitel(x)) {
↓ open down ↓ 11 lines elided ↑ open up ↑
 324  324              return  half+t;
 325  325          }
 326  326          if (x<=1.25L) {
 327  327              s = x-one;
 328  328              t = one-C1;
 329  329              return t - __poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1));
 330  330          } else if (x<=1.75L) {
 331  331              s = x-onehalf;
 332  332              t = one-C2;
 333  333              return t - __poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2));
 334      -        } 
      334 +        }
 335  335          if (x>=107.0L) return nearunfl*nearunfl;                /* underflow */
 336  336          else if (x >= L16_3) {
 337  337              y = __poly_libmq(x,15,R2);
 338  338              t = y/__poly_libmq(x,16,S2);
 339  339          } else {
 340  340              y = __poly_libmq(x,14,R1);
 341  341              t = y/__poly_libmq(x,15,S1);
 342  342          }
 343  343          /* see comment in ../Q/erfl.c */
 344  344          y = x;
 345  345          *(int*)&y = 0;
 346  346          t *= expl(-y*y)*expl(-(x-y)*(x+y));
 347  347          return  t;
 348  348  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX