/* specfunc/hermite.c

 * 

 * Copyright (C) 2011, 2012, 2013, 2014 Konrad Griessinger

 * (konradg(at)gmx.net)

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * This program 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

 * General Public License for more details.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

/*----------------------------------------------------------------------*

 * "The purpose of computing is insight, not numbers." - R.W. Hamming   *

 * Hermite polynomials, Hermite functions                               *

 * and their respective arbitrary derivatives                           *

 *----------------------------------------------------------------------*/

/* TODO:

 * - array functions for derivatives of Hermite functions

 * - asymptotic approximation for derivatives of Hermite functions

 * - refine existing asymptotic approximations, especially around x=sqrt(2*n+1) or x=sqrt(2*n+1)*sqrt(2), respectively

 */

#include <config.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_sf_airy.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_sf_pow_int.h>

#include <gsl/gsl_sf_gamma.h>

#include <gsl/gsl_sf_hermite.h>

#include "error.h"

#include "eval.h"

#define pow2(n) (gsl_sf_pow_int(2,n))

/* Evaluates the probabilists' Hermite polynomial of order n at position x using upward recurrence. */

static int 

gsl_sf_hermite_prob_iter_e(const int n, const double x, gsl_sf_result * result)

{

  result->val = 0.;

  result->err = 0.;

  if(n < 0) {

    DOMAIN_ERROR(result);

  }

  else if(n == 0) {

    result->val = 1.;

    result->err = 0.;

    return GSL_SUCCESS;

  }

  else if(n == 1) {

    result->val = x;

    result->err = 0.;

    return GSL_SUCCESS;

  }

  else if(x == 0.){

    if(GSL_IS_ODD(n)){

      result->val = 0.;

      result->err = 0.;

      return GSL_SUCCESS;

    }

    else{

      if(n < 301){

	/*

	double f;

	int j;

	f = (GSL_IS_ODD(n/2)?-1.:1.);

	for(j=1; j < n; j+=2) {

	  f*=j;

	}

	result->val = f;

	result->err = 0.;

	*/

	if(n < 297){

	  gsl_sf_doublefact_e(n-1, result);

	  (GSL_IS_ODD(n/2)?result->val = -result->val:1.);

	}

	else if (n == 298){

	  result->val = (GSL_IS_ODD(n/2)?-1.:1.)*1.25527562259930633890922678431e304;

	  result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);

	}

	else{

	  result->val = (GSL_IS_ODD(n/2)?-1.:1.)*3.7532741115719259533385880851e306;

	  result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);

	}

      }

      else{

	result->val = (GSL_IS_ODD(n/2)?GSL_NEGINF:GSL_POSINF);

	result->err = GSL_POSINF;

      }

      return GSL_SUCCESS;

    }

  }

/*

  else if(x*x < 4.0*n && n > 100000) {

    // asymptotic formula

    double f = 1.0;

    int j;

    if(GSL_IS_ODD(n)) {

      f=gsl_sf_fact((n-1)/2)*gsl_sf_pow_int(2,n/2)*M_SQRT2/M_SQRTPI;

    }

    else {

      for(j=1; j < n; j+=2) {

	f*=j;

      }

    }

    return f*exp(x*x/4)*cos(x*sqrt(n)-(n%4)*M_PI_2)/sqrt(sqrt(1-x*x/4.0/n));

    // return f*exp(x*x/4)*cos(x*sqrt(n)-n*M_PI_2)/sqrt(sqrt(1-x*x/4.0/n));

  }

*/

  else{

    /* upward recurrence: He_{n+1} = x He_n - n He_{n-1} */

    double p_n0 = 1.0;  /* He_0(x) */

    double p_n1 = x;    /* He_1(x) */

    double p_n = p_n1;

    double e_n0 = GSL_DBL_EPSILON;

    double e_n1 = fabs(x)*GSL_DBL_EPSILON;

    double e_n = e_n1;

    int j=0, c=0;

    for(j=1; j <= n-1; j++){

      if (gsl_isnan(p_n) == 1){

	break;

      }

      p_n  = x*p_n1-j*p_n0;

      p_n0 = p_n1;

      p_n1 = p_n;

      e_n  = (fabs(x)*e_n1+j*e_n0);

      e_n0 = e_n1;

      e_n1 = e_n;

      while(( GSL_MIN(fabs(p_n0),fabs(p_n1)) > 2.0*GSL_SQRT_DBL_MIN ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) > GSL_SQRT_DBL_MAX )){

	p_n0 *= 0.5;

	p_n1 *= 0.5;

	p_n = p_n1;

	e_n0 *= 0.5;

	e_n1 *= 0.5;

	e_n = e_n1;

	c++;

      }

      while(( ( ( fabs(p_n0) < GSL_SQRT_DBL_MIN ) && ( p_n0 != 0) ) || ( ( fabs(p_n1) < GSL_SQRT_DBL_MIN ) && ( p_n1 != 0) ) ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) < 0.5*GSL_SQRT_DBL_MAX )){

	p_n0 *= 2.0;

	p_n1 *= 2.0;

	p_n = p_n1;

	e_n0 *= 2.0;

	e_n1 *= 2.0;

	e_n = e_n1;

	c--;

      }

    }

    /*

    // check to see that the correct values are computed, even when overflow strikes in the end; works, thus very large results are accessible by determining mantissa and exponent separately

    double lg2 = 0.30102999566398119521467838;

    double ln10 = 2.3025850929940456840179914546843642076011014886;

    printf("res= %g\n", p_n*pow(10.,((lg2*c)-((long)(lg2*c)))) );

    printf("res= %g * 10^(%ld)\n", p_n*pow(10.,((lg2*c)-((long)(lg2*c))))/pow(10.,((long)(log(fabs(p_n*pow(10.,((lg2*c)-((long)(lg2*c))))))/ln10))), ((long)(log(fabs(p_n*pow(10.,((lg2*c)-((long)(lg2*c))))))/ln10))+((long)(lg2*c)) );

    */

    result->val = pow2(c)*p_n;

    result->err = pow2(c)*e_n + fabs(result->val)*GSL_DBL_EPSILON;

    /* result->err = e_n + n*fabs(p_n)*GSL_DBL_EPSILON;

       no idea, where the factor n came from => removed

     */

    if (gsl_isnan(result->val) != 1){

      return GSL_SUCCESS;

    }

    else{

      return GSL_ERANGE;

    }

  }

}

/* Approximatively evaluates the probabilists' Hermite polynomial of order n at position x.

 * An approximation depending on the x-range (see Szego, Gabor (1939, 1958, 1967), Orthogonal Polynomials, American Mathematical Society) is used. */

static int 

gsl_sf_hermite_prob_appr_e(const int n, const double x, gsl_sf_result * result)

{

  /* Plancherel-Rotach approximation (note: Szego defines the Airy function differently!) */

  const double aizero1 = -2.3381074104597670384891972524467; /* first zero of the Airy function Ai */

  double z = fabs(x)*M_SQRT1_2;

  double f = 1.;

  int j;

  for(j=1; j <= n; j++) {

    f*=sqrt(j);

  }

  if (z < sqrt(2*n+1.)+aizero1/M_SQRT2/pow(n,1/6.)){

    double phi = acos(z/sqrt(2*n+1.));

    result->val = f*(GSL_IS_ODD(n)&&(x<0.)?-1.:1.)*pow(2./n,0.25)/sqrt(M_SQRTPI*sin(phi))*sin(M_PI*0.75+(0.5*n+0.25)*(sin(2*phi)-2*phi))*exp(0.5*z*z);

    result->err = 2. * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_SUCCESS;

  }

  else if (z > sqrt(2*n+1.)-aizero1/M_SQRT2/pow(n,1/6.)){

    double phi = acosh(z/sqrt(2*n+1.));

    result->val = f*(GSL_IS_ODD(n)&&(x<0.)?-1.:1.)*pow(n,-0.25)/M_SQRT2/sqrt(M_SQRT2*M_SQRTPI*sinh(phi))*exp((0.5*n+0.25)*(2*phi-sinh(2*phi)))*exp(0.5*z*z);

    result->err = 2. * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_SUCCESS;

  }

  else{

    gsl_sf_result Ai;

    gsl_sf_airy_Ai_e((z-sqrt(2*n+1.))*M_SQRT2*pow(n,1/6.),0,&Ai;;

    result->val = f*(GSL_IS_ODD(n)&&(x<0.)?-1.:1.)*sqrt(M_SQRTPI*M_SQRT2)*pow(n,-1/12.)*Ai.val*exp(0.5*z*z);

    result->err = f*sqrt(M_SQRTPI*M_SQRT2)*pow(n,-1/12.)*Ai.err*exp(0.5*z*z) + GSL_DBL_EPSILON*fabs(result->val);

    return GSL_SUCCESS;

  }

}

/* Evaluates the probabilists' Hermite polynomial of order n at position x.

 * For small n upward recurrence is employed, while for large n and NaNs from the iteration an approximation depending on the x-range (see Szego, Gabor (1939, 1958, 1967), Orthogonal Polynomials, American Mathematical Society) is used. */

int 

gsl_sf_hermite_prob_e(const int n, const double x, gsl_sf_result * result)

{

  if( (x==0. || n<=100000) && (gsl_sf_hermite_prob_iter_e(n,x,result)==GSL_SUCCESS) ){

    return GSL_SUCCESS;

  }

  else{

    return gsl_sf_hermite_prob_appr_e(n,x,result);

  }

}

double gsl_sf_hermite_prob(const int n, const double x)

{

  EVAL_RESULT(gsl_sf_hermite_prob_e(n, x, &result));

}

/* Evaluates the m-th derivative of the probabilists' Hermite polynomial of order n at position x.

 * The direct formula He^{(m)}_n = n!/(n-m)!*He_{n-m}(x) (where He_j(x) is the j-th probabilists' Hermite polynomial and He^{(m)}_j(x) its m-th derivative) is employed. */

int

gsl_sf_hermite_prob_der_e(const int m, const int n, const double x, gsl_sf_result * result)

{

  if(n < 0 || m < 0) {

    DOMAIN_ERROR(result);

  }

  else if(n < m) {

    result->val = 0.;

    result->err = 0.;

    return GSL_SUCCESS;

  }

  else{

    double f = gsl_sf_choose(n,m)*gsl_sf_fact(m);

    gsl_sf_result He;

    gsl_sf_hermite_prob_e(n-m,x,&He;;

    result->val = He.val*f;

    result->err = He.err*f + GSL_DBL_EPSILON*fabs(result->val);

    return GSL_SUCCESS;

  }

}

double

gsl_sf_hermite_prob_der(const int m, const int n, const double x)

{

  EVAL_RESULT(gsl_sf_hermite_prob_der_e(m, n, x, &result));

}

/* Evaluates the physicists' Hermite polynomial of order n at position x.

 * For small n upward recurrence is employed, while for large n and NaNs from the iteration an approximation depending on the x-range (see Szego, Gabor (1939, 1958, 1967), Orthogonal Polynomials, American Mathematical Society) is used. */

int

gsl_sf_hermite_phys_e(const int n, const double x, gsl_sf_result * result)

{

  result->val = 0.;

  result->err = 0.;

  if(n < 0) {

    DOMAIN_ERROR(result);

  }

  else if(n == 0) {

    result->val = 1.;

    result->err = 0.;

    return GSL_SUCCESS;

  }

  else if(n == 1) {

    result->val = 2.0*x;

    result->err = 0.;

    return GSL_SUCCESS;

  }

  else if(x == 0.){

    if(GSL_IS_ODD(n)){

      result->val = 0.;

      result->err = 0.;

      return GSL_SUCCESS;

    }

    else{

      if(n < 269){

	double f = pow2(n/2);

	gsl_sf_doublefact_e(n-1, result);

	result->val *= f;

	result->err *= f;

	(GSL_IS_ODD(n/2)?result->val = -result->val:1.);

	/*

	double f;

	int j;

	f = (GSL_IS_ODD(n/2)?-1.:1.);

	for(j=1; j < n; j+=2) {

	  f*=2*j;

	}

	result->val = f;

	result->err = 0.;

	*/

      }

      else{

	result->val = (GSL_IS_ODD(n/2)?GSL_NEGINF:GSL_POSINF);

	result->err = GSL_POSINF;

      }

      return GSL_SUCCESS;

    }

  }

  /*

  else if(x*x < 2.0*n && n > 100000) {

    // asymptotic formula

    double f = 1.0;

    int j;

    if(GSL_IS_ODD(n)) {

      f=gsl_sf_fact((n-1)/2)*gsl_sf_pow_int(2,n)/M_SQRTPI;

    }

    else {

      for(j=1; j < n; j+=2) {

	f*=j;

      }

      f*=gsl_sf_pow_int(2,n/2);

    }

    return f*exp(x*x/2)*cos(x*sqrt(2.0*n)-(n%4)*M_PI_2)/sqrt(sqrt(1-x*x/2.0/n));

    // return f*exp(x*x/2)*cos(x*sqrt(2.0*n)-n*M_PI_2)/sqrt(sqrt(1-x*x/2.0/n));

  }

  */

  else if (n <= 100000){

    /* upward recurrence: H_{n+1} = 2x H_n - 2j H_{n-1} */

    double p_n0 = 1.0;    /* H_0(x) */

    double p_n1 = 2.0*x;  /* H_1(x) */

    double p_n = p_n1;

    double e_n0 = GSL_DBL_EPSILON;

    double e_n1 = 2.*fabs(x)*GSL_DBL_EPSILON;

    double e_n = e_n1;

    int j=0, c=0;

    for(j=1; j <= n-1; j++){

      if (gsl_isnan(p_n) == 1){

	break;

      }

      p_n  = 2.0*x*p_n1-2.0*j*p_n0;

      p_n0 = p_n1;

      p_n1 = p_n;

      e_n  = 2.*(fabs(x)*e_n1+j*e_n0);

      e_n0 = e_n1;

      e_n1 = e_n;

      while(( GSL_MIN(fabs(p_n0),fabs(p_n1)) > 2.0*GSL_SQRT_DBL_MIN ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) > GSL_SQRT_DBL_MAX )){

	p_n0 *= 0.5;

	p_n1 *= 0.5;

	p_n = p_n1;

	e_n0 *= 0.5;

	e_n1 *= 0.5;

	e_n = e_n1;

	c++;

      }

      while(( ( ( fabs(p_n0) < GSL_SQRT_DBL_MIN ) && ( p_n0 != 0) ) || ( ( fabs(p_n1) < GSL_SQRT_DBL_MIN ) && ( p_n1 != 0) ) ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) < 0.5*GSL_SQRT_DBL_MAX )){

	p_n0 *= 2.0;

	p_n1 *= 2.0;

	p_n = p_n1;

	e_n0 *= 2.0;

	e_n1 *= 2.0;

	e_n = e_n1;

	c--;

      }

    }

    result->val = pow2(c)*p_n;

    result->err = pow2(c)*e_n + fabs(result->val)*GSL_DBL_EPSILON;

    /* result->err = e_n + n*fabs(p_n)*GSL_DBL_EPSILON;

       no idea, where the factor n came from => removed */

    if (gsl_isnan(result->val) != 1){

      return GSL_SUCCESS;

    }

  }

  /* the following condition is implied by the logic above */

  {

    /* Plancherel-Rotach approximation (note: Szego defines the Airy function differently!) */

    const double aizero1 = -2.3381074104597670384891972524467; /* first zero of the Airy function Ai */

    double z = fabs(x);

    double f = 1.;

    int j;

    for(j=1; j <= n; j++) {

      f*=sqrt(j);

    }

    if (z < sqrt(2*n+1.)+aizero1/M_SQRT2/pow(n,1/6.)){

      double phi = acos(z/sqrt(2*n+1.));

      result->val = f*(GSL_IS_ODD(n)&&(x<0.)?-1.:1.)*(GSL_IS_ODD(n)?M_SQRT2:1.)*pow2(n/2)*pow(2./n,0.25)/sqrt(M_SQRTPI*sin(phi))*sin(M_PI*0.75+(0.5*n+0.25)*(sin(2*phi)-2*phi))*exp(0.5*z*z);

      result->err = 2. * GSL_DBL_EPSILON * fabs(result->val);

      return GSL_SUCCESS;

    }

    else if (z > sqrt(2*n+1.)-aizero1/M_SQRT2/pow(n,1/6.)){

      double phi = acosh(z/sqrt(2*n+1.));

      result->val = f*(GSL_IS_ODD(n)&&(x<0.)?-1.:1.)*(GSL_IS_ODD(n)?1.:M_SQRT1_2)*pow2(n/2)*pow(n,-0.25)/sqrt(M_SQRT2*M_SQRTPI*sinh(phi))*exp((0.5*n+0.25)*(2*phi-sinh(2*phi)))*exp(0.5*z*z);

      result->err = 2. * GSL_DBL_EPSILON * fabs(result->val);

      return GSL_SUCCESS;

    }

    else{

      gsl_sf_result Ai;

      gsl_sf_airy_Ai_e((z-sqrt(2*n+1.))*M_SQRT2*pow(n,1/6.),0,&Ai;;

      result->val = f*(GSL_IS_ODD(n)&&(x<0.)?-1.:1.)*(GSL_IS_ODD(n)?M_SQRT2:1.)*sqrt(M_SQRTPI*M_SQRT2)*pow2(n/2)*pow(n,-1/12.)*Ai.val*exp(0.5*z*z);

      result->err = f*(GSL_IS_ODD(n)?M_SQRT2:1.)*pow2(n/2)*sqrt(M_SQRTPI*M_SQRT2)*pow(n,-1/12.)*Ai.err*exp(0.5*z*z) + GSL_DBL_EPSILON*fabs(result->val);

      return GSL_SUCCESS;

    }

  }

}

double

gsl_sf_hermite_phys(const int n, const double x)

{

  EVAL_RESULT(gsl_sf_hermite_phys_e(n, x, &result));

}

/* Evaluates the m-th derivative of the physicists' Hermite polynomial of order n at position x.

 * The direct formula H^{(m)}_n = 2**m*n!/(n-m)!*H_{n-m}(x) (where H_j(x) is the j-th physicists' Hermite polynomial and H^{(m)}_j(x) its m-th derivative) is employed. */

int 

gsl_sf_hermite_phys_der_e(const int m, const int n, const double x, gsl_sf_result * result)

{

  if(n < 0 || m < 0) {

    DOMAIN_ERROR(result);

  }

  else if(n < m) {

    result->val = 0.;

    result->err = 0.;

    return GSL_SUCCESS;

  }

  else{

    double f = gsl_sf_choose(n,m)*gsl_sf_fact(m)*pow2(m);

    gsl_sf_result H;

    gsl_sf_hermite_phys_e(n-m,x,&H);

    result->val = H.val*f;

    result->err = H.err*f + GSL_DBL_EPSILON*fabs(result->val);

    return GSL_SUCCESS;

  }

}

double

gsl_sf_hermite_phys_der(const int m, const int n, const double x)

{

  EVAL_RESULT(gsl_sf_hermite_phys_der_e(m, n, x, &result));

}

/* Evaluates the Hermite function of order n at position x.

 * For large n an approximation depending on the x-range (see Szego, Gabor (1939, 1958, 1967), Orthogonal Polynomials, American Mathematical Society) is used, while for small n the direct formula via the probabilists' Hermite polynomial is applied. */

int

gsl_sf_hermite_func_e(const int n, const double x, gsl_sf_result * result)

{

  /*

  if (x*x < 2.0*n && n > 100000){

    // asymptotic formula

    double f = 1.0;

    int j;

    // return f*exp(x*x/4)*cos(x*sqrt(n)-n*M_PI_2)/sqrt(sqrt(1-x*x/4.0/n));

    return cos(x*sqrt(2.0*n)-(n%4)*M_PI_2)/sqrt(sqrt(0.5*n/M_PI*(1-0.5*x*x/n)))/M_PI;

  }

  */

  if (n < 0){

    DOMAIN_ERROR(result);

  }

  else if(n == 0 && x != 0.) {

    result->val = exp(-0.5*x*x)/sqrt(M_SQRTPI);

    result->err = GSL_DBL_EPSILON*fabs(result->val);

    return GSL_SUCCESS;

  }

  else if(n == 1 && x != 0.) {

    result->val = M_SQRT2*x*exp(-0.5*x*x)/sqrt(M_SQRTPI);

    result->err = GSL_DBL_EPSILON*fabs(result->val);

    return GSL_SUCCESS;

  }

  else if (x == 0.){

    if (GSL_IS_ODD(n)){

      result->val = 0.;

      result->err = 0.;

      return GSL_SUCCESS;

    }

    else{

      double f;

      int j;

      f = (GSL_IS_ODD(n/2)?-1.:1.);

      for(j=1; j < n; j+=2) {

	f*=sqrt(j/(j+1.));

      }

      result->val = f/sqrt(M_SQRTPI);

      result->err = GSL_DBL_EPSILON*fabs(result->val);

      return GSL_SUCCESS;

    }

  }

  else if (n <= 100000){

    double f = exp(-0.5*x*x)/sqrt(M_SQRTPI*gsl_sf_fact(n));

    gsl_sf_result He;

    gsl_sf_hermite_prob_iter_e(n,M_SQRT2*x,&He;;

    result->val = He.val*f;

    result->err = He.err*f + GSL_DBL_EPSILON*fabs(result->val);

    if (gsl_isnan(result->val) != 1 && f > GSL_DBL_MIN && gsl_finite(He.val) == 1){

      return GSL_SUCCESS;

    }

  }

  /* upward recurrence: Psi_{n+1} = sqrt(2/(n+1))*x Psi_n - sqrt(n/(n+1)) Psi_{n-1} */

  {

    double tw = exp(-x*x*0.5/n); /* "twiddle factor" (in the spirit of FFT) */

    double p_n0 = tw/sqrt(M_SQRTPI); /* Psi_0(x) */

    double p_n1 = p_n0*M_SQRT2*x;    /* Psi_1(x) */

    double p_n = p_n1;

    double e_n0 = p_n0*GSL_DBL_EPSILON;

    double e_n1 = p_n1*GSL_DBL_EPSILON;

    double e_n = e_n1;

    int j;

    int c = 0;

    for (j=1;j

      {

        if (gsl_isnan(p_n) == 1){

	  break;

        }

        p_n=tw*(M_SQRT2*x*p_n1-sqrt(j)*p_n0)/sqrt(j+1.);

        p_n0=tw*p_n1;

        p_n1=p_n;

        e_n=tw*(M_SQRT2*fabs(x)*e_n1+sqrt(j)*e_n0)/sqrt(j+1.);

        e_n0=tw*e_n1;

        e_n1=e_n;

        while(( GSL_MIN(fabs(p_n0),fabs(p_n1)) > 2.0*GSL_SQRT_DBL_MIN ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) > GSL_SQRT_DBL_MAX )){

  	p_n0 *= 0.5;

  	p_n1 *= 0.5;

	p_n = p_n1;

	e_n0 *= 0.5;

	e_n1 *= 0.5;

	e_n = e_n1;

	c++;

      }

	while(( ( ( fabs(p_n0) < GSL_SQRT_DBL_MIN ) && ( p_n0 != 0) ) || ( ( fabs(p_n1) < GSL_SQRT_DBL_MIN ) && ( p_n1 != 0) ) ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) < 0.5*GSL_SQRT_DBL_MAX )){

	p_n0 = p_n0*2;

	p_n1 = p_n1*2;

	p_n = p_n1;

	e_n0 = e_n0*2;

	e_n1 = e_n1*2;

	e_n = e_n1;

	c--;

      }

    }

  result->val = p_n*pow2(c);

  result->err = n*fabs(result->val)*GSL_DBL_EPSILON;

  if (gsl_isnan(result->val) != 1){

    return GSL_SUCCESS;

  }

  {

    /* Plancherel-Rotach approximation (note: Szego defines the Airy function differently!) */

    const double aizero1 = -2.3381074104597670384891972524467; /* first zero of the Airy function Ai */

    double z = fabs(x);

    if (z < sqrt(2*n+1.)+aizero1/M_SQRT2/pow(n,1/6.)){

      double phi = acos(z/sqrt(2*n+1.));

      result->val = (GSL_IS_ODD(n)&&(x<0.)?-1.:1.)*pow(2./n,0.25)/M_SQRTPI/sqrt(sin(phi))*sin(M_PI*0.75+(0.5*n+0.25)*(sin(2*phi)-2*phi));

      result->err = 2. * GSL_DBL_EPSILON * fabs(result->val);

      return GSL_SUCCESS;

    }

    else if (z > sqrt(2*n+1.)-aizero1/M_SQRT2/pow(n,1/6.)){

      double phi = acosh(z/sqrt(2*n+1.));

      result->val = (GSL_IS_ODD(n)&&(x<0.)?-1.:1.)*pow(n,-0.25)/

2/M_SQRTPI/sqrt(sinh(phi)/M_SQRT2)*exp((0.5*n+0.25)*(2*phi-sinh(2*phi)));

      result->err = 2. * GSL_DBL_EPSILON * fabs(result->val);

      return GSL_SUCCESS;

    }

    else{

      gsl_sf_result Ai;

      gsl_sf_airy_Ai_e((z-sqrt(2*n+1.))*M_SQRT2*pow(n,1/6.),0,&Ai;;

      result->val = (GSL_IS_ODD(n)&&(x<0.)?-1.:1.)*sqrt(M_SQRT2)*pow(n,-1/12.)*Ai.val;

      result->err = sqrt(M_SQRT2)*pow(n,-1/12.)*Ai.err +  GSL_DBL_EPSILON*fabs(result->val);

      return GSL_SUCCESS;

    }

  }

  }

}

double

gsl_sf_hermite_func(const int n, const double x)

{

  EVAL_RESULT(gsl_sf_hermite_func_e(n, x, &result));

}

/* Evaluates all probabilists' Hermite polynomials up to order nmax at position x. The results are stored in result_array.

 * Since all polynomial orders are needed, upward recurrence is employed. */

int

gsl_sf_hermite_prob_array(const int nmax, const double x, double * result_array)

{

  if(nmax < 0) {

    GSL_ERROR ("domain error", GSL_EDOM);

  }

  else if(nmax == 0) {

    result_array[0] = 1.0;

    return GSL_SUCCESS;

  }

  else if(nmax == 1) {

    result_array[0] = 1.0;

    result_array[1] = x;

    return GSL_SUCCESS;

  }

  else {

    /* upward recurrence: He_{n+1} = x He_n - n He_{n-1} */

    double p_n0 = 1.0;    /* He_0(x) */

    double p_n1 = x;      /* He_1(x) */

    double p_n = p_n1;

    int j=0, c=0;

    result_array[0] = 1.0;

    result_array[1] = x;

    for(j=1; j <= nmax-1; j++){

      p_n  = x*p_n1-j*p_n0;

      p_n0 = p_n1;

      p_n1 = p_n;

      while(( GSL_MIN(fabs(p_n0),fabs(p_n1)) > 2.0*GSL_SQRT_DBL_MIN ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) > GSL_SQRT_DBL_MAX )){

	p_n0 *= 0.5;

	p_n1 *= 0.5;

	p_n = p_n1;

	c++;

      }

      while(( ( ( fabs(p_n0) < GSL_SQRT_DBL_MIN ) && ( p_n0 != 0) ) || ( ( fabs(p_n1) < GSL_SQRT_DBL_MIN ) && ( p_n1 != 0) ) ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) < 0.5*GSL_SQRT_DBL_MAX )){

	p_n0 *= 2.0;

	p_n1 *= 2.0;

	p_n = p_n1;

	c--;

      }

      result_array[j+1] = pow2(c)*p_n;

    }

    return GSL_SUCCESS;

  }

}

/* Evaluates the m-th derivative of all probabilists' Hermite polynomials up to order nmax at position x. The results are stored in result_array.

 * Since all polynomial orders are needed, upward recurrence is employed. */

int

gsl_sf_hermite_prob_array_der(const int m, const int nmax, const double x, double * result_array)

{

  if(nmax < 0 || m < 0) {

    GSL_ERROR ("domain error", GSL_EDOM);

  }

  else if(m == 0) {

    gsl_sf_hermite_prob_array(nmax, x, result_array);

    return GSL_SUCCESS;

  }

  else if(nmax < m) {

    int j;

    for(j=0; j <= nmax; j++){

      result_array[j] = 0.0;

    }

    return GSL_SUCCESS;

  }

  else if(nmax == m) {

    int j;

    for(j=0; j < m; j++){

      result_array[j] = 0.0;

    }

    result_array[nmax] = gsl_sf_fact(m);

    return GSL_SUCCESS;

  }

  else if(nmax == m+1) {

    int j;

    for(j=0; j < m; j++){

      result_array[j] = 0.0;

    }

    result_array[nmax-1] = gsl_sf_fact(m);

    result_array[nmax] = result_array[nmax-1]*(m+1)*x;

    return GSL_SUCCESS;

  }

  else {

    /* upward recurrence: He^{(m)}_{n+1} = (n+1)/(n-m+1)*(x He^{(m)}_n - n He^{(m)}_{n-1}) */

    double p_n0 = gsl_sf_fact(m);    /* He^{(m)}_{m}(x) */

    double p_n1 = p_n0*(m+1)*x;      /* He^{(m)}_{m+1}(x) */

    double p_n = p_n1;

    int j=0, c=0;

    for(j=0; j < m; j++){

      result_array[j] = 0.0;

    }

    result_array[m] = p_n0;

    result_array[m+1] = p_n1;

    for(j=m+1; j <= nmax-1; j++){

      p_n  = (x*p_n1-j*p_n0)*(j+1)/(j-m+1);

      p_n0 = p_n1;

      p_n1 = p_n;

      while(( GSL_MIN(fabs(p_n0),fabs(p_n1)) > 2.0*GSL_SQRT_DBL_MIN ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) > GSL_SQRT_DBL_MAX )){

	p_n0 *= 0.5;

	p_n1 *= 0.5;

	p_n = p_n1;

	c++;

      }

      while(( ( ( fabs(p_n0) < GSL_SQRT_DBL_MIN ) && ( p_n0 != 0) ) || ( ( fabs(p_n1) < GSL_SQRT_DBL_MIN ) && ( p_n1 != 0) ) ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) < 0.5*GSL_SQRT_DBL_MAX )){

	p_n0 *= 2.0;

	p_n1 *= 2.0;

	p_n = p_n1;

	c--;

      }

      result_array[j+1] = pow2(c)*p_n;

    }

    return GSL_SUCCESS;

  }

}

/* Evaluates all derivatives (starting from 0) up to the mmax-th derivative of the probabilists' Hermite polynomial of order n at position x. The results are stored in result_array.

 * Since all polynomial orders are needed, upward recurrence is employed. */

int

gsl_sf_hermite_prob_der_array(const int mmax, const int n, const double x, double * result_array)

{

  if(n < 0 || mmax < 0) {

    GSL_ERROR ("domain error", GSL_EDOM);

  }

  else if(n == 0) {

    int j;

    result_array[0] = 1.0;

    for(j=1; j <= mmax; j++){

      result_array[j] = 0.0;

    }

    return GSL_SUCCESS;

  }

  else if(n == 1 && mmax > 0) {

    int j;

    result_array[0] = x;

    result_array[1] = 1.0;

    for(j=2; j <= mmax; j++){

      result_array[j] = 0.0;

    }

    return GSL_SUCCESS;

  }

  else if( mmax == 0) {

    result_array[0] = gsl_sf_hermite_prob(n,x);

    return GSL_SUCCESS;

  }

  else if( mmax == 1) {

    result_array[0] = gsl_sf_hermite_prob(n,x);

    result_array[1] = n*gsl_sf_hermite_prob(n-1,x);

    return GSL_SUCCESS;

  }

  else {

    /* upward recurrence */

    int k = GSL_MAX_INT(0,n-mmax);

    /* Getting a bit lazy here... */

    double p_n0 = gsl_sf_hermite_prob(k,x);        /* He_k(x) */

    double p_n1 = gsl_sf_hermite_prob(k+1,x);      /* He_{k+1}(x) */

    double p_n  = p_n1;

    int j=0, c=0;

    for(j=n+1; j <= mmax; j++){

      result_array[j] = 0.0;

    }

    result_array[GSL_MIN_INT(n,mmax)] = p_n0;

    result_array[GSL_MIN_INT(n,mmax)-1] = p_n1;

    for(j=GSL_MIN_INT(mmax,n)-1; j > 0; j--){

      k++;

      p_n  = x*p_n1-k*p_n0;

      p_n0 = p_n1;

      p_n1 = p_n;

      while(( GSL_MIN(fabs(p_n0),fabs(p_n1)) > 2.0*GSL_SQRT_DBL_MIN ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) > GSL_SQRT_DBL_MAX )){

	p_n0 *= 0.5;

	p_n1 *= 0.5;

	p_n = p_n1;

	c++;

      }

      while(( ( ( fabs(p_n0) < GSL_SQRT_DBL_MIN ) && ( p_n0 != 0) ) || ( ( fabs(p_n1) < GSL_SQRT_DBL_MIN ) && ( p_n1 != 0) ) ) && ( GSL_MAX(fabs(p_n0),fabs(p_n1)) < 0.5*GSL_SQRT_DBL_MAX )){

	p_n0 *= 2.0;

	p_n1 *= 2.0;

	p_n = p_n1;

	c--;

      }

      result_array[j-1] = pow2(c)*p_n;

    }

    p_n = 1.0;

    for(j=1; j <= GSL_MIN_INT(n,mmax); j++){

      p_n  = p_n*(n-j+1);

      result_array[j] = p_n*result_array[j];

    }

    return GSL_SUCCESS;

  }

}

/* Evaluates the series sum_{j=0}^n a_j*He_j(x) with He_j being the j-th probabilists' Hermite polynomial.

 * For improved numerical stability the Clenshaw algorithm (Clenshaw, C. W. (July 1955). "A note on the summation of Chebyshev series". Mathematical Tables and other Aids to Computation 9 (51): 118–110.) adapted to probabilists' Hermite polynomials is used. */

int

gsl_sf_hermite_prob_series_e(const int n, const double x, const double * a, gsl_sf_result * result)

{

  if(n < 0) {

    DOMAIN_ERROR(result);

  }

  else if(n == 0) {

    result->val = a[0];

    result->err = 0.;

    return GSL_SUCCESS;

  }

  else if(n == 1) {

    result->val = a[0]+a[1]*x;

    result->err = 2.*GSL_DBL_EPSILON * (fabs(a[0]) + fabs(a[1]*x)) ;

    return GSL_SUCCESS;

  }

  else {

    /* downward recurrence: b_n = a_n + x b_{n+1} - (n+1) b_{n+2} */

    double b0 = 0.;

    double b1 = 0.;

    double btmp = 0.;

    double e0 = 0.;

    double e1 = 0.;

    double etmp = e1;

    int j;

    for(j=n; j >= 0; j--){

      btmp = b0;

      b0  = a[j]+x*b0-(j+1)*b1;

      b1 = btmp;

      etmp = e0;

      e0  = (GSL_DBL_EPSILON*fabs(a[j])+fabs(x)*e0+(j+1)*e1);

      e1 = etmp;

    }

    result->val = b0;

    result->err = e0 + fabs(b0)*GSL_DBL_EPSILON;

    return GSL_SUCCESS;

  }

}

double

gsl_sf_hermite_prob_series(const int n, const double x, const double * a)

{

  EVAL_RESULT(gsl_sf_hermite_prob_series_e(n, x, a, &result));

}

/* Evaluates all physicists' Hermite polynomials up to order nmax at position x. The results are stored in result_array.