/*
 * special_functions - special functions (e.g.: gamma, zeta, psi)
 *
 * Copyright (C) 2013 Christoph Zurnieden
 *
 * special_functions is open software; you can redistribute it and/or modify
 * it under the terms of the version 2.1 of the GNU Lesser General Public
 * License as published by the Free Software Foundation.
 *
 * special_functions is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
 * General Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 *
 * @(#) $Revision: 30.4 $
 * @(#) $Id: specialfunctions.cal,v 30.4 2013/08/11 08:41:38 chongo Exp $
 * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/specialfunctions.cal,v $
 *
 * Under source code control:	2013/08/11 01:31:28
 * File existed as early as:	2013
 */


/*
 * hide internal function from resource debugging
 */
static resource_debug_level;
resource_debug_level = config("resource_debug", 0);


/*
  get dependencies
*/
read -once zeta2;


/*
  zeta2(x,s) is in the extra file "zeta2.cal" because of a different license:
  GPL instead of LGPL.
*/
define zeta(s)
{
  /* Not the best way, I'm afraid, but a way. */
  return hurwitzzeta(s,1);
}

define psi0(z){
  local i k m x y eps_digits eps ret;

  /*
    One can use the Stirling series, too, which might be faster for some
    values. The series used here converges very fast but needs a lot of
    bernoulli numbers which are quite expensive to compute.
  */

  eps = epsilon();
  epsilon(eps*1e-3);
  if(isint(z) && z<=0)return  newerror("psi(z); z is a negative integer or 0");
  if(re(z) < 0){
   return ( pi() * cot(pi() * (0-z) ) ) + psi0(1-z);
  }
  eps_digits = digits(1/epsilon());
  /*
    R.W. Gosper has r = .41 as the relation, empirical tests showed that
    for d<100 r = .375 is sufficient and r = .396 for d<200.
    It does not save much time but for a long series even a little
    can grow large.
   */
  if(eps_digits <  100)
    k = 2 * ceil(.375 * eps_digits);
  else if(eps_digits <  200)
    k = 2 * ceil(.395 * eps_digits);
  else
    k = 2 * ceil(11/27 * eps_digits);
  m = 0;
  y = (z+k)^2;
  x = 0.0;
  /*
    There is a chance to speed up the first partial sum with binary splitting.
    The second partial sum is dominated by the calculation of the Bernoulli
    numbers but can profit from binary splitting when the Bernoulli numbers are
    already cached.
   */
  for(i = 1; i <= (k/2); i++){
    m = 1 / ( z + 2*i - 1) + 1 / ( z + 2*i -2 )  + m;
    x = ( x + bernoulli(k-2*i+2) / (k-2*i+2) ) / y ;
  }
  ret = ln(z+k) - 1 / ( 2 * (z+k) ) - x - m ;
  epsilon(eps);
  return ret;
}

define psi(z)
{
  return psi0(z);
}

define polygamma(m,z)
{
  /*
    TODO:
    http://functions.wolfram.com/GammaBetaErf/PolyGamma2/16/01/01/0002/
  */
  if(!isint(m))return newerror("polygamma(m,z); m not an integer");
  if( m<0 )return newerror("polygamma(m,z); m is < 0");
  /*
    Reflection formula not implemented yet, needs cot-differentiation
    http://functions.wolfram.com/ElementaryFunctions/Cot/20/02/0003/
    which is not implemented yet.
   */
  if( m == 0){
    return psi0(z);
  }
  /*
    use factorial for m a (small) integer?
    use lngamma for m large?
  */
  if(isodd(m+1)){
    return (-1)* gamma(m+1) * hurwitzzeta(m+1,z)
  }
  return gamma(m+1) * hurwitzzeta(m+1,z);
}

/*
  Cache for the variable independent coefficients in the sum for the
  Gamma-function.
*/
static __CZ__Ck;
/*
  Log-Gamma function for Re(z) > 0.
*/
define __CZ__lngammarp(z)
{
  local epsilon accuracy a  factrl  sum n ret holds_enough term;

  epsilon  = epsilon();
  accuracy = digits(int(1/epsilon)) + 3;

  epsilon(1e-18);
  a = ceil(1.252850440912568095810522965 * accuracy);

  epsilon(epsilon*10^(-(digits(1/epsilon)//2)));

  holds_enough = 0;

  if(size( __CZ__Ck) != a) {
     __CZ__Ck = mat[a];
     holds_enough = 1;
  }

  factrl = 1.0;

  __CZ__Ck[0] = sqrt(2*pi()); /* c_0*/
  for(n = 1; n < a; n++){
    if(holds_enough == 1){
      __CZ__Ck[n] = (a - n)^(n - 0.5)  * exp(a - n)  /factrl;
    factrl *= -n
    }
  }
  sum = __CZ__Ck[0];
  for (n = 1; n < a; n++){
    sum += __CZ__Ck[n]/(z+n);
  }

  ret = ln(sum)+(-(z+a)) + ln(z+a)*( z+0.5);
  ret = ret-ln(z);

   /*
      Will take some time for large im(z) but almost all time is spend above
      in that case.
   */
   if(im(ret))
     ret = re(ret) + ln( exp( im(ret) *1i ) );

  epsilon(epsilon);
  return ret;
}

/* Simple lngamma with low precision*/
define __CZ__lngamma_lanczos(z){
   local a k g zghalf lanczos;
   mat lanczos[15] = {
                  9.9999999999999709182e-1,
                  5.7156235665862923516e1,
                 -5.9597960355475491248e1,
                  1.4136097974741747173e1,
                 -4.9191381609762019978e-1,
                  3.3994649984811888638e-5,
                  4.6523628927048576010e-5,
                 -9.8374475304879566105e-5,
                  1.5808870322491249322e-4,
                 -2.1026444172410489480e-4,
                  2.1743961811521265523e-4,
                 -1.6431810653676390482e-4,
                  8.4418223983852751308e-5,
                 -2.6190838401581411237e-5,
                  3.6899182659531626821e-6
                     };
   g = 607/128;
   z = z-1;
   a = 0;
   for(k = 12; k >= 1; k--){
     a += lanczos[k]/(z+k);
   }
   a += lanczos[0];
   zghalf = z + g + .5;
   return ( ln(sqrt(2*pi())) + ln(a)  -zghalf  ) + (z+.5)*ln( zghalf );
}

/* Prints the Spouge coefficients actually in use. */
define __CZ__print__CZ__Ck(){
  local k;
  if(size(__CZ__Ck) <=1){
    __CZ__lngammarp(2.2-2.2i);
  }
  for(k=0;k<size(__CZ__Ck);k++){
    print __CZ__Ck[k];
  }
}

/*Kronecker delta function */
define kroneckerdelta(i,j){
    if(isnull(j)){
      if(i==0) return 1;
      else return 0;
    }
    if(i!=j) return 0;
    else return 1;
}

/* (Pre-)Computed terms of the Stirling series */
static __CZ__precomp_stirling;
/* Number of terms in mat[0,0], precision in mat[0,1] with |z| >= mat[0,2]*/
static __CZ__stirling_params;

define __CZ__lngstirling(z,n){
  local k head sum z2 bernterm zz;
  head = (z-1/2)*ln(z)-z+(ln(2*pi())/2);
  sum = 0;
  bernterm=0;
  zz = z;
  z2 = z^2;

  if(size(__CZ__precomp_stirling)<n){
    __CZ__precomp_stirling = mat[n];
    for(k=1;k<=n;k++){
      bernterm = bernoulli(2*k)/(   2*k*(2*k-1 ) );
      __CZ__precomp_stirling[k-1] = bernterm;
      sum +=  bernterm /zz;
      zz *= z2;
    }
  }
  else{
    for(k=1;k<=n;k++){
      bernterm = __CZ__precomp_stirling[k-1];
      sum +=  bernterm /zz;
      zz *= z2;
    }
  }
  return head + sum;
}

/* Compute error for Stirling series for "z" with "k" terms*/
define R(z,k){
  local a b ab stirlingterm;
  stirlingterm = bernoulli(2*k)/( 2*k*(2*k-1)*z^(2*k));
  a = abs( stirlingterm );
  b = (cos(.5*arg(z)^(2*k)));
  ab = a/b;
  /* a/b might round to zero */
  if( abs(ab) < epsilon()){
    return digits(1/epsilon());
  }
  return digits(1/(a/b));
  ##return a/b;
}

/*Low precision lngamma(z) for testing the branch correction */
define lngammasmall(z){
  local ret eps flag corr;
  flag = 0;
  if(isint(z)){
    if(z <=0)
      return newerror("gamma(z): z is a negative integer");
    else{
      /* may hold up accuracy a bit longer, but YMMV */
      if(z < 20)
        return ln((z-1)!);
      else
        return __CZ__lngamma_lanczos(z);
    }
  }
  else{
    if(re(z)<0){
      if(im(z) && im(z)<0){
        z = conj(z);
        flag = 1;
      }

      ret = ln( pi() / sin(pi() *z ) ) -__CZ__lngamma_lanczos(1-z);

      if(!im(z)){
        if(abs(z) <= 1/2)
          ret = re(ret) - pi()*1i;
        else if( isodd(floor(abs(re(z)))) ){
          ret = re(ret) + ( ceil(abs(z))   * pi() * 1i);
        }
        else if( iseven(floor(abs(re(z)))) ){
          /* < n+1/2 */
          if(iseven(floor(abs(z)))){
            ret = re(ret) + ( ceil(abs(z)-1/2 -1 )   * pi() * 1i);
          }
          else{
            ret = re(ret) + ( ceil(abs(z) -1/2 )   * pi() * 1i);
          }
        }
      }
      else{
        corr = ceil( re(z)/2 -3/4 - kroneckerdelta(im(z))/4);
        ret = ret +2*(corr *pi() )*1i;
      }
      if(flag == 1)
         ret = conj(ret);
      return ret;
    }
    ret = (__CZ__lngamma_lanczos(z));
    return ret;
  }
}


/*
  The logarithmic gamma function lngamma(z)
  It has a different principal branch than log(gamma(z))
*/
define lngamma(z)
{
  local ret eps  flag  increasedby Z k tmp tmp2 decdigits;
  local corr d37 d52 termcount;
  /* z \in \mathbf{Z}*/
  if(isint(z)){
    /*The gamma-function has poles at zero and the negative integers*/
    if(z <=0)
      return newerror("lngamma(z): z is a negative integer");
    else{
      /* may hold up accuracy a bit longer, but YMMV */
      if(z < 20)
        return ln( (z-1)! );
      else
        return  __CZ__lngammarp(z);
    }
  }
  else{/*if(isint(z))*/
    if(re(z) > 0){
      flag = 0;
      eps = epsilon();
      epsilon(eps*1e-3);

      /* Compute values on the real line with Spouge's algorithm*/
      if(!im(z)){#
        ret = __CZ__lngammarp(z);
        epsilon(eps);
        return ret;
      }
      /* Do the rest with the Stirling series.*/
      /* This code repeats down under. Booh! */
      /* Make it a positive im(z) */
      if(im(z)<0){
        z = conj(z);
        flag = 1;
      }
      /* Evaluate the number of terms needed */
      decdigits = floor( digits(1/eps) );
      /* 20 dec. digits is the default in calc(?)*/
      if(decdigits <= 21){
        /* set 20 as the minimum */
        epsilon(1e-22);
        increasedby = 0;
        /* inflate z */
        Z=z;
        while(abs(z) < 10){
          z++;
          increasedby++;
        }

        ret = __CZ__lngstirling(z,11);
        /* deflate z */
        if(increasedby > 1){
          for(k=0;k<increasedby;k++){
            ret = ret - ln(Z+(k));
          }
        }
        ## ret = ret + 2*pi()*1i*ceil(re(x+1)/2-3/4 );
        /* Undo conjugate if one has been applied */
        if(flag == 1)
          ret = conj(ret);
        epsilon(eps);
        return ret;
      }/* if(decdigits <= 20) */
      else{
        /*
           Compute the number of terms for the Stirling-sum first.
           The error is bound by
                                            |T(k)|
              |R(k,z)| <=  ----------------------------------------
                                             2 k                2 k
                              2 k (2 k - 1) z     cos(0.5 arg(z)   )
             with
                                  1 - 2 k
                          B(2 k) z
              T(k) =  ----------------------
                          2 k (2 k - 1)
            the terms of the Stirling-sum and B(n) the n-th Bernoulli number.

            We do know the current precision and the result of the helper
            function __CZ__R(z,k) returns the same number as "decdigits" if
            fed with the correct "z" and "n".
            The largest error is near the real line, so we try Z = re(1-z).
            If Z is small, say < 20 we increment it until it is.
            Computing Bernoulli numbers is expensive but they are small on
            average, get cached and because of the method used, computing
            any Bernoulli number > 0 produces the rest as a by-product.
            Nevertheless we could use a faster asymptotic approximation.

            We increment the number of terms "k" in R(z,k) until the result
            quits incrementing (it may even get smaller!). If the result is
            still smaller than the current precision we increment "z" with
            fixed "k" untill the result quits incrementing.
            The results, the current precision, abs(re(z)) and "k" are kept.

            BTW: incrementing the number of terms might be more costly than
            incrementing "z" -- computing large Bernoulli numbers vs.
            computing a large number of complex logarithms is a fight with
            a hard to know result -- and that the series isn't convergent
            is of not much help either. E.g:
              R(25,68)   = 71 max
              R(50,55)   = 101
              R(50,145)  = 140 max
              R(60,170)  = 167 max
              R(70,209)  = 195 max
              R(75,173)  = 200 max
              R(80,147)  = 200 max
              R(90,124)  = 200 max
              R(100,111) = 200 max
            Bernoulli(222) has a denominator of 9388 with a 254 digit
            numerator. Computing up to 100 complex logarithms on the
            other side ...

            D.E.G. Hare has found the bounds
               |im(z)| > .37d or re(z) >= 0 and |z| > .52d
            to be usefull to compute "z" to d digits accuracy. The numbers
            correspond to the table above.

            To avoid repeated expensive computation, the result is cached
            together with the current precision. It might be a good idea
            to keep it more permanently in a config-file?
         */
        d37 = decdigits * .37;
        d52 = decdigits * .52;
        termcount = ceil(d52);
        if(abs(z) >= d52){
         if(abs(im(z))>= d37 )
           termcount = ceil(d37);
         else
           termcount = ceil(d52);
        }

        Z=z;
        increasedby = 0;
        /* inflate z */
        if( abs(im(z))>= d37){
          while(abs(z) < d52+1){
            z++;
            increasedby++;
          }
        }
        else{
          tmp = R(z,termcount);
          tmp2 = tmp;
          while(tmp2 < decdigits){
            z++;
            increasedby++;
            tmp2 = R(z,termcount);
            if(tmp2 < tmp)
              return newerror("lngamma(1): something happend that "
	      		      "should not have happend");
          }
        }

        corr = ceil( re(z)/2 -3/4 - kroneckerdelta(im(z))/4);

        ret = __CZ__lngstirling(z,termcount);

        /* deflate z */
        if(increasedby > 1){
          for(k=0;k<increasedby;k++){
            ret = ret - ln(Z+(k));
          }
        }
        /* Undo conjugate if one has been applied */
        if(flag == 1)
          ret = conj(ret);
        epsilon(eps);
        return ret;
      }/* if(decdigits <= 20){...} else */
    }/* if(re(z) > 0) */
    else{/* re(z)<0 */
      eps = epsilon();
      epsilon(eps*1e-3);

      /* Use Spouge's approximation on the real line */
      if(!im(z)){
        /* reflection */
        ret = ln( pi() / sin(pi() *z ) ) -  __CZ__lngammarp(1-z);
        /* it is log(gamma) and im(log(even(-x))) = k\pi, therefore: */
        if(abs(z) <= 1/2)
          ret = re(ret) - pi()*1i;
        else if( isodd(floor(abs(re(z)))) ){
          ret = re(ret) + ( ceil(abs(z))   * pi() * 1i);
        }
        else if( iseven(floor(abs(re(z)))) ){
          /* < n+1/2 */
          if(iseven(floor(abs(z)))){
            ret = re(ret) + ( ceil(abs(z)-1/2 -1 )   * pi() * 1i);
          }
          else{
            ret = re(ret) + ( ceil(abs(z) -1/2 )   * pi() * 1i);
          }
        }
        epsilon(eps);
        return ret;
      }/*if(!im(z))*/
      /* Use Stirlinsg approximation for the rest of the complex plane */
      else{
        /* Make it a positive im(z) */
        if(im(z)<0){
          z = conj(z);
          flag = 1;
        }
        /* Evaluate the number of terms needed */
        decdigits = floor( digits(1/eps) );
        /*
          Evaluate the correction term for the imaginary part needed because
          of the reflection.
          See
             D. E. G. Hare, "Computing the Principal Branch of log-Gamma",
             Journal of Algorithms 25(2):221-236 (1997)
             http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.84.2063
         */

        /* 20 dec. digits is the default in calc(?)*/
        if(decdigits <= 21){
          /* set 20 as the minimum */
          epsilon(1e-22);
          termcount = 11;
          increasedby = 0;
          Z=z;
          /* inflate z */
          if( im(z)>= digits(1/epsilon()) * .37){
            while(abs(1-z) < 10){
              /* making z more negative makes 1-z more positive. */
              z--;
              increasedby++;
            }
          }
          else{
            tmp = R(1-z,termcount);
            tmp2 = tmp;
            while(tmp2 < 21){
              z--;
              increasedby++;
              tmp2 = R(1-z,termcount);
              if(tmp2 < tmp)
                return newerror("lngamma(1): something happend "
				"that should not have happend");
            }
          }

          corr = ceil( re(z)/2 -3/4 - kroneckerdelta(im(z))/4);

          /* reflection */
          ret = ln( pi() / sin(pi() * z ) ) -  __CZ__lngstirling(1-z,termcount);

          /* deflate z */
          if(increasedby > 0){
            for(k=0;k<increasedby;k++){
              ret = ret + ln(z+(k));
            }
          }
          /* Apply correction term */
          ret = ret +2*(corr *pi() )*1i;
          /* Undo conjugate if one has been applied */
          if(flag == 1)
            ret = conj(ret);

          epsilon(eps);
          return ret;
        }/* if(decdigits <= 20) */
        else{
          d37 = decdigits * .37;
          d52 = decdigits * .52;
          termcount = ceil(d52);
          if(abs(z) >= d52){
             if(abs(im(z))>= d37 )
               termcount = ceil(d37);
             else
               termcount = ceil(d52);
          }
          increasedby = 0;
          ##print "Z 1: ",z;
          Z=z;
          /* inflate z */
          if( abs(im(z))>= d37){
            while(abs(1-z) < d52+1){
              /* making z more negative makes 1-z more positive. */
              z--;
              increasedby++;
            }
          }
          else{
            tmp = R(1-z,termcount);
            tmp2 = tmp;
            while(tmp2 < decdigits){
              z--;
              increasedby++;
              tmp2 = R(1-z,termcount);
              if(tmp2 < tmp)
                return newerror("lngamma(1): something happend that "
				"should not have happend");
            }
          }
          corr = ceil( re(z)/2 -3/4 - kroneckerdelta(im(z))/4);
          /* reflection */
          ret = ln( pi() / sin(pi() *(z) ) ) -
	        __CZ__lngstirling(1-z,termcount);
          /* deflate z */
          if(increasedby > 0){
            for(k=0;k<increasedby;k++){
              ret = ret + ln(z+(k));
            }
          }
          /* Apply correction term */
          ret = ret +2*((corr) *pi() )*1i;
          /* Undo conjugate if one has been applied */
          if(flag == 1)
            ret = conj(ret);
          epsilon(eps);
          return ret;
        }/* if(decdigits <= 20){...} else */
      }/*if(!im(z)){...} else*/
    }/* if(re(z) > 0){...} else */
  }/*if(isint(z)){...} else*/
}

/* Warning about large values? */
define gamma(z)
{

  /* exp(log(gamma(z))) = exp(lngamma(z)), so use Spouge here?*/
  local ret eps;
  if(isint(z)){
    if(z <=0)
      return newerror("gamma(z): z is a negative integer");
    else{
      /* may hold up accuracy a bit longer, but YMMV */
      if(z < 20)
        return (z-1)!*1.0;
      else{
        eps = epsilon(epsilon()*1e-2);
        ret = exp(lngamma(z));
        epsilon(eps);
        return ret;
      }
    }
  }
  else{
    eps = epsilon(epsilon()*1e-2);
    ret = exp(lngamma(z));
    epsilon(eps);
    return ret;
  }
}

define __CZ__harmonicF( a, b ,s)
{
  local c;
  if( b == a) return s ;
  if( b-a > 1){
    c= (b + a) >> 1;
    return( __CZ__harmonicF(a, c,1/a) + __CZ__harmonicF(c+1, b,1/b));
  }
  return  (1/a+1/b);
}

define harmonic(limit)
{
    if( !isint(limit) )
      return newerror("harmonic(limit): limit is not an integer");
    if(  limit <= 0   )
      return newerror("harmonic(limit): limit is <=0");
    /* The binary splitting algorithm returns 0 for f(1,1,0) */
    if( limit == 1    ) return 1;
    return __CZ__harmonicF( 1, limit ,0);
}

/*  lower incomplete gamma function */

/* lower
   for z <= 1.1
 */
define __CZ__gammainc_series_lower(a,z){
  local k ret tmp eps fact;
  ret = 0;
  k=0;
  tmp=1;
  fact = 1;
  eps = epsilon();
  while(abs(tmp-ret) > eps){
    tmp = ret;
    ret += (z^(k+a))/(  (a+k)*fact );
    k++;
    fact *= -k;
  }
  return gamma(a)-ret;
}

/* lower
   for z > 1.1
 */
define __CZ__gammainc_cf(a,z,n){
  local ret k  num1 denom1  num2 denom2 ;
  ret = 0;
  for(k=n+1;k>1;k--){
    ret = ((1-k)*(k-1-a))/(2*k-1+z-a+ret);
  }
  return ((z^a*exp(-z))/(1+z-a+ret));
}

/* G(n,z) lower*/
define __CZ__gammainc_integer_a(a,z){
  local k sum fact zz;
  for(k=0;k<a;k++){
    sum += (z^k)/(k!);
  }
  return  (a-1)!*exp(-z)*sum;
}

/*
  P = precision in dec digits
  n = 1,2,3...
  a,z => G(a,z)
*/
define __CZ__endcf(n,a,z,P){
  local ret;

  ret = P*ln(10)+ln(4*pi()*sqrt(n)) + re(z+(3/2-a)*ln(z)-lngamma(1-a));
  ret = ret /( sqrt( 8*(abs(z)+re(z)) )  );
  return ret^2;

}

/* lower incomplete gamma function */ define gammainc(a,z){
  local ret nterms eps epsilon tmp_before tmp_after n A B C sum k;
  if(z == 0)
    return 1;
  if(isint(a)){
    if(a>0){
      if(a==1) return exp(-z); return __CZ__gammainc_integer_a(a,z);
    } else{
      if(a==0){
        return -expoint(-z)+1/2*( ln(-z) - ln(-1/z) ) -ln(z);
      } else if(a==-1){
        return  expoint(-z)+1/2*(ln(-1/z)-ln(-z))+ln(z)+(exp(-z)/z);
      } else{
        A = (-1)^((-a)-1)/((-a)!); B =
        expoint(-z)-1/2*(ln(-z)-ln(1/-z))+ln(z); C = exp(-z); sum =0;
        for(k=1;k<-a;k++){
          sum += (z^(k-a-1))/( fallingfactorial(-a,k)  );
        } return A * B - C *sum;
      }
    }
  } if(re(z)<=1.1|| re(z) < a+1){##print "series";
    eps = epsilon(epsilon()*1e-10); ret =
    __CZ__gammainc_series_lower(a,z); epsilon(eps); return ret;
  } else{##print "cf";
    eps = epsilon(epsilon()*1e-10); if(abs(exp(-z)) <= eps) return 0;
    tmp_before = 0; tmp_after  = 1; n = 1; while(ceil(tmp_before) !=
    ceil(tmp_after)){
      tmp_before = tmp_after; tmp_after =
      __CZ__endcf(n++,a,z,digits(1/eps)); /* a still quite arbitrary
      limit */ if(n > 10){
	return newerror("gammainc: evaluating limit for continued "
			"fraction does not converge");
      }
    } ret  = __CZ__gammainc_cf(a,z,ceil(tmp_after)); epsilon(eps);
    return  ret;
  }
}

define heavisidestep(x){
   return (1+sign(x))/2;
}

define NUMBER_POSITIVE_INFINITY(){ return 1/epsilon();}

define NUMBER_NEGATIVE_INFINITY(){ return -(1/epsilon());}

static TRUE = 1;
static FALSE = 0;

define g(prec){
  local eps ret;
  if(!isnull(prec)){
    eps = epsilon(prec);
    ret = -psi(1);
    epsilon(eps);
    return ret;
  }
  return -psi(1);
}

define __CZ__series_converged(new,old,max){
  local eps;
  if(isnull(max))
    eps = epsilon();
  else
    eps = max;
  if(    abs(re(new - old)) <= eps * abs(re(new))
      && abs(im(new - old)) <= eps * abs(im(new)))
    return TRUE;
  return FALSE;
}

define __CZ__ei_power(z){
  local tmp ei k old;
  ei = g() + ln(abs(z)) + sgn(im(z)) * 1i * abs(arg(z));
  ##ei = g() + ln(z) -1i*pi()*floor( (arg(z)+pi()) / (2*pi()) );
  tmp = 1;
  k = 1;
  while(k){
    tmp *= z / k;
    old = ei;
    ei += tmp / k;
    if (__CZ__series_converged(ei,old)) break;
    k++;
  }
  return ei;
}

define __CZ__ei_asymp(z){
  local ei old tmp k;
  ei = sgn(im(z)) * 1i * pi();
  tmp = exp(z) / z;
  for(k=1; k<=floor(abs(z))+1; k++){
    old = ei;
    ei += tmp;
    if (__CZ__series_converged(ei, old)) return ei;
    tmp *= k / z;
  }
  return newerror("expoint: asymptotic series does not converge");
}

define __CZ__ei_cf(z){
  local ei c d k old;
  ei = sgn(im(z)) * 1i * pi();
  if(ei != 0){
    c = 1 / ei;
    d = 0;
    c = 1 / (1 - z - exp(z) * c);
    d = 1 / (1 - z - exp(z) * d);
    ei *= d / c;
  }
  else{
    c = NUMBER_POSITIVE_INFINITY();
    d = 0;
    c = 0;
    d = 1 / (1 - z - exp(z) * d);
    ei = d * (- exp(z));
  }
  k = 1;
  while(1){
     c = 1 / (2 * k + 1 - z - k * k * c);
     d = 1 / (2 * k + 1 - z - k * k * d);
     old = ei;
     ei *= d / c;
     if (__CZ__series_converged(ei, old)) break;
     k++;
  }
  return ei;
}

define expoint(z){
  local ei eps ret;
  eps=epsilon(epsilon()*1e-5);
  if(abs(z) >= NUMBER_POSITIVE_INFINITY()){
    if (config("user_debug") > 0) {
     print "expoint: abs(z) > +inf";
    }
    ret = sgn(im(z)) * 1i * pi() + exp(z) / z;
    epsilon(eps);
    return ret;
  }
  if(abs(z) > 2 - 1.035 * log(epsilon())){
    if (config("user_debug") > 0) {
     print "expoint: using asymptotic series";
    }
    ei = __CZ__ei_asymp(z);
    if (!iserror(ei)){
        ret = ei;
        epsilon(eps);
        return ret;
    }
  }
  if(abs(z) > 1 && (re(z) < 0 || abs(im(z)) > 1)){
    if (config("user_debug") > 0) {
      print "expoint: using continued fraction";
    }
    ret =  __CZ__ei_cf(z);
    epsilon(eps);
    return ret;
  }
  if(abs(z) > 0){
    if (config("user_debug") > 0) {
      print "expoint: using power series";
    }
    ret =  __CZ__ei_power(z);
    epsilon(eps);
    return ret;
  }
  if(abs(z) == 0){
    if (config("user_debug") > 0) {
     print "expoint: abs(z) = zero ";
    }
    epsilon(eps);
    return  NUMBER_NEGATIVE_INFINITY();
  }
}

define erf(z){
  return sqrt(z^2)/z * ( 1-1/sqrt(pi()) *gammainc(1/2,z^2)   );
}

define erfc(z){
  return 1-erf(z);
}

define erfi(z){
  return -1i*erf(1i*z);
}

define faddeeva(z){
  return exp(-z^2)*erfc(-1i*z);
}

define gammap(a,z){
  return gammainc(a,z)/gamma(a);
}

define gammaq(a,z){
  return 1-gammap(a,z);
}

define lnbeta(a,b){
  local ret eps;
  eps=epsilon(epsilon()*1e-3);
  ret = (lngamma(a)+lngamma(b))-lngamma(a+b);
  epsilon(eps);
  return ret;
}

define beta(a,b){
  return exp(lnbeta(a,b));
}

define __CZ__ibetacf_a(a,b,z,n){
  local A B m places;
  if(n==1) return 1;
  m=n-1;
  places = highbit(1 + int(1/epsilon())) + 1;
  A = bround((a+m-1) * (a+b+m-1) * m * (b-m) * z^2,places++);
  B = bround((a+2*(m)-1)^2,places++);
  return A/B;
}

define __CZ__ibetacf_b(a,b,z,n){
  local A B m places;
  places = highbit(1 + int(1/epsilon())) + 1;
  m=n-1;
  A = bround((m*(b-m)*z)/(a+2*m-1),places++);
  B = bround(( (a+m) * (a-(a+b)*z+1+m*(2-z)) )/(a+2*m+1),places++);
  return m+A+B;
}

/* Didonato-Morris */
define __CZ__ibeta_cf_var_dm(a,b,z,max){
  local m  f c d  check  del h qab qam qap eps places;

  eps= epsilon();

  if(isnull(max)) max = 100;
  places = highbit(1 + int(1/epsilon())) + 1;
  f = eps;
  c = f;
  d = 0;
  for(m=1;m<=max;m++){
    d =bround( __CZ__ibetacf_b(a,b,z,m)+__CZ__ibetacf_a(a,b,z,m)*d,places++);
    if(abs(d)<eps ) d=eps;
    c =bround( __CZ__ibetacf_b(a,b,z,m)+__CZ__ibetacf_a(a,b,z,m)/c,places++);
    if(abs(c)<eps ) c=eps;
    d=1/d;
    check = c*d;
    f = f*check;
    if(abs(check-1)<eps)break;
  }
  if(m > max) return newerror("ibeta: continous fraction does not converge");
  return f;
}

define betainc_complex(z,a,b){
  local factor ret eps cf sum k N places tmp tmp2;

  if(z == 0){
    if(re(a) > 0) return 0;
    if(re(a) < 0) return newerror("betainc_complex: z == 0 and re(a) < 0");
  }
  if(z == 1){
     if(re(b)>0) return 1;
     else return newerror("betainc_complex: z == 1 and re(b) < 0");
  }
  if(b<=0){
    if(isint(b)) return 0;
    else return newerror("betainc_complex: b <= 0");
  }
  if(z==1/2 && (a==b)){return 1/2;
  }
  ##if(2==1){
  if(isint(a) && isint(b)){
    eps=epsilon(epsilon()*1e-10);
    N = a+b-1;
    sum = 0;
    for(k=a;k<=N;k++){
      tmp  =  ln(z)*k+ln(1-z)*(N-k);
      tmp2 = exp(ln(comb(N,k))+tmp);
      sum += tmp2;
    }
    epsilon(eps);
    return sum
  }
  else if(re(z) <= re((a+1)/(a+b+2))){
    eps=epsilon(epsilon()*1e-10);
    places = highbit(1 + int(1/epsilon())) + 1;
    factor = bround(( ln(z^a*(1-z)^b ) - lnbeta(a,b) ),places);
    cf =bround( __CZ__ibeta_cf_var_dm(a,b,z),places);
    ret = factor + ln(cf);
    if(abs(ret//ln(2)) >= places)
      ret = 0;
    else
      ret = bround(exp(factor + ln(cf)),places);
    epsilon(eps);
    return ret;
  }
  else if( re(z) > re( (a+1)/(a+b+2) ) || re(1-z) < re( (b+1)/(a+b+2) ) ){
    ret = 1 - betainc_complex(1-z,b,a);
  }
  return ret;
}


/******************************************************************************/
/*
  Purpose:

    __CZ__ibetaas63 computes the incomplete Beta function ratio.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    2013-08-03 20:52:05 +0000

  Author:

    Original FORTRAN77 version by KL Majumder, GP Bhattacharjee.
    C version by John Burkardt
    Calc version by Christoph Zurnieden

  Reference:

    KL Majumder, GP Bhattacharjee,
    Algorithm AS 63:
    The incomplete Beta Integral,
    Applied Statistics,
    Volume 22, Number 3, 1973, pages 409-411.

  Parameters:

    Input, x, the argument, between 0 and 1.

    Input, a, b, the parameters, which
    must be positive.


    Output, the value of the incomplete
    Beta function ratio.
*/
define __CZ__ibetaas63(x, a, b,beta){
  local ai betain cx indx ns aa asb bb rx temp term value xx acu places;
  acu = epsilon();

  value = x;
  /* inverse incbeta calculates it already */
  if(isnull(beta))
    beta = lnbeta(a,b);

  if ( a <= 0.0 || b <= 0.0 ){
    return newerror("betainc: domain error: a < 0 and/or b < 0");
  }
  if ( x < 0.0 || 1.0 < x ){
    return newerror("betainc: domain error: x<0 or x>1");
  }
  if ( x == 0.0 || x == 1.0 ){
    return value;
  }
  asb = a + b;
  cx = 1.0 - x;

  if ( a < asb * x ){
    xx = cx;
    cx = x;
    aa = b;
    bb = a;
    indx = 1;
  }
  else{
    xx = x;
    aa = a;
    bb = b;
    indx = 0;
  }

  term = 1.0;
  ai = 1.0;
  value = 1.0;
  ns = floor( bb + cx * asb );

  rx = xx / cx;
  temp = bb - ai;
  if ( ns == 0 ){
    rx = xx;
  }
  places = highbit(1 + int(1/acu)) + 1;
  while(1){
    term = bround(term * temp * rx / ( aa + ai ),places++);
    value = value + term;;
    temp = abs ( term );

    if ( temp <= acu && temp <= abs(acu * value) ){
      value = value * exp ( aa * ln ( xx )
      + ( bb - 1.0 ) * ln ( cx ) - beta ) / aa;

      if ( indx ){
        value = 1.0 - value;
      }
      break;
    }

    ai = ai + 1.0;
    ns = ns - 1;

    if ( 0 <= ns ) {
      temp = bb - ai;
      if ( ns == 0 ) {
        rx = xx;
      }
    }
    else {
      temp = asb;
      asb = asb + 1.0;
    }
  }
  epsilon(acu);
  return value;
}

/*
                    z
                  /
                  [         b - 1  a - 1
   1/beta(a,b) *  I  (1 - t)      t      dt
                  ]
                  /
                 0

*/

define betainc(z,a,b){
  local factor ret eps cf sum k N places tmp tmp2;

  if(im(z) || im(a) || im(b))
   return betainc_complex(z,a,b);

  if(z == 0){
    if(re(a) > 0) return 0;
    if(re(a) < 0) return newerror("betainc: z == 0 and re(a) < 0");
  }
  if(z == 1){
     if(re(b)>0) return 1;
     else return newerror("betainc: z == 1 and re(b) < 0");
  }
  if(b<=0){
    if(isint(b)) return 0;
    else return newerror("betainc: b <= 0");
  }
  if(z==1/2 && a==b){
    return 1/2;
  }
  return __CZ__ibetaas63(z,a,b);

}

define __CZ__erfinvapprox(x){
  local a;
  a =0.147;
  return sgn(x)*sqrt(sqrt((2/(pi()*a)+(ln(1-x^2))/2)^2-(ln(1-x^2))/a)-
    (2/(pi()*a)+(ln(1-x^2))/2));
}

/* complementary inverse errror function, faster at about x < 1-.91
   Henry E. Fettis. "A stable algorithm for computing the inverse error function
   in the 'tail-end' region" Math. Comp., 28:585-587, 1974.
*/
define __CZ__inverffettis(x,n){
   local y sqrtpi oldy k places;
   if (isnull(n))
    n = 205;
   y = erfinvapprox(1-x);
   places = highbit(1 + int(1/epsilon())) + 1;
   sqrtpi = sqrt(pi());
   do
   {
      oldy = y;k++;
      y = bround((ln( __CZ__fettiscf(y,n) / (sqrtpi * x)))^(1/2),places);
   } while( abs(y - oldy)/y > epsilon());
   return y;
}

/* cf for erfc() */
define __CZ__fettiscf(y,n){
  local k t tt r a b ;
  t = 1/y;
  tt = t^2/2;
  for (k=n;k> 0;k--){
 		a = 1;
		b = k*tt;
		r = b / (a + r);
  }
  return t / (1+r);
}

/* inverse errror function, faster at about x<=.91*/
define __CZ__inverfbin(x){
  local places approx flow fhigh eps high low mid fmid epsilon;
  approx = erfinvapprox(x);
  epsilon = epsilon();
  high = approx + 1e-4;
  low  = -1;
  places = highbit(1 + int(1/epsilon)) + 1;
  fhigh = x-erf(high);
  flow  = x-erf(low);
  while(1){
    mid = bround(high - fhigh * (high - low) / (fhigh - flow), places);
    if ((mid == low) || (mid == high))
      places++;
    fmid = x-erf(mid);
    if (abs(fmid) < epsilon)
      return mid;
    if (sgn(fmid) == sgn(flow)) {
      low = mid;
      flow = fmid;
    }
    else {
      high = mid;
      fhigh = fmid;
    }
  }
}

define erfinv(x){
  local ret approx a eps y old places errfunc sqrtpihalf flag k;
  if(x<-1 || x > 1) return newerror("erfinv: input out of domain (-1<=x<=1)");
  if(x == 0) return 0;
  if(x == -1) return NUMBER_NEGATIVE_INFINITY();
  if(x == +1) return NUMBER_POSITIVE_INFINITY();

  if(x<0){
    x = -x;
    flag = 1;
  }
  /* No need for full pecision */
  eps=epsilon(1e-20);
  if(eps >= 1e-40){
    /* Winitzki, Sergei (6 February 2008). "A handy approximation for the error
                                          function and its inverse"*/
    a = 0.147;
    y = sgn(x)*sqrt(sqrt((2/(pi()*a)
                  +(ln(1-x^2))/2)^2
                  -(ln(1-x^2))/a)
                  -(2/(pi()*a)+(ln(1-x^2))/2));

  }
  else {
    /* 20 digits instead of 5 */
    if(x <= .91)
      y = __CZ__inverfbin(x);
    else
      y = __CZ__inverffettis(1-x);

    if(eps <= 1e-20){
      epsilon(eps);
      return y;
    }
  }
  epsilon(eps);
  /* binary digits in number (here: number = epsilon()) */
  places = highbit(1 + int(1/eps)) + 1;
  sqrtpihalf = 2/sqrt(pi());
  k = 0;
  /*
     Do some Newton-Raphson steps to reach final accuracy.
     Only a couple of steps are necessary but calculating the error function at
     higher precision is quite costly;
   */
  do{
    old = y;
    errfunc = bround( erf(y),places);
    if( abs(errfunc-x) <= eps ) break;
    y = bround(y-( errfunc  -x) / ( sqrtpihalf * exp(-y^2)),places);
    k++;
  }while(1);
  /*
    This is not really necessary but e.g:
    ; epsilon(1e-50)
	    0.00000000000000000000000000000000000000000000000001
    ; erfinv(.9999999999999999999999999999999)
	    8.28769266865549025938
    ; erfinv(.999999999999999999999999999999)
	    8.14861622316986460738453487549552168842204512959346
    ; erf(8.28769266865549025938)
	    0.99999999999999999999999999999990000000000000000000
    ; erf(8.14861622316986460738453487549552168842204512959346)
	    0.99999999999999999999999999999900000000000000000000
    The precision "looks too short".
  */
  if(k == 0)
    y = bround(y-( errfunc  -x) / ( sqrtpihalf * exp(-y^2)),places);
  if(flag == 1)
   y = -y;
  return y;
}


/*
 * restore internal function from resource debugging
 */
config("resource_debug", resource_debug_level),;
if (config("resource_debug") & 3) {
    print "zeta(z)";
    print "psi(z)";
    print "polygamma(m,z)";
    print "lngamma(z)";
    print "gamma(z)";
    print "harmonic(limit)";
    print "gammainc(a,z)";
    print "heavisidestep(x)";
    print "expoint(z)";
    print "erf(z)";
    print "erfinv(x)";
    print "erfc(z)";
    print "erfi(z)";
    print "erfinv(x)";
    print "faddeeva(z)";
    print "gammap(a,z)";
    print "gammaq(a,z)";
    print "beta(a,b)";
    print "lnbeta(a,b)";
    print "betainc(z,a,b)";
}