/*							i0e.c
 *
 *	Modified Bessel function of order zero,
 *	exponentially scaled
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, i0e();
 *
 * y = i0e( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns exponentially scaled modified Bessel function
 * of order zero of the argument.
 *
 * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0,30        30000       5.4e-16     1.2e-16
 * See i0().
 *
 */

/*							i0.c		*/
#include <stdint.h>
#include <math.h>

/* Chebyshev coefficients for exp(-x) I0(x)
 * in the interval [0,8].
 *
 * lim(x->0){ exp(-x) I0(x) } = 1.
 */
static const double P0[30] = {
   0.20700192122398670,
  -1.13004326886205e-01,
   9.45156285399510e-02,
  -8.90434413607596e-02,
   8.73434161755231e-02,
  -8.49633809161631e-02,
   7.92752259071635e-02,
  -6.95411956171938e-02,
   5.67770320199169e-02,
  -4.29859494547961e-02,
   3.01779120229784e-02,
  -1.96826546905923e-02,
   1.19600779070262e-02,
  -6.79239987307317e-03,
   3.61712644889704e-03,
  -1.81190575768570e-03,
   8.56172218291950e-04,
  -3.82864168004715e-04,
   1.62548478164061e-04,
  -6.45212783972906e-05,
   2.49162142005162e-05,
  -1.02200523353782e-05,
   3.37239561792087e-06,
  -1.14523064092414e-06,
   2.24527757375167e-07,
  -9.48773984451770e-08,
   2.62653631519701e-08,
  -7.72309455684633e-09,
   4.47051672751174e-09,
  -1.18523424826847e-09
};

static const double A[30] =
{
 6.76795274409476084995E-01,
-3.04682672343198398683E-01,
 1.71620901522208775349E-01,
-9.49010970480476444210E-02,
 4.93052842396707084878E-02,
-2.37374148058994688156E-02,
 1.05464603945949983183E-02,
-4.32430999505057594430E-03,
 1.63947561694133579842E-03,
-5.76375574538582365885E-04,
 1.88502885095841655729E-04,
-5.75419501008210370398E-05,
 1.64484480707288970893E-05,
-4.41673835845875056359E-06,
 1.11738753912010371815E-06,
-2.67079385394061173391E-07,
 6.04699502254191894932E-08,
-1.30002500998624804212E-08,
 2.65982372468238665035E-09,
-5.18979560163526290666E-10,
 9.67580903537323691224E-11,
-1.72682629144155570723E-11,
 2.95505266312963983461E-12,
-4.85644678311192946090E-13,
 7.67618549860493561688E-14,
-1.16853328779934516808E-14,
 1.71539128555513303061E-15,
-2.43127984654795469359E-16,
 3.33079451882223809783E-17,
-4.41534164647933937950E-18
};

/* Chebyshev coefficients for exp(-x) sqrt(x) I0(x)
 * in the inverted interval [8,infinity].
 *
 * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
 */
/* Polynomial coefficients for exp(-x) sqrt(x) I0(x)
 * in the inverted interval [7.75,Inf].
 * (Only used over [8,Inf]).
 *
 * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
 *
 * Coefficients and technique from:
 * https://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision/
 */
static const double P2[23] =
{
 3.9894228040143265335649948e-01,
 4.9867785050353992900698488e-02,
 2.8050628884163787533196746e-02,
 2.9219501690198775910219311e-02,
 4.4718622769244715693031735e-02,
 9.4085204199017869159183831e-02,
-1.0699095472110916094973951e-01,
 2.2725199603010833194037016e+01,
-1.0026890180180668595066918e+03,
 3.1275740782277570164423916e+04,
-5.9355022509673600842060002e+05,
 2.6092888649549172879282592e+06,
 2.3518420447411254516178388e+08,
-8.9270060370015930749184222e+09,
 1.8592340458074104721496236e+11,
-2.6632742974569782078420204e+12,
 2.7752144774934763122129261e+13,
-2.1323049786724612220362154e+14,
 1.1989242681178569338129044e+15,
-4.8049082153027457378879746e+15,
 1.3012646806421079076251950e+16,
-2.1363029690365351606041265e+16,
 1.6069467093441596329340754e+16
};

static inline double
cheb_eval(const double *cs, int order, double y)
{
  double d = 0.0, dd = 0.0, y2 = y+y;
  int j;
  order--;
  for(j = order; j>=1; j--) {
    double temp = d;
    d = y2*d - dd + cs[j];
    dd = temp;
  }
  d = y*d + (0.5*cs[0] - dd);
  return d;
}

double I0_scaled(double x)
{
    double val = 0.0;
    int i;
    x = fabs(x);
    if(x <= 8.0) {
	    double z = (0.25*x - 1.0);
#if 0
        val = P0[29];
        for(i=30; i>0; i--) val = P0[i-1] + z * val;
#else
        val = cheb_eval(A, 29, z);
#endif
	} else {
        //double z = (32.0/x - 2.0);
        double z = 1.0/x;
        val = P2[29];
        for(i=29; i>0; i--) val = P2[i-1] + z * val;
        val = (val / sqrt(x));
    }
    return val;
}

#ifdef TEST

#include <stdio.h>

static const double I0tests[21][2] = {
    {   1e-10, 0.99999999990000000001 },
    {     0.1, 0.90710092578230109640 },
    {       2, 0.30850832255367103953 },
    {   100.0, 0.03994437929909668265 },
    { 65536.0, 0.0015583712551952223537 },
    { 0.177219114266335964202880859375e-2, 0.1000000785165515654790765567976560128368e1 }, 
    { 0.22177286446094512939453125e-2, 0.1000001229580463247166647265657459686513e1 }, 
    { 0.7444499991834163665771484375e-2, 0.1000013855193023439561602999309818071557e1 }, 
    { 0.1433600485324859619140625e-1, 0.100005138091877460790808253311381371617e1 }, 
    { 0.1760916970670223236083984375e-1, 0.1000077522216818092267597684468025144708e1 }, 
    { 0.6152711808681488037109375e-1, 0.1000946620505179168127565717529250425225e1 }, 
    { 0.11958599090576171875e0, 0.1003578399092795290959506558979668000986e1 }, 
    { 0.15262925624847412109375e0, 0.1005832407473077490942708958509537944332e1 }, 
    { 0.408089816570281982421875e0, 0.1042069688780878266908906813963061201236e1 }, 
    { 0.6540834903717041015625e0, 0.1109850431257999241542412962070021625239e1 }, 
    { 0.1097540378570556640625e1, 0.1324594459649893673687924138468203300537e1 }, 
    { 0.30944411754608154296875e1, 0.5270503008520823311035815313268274322688e1 }, 
    { 0.51139926910400390625e1, 0.3016460113324413124385977969424656719006e2 }, 
    { 0.95070552825927734375e1, 0.1765221230864037183559483071058505259281e4 }, 
    { 0.24750102996826171875e2, 0.4520582716299093122418957895745838047526e10 }, 
    { 0.637722015380859375e2, 0.248523599457696979135532431768725854655e27 }
};

static double sf_frac_diff(double x1, double x2)
{
  if (x2 == 0.0) {
    return fabs(x1);
  } else if(x1 <= 1e308 && x2 <= 1e308 && (x1 + x2 != 0.0)) {
    return fabs((x1-x2)/(x1+x2));
  } else {
    return 1.0;
  }
}

static __attribute__((always_inline)) inline void disable_denormals(void)
{
    unsigned int a = 0x9f80;
    __asm__ __volatile__ ("ldmxcsr %0" : "+m" (a));
}

int main(void)
{
  double diff;
  unsigned int i;

  disable_denormals();
  setvbuf(stdout, NULL, _IONBF, 0);

  for (i = 0; i < 21; i++) {
    double x = I0tests[i][0];
    double Jnx_ans = I0tests[i][1];
    double Jnx = I0_scaled(x);
    if (i >= 5) Jnx_ans *= exp(-x);
    if ((diff = sf_frac_diff(Jnx, Jnx_ans)) > 0x1p-51)
      printf("%u: I0(%#.7g): %#.10g != %#.10g (diff: %#.10g / %a)\n", i, x, Jnx, Jnx_ans, diff, diff);
  }

  return 0;
}

#endif