erf.inl

//
//
/*-
 * Copyright (c) 1992, 1993
 *      The Regents of the University of California.  All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */
 
/* Modified Nov 30, 1992 P. McILROY:
 *      Replaced expansions for x >= 1.25 (error 1.7ulp vs ~6ulp)
 * Replaced even+odd with direct calculation for x < .84375,
 * to avoid destructive cancellation.
 *
 * Performance of erfc(x):
 * In 300000 trials in the range [.83, .84375] the
 * maximum observed error was 3.6ulp.
 *
 * In [.84735,1.25] the maximum observed error was <2.5ulp in
 * 100000 runs in the range [1.2, 1.25].
 *
 * In [1.25,26] (Not including subnormal results)
 * the error is < 1.7ulp.
 */
 
/* double erf(double x)
 * double erfc(double x)
 *                           x
 *                    2      |\
 *     erf(x)  =  ---------  | exp(-t*t)dt
 *                 sqrt(pi) \|
 *                           0
 *
 *     erfc(x) =  1-erf(x)
 *
 * Method:
 *      1. Reduce x to |x| by erf(-x) = -erf(x)
 *      2. For x in [0, 0.84375]
 *          erf(x)  = x + x*P(x^2)
 *          erfc(x) = 1 - erf(x)           if x<=0.25
 *                  = 0.5 + ((0.5-x)-x*P)  if x in [0.25,0.84375]
 *         where
 *                      2                2        4               20
 *              P =  P(x ) = (p0 + p1 * x + p2 * x + ... + p10 * x  )
 *         is an approximation to (erf(x)-x)/x with precision
 *
 *                                               -56.45
 *                      | P - (erf(x)-x)/x | <= 2
 *
 *
 *         Remark. The formula is derived by noting
 *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
 *         and that
 *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
 *         is close to one. The interval is chosen because the fixed
 *         point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
 *         near 0.6174), and by some experiment, 0.84375 is chosen to
 *         guarantee the error is less than one ulp for erf.
 *
 *      3. For x in [0.84375,1.25], let s = x - 1, and
 *         c = 0.84506291151 rounded to single (24 bits)
 *              erf(x)  = c  + P1(s)/Q1(s)
 *              erfc(x) = (1-c)  - P1(s)/Q1(s)
 *              |P1/Q1 - (erf(x)-c)| <= 2**-59.06
 *         Remark: here we use the taylor series expansion at x=1.
 *              erf(1+s) = erf(1) + s*Poly(s)
 *                       = 0.845.. + P1(s)/Q1(s)
 *         That is, we use rational approximation to approximate
 *                      erf(1+s) - (c = (single)0.84506291151)
 *         Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
 *         where
 *              P1(s) = degree 6 poly in s
 *              Q1(s) = degree 6 poly in s
 *
 *      4. For x in [1.25, 2]; [2, 4]
 *              erf(x)  = 1.0 - tiny
 *              erfc(x) = (1/x)exp(-x*x-(.5*log(pi) -.5z + R(z)/S(z))
 *
 *      Where z = 1/(x*x), R is degree 9, and S is degree 3;
 *
 *      5. For x in [4,28]
 *              erf(x)  = 1.0 - tiny
 *              erfc(x) = (1/x)exp(-x*x-(.5*log(pi)+eps + zP(z))
 *
 *      Where P is degree 14 polynomial in 1/(x*x).
 *
 *      Notes:
 *         Here 4 and 5 make use of the asymptotic series
 *                        exp(-x*x)
 *              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) );
 *                        x*sqrt(pi)
 *
 *              where for z = 1/(x*x)
 *              P(z) ~ z/2*(-1 + z*3/2*(1 + z*5/2*(-1 + z*7/2*(1 +...))))
 *
 *         Thus we use rational approximation to approximate
 *              erfc*x*exp(x*x) ~ 1/sqrt(pi);
 *
 *              The error bound for the target function, G(z) for
 *              the interval
 *              [4, 28]:
 *              |eps + 1/(z)P(z) - G(z)| < 2**(-56.61)
 *              for [2, 4]:
 *              |R(z)/S(z) - G(z)|       < 2**(-58.24)
 *              for [1.25, 2]:
 *              |R(z)/S(z) - G(z)|       < 2**(-58.12)
 *
 *      6. For inf > x >= 28
 *              erf(x)  = 1 - tiny  (raise inexact)
 *              erfc(x) = tiny*tiny (raise underflow)
 *
 *      7. Special cases:
 *              erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
 *              erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
 *              erfc/erf(NaN) is NaN
 */
 
/* Modified for STIR library Jun 05, 2005 Ch Tsoumpas and K Thielemans.
 
   WARNING:
   There are some assumptions here about IEEE arithmetic etc.
   On non-IEEE machines you're bound to have problems.
 
   We put in some work-arounds by using (hopefully) portable STIR_isnan,
   and STIR_finite().
*/
 
#ifdef _IEEE_LIBM
/*
 * redefining "___function" to "function" in _IEEE_LIBM mode
 */
#  include "ieee_libm.h"
#endif
 
#if defined(__vax__) || defined(tahoe)
// non-IEEE machines
#  define _IEEE 0
#  define TRUNC(x) (double)(x) = (float)(x)
#else
// assume everything uses IEEE floating point arithmetic
#  define _IEEE 1
// warning: strange definition of TRUNC that will go very weird on non-IEEE machines
#  define TRUNC(x) *(((int*)&x) + 1) &= 0xf8000000
#endif
 
#include <cmath>
#include "stir/numerics/ieeefp.h"
 
START_NAMESPACE_STIR
 
static const double tiny = 1e-300, half = 0.5, one = 1.0, two = 2.0, c = 8.45062911510467529297e-01, /* (float)0.84506291151 */
    /*
     * Coefficients for approximation to erf in [0,0.84375]
     */
    p0t8 = 1.02703333676410051049867154944018394163280, p0 = 1.283791670955125638123339436800229927041e-0001,
                    p1 = -3.761263890318340796574473028946097022260e-0001, p2 = 1.128379167093567004871858633779992337238e-0001,
                    p3 = -2.686617064084433642889526516177508374437e-0002, p4 = 5.223977576966219409445780927846432273191e-0003,
                    p5 = -8.548323822001639515038738961618255438422e-0004, p6 = 1.205520092530505090384383082516403772317e-0004,
                    p7 = -1.492214100762529635365672665955239554276e-0005, p8 = 1.640186161764254363152286358441771740838e-0006,
                    p9 = -1.571599331700515057841960987689515895479e-0007, p10 = 1.073087585213621540635426191486561494058e-0008;
/*
 * Coefficients for approximation to erf in [0.84375,1.25]
 */
static const double pa0 = -2.362118560752659485957248365514511540287e-0003, pa1 = 4.148561186837483359654781492060070469522e-0001,
                    pa2 = -3.722078760357013107593507594535478633044e-0001, pa3 = 3.183466199011617316853636418691420262160e-0001,
                    pa4 = -1.108946942823966771253985510891237782544e-0001, pa5 = 3.547830432561823343969797140537411825179e-0002,
                    pa6 = -2.166375594868790886906539848893221184820e-0003, qa1 = 1.064208804008442270765369280952419863524e-0001,
                    qa2 = 5.403979177021710663441167681878575087235e-0001, qa3 = 7.182865441419627066207655332170665812023e-0002,
                    qa4 = 1.261712198087616469108438860983447773726e-0001, qa5 = 1.363708391202905087876983523620537833157e-0002,
                    qa6 = 1.198449984679910764099772682882189711364e-0002;
/*
 * log(sqrt(pi)) for large x expansions.
 * The tail (lsqrtPI_lo) is included in the rational
 * approximations.
 */
static const double lsqrtPI_hi = .5723649429247000819387380943226;
/*
 * lsqrtPI_lo = .000000000000000005132975581353913;
 *
 * Coefficients for approximation to erfc in [2, 4]
 */
static const double rb0 = -1.5306508387410807582e-010, /* includes lsqrtPI_lo */
    rb1 = 2.15592846101742183841910806188e-008, rb2 = 6.24998557732436510470108714799e-001,
                    rb3 = 8.24849222231141787631258921465e+000, rb4 = 2.63974967372233173534823436057e+001,
                    rb5 = 9.86383092541570505318304640241e+000, rb6 = -7.28024154841991322228977878694e+000,
                    rb7 = 5.96303287280680116566600190708e+000, rb8 = -4.40070358507372993983608466806e+000,
                    rb9 = 2.39923700182518073731330332521e+000, rb10 = -6.89257464785841156285073338950e-001,
                    sb1 = 1.56641558965626774835300238919e+001, sb2 = 7.20522741000949622502957936376e+001,
                    sb3 = 9.60121069770492994166488642804e+001;
/*
 * Coefficients for approximation to erfc in [1.25, 2]
 */
static const double rc0 = -2.47925334685189288817e-007, /* includes lsqrtPI_lo */
    rc1 = 1.28735722546372485255126993930e-005, rc2 = 6.24664954087883916855616917019e-001,
                    rc3 = 4.69798884785807402408863708843e+000, rc4 = 7.61618295853929705430118701770e+000,
                    rc5 = 9.15640208659364240872946538730e-001, rc6 = -3.59753040425048631334448145935e-001,
                    rc7 = 1.42862267989304403403849619281e-001, rc8 = -4.74392758811439801958087514322e-002,
                    rc9 = 1.09964787987580810135757047874e-002, rc10 = -1.28856240494889325194638463046e-003,
                    sc1 = 9.97395106984001955652274773456e+000, sc2 = 2.80952153365721279953959310660e+001,
                    sc3 = 2.19826478142545234106819407316e+001;
/*
 * Coefficients for approximation to  erfc in [4,28]
 */
static const double rd0 = -2.1491361969012978677e-016, /* includes lsqrtPI_lo */
    rd1 = -4.99999999999640086151350330820e-001, rd2 = 6.24999999772906433825880867516e-001,
                    rd3 = -1.54166659428052432723177389562e+000, rd4 = 5.51561147405411844601985649206e+000,
                    rd5 = -2.55046307982949826964613748714e+001, rd6 = 1.43631424382843846387913799845e+002,
                    rd7 = -9.45789244999420134263345971704e+002, rd8 = 6.94834146607051206956384703517e+003,
                    rd9 = -5.27176414235983393155038356781e+004, rd10 = 3.68530281128672766499221324921e+005,
                    rd11 = -2.06466642800404317677021026611e+006, rd12 = 7.78293889471135381609201431274e+006,
                    rd13 = -1.42821001129434127360582351685e+007;
 
inline double
erf(double x)
{
  double R, S, P, Q, ax, s, y, z, r;
  if (!STIR_finite(x))
    { /* erf(nan)=nan */
      if (STIR_isnan(x))
        return (x);
      return (x > 0 ? one : -one); /* erf(+/-inf)= +/-1 */
    }
  if ((ax = x) < 0)
    ax = -ax;
  if (ax < .84375)
    {
      if (ax < 3.7e-09)
        {
          if (ax < 1.0e-308)
            return 0.125 * (8.0 * x + p0t8 * x); /*avoid underflow */
          return x + p0 * x;
        }
      y = x * x;
      r = y * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * (p7 + y * (p8 + y * (p9 + y * p10)))))))));
      return x + x * (p0 + r);
    }
  if (ax < 1.25)
    { /* 0.84375 <= |x| < 1.25 */
      s = std::fabs(x) - one;
      P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 + s * (pa5 + s * pa6)))));
      Q = one + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 + s * (qa5 + s * qa6)))));
      if (x >= 0)
        return (c + P / Q);
      else
        return (-c - P / Q);
    }
  if (ax >= 6.0)
    { /* inf>|x|>=6 */
      if (x >= 0.0)
        return (one - tiny);
      else
        return (tiny - one);
    }
  /* 1.25 <= |x| < 6 */
  z = -ax * ax;
  s = -one / z;
  if (ax < 2.0)
    {
      R = rc0
          + s * (rc1 + s * (rc2 + s * (rc3 + s * (rc4 + s * (rc5 + s * (rc6 + s * (rc7 + s * (rc8 + s * (rc9 + s * rc10)))))))));
      S = one + s * (sc1 + s * (sc2 + s * sc3));
    }
  else
    {
      R = rb0
          + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 + s * (rb5 + s * (rb6 + s * (rb7 + s * (rb8 + s * (rb9 + s * rb10)))))))));
      S = one + s * (sb1 + s * (sb2 + s * sb3));
    }
  y = (R / S - .5 * s) - lsqrtPI_hi;
  z += y;
  z = std::exp(z) / ax;
  if (x >= 0)
    return (one - z);
  else
    return (z - one);
}
 
inline double
erfc(double x)
{
  double R, S, P, Q, s, ax, y, z, r;
  if (!STIR_finite(x))
    {
      if (STIR_isnan(x)) /* erfc(NaN) = NaN */
        return (x);
      else if (x > 0) /* erfc(+-inf)=0,2 */
        return 0.0;
      else
        return 2.0;
    }
  if ((ax = x) < 0)
    ax = -ax;
  if (ax < .84375)
    {                                   /* |x|<0.84375 */
      if (ax < 1.38777878078144568e-17) /* |x|<2**-56 */
        return one - x;
      y = x * x;
      r = y * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * (p7 + y * (p8 + y * (p9 + y * p10)))))))));
      if (ax < .0625)
        { /* |x|<2**-4 */
          return (one - (x + x * (p0 + r)));
        }
      else
        {
          r = x * (p0 + r);
          r += (x - half);
          return (half - r);
        }
    }
  if (ax < 1.25)
    { /* 0.84375 <= |x| < 1.25 */
      s = ax - one;
      P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 + s * (pa5 + s * pa6)))));
      Q = one + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 + s * (qa5 + s * qa6)))));
      if (x >= 0)
        {
          z = one - c;
          return z - P / Q;
        }
      else
        {
          z = c + P / Q;
          return one + z;
        }
    }
  if (ax >= 28)
    { /* Out of range */
      if (x > 0)
        return (tiny * tiny);
      else
        return (two - tiny);
    }
  z = ax;
  TRUNC(z);
  y = z - ax;
  y *= (ax + z);
  z *= -z;            /* Here z + y = -x^2 */
  s = one / (-z - y); /* 1/(x*x) */
  if (ax >= 4)
    { /* 6 <= ax */
      R = s
          * (rd1
             + s
                   * (rd2
                      + s
                            * (rd3
                               + s
                                     * (rd4
                                        + s
                                              * (rd5
                                                 + s
                                                       * (rd6
                                                          + s
                                                                * (rd7
                                                                   + s
                                                                         * (rd8
                                                                            + s
                                                                                  * (rd9
                                                                                     + s
                                                                                           * (rd10
                                                                                              + s
                                                                                                    * (rd11
                                                                                                       + s
                                                                                                             * (rd12
                                                                                                                + s * rd13))))))))))));
      y += rd0;
    }
  else if (ax >= 2)
    {
      R = rb0
          + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 + s * (rb5 + s * (rb6 + s * (rb7 + s * (rb8 + s * (rb9 + s * rb10)))))))));
      S = one + s * (sb1 + s * (sb2 + s * sb3));
      y += R / S;
      R = -.5 * s;
    }
  else
    {
      R = rc0
          + s * (rc1 + s * (rc2 + s * (rc3 + s * (rc4 + s * (rc5 + s * (rc6 + s * (rc7 + s * (rc8 + s * (rc9 + s * rc10)))))))));
      S = one + s * (sc1 + s * (sc2 + s * sc3));
      y += R / S;
      R = -.5 * s;
    }
  /* return std::exp(-x^2 - lsqrtPI_hi + R + y)/x;      */
  //    s = ((R + y) - lsqrtPI_hi) + z;
  //    y = (((z-s) - lsqrtPI_hi) + R) + y;
  //    r = std::exp(-x + s/x);
  // r = __exp__D(s, y)/x;
  if (x > 0)
    return (std::exp(-x * x - lsqrtPI_hi + R + y)) / x;
  // return r;
  else
    return 2. - (std::exp(-x * x - lsqrtPI_hi + R + y)) / x;
  //    return two-r;
}
 
#undef TRUNC
#undef _IEEE
END_NAMESPACE_STIR