/**
 * A math namespace - all functions can be exported from base/math.mjs.
 * Also all these functions can be used with TFormula calcualtions
 * @namespace Math
 */

/* eslint-disable no-loss-of-precision */
/* eslint-disable space-in-parens */
/* eslint-disable curly */
/* eslint-disable operator-linebreak */
/* eslint-disable no-floating-decimal */
/* eslint-disable brace-style */
/* eslint-disable comma-spacing */

// this can be improved later

/* eslint-disable no-unreachable-loop */
/* eslint-disable eqeqeq */

const kMACHEP = 1.11022302462515654042363166809e-16,
      kMINLOG = -708.396418532264078748994506896,
      kMAXLOG = 709.782712893383973096206318587,
      kMAXSTIR = 108.116855767857671821730036754,
      kBig = 4.503599627370496e15,
      kBiginv = 2.22044604925031308085e-16,
      kSqrt2 = 1.41421356237309515,
      M_PI = 3.14159265358979323846264338328;

/** @summary Polynomialeval function
  * @desc calculates a value of a polynomial of the form:
  * a[0]x^N+a[1]x^(N-1) + ... + a[N]
  * @memberof Math */
function Polynomialeval(x, a, N) {
   if (!N) return a[0];

   let pom = a[0];
   for (let i = 1; i <= N; ++i)
      pom = pom *x + a[i];
   return pom;
}

/** @summary Polynomial1eval function
  * @desc calculates a value of a polynomial of the form:
  * x^N+a[0]x^(N-1) + ... + a[N-1]
  * @memberof Math */
function Polynomial1eval(x, a, N) {
   if (!N) return a[0];

   let pom = x + a[0];
   for (let i = 1; i < N; ++i)
      pom = pom *x + a[i];
   return pom;
}

/** @summary lgam function, logarithm from gamma
  * @memberof Math */
function lgam(x) {
   let p, q, u, w, z;
   const kMAXLGM = 2.556348e305,
         LS2PI = 0.91893853320467274178,
   A = [
      8.11614167470508450300E-4,
      -5.95061904284301438324E-4,
      7.93650340457716943945E-4,
      -2.77777777730099687205E-3,
      8.33333333333331927722E-2
   ], B = [
      -1.37825152569120859100E3,
      -3.88016315134637840924E4,
      -3.31612992738871184744E5,
      -1.16237097492762307383E6,
      -1.72173700820839662146E6,
      -8.53555664245765465627E5
   ], C = [
   /* 1.00000000000000000000E0, */
      -3.51815701436523470549E2,
      -1.70642106651881159223E4,
      -2.20528590553854454839E5,
      -1.13933444367982507207E6,
      -2.53252307177582951285E6,
      -2.01889141433532773231E6
   ];

   if ((x >= Number.MAX_VALUE) || (x == Number.POSITIVE_INFINITY))
      return Number.POSITIVE_INFINITY;

   if ( x < -34.0 ) {
      q = -x;
      w = lgam(q);
      p = Math.floor(q);
      if ( p === q ) // _unur_FP_same(p,q)
         return Number.POSITIVE_INFINITY;
      z = q - p;
      if ( z > 0.5 ) {
         p += 1.0;
         z = p - q;
      }
      z = q * Math.sin( Math.PI * z );
      if ( z < 1e-300 )
         return Number.POSITIVE_INFINITY;
      z = Math.log(Math.PI) - Math.log( z ) - w;
      return z;
   }
   if ( x < 13.0 ) {
      z = 1.0;
      p = 0.0;
      u = x;
      while ( u >= 3.0 ) {
         p -= 1.0;
         u = x + p;
         z *= u;
      }
      while ( u < 2.0 ) {
         if ( u < 1e-300 )
            return Number.POSITIVE_INFINITY;
         z /= u;
         p += 1.0;
         u = x + p;
      }
      if ( z < 0.0 ) {
         z = -z;
      }
      if ( u === 2.0 )
         return Math.log(z);
      p -= 2.0;
      x = x + p;
      p = x * Polynomialeval(x, B, 5 ) / Polynomial1eval( x, C, 6);
      return Math.log(z) + p;
   }
   if ( x > kMAXLGM )
      return Number.POSITIVE_INFINITY;

   q = ( x - 0.5 ) * Math.log(x) - x + LS2PI;
   if ( x > 1.0e8 )
      return q;

   p = 1.0/(x*x);
   if ( x >= 1000.0 )
      q += ((7.9365079365079365079365e-4 * p
            - 2.7777777777777777777778e-3) *p
            + 0.0833333333333333333333) / x;
   else
      q += Polynomialeval( p, A, 4 ) / x;
   return q;
}

/** @summary Stirling formula for the gamma function
  * @memberof Math */
function stirf(x) {
   let y, w, v;

   const STIR = [
      7.87311395793093628397E-4,
      -2.29549961613378126380E-4,
      -2.68132617805781232825E-3,
      3.47222221605458667310E-3,
      8.33333333333482257126E-2
   ], SQTPI = Math.sqrt(2*Math.PI);

   w = 1.0/x;
   w = 1.0 + w * Polynomialeval( w, STIR, 4 );
   y = Math.exp(x);

/*   #define kMAXSTIR kMAXLOG/log(kMAXLOG)  */

   if ( x > kMAXSTIR )
   { /* Avoid overflow in pow() */
      v = Math.pow( x, 0.5 * x - 0.25 );
      y = v * (v / y);
   }
   else
   {
      y = Math.pow( x, x - 0.5 ) / y;
   }
   y = SQTPI * y * w;
   return y;
}

/** @summary complementary error function
  * @memberof Math */
function erfc(a) {
   const erfP = [
      2.46196981473530512524E-10,
      5.64189564831068821977E-1,
      7.46321056442269912687E0,
      4.86371970985681366614E1,
      1.96520832956077098242E2,
      5.26445194995477358631E2,
      9.34528527171957607540E2,
      1.02755188689515710272E3,
      5.57535335369399327526E2
   ], erfQ = [
      1.32281951154744992508E1,
      8.67072140885989742329E1,
      3.54937778887819891062E2,
      9.75708501743205489753E2,
      1.82390916687909736289E3,
      2.24633760818710981792E3,
      1.65666309194161350182E3,
      5.57535340817727675546E2
   ], erfR = [
      5.64189583547755073984E-1,
      1.27536670759978104416E0,
      5.01905042251180477414E0,
      6.16021097993053585195E0,
      7.40974269950448939160E0,
      2.97886665372100240670E0
   ], erfS = [
      2.26052863220117276590E0,
      9.39603524938001434673E0,
      1.20489539808096656605E1,
      1.70814450747565897222E1,
      9.60896809063285878198E0,
      3.36907645100081516050E0
   ];

   let p,q,x,y,z;

   if ( a < 0.0 )
      x = -a;
   else
      x = a;

   if ( x < 1.0 )
      return 1.0 - erf(a);

   z = -a * a;

   if (z < -kMAXLOG)
      return (a < 0) ? 2.0 : 0.0;

   z = Math.exp(z);

   if ( x < 8.0 ) {
      p = Polynomialeval( x, erfP, 8 );
      q = Polynomial1eval( x, erfQ, 8 );
   } else {
      p = Polynomialeval( x, erfR, 5 );
      q = Polynomial1eval( x, erfS, 6 );
   }
   y = (z * p)/q;

   if (a < 0)
      y = 2.0 - y;

   if (y == 0)
      return (a < 0) ? 2.0 : 0.0;

   return y;
}

/** @summary error function
  * @memberof Math */
function erf(x) {
   if (Math.abs(x) > 1.0)
      return 1.0 - erfc(x);

   const erfT = [
      9.60497373987051638749E0,
      9.00260197203842689217E1,
      2.23200534594684319226E3,
      7.00332514112805075473E3,
      5.55923013010394962768E4
   ], erfU = [
      3.35617141647503099647E1,
      5.21357949780152679795E2,
      4.59432382970980127987E3,
      2.26290000613890934246E4,
      4.92673942608635921086E4
   ],

   z = x * x;

   return x * Polynomialeval(z, erfT, 4) / Polynomial1eval(z, erfU, 5);
}

/** @summary lognormal_cdf_c function
  * @memberof Math */
function lognormal_cdf_c(x, m, s, x0) {
   if (x0 === undefined) x0 = 0;
   const z = (Math.log((x-x0))-m)/(s*kSqrt2);
   if (z > 1.) return 0.5*erfc(z);
   else return 0.5*(1.0 - erf(z));
}

/** @summary lognormal_cdf_c function
  * @memberof Math */
function lognormal_cdf(x, m, s, x0 = 0) {
   const z = (Math.log((x-x0))-m)/(s*kSqrt2);
   if (z < -1.) return 0.5*erfc(-z);
   else return 0.5*(1.0 + erf(z));
}

/** @summary normal_cdf_c function
  * @memberof Math */
function normal_cdf_c(x, sigma, x0 = 0) {
   const z = (x-x0)/(sigma*kSqrt2);
   if (z > 1.) return 0.5*erfc(z);
   else return 0.5*(1.-erf(z));
}

/** @summary normal_cdf function
  * @memberof Math */
function normal_cdf(x, sigma, x0 = 0) {
   const z = (x-x0)/(sigma*kSqrt2);
   if (z < -1.) return 0.5*erfc(-z);
   else return 0.5*(1.0 + erf(z));
}

/** @summary log normal pdf
  * @memberof Math */
function lognormal_pdf(x, m, s, x0 = 0) {
   if ((x-x0) <= 0)
      return 0.0;
   const tmp = (Math.log((x-x0)) - m)/s;
   return 1.0 / ((x-x0) * Math.abs(s) * Math.sqrt(2 * M_PI)) * Math.exp(-(tmp * tmp) /2);
}

/** @summary normal pdf
  * @memberof Math */
function normal_pdf(x, sigma = 1, x0 = 0) {
   const tmp = (x-x0)/sigma;
   return (1.0/(Math.sqrt(2 * M_PI) * Math.abs(sigma))) * Math.exp(-tmp*tmp/2);
}

/** @summary gamma calculation
  * @memberof Math */
function gamma(x) {
   let p, q, z, i, sgngam = 1;

   if (x >= Number.MAX_VALUE)
      return x;

   q = Math.abs(x);

   if ( q > 33.0 )
   {
      if ( x < 0.0 )
      {
         p = Math.floor(q);
         if ( p == q )
            return Number.POSITIVE_INFINITY;
         i = Math.round(p);
         if ( (i & 1) == 0 )
            sgngam = -1;
         z = q - p;
         if ( z > 0.5 )
         {
            p += 1.0;
            z = q - p;
         }
         z = q * Math.sin(Math.PI * z);
         if ( z == 0 )
         {
            return sgngam > 0 ? Number.POSITIVE_INFINITY : Number.NEGATIVE_INFINITY;
         }
         z = Math.abs(z);
         z = Math.PI / (z * stirf(q));
      }
      else
      {
         z = stirf(x);
      }
      return sgngam * z;
   }

   z = 1.0;
   while ( x >= 3.0 ) {
      x -= 1.0;
      z *= x;
   }

  let small = false;

   while (( x < 0.0 ) && !small) {
      if ( x > -1.E-9 )
         small = true;
      else {
         z /= x;
         x += 1.0;
      }
   }

   while (( x < 2.0 ) && !small) {
      if ( x < 1.e-9 )
         small = true;
      else {
         z /= x;
         x += 1.0;
      }
   }

   if (small) {
      if ( x == 0 )
         return Number.POSITIVE_INFINITY;
      else
         return z/((1.0 + 0.5772156649015329 * x) * x);
   }

   if ( x == 2.0 )
      return z;

   const P = [
      1.60119522476751861407E-4,
      1.19135147006586384913E-3,
      1.04213797561761569935E-2,
      4.76367800457137231464E-2,
      2.07448227648435975150E-1,
      4.94214826801497100753E-1,
      9.99999999999999996796E-1
   ], Q = [
      -2.31581873324120129819E-5,
      5.39605580493303397842E-4,
      -4.45641913851797240494E-3,
      1.18139785222060435552E-2,
      3.58236398605498653373E-2,
      -2.34591795718243348568E-1,
      7.14304917030273074085E-2,
      1.00000000000000000320E0];

   x -= 2.0;
   p = Polynomialeval( x, P, 6 );
   q = Polynomialeval( x, Q, 7 );
   return z * p / q;
}

/** @summary ndtri function
  * @memberof Math */
function ndtri(y0) {
   if ( y0 <= 0.0 )
      return Number.NEGATIVE_INFINITY;
   if ( y0 >= 1.0 )
      return Number.POSITIVE_INFINITY;

   const P0 = [
        -5.99633501014107895267E1,
         9.80010754185999661536E1,
        -5.66762857469070293439E1,
         1.39312609387279679503E1,
        -1.23916583867381258016E0
   ], Q0 = [
         1.95448858338141759834E0,
         4.67627912898881538453E0,
         8.63602421390890590575E1,
        -2.25462687854119370527E2,
         2.00260212380060660359E2,
        -8.20372256168333339912E1,
         1.59056225126211695515E1,
        -1.18331621121330003142E0
   ], P1 = [
         4.05544892305962419923E0,
         3.15251094599893866154E1,
         5.71628192246421288162E1,
         4.40805073893200834700E1,
         1.46849561928858024014E1,
         2.18663306850790267539E0,
        -1.40256079171354495875E-1,
        -3.50424626827848203418E-2,
        -8.57456785154685413611E-4
   ], Q1 = [
         1.57799883256466749731E1,
         4.53907635128879210584E1,
         4.13172038254672030440E1,
         1.50425385692907503408E1,
         2.50464946208309415979E0,
        -1.42182922854787788574E-1,
        -3.80806407691578277194E-2,
        -9.33259480895457427372E-4
   ], P2 = [
         3.23774891776946035970E0,
         6.91522889068984211695E0,
         3.93881025292474443415E0,
         1.33303460815807542389E0,
         2.01485389549179081538E-1,
         1.23716634817820021358E-2,
         3.01581553508235416007E-4,
         2.65806974686737550832E-6,
         6.23974539184983293730E-9
   ], Q2 = [
         6.02427039364742014255E0,
         3.67983563856160859403E0,
         1.37702099489081330271E0,
         2.16236993594496635890E-1,
         1.34204006088543189037E-2,
         3.28014464682127739104E-4,
         2.89247864745380683936E-6,
         6.79019408009981274425E-9
   ], s2pi = 2.50662827463100050242e0, dd = 0.13533528323661269189;

   let code = 1, y = y0, x, y2, x1;

   if (y > (1.0 - dd)) {
      y = 1.0 - y;
      code = 0;
   }
   if ( y > dd ) {
      y = y - 0.5;
      y2 = y * y;
      x = y + y * (y2 * Polynomialeval( y2, P0, 4)/ Polynomial1eval( y2, Q0, 8 ));
      x = x * s2pi;
      return x;
   }
   x = Math.sqrt(-2.0 * Math.log(y));
   const x0 = x - Math.log(x)/x,
         z = 1.0/x;
   if ( x < 8.0 )
      x1 = z * Polynomialeval( z, P1, 8 )/ Polynomial1eval( z, Q1, 8 );
   else
      x1 = z * Polynomialeval( z, P2, 8 )/ Polynomial1eval( z, Q2, 8 );
   x = x0 - x1;
   if ( code != 0 )
      x = -x;
   return x;
}

/** @summary normal_quantile function
  * @memberof Math */
function normal_quantile(z, sigma) {
   return sigma * ndtri(z);
}

/** @summary normal_quantile_c function
  * @memberof Math */
function normal_quantile_c(z, sigma) {
   return -sigma * ndtri(z);
}

/** @summary igamc function
  * @memberof Math */
function igamc(a,x) {
   // LM: for negative values returns 0.0
   // This is correct if a is a negative integer since Gamma(-n) = +/- inf
   if (a <= 0) return 0.0;

   if (x <= 0) return 1.0;

   if ((x < 1.0) || (x < a))
      return (1.0 - igam(a,x));

   let ax = a * Math.log(x) - x - lgam(a);
   if ( ax < -kMAXLOG )
      return 0.0;

   ax = Math.exp(ax);

   /* continued fraction */
   let y = 1.0 - a,
       z = x + y + 1.0,
       c = 0.0,
       pkm2 = 1.0,
       qkm2 = x,
       pkm1 = x + 1.0,
       qkm1 = z * x,
       ans = pkm1/qkm1,
       yc, r, t, pk, qk;

   do {
      c += 1.0;
      y += 1.0;
      z += 2.0;
      yc = y * c;
      pk = pkm1 * z - pkm2 * yc;
      qk = qkm1 * z - qkm2 * yc;
      if (qk)
      {
         r = pk/qk;
         t = Math.abs( (ans - r)/r );
         ans = r;
      }
      else
         t = 1.0;
      pkm2 = pkm1;
      pkm1 = pk;
      qkm2 = qkm1;
      qkm1 = qk;
      if ( Math.abs(pk) > kBig )
      {
         pkm2 *= kBiginv;
         pkm1 *= kBiginv;
         qkm2 *= kBiginv;
         qkm1 *= kBiginv;
      }
   } while ( t > kMACHEP );

   return ans * ax;
}

/** @summary igam function
  * @memberof Math */
function igam(a, x) {
   // LM: for negative values returns 1.0 instead of zero
   // This is correct if a is a negative integer since Gamma(-n) = +/- inf
   if (a <= 0) return 1.0;

   if (x <= 0) return 0.0;

   if ((x > 1.0) && (x > a))
      return 1.0 - igamc(a,x);

   /* Compute  x**a * exp(-x) / gamma(a)  */
   let ax = a * Math.log(x) - x - lgam(a);
   if ( ax < -kMAXLOG )
      return 0.0;

   ax = Math.exp(ax);

   /* power series */
   let r = a, c = 1.0, ans = 1.0;

   do {
      r += 1.0;
      c *= x/r;
      ans += c;
   } while ( c/ans > kMACHEP );

   return ans * ax/a;
}


/** @summary igami function
  * @memberof Math */
function igami(a, y0) {
   // check the domain
   if (a <= 0) {
      console.error(`igami : Wrong domain for parameter a = ${a} (must be > 0)`);
      return 0;
   }
   if (y0 <= 0) {
      return Number.POSITIVE_INFINITY;
   }
   if (y0 >= 1) {
      return 0;
   }
   const kMAXNUM = Number.MAX_VALUE, dithresh = 5.0 * kMACHEP;
   let x0 = kMAXNUM, x1 = 0, x, yl = 0, yh = 1, y, d, lgm, i, dir;

   /* approximation to inverse function */
   d = 1.0/(9.0*a);
   y = 1.0 - d - ndtri(y0) * Math.sqrt(d);
   x = a * y * y * y;

   lgm = lgam(a);

   for ( i=0; i<10; ++i ) {
      if ( x > x0 || x < x1 )
         break;
      y = igamc(a,x);
      if ( y < yl || y > yh )
         break;
      if ( y < y0 ) {
         x0 = x;
         yl = y;
      }
      else {
         x1 = x;
         yh = y;
      }
      /* compute the derivative of the function at this point */
      d = (a - 1.0) * Math.log(x) - x - lgm;
      if ( d < -kMAXLOG )
         break;
      d = -Math.exp(d);
      /* compute the step to the next approximation of x */
      d = (y - y0)/d;
      if ( Math.abs(d/x) < kMACHEP )
         return x;
      x = x - d;
   }
   /* Resort to interval halving if Newton iteration did not converge. */
   d = 0.0625;
   if ( x0 == kMAXNUM ) {
      if ( x <= 0.0 )
         x = 1.0;
      while ( x0 == kMAXNUM ) {
         x = (1.0 + d) * x;
         y = igamc( a, x );
         if ( y < y0 ) {
            x0 = x;
            yl = y;
            break;
         }
         d = d + d;
      }
   }
   d = 0.5;
   dir = 0;

   for ( i=0; i<400; ++i ) {
      x = x1 + d * (x0 - x1);
      y = igamc( a, x );
      lgm = (x0 - x1)/(x1 + x0);
      if ( Math.abs(lgm) < dithresh )
         break;
      lgm = (y - y0)/y0;
      if ( Math.abs(lgm) < dithresh )
         break;
      if ( x <= 0.0 )
         break;
      if ( y >= y0 ) {
         x1 = x;
         yh = y;
         if ( dir < 0 ) {
            dir = 0;
            d = 0.5;
         }
         else if ( dir > 1 )
            d = 0.5 * d + 0.5;
         else
            d = (y0 - yl)/(yh - yl);
         dir += 1;
      }
      else {
         x0 = x;
         yl = y;
         if ( dir > 0 ) {
            dir = 0;
            d = 0.5;
         }
         else if ( dir < -1 )
            d = 0.5 * d;
         else
            d = (y0 - yl)/(yh - yl);
         dir -= 1;
      }
   }
   return x;
}

/** @summary landau_pdf function
  * @desc LANDAU pdf : algorithm from CERNLIB G110 denlan
  *  same algorithm is used in GSL
  * @memberof Math */
function landau_pdf(x, xi, x0 = 0) {
   if (xi <= 0) return 0;
   const v = (x - x0)/xi;
   let u, ue, us, denlan;
   const p1 = [0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253],
         q1 = [1.0 ,-0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063],
         p2 = [0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211],
         q2 = [1.0 , 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714],
         p3 = [0.1788544503, 0.09359161662,0.006325387654, 0.00006611667319,-0.000002031049101],
         q3 = [1.0 , 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675],
         p4 = [0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186],
         q4 = [1.0 , 106.8615961, 337.6496214, 2016.712389, 1597.063511],
         p5 = [1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910],
         q5 = [1.0 , 156.9424537, 3745.310488, 9834.698876, 66924.28357],
         p6 = [1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109],
         q6 = [1.0 , 651.4101098, 56974.73333, 165917.4725, -2815759.939],
         a1 = [0.04166666667,-0.01996527778, 0.02709538966],
         a2 = [-1.845568670,-4.284640743];

   if (v < -5.5) {
      u = Math.exp(v+1.0);
      if (u < 1e-10) return 0.0;
      ue = Math.exp(-1/u);
      us = Math.sqrt(u);
      denlan = 0.3989422803*(ue/us)*(1+(a1[0]+(a1[1]+a1[2]*u)*u)*u);
   } else if (v < -1) {
      u = Math.exp(-v-1);
      denlan = Math.exp(-u)*Math.sqrt(u)*
         (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
         (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
   } else if (v < 1) {
      denlan = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*v)*v)*v)*v)/
         (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*v)*v)*v)*v);
   } else if (v < 5) {
      denlan = (p3[0]+(p3[1]+(p3[2]+(p3[3]+p3[4]*v)*v)*v)*v)/
         (q3[0]+(q3[1]+(q3[2]+(q3[3]+q3[4]*v)*v)*v)*v);
   } else if (v < 12) {
      u = 1/v;
      denlan = u*u*(p4[0]+(p4[1]+(p4[2]+(p4[3]+p4[4]*u)*u)*u)*u)/
         (q4[0]+(q4[1]+(q4[2]+(q4[3]+q4[4]*u)*u)*u)*u);
   } else if (v < 50) {
      u = 1/v;
      denlan = u*u*(p5[0]+(p5[1]+(p5[2]+(p5[3]+p5[4]*u)*u)*u)*u)/
         (q5[0]+(q5[1]+(q5[2]+(q5[3]+q5[4]*u)*u)*u)*u);
   } else if (v < 300) {
      u = 1/v;
      denlan = u*u*(p6[0]+(p6[1]+(p6[2]+(p6[3]+p6[4]*u)*u)*u)*u)/
         (q6[0]+(q6[1]+(q6[2]+(q6[3]+q6[4]*u)*u)*u)*u);
   } else {
      u = 1/(v-v*Math.log(v)/(v+1));
      denlan = u*u*(1+(a2[0]+a2[1]*u)*u);
   }
   return denlan/xi;
}

/** @summary Landau function
  * @memberof Math */
function Landau(x, mpv, sigma, norm) {
   if (sigma <= 0) return 0;
   const den = landau_pdf((x - mpv) / sigma, 1, 0);
   if (!norm) return den;
   return den/sigma;
}

/** @summary inc_gamma_c
  * @memberof Math */
function inc_gamma_c(a,x) {
   return igamc(a,x);
}

/** @summary inc_gamma
  * @memberof Math */
function inc_gamma(a,x) {
   return igam(a,x);
}

/** @summary lgamma
  * @memberof Math */
function lgamma(z) {
   return lgam(z);
}

/** @summary Probability density function of the beta distribution.
  * @memberof Math */
function beta_pdf(x, a, b) {
  if (x < 0 || x > 1.0) return 0;
  if (x == 0 ) {
     if (a < 1) return Number.POSITIVE_INFINITY;
     else if (a > 1) return 0;
     else if ( a == 1) return b; // to avoid a nan from log(0)*0
   }
   if (x == 1 ) {
      if (b < 1) return Number.POSITIVE_INFINITY;
      else if (b > 1) return 0;
      else if ( b == 1) return a; // to avoid a nan from log(0)*0
   }
   return Math.exp(lgamma(a + b) - lgamma(a) - lgamma(b) +
                    Math.log(x) * (a -1.) + Math.log1p(-x) * (b - 1.));
}

/** @summary beta
  * @memberof Math */
function beta(x,y) {
   return Math.exp(lgamma(x)+lgamma(y)-lgamma(x+y));
}

/** @summary chisquared_cdf_c
  * @memberof Math */
function chisquared_cdf_c(x,r,x0 = 0) {
   return inc_gamma_c(0.5 * r, 0.5*(x-x0));
}

/** @summary Continued fraction expansion #1 for incomplete beta integral
  * @memberof Math */
function incbcf(a,b,x) {
   let xk, pk, pkm1, pkm2, qk, qkm1, qkm2,
       k1, k2, k3, k4, k5, k6, k7, k8,
       r, t, ans, n;
   const thresh = 3.0 * kMACHEP;

   k1 = a;
   k2 = a + b;
   k3 = a;
   k4 = a + 1.0;
   k5 = 1.0;
   k6 = b - 1.0;
   k7 = k4;
   k8 = a + 2.0;

   pkm2 = 0.0;
   qkm2 = 1.0;
   pkm1 = 1.0;
   qkm1 = 1.0;
   ans = 1.0;
   r = 1.0;
   n = 0;

   do {
      xk = -( x * k1 * k2 )/( k3 * k4 );
      pk = pkm1 + pkm2 * xk;
      qk = qkm1 + qkm2 * xk;
      pkm2 = pkm1;
      pkm1 = pk;
      qkm2 = qkm1;
      qkm1 = qk;

      xk = ( x * k5 * k6 )/( k7 * k8 );
      pk = pkm1 + pkm2 * xk;
      qk = qkm1 + qkm2 * xk;
      pkm2 = pkm1;
      pkm1 = pk;
      qkm2 = qkm1;
      qkm1 = qk;

      if ( qk !=0 )
         r = pk/qk;
      if ( r != 0 )
      {
         t = Math.abs( (ans - r)/r );
         ans = r;
      }
      else
         t = 1.0;

      if ( t < thresh )
         break; // goto cdone;

      k1 += 1.0;
      k2 += 1.0;
      k3 += 2.0;
      k4 += 2.0;
      k5 += 1.0;
      k6 -= 1.0;
      k7 += 2.0;
      k8 += 2.0;

      if ((Math.abs(qk) + Math.abs(pk)) > kBig) {
         pkm2 *= kBiginv;
         pkm1 *= kBiginv;
         qkm2 *= kBiginv;
         qkm1 *= kBiginv;
      }
      if ((Math.abs(qk) < kBiginv) || (Math.abs(pk) < kBiginv)) {
         pkm2 *= kBig;
         pkm1 *= kBig;
         qkm2 *= kBig;
         qkm1 *= kBig;
      }
   }
   while ( ++n < 300 );

// cdone:
   return ans;
}

/** @summary Continued fraction expansion #2 for incomplete beta integral
  * @memberof Math */
function incbd(a,b,x) {
   const z = x / (1.0-x),
         thresh = 3.0 * kMACHEP;

   let xk, pk, pkm1, pkm2, qk, qkm1, qkm2,
       k1, k2, k3, k4, k5, k6, k7, k8,
       r, t, ans, n;

   k1 = a;
   k2 = b - 1.0;
   k3 = a;
   k4 = a + 1.0;
   k5 = 1.0;
   k6 = a + b;
   k7 = a + 1.0;
   k8 = a + 2.0;

   pkm2 = 0.0;
   qkm2 = 1.0;
   pkm1 = 1.0;
   qkm1 = 1.0;
   ans = 1.0;
   r = 1.0;
   n = 0;
   do {
      xk = -( z * k1 * k2 )/( k3 * k4 );
      pk = pkm1 + pkm2 * xk;
      qk = qkm1 + qkm2 * xk;
      pkm2 = pkm1;
      pkm1 = pk;
      qkm2 = qkm1;
      qkm1 = qk;

      xk = ( z * k5 * k6 )/( k7 * k8 );
      pk = pkm1 + pkm2 * xk;
      qk = qkm1 + qkm2 * xk;
      pkm2 = pkm1;
      pkm1 = pk;
      qkm2 = qkm1;
      qkm1 = qk;

      if ( qk != 0 )
         r = pk/qk;
      if ( r != 0 )
      {
         t = Math.abs( (ans - r)/r );
         ans = r;
      }
      else
         t = 1.0;

      if ( t < thresh )
         break; // goto cdone;

      k1 += 1.0;
      k2 -= 1.0;
      k3 += 2.0;
      k4 += 2.0;
      k5 += 1.0;
      k6 += 1.0;
      k7 += 2.0;
      k8 += 2.0;

      if ((Math.abs(qk) + Math.abs(pk)) > kBig) {
         pkm2 *= kBiginv;
         pkm1 *= kBiginv;
         qkm2 *= kBiginv;
         qkm1 *= kBiginv;
      }
      if ((Math.abs(qk) < kBiginv) || (Math.abs(pk) < kBiginv)) {
         pkm2 *= kBig;
         pkm1 *= kBig;
         qkm2 *= kBig;
         qkm1 *= kBig;
      }
   }
   while ( ++n < 300 );
// cdone:
   return ans;
}

/** @summary ROOT::Math::Cephes::pseries
  * @memberof Math */
function pseries(a,b,x) {
   let s, t, u, v, n;

   const ai = 1.0 / a;
   u = (1.0 - b) * x;
   v = u / (a + 1.0);
   const t1 = v;
   t = u;
   n = 2.0;
   s = 0.0;
   const z = kMACHEP * ai;
   while ( Math.abs(v) > z ) {
      u = (n - b) * x / n;
      t *= u;
      v = t / (a + n);
      s += v;
      n += 1.0;
   }
   s += t1;
   s += ai;

   u = a * Math.log(x);
   if ( (a+b) < kMAXSTIR && Math.abs(u) < kMAXLOG )
   {
      t = gamma(a+b) / (gamma(a)*gamma(b));
      s = s * t * Math.pow(x,a);
   }
   else
   {
      t = lgam(a+b) - lgam(a) - lgam(b) + u + Math.log(s);
      if ( t < kMINLOG )
         s = 0.0;
      else
         s = Math.exp(t);
   }
   return s;
}

/** @summary ROOT::Math::Cephes::incbet
  * @memberof Math */
function incbet(aa,bb,xx) {
   let a, b, t, x, xc, w, y, flag;

   if ( aa <= 0.0 || bb <= 0.0 )
      return 0.0;

   // LM: changed: for X > 1 return 1.
   if (xx <= 0.0) return 0.0;
   if ( xx >= 1.0) return 1.0;

   flag = 0;

/* - to test if that way is better for large b/  (comment out from Cephes version)
   if ( (bb * xx) <= 1.0 && xx <= 0.95)
   {
   t = pseries(aa, bb, xx);
   goto done;
   }

**/
   w = 1.0 - xx;

/* Reverse a and b if x is greater than the mean. */
/* aa,bb > 1 -> sharp rise at x=aa/(aa+bb) */
   if (xx > (aa/(aa+bb)))
   {
      flag = 1;
      a = bb;
      b = aa;
      xc = xx;
      x = w;
   }
   else
   {
      a = aa;
      b = bb;
      xc = w;
      x = xx;
   }

   if ( flag == 1 && (b * x) <= 1.0 && x <= 0.95) {
      t = pseries(a, b, x);
      // goto done;
   } else {
      /* Choose expansion for better convergence. */
      y = x * (a+b-2.0) - (a-1.0);
      if ( y < 0.0 )
         w = incbcf( a, b, x );
      else
         w = incbd( a, b, x ) / xc;

   /* Multiply w by the factor
      a      b   _             _     _
      x  (1-x)   | (a+b) / (a | (a) | (b)) .   */

      y = a * Math.log(x);
      t = b * Math.log(xc);
      if ( (a+b) < kMAXSTIR && Math.abs(y) < kMAXLOG && Math.abs(t) < kMAXLOG )
      {
         t = Math.pow(xc,b);
         t *= Math.pow(x,a);
         t /= a;
         t *= w;
         t *= gamma(a+b) / (gamma(a) * gamma(b));
         // goto done;
      } else {
      /* Resort to logarithms.  */
         y += t + lgam(a+b) - lgam(a) - lgam(b);
         y += Math.log(w/a);
         if ( y < kMINLOG )
            t = 0.0;
         else
            t = Math.exp(y);
      }
   }

// done:

   if (flag == 1) {
      if ( t <= kMACHEP )
         t = 1.0 - kMACHEP;
      else
         t = 1.0 - t;
   }
   return t;
}

/** @summary copy of ROOT::Math::Cephes::incbi
  * @memberof Math */
function incbi(aa,bb,yy0) {
   let a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt,
       i, rflg, dir, nflg, ihalve = true;

   // check the domain
   if (aa <= 0) {
      // MATH_ERROR_MSG('Cephes::incbi','Wrong domain for parameter a (must be > 0)');
      return 0;
   }
   if (bb <= 0) {
      // MATH_ERROR_MSG('Cephes::incbi','Wrong domain for parameter b (must be > 0)');
      return 0;
   }

   const process_done = () => {
      if ( rflg ) {
         if ( x <= kMACHEP )
            x = 1.0 - kMACHEP;
         else
            x = 1.0 - x;
      }
      return x;
   };

   i = 0;
   if ( yy0 <= 0 )
      return 0.0;
   if ( yy0 >= 1.0 )
      return 1.0;
   x0 = 0.0;
   yl = 0.0;
   x1 = 1.0;
   yh = 1.0;
   nflg = 0;

   if ( aa <= 1.0 || bb <= 1.0 )
   {
      dithresh = 1.0e-6;
      rflg = 0;
      a = aa;
      b = bb;
      y0 = yy0;
      x = a/(a+b);
      y = incbet( a, b, x );
      // goto ihalve; // will start
   }
   else
   {
      dithresh = 1.0e-4;
/* approximation to inverse function */

      yp = -ndtri(yy0);

      if ( yy0 > 0.5 )
      {
         rflg = 1;
         a = bb;
         b = aa;
         y0 = 1.0 - yy0;
         yp = -yp;
      }
      else
      {
         rflg = 0;
         a = aa;
         b = bb;
         y0 = yy0;
      }

      lgm = (yp * yp - 3.0)/6.0;
      x = 2.0/(1.0/(2.0*a-1.0) + 1.0/(2.0*b-1.0));
      d = yp * Math.sqrt( x + lgm ) / x
         - (1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0))
         * (lgm + 5.0/6.0 - 2.0/(3.0*x));
      d = 2.0 * d;
      if ( d < kMINLOG ) {
         // x = 1.0;
         // goto under;
         x = 0.0;
         return process_done();
      }
      x = a/(a + b * Math.exp(d));
      y = incbet( a, b, x );
      yp = (y - y0)/y0;
      if ( Math.abs(yp) < 0.2 )
         ihalve = false; // instead goto newt; exclude ihalve for the first time
   }

  let mainloop = 1000;

  // endless loop until coverage
  while (mainloop-- > 0) {
   /* Resort to interval halving if not close enough. */
   // ihalve:
      if (ihalve) {
         dir = 0;
         di = 0.5;
         for ( i=0; i<100; i++ )
         {
            if ( i != 0 )
            {
               x = x0 + di * (x1 - x0);
               if ( x == 1.0 )
                  x = 1.0 - kMACHEP;
               if ( x == 0.0 )
               {
                  di = 0.5;
                  x = x0 + di * (x1 - x0);
                  if ( x == 0.0 )
                     return process_done(); // goto under;
               }
               y = incbet( a, b, x );
               yp = (x1 - x0)/(x1 + x0);
               if ( Math.abs(yp) < dithresh )
                  break; // goto newt;
               yp = (y-y0)/y0;
               if ( Math.abs(yp) < dithresh )
                  break; // goto newt;
            }
            if ( y < y0 )
            {
               x0 = x;
               yl = y;
               if ( dir < 0 )
               {
                  dir = 0;
                  di = 0.5;
               }
               else if ( dir > 3 )
                  di = 1.0 - (1.0 - di) * (1.0 - di);
               else if ( dir > 1 )
                  di = 0.5 * di + 0.5;
               else
                  di = (y0 - y)/(yh - yl);
               dir += 1;
               if ( x0 > 0.75 )
               {
                  if ( rflg == 1 )
                  {
                     rflg = 0;
                     a = aa;
                     b = bb;
                     y0 = yy0;
                  }
                  else
                  {
                     rflg = 1;
                     a = bb;
                     b = aa;
                     y0 = 1.0 - yy0;
                  }
                  x = 1.0 - x;
                  y = incbet( a, b, x );
                  x0 = 0.0;
                  yl = 0.0;
                  x1 = 1.0;
                  yh = 1.0;
                  continue; // goto ihalve;
               }
            }
            else
            {
               x1 = x;
               if ( rflg == 1 && x1 < kMACHEP )
               {
                  x = 0.0;
                  return process_done(); // goto done;
               }
               yh = y;
               if ( dir > 0 )
               {
                  dir = 0;
                  di = 0.5;
               }
               else if ( dir < -3 )
                  di = di * di;
               else if ( dir < -1 )
                  di = 0.5 * di;
               else
                  di = (y - y0)/(yh - yl);
               dir -= 1;
            }
         }
         // math_error( 'incbi', PLOSS );
         if ( x0 >= 1.0 ) {
            x = 1.0 - kMACHEP;
            return process_done(); // goto done;
         }
         if ( x <= 0.0 ) {
            // math_error( 'incbi', UNDERFLOW );
            x = 0.0;
            return process_done(); // goto done;
         }
         break; // if here, break ihalve
      } // end of ihalve

      ihalve = true; // enter loop next time

   // newt:

      if ( nflg )
         return process_done(); // goto done;
      nflg = 1;
      lgm = lgam(a+b) - lgam(a) - lgam(b);

      for ( i=0; i<8; i++ )
      {
         /* Compute the function at this point. */
         if ( i != 0 )
            y = incbet(a,b,x);
         if ( y < yl )
         {
            x = x0;
            y = yl;
         }
         else if ( y > yh )
         {
            x = x1;
            y = yh;
         }
         else if ( y < y0 )
         {
            x0 = x;
            yl = y;
         }
         else
         {
            x1 = x;
            yh = y;
         }
         if ( x == 1.0 || x == 0.0 )
            break;
         /* Compute the derivative of the function at this point. */
         d = (a - 1.0) * Math.log(x) + (b - 1.0) * Math.log(1.0-x) + lgm;
         if ( d < kMINLOG )
            return process_done(); // goto done;
         if ( d > kMAXLOG )
            break;
         d = Math.exp(d);
         /* Compute the step to the next approximation of x. */
         d = (y - y0)/d;
         xt = x - d;
         if ( xt <= x0 )
         {
            y = (x - x0) / (x1 - x0);
            xt = x0 + 0.5 * y * (x - x0);
            if ( xt <= 0.0 )
               break;
         }
         if ( xt >= x1 )
         {
            y = (x1 - x) / (x1 - x0);
            xt = x1 - 0.5 * y * (x1 - x);
            if ( xt >= 1.0 )
               break;
         }
         x = xt;
         if ( Math.abs(d/x) < 128.0 * kMACHEP )
            return process_done(); // goto done;
      }
      /* Did not converge.  */
      dithresh = 256.0 * kMACHEP;
   } // endless loop instead of // goto ihalve;

// done:

   return process_done();
}

/** @summary Calculates the normalized (regularized) incomplete beta function.
  * @memberof Math */
function inc_beta(x,a,b) {
   return incbet(a,b,x);
}

const BetaIncomplete = inc_beta;

/** @summary ROOT::Math::beta_quantile
  * @memberof Math */
function beta_quantile(z,a,b) {
   return incbi(a,b,z);
}

/** @summary Complement of the cumulative distribution function of the beta distribution.
  * @memberof Math */
function beta_cdf_c(x,a,b) {
   return inc_beta(1-x, b, a);
}

/** @summary chisquared_cdf
  * @memberof Math */
function chisquared_cdf(x,r,x0=0) {
   return inc_gamma(0.5 * r, 0.5*(x-x0));
}

/** @summary gamma_quantile_c function
  * @memberof Math */
function gamma_quantile_c(z, alpha, theta) {
   return theta * igami( alpha, z);
}

/** @summary gamma_quantile function
  * @memberof Math */
function gamma_quantile(z, alpha, theta) {
   return theta * igami( alpha, 1.- z);
}

/** @summary breitwigner_cdf_c function
  * @memberof Math */
function breitwigner_cdf_c(x,gamma, x0 = 0) {
   return 0.5 - Math.atan(2.0 * (x-x0) / gamma) / M_PI;
}

/** @summary breitwigner_cdf function
  * @memberof Math */
function breitwigner_cdf(x, gamma, x0 = 0) {
   return 0.5 + Math.atan(2.0 * (x-x0) / gamma) / M_PI;
}

/** @summary cauchy_cdf_c function
  * @memberof Math */
function cauchy_cdf_c(x, b, x0 = 0) {
   return 0.5 - Math.atan( (x-x0) / b) / M_PI;
}

/** @summary cauchy_cdf function
  * @memberof Math */
function cauchy_cdf(x, b, x0 = 0) {
   return 0.5 + Math.atan( (x-x0) / b) / M_PI;
}

/** @summary cauchy_pdf function
  * @memberof Math */
function cauchy_pdf(x, b = 1, x0 = 0) {
   return b/(M_PI * ((x-x0)*(x-x0) + b*b));
}

/** @summary gaussian_pdf function
  * @memberof Math */
function gaussian_pdf(x, sigma = 1, x0 = 0) {
   const tmp = (x-x0)/sigma;
   return (1.0/(Math.sqrt(2 * M_PI) * Math.abs(sigma))) * Math.exp(-tmp*tmp/2);
}

/** @summary gamma_pdf function
  * @memberof Math */
function gamma_pdf(x, alpha, theta, x0 = 0) {
   if ((x - x0) < 0) {
      return 0.0;
   } else if ((x - x0) == 0) {
      return (alpha == 1) ? 1.0 / theta : 0;
   } else if (alpha == 1) {
      return Math.exp(-(x - x0) / theta) / theta;
   }
   return Math.exp((alpha - 1) * Math.log((x - x0) / theta) - (x - x0) / theta - lgamma(alpha)) / theta;
}

/** @summary tdistribution_cdf_c function
  * @memberof Math */
function tdistribution_cdf_c(x, r, x0 = 0) {
   const p = x - x0,
         sign = (p > 0) ? 1. : -1;
   return .5 - .5*inc_beta(p*p/(r + p*p), .5, .5*r)*sign;
}

/** @summary tdistribution_cdf function
  * @memberof Math */
function tdistribution_cdf(x, r, x0 = 0) {
   const p = x - x0,
         sign = (p > 0) ? 1. : -1;
   return .5 + .5*inc_beta(p*p/(r + p*p), .5, .5*r)*sign;
}

/** @summary tdistribution_pdf function
  * @memberof Math */
function tdistribution_pdf(x, r, x0 = 0) {
   return (Math.exp(lgamma((r + 1.0)/2.0) - lgamma(r/2.0)) / Math.sqrt(M_PI * r))
          * Math.pow((1.0 + (x-x0)*(x-x0)/r), -(r + 1.0)/2.0);
}

/** @summary exponential_cdf_c function
  * @memberof Math */
function exponential_cdf_c(x, lambda, x0 = 0) {
   return ((x-x0) < 0) ? 1.0 : Math.exp(-lambda * (x-x0));
}

/** @summary exponential_cdf function
  * @memberof Math */
function exponential_cdf(x, lambda, x0 = 0) {
   return ((x-x0) < 0) ? 0.0 : -Math.expm1(-lambda * (x-x0));
}

/** @summary chisquared_pdf
  * @memberof Math */
function chisquared_pdf(x, r, x0 = 0) {
   if ((x-x0) < 0) return 0.0;
   const a = r/2 -1.;
   // let return inf for case x  = x0 and treat special case of r = 2 otherwise will return nan
   if (x == x0 && a == 0) return 0.5;

   return Math.exp((r/2 - 1) * Math.log((x-x0)/2) - (x-x0)/2 - lgamma(r/2))/2;
}

/** @summary Probability density function of the F-distribution.
  * @memberof Math */
function fdistribution_pdf(x, n, m, x0 = 0) {
   if (n < 0 || m < 0)
      return Number.NaN;
   if ((x-x0) < 0)
      return 0.0;

   return Math.exp((n/2) * Math.log(n) + (m/2) * Math.log(m) + lgamma((n+m)/2) - lgamma(n/2) - lgamma(m/2)
                 + (n/2 -1) * Math.log(x-x0) - ((n+m)/2) * Math.log(m + n*(x-x0)));
}

/** @summary fdistribution_cdf_c function
  * @memberof Math */
function fdistribution_cdf_c(x, n, m, x0 = 0) {
   if (n < 0 || m < 0) return Number.NaN;

   const z = m / (m + n * (x - x0));
   // fox z->1 and large a and b IB looses precision use complement function
   if (z > 0.9 && n > 1 && m > 1) return 1. - fdistribution_cdf(x, n, m, x0);

   // for the complement use the fact that IB(x,a,b) = 1. - IB(1-x,b,a)
   return inc_beta(m / (m + n * (x - x0)), .5 * m, .5 * n);
}

/** @summary fdistribution_cdf function
  * @memberof Math */
function fdistribution_cdf(x, n, m, x0 = 0) {
   if (n < 0 || m < 0) return Number.NaN;

   const z = n * (x - x0) / (m + n * (x - x0));
   // fox z->1 and large a and b IB looses precision use complement function
   if (z > 0.9 && n > 1 && m > 1)
      return 1. - fdistribution_cdf_c(x, n, m, x0);

   return inc_beta(z, .5 * n, .5 * m);
}

/** @summary Prob function
  * @memberof Math */
function Prob(chi2, ndf) {
   if (ndf <= 0) return 0; // Set CL to zero in case ndf <= 0

   if (chi2 <= 0) {
      if (chi2 < 0) return 0;
      else return 1;
   }

   return chisquared_cdf_c(chi2,ndf,0);
}

/** @summary Square function
  * @memberof Math */
function Sq(x) {
   return x * x;
}

/** @summary Gaus function
  * @memberof Math */
function Gaus(x, mean, sigma, norm) {
   if (!sigma) return 1e30;
   const arg = (x - mean) / sigma;
   if (arg < -39 || arg > 39) return 0;
   const res = Math.exp(-0.5*arg*arg);
   return norm ? res/(2.50662827463100024*sigma) : res; // sqrt(2*Pi)=2.50662827463100024
}

/** @summary BreitWigner function
  * @memberof Math */
function BreitWigner(x, mean, gamma) {
   return gamma/((x-mean)*(x-mean) + gamma*gamma/4) / 2 / Math.PI;
}

/** @summary Calculates Beta-function Gamma(p)*Gamma(q)/Gamma(p+q).
  * @memberof Math */
function Beta(x,y) {
   return Math.exp(lgamma(x) + lgamma(y) - lgamma(x+y));
}

/** @summary GammaDist function
  * @memberof Math */
function GammaDist(x, gamma, mu = 0, beta = 1) {
   if ((x < mu) || (gamma <= 0) || (beta <= 0)) return 0;
   return gamma_pdf(x, gamma, beta, mu);
}

/** @summary probability density function of Laplace distribution
  * @memberof Math */
function LaplaceDist(x, alpha = 0, beta = 1) {
   return Math.exp(-Math.abs((x-alpha)/beta)) / (2.*beta);
}

/** @summary distribution function of Laplace distribution
  * @memberof Math */
function LaplaceDistI(x, alpha = 0, beta = 1) {
   return (x <= alpha) ? 0.5*Math.exp(-Math.abs((x-alpha)/beta)) : 1 - 0.5*Math.exp(-Math.abs((x-alpha)/beta));
}

/** @summary density function for Student's t- distribution
  * @memberof Math */
function Student(T, ndf) {
   if (ndf < 1) return 0;

   const r = ndf,
         rh = 0.5*r,
         rh1 = rh + 0.5,
         denom = Math.sqrt(r*Math.PI)*gamma(rh)*Math.pow(1+T*T/r, rh1);
   return gamma(rh1)/denom;
}

/** @summary cumulative distribution function of Student's
  * @memberof Math */
function StudentI(T, ndf) {
   const r = ndf;

   return (T > 0)
     ? (1 - 0.5*BetaIncomplete((r/(r + T*T)), r*0.5, 0.5))
     : 0.5*BetaIncomplete((r/(r + T*T)), r*0.5, 0.5);
}

/** @summary LogNormal function
  * @memberof Math */
function LogNormal(x, sigma, theta = 0, m = 1) {
   if ((x < theta) || (sigma <= 0) || (m <= 0)) return 0;
   return lognormal_pdf(x, Math.log(m), sigma, theta);
}

/** @summary Computes the probability density function of the Beta distribution
  * @memberof Math */
function BetaDist(x, p, q) {
   if ((x < 0) || (x > 1) || (p <= 0) || (q <= 0))
     return 0;
   const beta = Beta(p, q);
   return Math.pow(x, p-1) * Math.pow(1-x, q-1) / beta;
}

/** @summary Computes the distribution function of the Beta distribution.
  * @memberof Math */
function BetaDistI(x, p, q) {
   if ((x < 0) || (x > 1) || (p <= 0) || (q <= 0)) return 0;
   return BetaIncomplete(x, p, q);
}

/** @summary gaus function for TFormula
  * @memberof Math */
function gaus(f, x, i) {
   return f.GetParValue(i+0) * Math.exp(-0.5 * Math.pow((x-f.GetParValue(i+1)) / f.GetParValue(i+2), 2));
}

/** @summary gausn function for TFormula
  * @memberof Math */
function gausn(f, x, i) {
   return gaus(f, x, i)/(Math.sqrt(2 * Math.PI) * f.GetParValue(i+2));
}

/** @summary gausxy function for TFormula
  * @memberof Math */
function gausxy(f, x, y, i) {
   return f.GetParValue(i+0) * Math.exp(-0.5 * Math.pow((x-f.GetParValue(i+1)) / f.GetParValue(i+2), 2))
                             * Math.exp(-0.5 * Math.pow((y-f.GetParValue(i+3)) / f.GetParValue(i+4), 2));
}

/** @summary expo function for TFormula
  * @memberof Math */
function expo(f, x, i) {
   return Math.exp(f.GetParValue(i+0) + f.GetParValue(i+1) * x);
}

/** @summary landau function for TFormula
  * @memberof Math */
function landau(f, x, i) {
   return Landau(x, f.GetParValue(i+1),f.GetParValue(i+2), false);
}

/** @summary landaun function for TFormula
  * @memberof Math */
function landaun(f, x, i) {
   return Landau(x, f.GetParValue(i+1),f.GetParValue(i+2), true);
}

/** @summary Crystal ball function
  * @memberof Math */
function crystalball_function(x, alpha, n, sigma, mean = 0) {
   if (sigma < 0.) return 0.;
   let z = (x - mean)/sigma;
   if (alpha < 0) z = -z;
   const abs_alpha = Math.abs(alpha);
   if (z > -abs_alpha)
      return Math.exp(-0.5 * z * z);
   const nDivAlpha = n/abs_alpha,
         AA = Math.exp(-0.5*abs_alpha*abs_alpha),
         B = nDivAlpha - abs_alpha,
         arg = nDivAlpha/(B-z);
  return AA * Math.pow(arg,n);
}

/** @summary pdf definition of the crystal_ball which is defined only for n > 1 otherwise integral is diverging
  * @memberof Math */
function crystalball_pdf(x, alpha, n, sigma, mean = 0) {
   if (sigma < 0.) return 0.;
   if (n <= 1) return Number.NaN;  // pdf is not normalized for n <=1
   const abs_alpha = Math.abs(alpha),
         C = n/abs_alpha * 1./(n-1.) * Math.exp(-alpha*alpha/2.),
         D = Math.sqrt(M_PI/2.)*(1.+erf(abs_alpha/Math.sqrt(2.))),
         N = 1./(sigma*(C+D));
   return N * crystalball_function(x,alpha,n,sigma,mean);
}

/** @summary compute the integral of the crystal ball function
  * @memberof Math */
function crystalball_integral(x, alpha, n, sigma, mean = 0) {
   if (sigma == 0) return 0;
   if (alpha == 0) return 0.;
   const useLog = (n == 1.0),
         abs_alpha = Math.abs(alpha);

   let z = (x-mean)/sigma, intgaus = 0., intpow = 0.;
   if (alpha < 0 ) z = -z;

   const sqrtpiover2 = Math.sqrt(M_PI/2.),
         sqrt2pi = Math.sqrt( 2.*M_PI),
         oneoversqrt2 = 1./Math.sqrt(2.);
   if (z <= -abs_alpha) {
      const A = Math.pow(n/abs_alpha,n) * Math.exp(-0.5 * alpha*alpha),
            B = n/abs_alpha - abs_alpha;

      if (!useLog) {
         const C = (n/abs_alpha) * (1./(n-1)) * Math.exp(-alpha*alpha/2.);
         intpow = C - A /(n-1.) * Math.pow(B-z,-n+1);
      }
      else {
         // for n=1 the primitive of 1/x is log(x)
         intpow = -A * Math.log( n / abs_alpha ) + A * Math.log(B - z);
      }
      intgaus = sqrtpiover2*(1. + erf(abs_alpha*oneoversqrt2));
   } else {
      intgaus = normal_cdf_c(z, 1);
      intgaus *= sqrt2pi;
      intpow = 0;
   }
   return sigma * (intgaus + intpow);
}

/** @summary crystalball_cdf function
  * @memberof Math */
function crystalball_cdf(x, alpha, n, sigma, mean = 0) {
   if (n <= 1.)
      return Number.NaN;

   const abs_alpha = Math.abs(alpha),
         C = n/abs_alpha * 1./(n-1.) * Math.exp(-alpha*alpha/2.),
         D = Math.sqrt(M_PI/2.)*(1. + erf(abs_alpha/Math.sqrt(2.))),
         totIntegral = sigma*(C+D),
         integral = crystalball_integral(x,alpha,n,sigma,mean);

   return (alpha > 0) ? 1. - integral/totIntegral : integral/totIntegral;
}

/** @summary crystalball_cdf_c function
  * @memberof Math */
function crystalball_cdf_c(x, alpha, n, sigma, mean = 0) {
   if (n <= 1.)
      return Number.NaN;

   const abs_alpha = Math.abs(alpha),
         C = n/abs_alpha * 1./(n-1.) * Math.exp(-alpha*alpha/2.),
         D = Math.sqrt(M_PI/2.)*(1. + erf(abs_alpha/Math.sqrt(2.))),
         totIntegral = sigma*(C+D),
         integral = crystalball_integral(x,alpha,n,sigma,mean);

   return (alpha > 0) ? integral/totIntegral : 1. - (integral/totIntegral);
}

/** @summary ChebyshevN function
  * @memberof Math */
function ChebyshevN(n, x, c) {
   let d1 = 0.0, d2 = 0.0;
   const y2 = 2.0 * x;

   for (let i = n; i >= 1; i--) {
      const temp = d1;
      d1 = y2 * d1 - d2 + c[i];
      d2 = temp;
   }

   return x * d1 - d2 + c[0];
}

/** @summary Chebyshev0 function
  * @memberof Math */
function Chebyshev0(x, c0) {
   return c0;
}

/** @summary Chebyshev1 function
  * @memberof Math */
function Chebyshev1(x, c0, c1) {
   return c0 + c1*x;
}

/** @summary Chebyshev2 function
  * @memberof Math */
function Chebyshev2(x, c0, c1, c2) {
   return c0 + c1*x + c2*(2.0*x*x - 1.0);
}

/** @summary Chebyshev3 function
  * @memberof Math */
function Chebyshev3(x, ...args) {
   return ChebyshevN(3, x, args);
}

/** @summary Chebyshev4 function
  * @memberof Math */
function Chebyshev4(x, ...args) {
   return ChebyshevN(4, x, args);
}

/** @summary Chebyshev5 function
  * @memberof Math */
function Chebyshev5(x, ...args) {
   return ChebyshevN(5, x, args);
}

/** @summary Chebyshev6 function
  * @memberof Math */
function Chebyshev6(x, ...args) {
   return ChebyshevN(6, x, args);
}

/** @summary Chebyshev7 function
  * @memberof Math */
function Chebyshev7(x, ...args) {
   return ChebyshevN(7, x, args);
}

/** @summary Chebyshev8 function
  * @memberof Math */
function Chebyshev8(x, ...args) {
   return ChebyshevN(8, x, args);
}

/** @summary Chebyshev9 function
  * @memberof Math */
function Chebyshev9(x, ...args) {
   return ChebyshevN(9, x, args);
}

/** @summary Chebyshev10 function
  * @memberof Math */
function Chebyshev10(x, ...args) {
   return ChebyshevN(10, x, args);
}

// =========================================================================

/** @summary Caluclate ClopperPearson
  * @memberof Math */
function eff_ClopperPearson(total,passed,level,bUpper) {
   const alpha = (1.0 - level) / 2;
   if (bUpper)
      return ((passed == total) ? 1.0 : beta_quantile(1 - alpha,passed + 1,total-passed));

   return ((passed == 0) ? 0.0 : beta_quantile(alpha,passed,total-passed+1.0));
}

/** @summary Caluclate normal
  * @memberof Math */
function eff_Normal(total,passed,level,bUpper) {
   if (total == 0) return bUpper ? 1 : 0;

   const alpha = (1.0 - level)/2,
         average = passed / total,
         sigma = Math.sqrt(average * (1 - average) / total),
         delta = normal_quantile(1 - alpha, sigma);

   if (bUpper)
      return ((average + delta) > 1) ? 1.0 : (average + delta);

   return ((average - delta) < 0) ? 0.0 : (average - delta);
}

/** @summary Calculates the boundaries for the frequentist Wilson interval
  * @memberof Math */
function eff_Wilson(total,passed,level,bUpper) {
   const alpha = (1.0 - level)/2;
   if (total == 0) return bUpper ? 1 : 0;
   const average = passed / total,
         kappa = normal_quantile(1 - alpha,1),
         mode = (passed + 0.5 * kappa * kappa) / (total + kappa * kappa),
         delta = kappa / (total + kappa*kappa) * Math.sqrt(total * average * (1 - average) + kappa * kappa / 4);

   if (bUpper)
      return ((mode + delta) > 1) ? 1.0 : (mode + delta);

   return ((mode - delta) < 0) ? 0.0 : (mode - delta);
}

/** @summary Calculates the boundaries for the frequentist Agresti-Coull interval
  * @memberof Math */
function eff_AgrestiCoull(total,passed,level,bUpper) {
   const alpha = (1.0 - level)/2,
         kappa = normal_quantile(1 - alpha,1),
         mode = (passed + 0.5 * kappa * kappa) / (total + kappa * kappa),
         delta = kappa * Math.sqrt(mode * (1 - mode) / (total + kappa * kappa));

  if (bUpper)
     return ((mode + delta) > 1) ? 1.0 : (mode + delta);

  return ((mode - delta) < 0) ? 0.0 : (mode - delta);
}

/** @summary Calculates the boundaries using the  mid-P binomial
  * @memberof Math */
function eff_MidPInterval(total,passed,level,bUpper) {
   const alpha = 1. - level, equal_tailed = true, alpha_min = equal_tailed ? alpha/2 : alpha, tol = 1e-9; // tolerance
   let pmin = 0, pmax = 1, p = 0;

   // treat special case for 0<passed<1
   // do a linear interpolation of the upper limit values
   if (passed > 0 && passed < 1) {
      const p0 = eff_MidPInterval(total, 0.0, level, bUpper),
            p1 = eff_MidPInterval(total, 1.0, level, bUpper);
      p = (p1 - p0) * passed + p0;
      return p;
   }

   while (Math.abs(pmax - pmin) > tol) {
      p = (pmin + pmax)/2;
      // double v = 0.5 * ROOT::Math::binomial_pdf(int(passed), p, int(total));
      // make it work for non integer using the binomial - beta relationship
      let v = 0.5 * beta_pdf(p, passed+1., total-passed+1)/(total+1);
      // if (passed > 0) v += ROOT::Math::binomial_cdf(int(passed - 1), p, int(total));
      // compute the binomial cdf at passed -1
      if ( (passed-1) >= 0) v += beta_cdf_c(p, passed, total-passed+1);

      const vmin = bUpper ? alpha_min : 1.- alpha_min;
      if (v > vmin)
         pmin = p;
      else
         pmax = p;
   }

   return p;
}

/** @summary for a central confidence interval for a Beta distribution
  * @memberof Math */
function eff_Bayesian(total,passed,level,bUpper,alpha,beta) {
   const a = passed + alpha,
          b = total - passed + beta;
   if (bUpper) {
      if ((a > 0) && (b > 0))
         return beta_quantile((1+level)/2,a,b);
      else
         return 1;
   } else {
      if ((a > 0) && (b > 0))
         return beta_quantile((1-level)/2,a,b);
      else
         return 0;
   }
}

/** @summary Return function to calculate boundary of TEfficiency
  * @memberof Math */
function getTEfficiencyBoundaryFunc(option, isbayessian) {
   const kFCP = 0,       // Clopper-Pearson interval (recommended by PDG)
          kFNormal = 1,   // Normal approximation
          kFWilson = 2,   // Wilson interval
          kFAC = 3,       // Agresti-Coull interval
          kFFC = 4,       // Feldman-Cousins interval, too complicated for JavaScript
          // kBJeffrey = 5,  // Jeffrey interval (Prior ~ Beta(0.5,0.5)
          // kBUniform = 6,  // Prior ~ Uniform = Beta(1,1)
          // kBBayesian = 7, // User specified Prior ~ Beta(fBeta_alpha,fBeta_beta)
          kMidP = 8;      // Mid-P Lancaster interval

   if (isbayessian)
      return eff_Bayesian;

   switch (option) {
      case kFCP: return eff_ClopperPearson;
      case kFNormal: return eff_Normal;
      case kFWilson: return eff_Wilson;
      case kFAC: return eff_AgrestiCoull;
      case kFFC: console.log('Feldman-Cousins interval kFFC not supported; using kFCP'); return eff_ClopperPearson;
      case kMidP: return eff_MidPInterval;
      // case kBJeffrey:
      // case kBUniform:
      // case kBBayesian: return eff_ClopperPearson;
   }
   console.log(`Not recognized stat option ${option}, using kFCP`);
   return eff_ClopperPearson;
}

export { gamma, gamma as tgamma, gamma as Gamma,
   Polynomialeval, Polynomial1eval, stirf,
   gamma_pdf, ndtri, normal_quantile, normal_quantile_c, lognormal_cdf_c, lognormal_cdf,
   igami, igamc, igam, lgam, lgamma, erfc, erf,
   beta_pdf, inc_beta, BetaIncomplete,
   pseries, incbet, incbi, beta_quantile, chisquared_cdf_c,
   beta, inc_gamma, inc_gamma_c, landau_pdf, beta_cdf_c, Landau,
   fdistribution_pdf, fdistribution_pdf as FDist,
   fdistribution_cdf, fdistribution_cdf as FDistI,
   fdistribution_cdf_c,
   normal_cdf_c, normal_cdf_c as gaussian_cdf_c,
   normal_cdf, normal_cdf as gaussian_cdf,
   lognormal_pdf, normal_pdf, crystalball_function, crystalball_pdf, crystalball_cdf, crystalball_cdf_c,

   chisquared_cdf, gamma_quantile_c, gamma_quantile, breitwigner_cdf_c, breitwigner_cdf,
   cauchy_cdf_c, cauchy_cdf, cauchy_pdf, gaussian_pdf,
   tdistribution_cdf_c, tdistribution_cdf, tdistribution_pdf, exponential_cdf_c, exponential_cdf, chisquared_pdf,
   Beta, GammaDist, LaplaceDist, LaplaceDistI, LogNormal, Student, StudentI,
   gaus, gausn, gausxy, expo,
   Prob, Sq, Gaus, BreitWigner, BetaDist, BetaDistI, landau, landaun,

   ChebyshevN, Chebyshev0, Chebyshev1, Chebyshev2, Chebyshev3, Chebyshev4,
   Chebyshev5, Chebyshev6, Chebyshev7, Chebyshev8, Chebyshev9, Chebyshev10,

   getTEfficiencyBoundaryFunc };