Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/Makefile (revision 37532) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/Makefile (revision 37533) @@ -35,4 +35,6 @@ $(SRC)/time.$(ARCH).o \ $(SRC)/gaussj.$(ARCH).o \ +$(SRC)/legendre.$(ARCH).o \ +$(SRC)/gsl_utils.$(ARCH).o \ $(SRC)/vstats.$(ARCH).o \ $(SRC)/config.$(ARCH).o \ Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/gsl_utils.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/gsl_utils.h (revision 37533) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/gsl_utils.h (revision 37533) @@ -0,0 +1,9 @@ +# define GSL_DBL_EPSILON 2.2204460492503131e-16 +# define GSL_LOG_DBL_MIN (-7.0839641853226408e+02) + +# define GSL_ROOT6_DBL_EPSILON 2.4607833005759251e-03 + +# define GSL_SF_GAMMA_XMAX 171.0 + +# define GSL_IS_ODD(n) ((n) & 1) + Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h (revision 37532) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/include/ohana.h (revision 37533) @@ -435,3 +435,15 @@ */ -# endif + +/** gsl-based functions for legendre and related polynomials **/ + +double ohana_gsl_sf_log_1plusx (double x, double *err); +double ohana_gsl_sf_lnpoch(const double a, const double x, double *err); +double ohana_gsl_lnpoch_pos (const double a, const double x, double *err); +double ohana_gsl_pochrel_smallx(const double a, const double x, double *err); + +double legendre_Pl (int l, double x, double *err); +double legendre_Plm (int l, int m, double x, double *err); +double legendre_Plm_sphere(int l, int m, double x, double *err); + +# endif Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/gsl_utils.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/gsl_utils.c (revision 37533) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/gsl_utils.c (revision 37533) @@ -0,0 +1,395 @@ +# include +# include + +typedef struct { + double *c; /* coefficients */ + int order; /* order of expansion */ + double a; /* lower interval point */ + double b; /* upper interval point */ + int order_sp; /* effective single precision order */ +} cheb_series; + +/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/ + +/* Chebyshev expansion for log(1 + x(t))/x(t) + * + * x(t) = (4t-1)/(2(4-t)) + * t(x) = (8x+1)/(2(x+2)) + * -1/2 < x < 1/2 + * -1 < t < 1 + */ +static double log_1plusx_data[21] = { + 2.16647910664395270521272590407, + -0.28565398551049742084877469679, + 0.01517767255690553732382488171, + -0.00200215904941415466274422081, + 0.00019211375164056698287947962, + -0.00002553258886105542567601400, + 2.9004512660400621301999384544e-06, + -3.8873813517057343800270917900e-07, + 4.7743678729400456026672697926e-08, + -6.4501969776090319441714445454e-09, + 8.2751976628812389601561347296e-10, + -1.1260499376492049411710290413e-10, + 1.4844576692270934446023686322e-11, + -2.0328515972462118942821556033e-12, + 2.7291231220549214896095654769e-13, + -3.7581977830387938294437434651e-14, + 5.1107345870861673561462339876e-15, + -7.0722150011433276578323272272e-16, + 9.7089758328248469219003866867e-17, + -1.3492637457521938883731579510e-17, + 1.8657327910677296608121390705e-18 +}; + +static cheb_series log_1plusx_cs = { + log_1plusx_data, + 20, + -1, 1, + 10 +}; + +static inline double ohana_gsl_cheb_eval(const cheb_series *cs, const double x, double *err) { + + double d = 0.0; + double dd = 0.0; + + double y = (2.0*x - cs->a - cs->b) / (cs->b - cs->a); + double y2 = 2.0 * y; + + double e = 0.0; + + int j; + for (j = cs->order; j >= 1; j--) { + double temp = d; + d = y2*d - dd + cs->c[j]; + e += fabs(y2*temp) + fabs(dd) + fabs(cs->c[j]); + dd = temp; + } + + { + double temp = d; + d = y*d - dd + 0.5 * cs->c[0]; + e += fabs(y*temp) + fabs(dd) + 0.5 * fabs(cs->c[0]); + } + + double value = d; + if (err) *err = GSL_DBL_EPSILON * e + fabs(cs->c[cs->order]); + + return value; +} + +/*** end chebyshev approximation code *******************************************/ + +double ohana_gsl_sf_log_1plusx (double x, double *err) { + + if (x <= -1.0) return NAN; + + if (fabs(x) < GSL_ROOT6_DBL_EPSILON) { + const double c1 = -0.5; + const double c2 = 1.0/3.0; + const double c3 = -1.0/4.0; + const double c4 = 1.0/5.0; + const double c5 = -1.0/6.0; + const double c6 = 1.0/7.0; + const double c7 = -1.0/8.0; + const double c8 = 1.0/9.0; + const double c9 = -1.0/10.0; + const double t = c5 + x*(c6 + x*(c7 + x*(c8 + x*c9))); + double value = x * (1.0 + x*(c1 + x*(c2 + x*(c3 + x*(c4 + x*t))))); + if (err) *err = GSL_DBL_EPSILON * fabs(value); + return value; + } + + if (fabs(x) < 0.5) { + double t = 0.5*(8.0*x + 1.0)/(x+2.0); + + // log_1plusx_cs is a static cheb_series structure + double c = ohana_gsl_cheb_eval (&log_1plusx_cs, t, err); + + double value = x * c; + + if (err) *err = fabs(x * *err); + + return value; + } + + double value = log(1.0 + x); + + if (err) *err = GSL_DBL_EPSILON * fabs(value); + + return value; +} + +double ohana_gsl_sf_lnpoch(const double a, const double x, double *err) { + + if (a <= 0.0) return NAN; + if (a+x <= 0.0) return NAN; + + if (x == 0.0) { + double value = 0.0; + if (err) *err = 0.0; + return value; + } + + double value = ohana_gsl_lnpoch_pos (a, x, err); + return value; +} + +/* Assumes a>0 and a+x>0. + */ +double ohana_gsl_lnpoch_pos (const double a, const double x, double *err) { + + double absx = fabs(x); + + if (absx > 0.1*a || absx*log(MAX(a,2.0)) > 0.1) { + if(a < GSL_SF_GAMMA_XMAX && a+x < GSL_SF_GAMMA_XMAX) { + /* If we can do it by calculating the gamma functions + * directly, then that will be more accurate than + * doing the subtraction of the logs. + */ + + double g1_error = NAN; + double g2_error = NAN; + double *ptrerr = NULL; + + if (err) ptrerr = &g1_error; + double g1_value = ohana_gsl_sf_gammainv(a, ptrerr); + + if (err) ptrerr = &g2_error; + double g2_value = ohana_gsl_sf_gammainv(a+x, ptrerr); + + value = -log(g2_value / g1_value); + + if (err) { + double tmp_err = (g1_error / fabs(g1_value)) + (g2_error / fabs(g2_value)); + tmp_err += 2.0 * GSL_DBL_EPSILON * fabs(value); + *err = tmp_err; + } + return value; + } + + /* Otherwise we must do the subtraction. + */ + double lg1_error = NAN; + double lg2_error = NAN; + double *ptrerr = NULL; + + if (err) ptrerr = &lg1_error; + double lg1_value = ohana_gsl_sf_lngamma(a, ptrerr); + + if (err) ptrerr = &lg2_error; + double lg2_value = ohana_gsl_sf_lngamma(a+x, ptrerr); + + double value = lg2_value - lg1_value; + + if (err) { + double tmp_err = lg2_error + lg1_error; + tmp_err += 2.0 * GSL_DBL_EPSILON * fabs(value); + *err = tmp_err; + } + return value; + } + + if ((absx < 0.1*a) && (a > 15.0)) { + /* Be careful about the implied subtraction. + * Note that both a+x and and a must be + * large here since a is not small + * and x is not relatively large. + * So we calculate using Stirling for Log[Gamma(z)]. + * + * Log[Gamma(a+x)/Gamma(a)] = x(Log[a]-1) + (x+a-1/2)Log[1+x/a] + * + (1/(1+eps) - 1) / (12 a) + * - (1/(1+eps)^3 - 1) / (360 a^3) + * + (1/(1+eps)^5 - 1) / (1260 a^5) + * - (1/(1+eps)^7 - 1) / (1680 a^7) + * + ... + */ + const double eps = x/a; + const double den = 1.0 + eps; + const double d3 = den*den*den; + const double d5 = d3*den*den; + const double d7 = d5*den*den; + const double c1 = -eps/den; + const double c3 = -eps*(3.0+eps*(3.0+eps))/d3; + const double c5 = -eps*(5.0+eps*(10.0+eps*(10.0+eps*(5.0+eps))))/d5; + const double c7 = -eps*(7.0+eps*(21.0+eps*(35.0+eps*(35.0+eps*(21.0+eps*(7.0+eps))))))/d7; + + const double p8 = ohana_gsl_sf_pow_int(1.0+eps,8); + const double c8 = 1.0/p8 - 1.0; /* these need not */ + const double c9 = 1.0/(p8*(1.0+eps)) - 1.0; /* be very accurate */ + const double a4 = a*a*a*a; + const double a6 = a4*a*a; + const double ser_1 = c1 + c3/(30.0*a*a) + c5/(105.0*a4) + c7/(140.0*a6); + const double ser_2 = c8/(99.0*a6*a*a) - 691.0/360360.0 * c9/(a6*a4); + const double ser = (ser_1 + ser_2)/ (12.0*a); + + double term1 = x * log(a/M_E); + double term2; + + double ln_1peps_error = NAN; + double *ptr = err ? &ln_1peps_error : NULL; + + double ln_1peps_value = ohana_gsl_sf_log_1plusx(eps, ptr); /* log(1 + x/a) */ + term2 = (x + a - 0.5) * ln_1peps_value; + + double value = term1 + term2 + ser; + + if (err) { + double tmp_err = GSL_DBL_EPSILON*fabs(term1); + tmp_err += fabs((x + a - 0.5)*ln_1peps_error); + tmp_err += fabs(ln_1peps_value) * GSL_DBL_EPSILON * (fabs(x) + fabs(a) + 0.5); + tmp_err += 2.0 * GSL_DBL_EPSILON * fabs(value); + *err = tmp_err; + } + return value; + } + + double poch_rel_error = NAN; + double *ptr = err ? &poch_rel_error : NULL; + double poch_rel_value = ohana_gsl_pochrel_smallx (a, x, ptr); + + double eps = x*poch_rel_value; + + double value = ohana_gsl_sf_log_1plusx(eps, NULL); + + if (err) { + double tmp_err = 2.0 * fabs(x * poch_rel_error / (1.0 + eps)); + tmp_err += 2.0 * GSL_DBL_EPSILON * fabs(value); + *err = tmp_err; + } + + return value; +} + +/* ((a)_x - 1)/x in the "small x" region where + * cancellation must be controlled. + * + * Based on SLATEC DPOCH1(). + */ +/* +C When ABS(X) is so small that substantial cancellation will occur if +C the straightforward formula is used, we use an expansion due +C to Fields and discussed by Y. L. Luke, The Special Functions and Their +C Approximations, Vol. 1, Academic Press, 1969, page 34. +C +C The ratio POCH(A,X) = GAMMA(A+X)/GAMMA(A) is written by Luke as +C (A+(X-1)/2)**X * polynomial in (A+(X-1)/2)**(-2) . +C In order to maintain significance in POCH1, we write for positive a +C (A+(X-1)/2)**X = EXP(X*LOG(A+(X-1)/2)) = EXP(Q) +C = 1.0 + Q*EXPREL(Q) . +C Likewise the polynomial is written +C POLY = 1.0 + X*POLY1(A,X) . +C Thus, +C POCH1(A,X) = (POCH(A,X) - 1) / X +C = EXPREL(Q)*(Q/X + Q*POLY1(A,X)) + POLY1(A,X) +C +*/ +double ohana_gsl_pochrel_smallx(const double a, const double x, double *err) { + /* + SQTBIG = 1.0D0/SQRT(24.0D0*D1MACH(1)) + ALNEPS = LOG(D1MACH(3)) + */ + const double SQTBIG = 1.0/(2.0*M_SQRT2*M_SQRT3*GSL_SQRT_DBL_MIN); + const double ALNEPS = GSL_LOG_DBL_EPSILON - M_LN2; + + if (x == 0.0) { + double value = ohana_gsl_sf_psi(a, err); + return value; + } + + const double bp = ( (a < -0.5) ? 1.0-a-x : a ); + const int incr = ( (bp < 10.0) ? 11.0-bp : 0 ); + const double b = bp + incr; + + double var = b + 0.5*(x-1.0); + double alnvar = log(var); + double q = x*alnvar; + + double poly1 = 0.0; + + if (var < SQTBIG) { + const int nterms = (int)(-0.5*ALNEPS/alnvar + 1.0); + const double var2 = (1.0/var)/var; + const double rho = 0.5 * (x + 1.0); + double term = var2; + double gbern[24]; + int k, j; + + gbern[1] = 1.0; + gbern[2] = -rho/12.0; + poly1 = gbern[2] * term; + + if (nterms > 20) { + /* NTERMS IS TOO BIG, MAYBE D1MACH(3) IS BAD */ + /* nterms = 20; */ + value = NAN; + if (err) *err = NAN; + return value; + } + + for (k = 2; k <= nterms; k++) { + double gbk = 0.0; + for (j = 1; j <= k; j++) { + gbk += bern[k-j+1]*gbern[j]; + } + gbern[k+1] = -rho*gbk/k; + + term *= (2*k-2-x)*(2*k-1-x)*var2; + poly1 += gbern[k+1]*term; + } + } + + double dexprl = ohana_gsl_sf_expm1(q); + if (!isfinite(dexprl)) { + double value = NAN; + if (err) *err = NAN; + return value; + } + + dexprl = dexprl / q; + poly1 *= (x - 1.0); + double dpoch1 = dexprl * (alnvar + q * poly1) + poly1; + + int i; + for (i = incr-1; i >= 0; i--) { + /* + C WE HAVE DPOCH1(B,X), BUT BP IS SMALL, SO WE USE BACKWARDS RECURSION + C TO OBTAIN DPOCH1(BP,X). + */ + double binv = 1.0/(bp+i); + dpoch1 = (dpoch1 - binv) / (1.0 + x*binv); + } + + if (bp == a) { + double value = dpoch1; + if (err) *err = 2.0 * GSL_DBL_EPSILON * (fabs(incr) + 1.0) * fabs(value); + return value; + } + + /* + C WE HAVE DPOCH1(BP,X), BUT A IS LT -0.5. WE THEREFORE USE A + C REFLECTION FORMULA TO OBTAIN DPOCH1(A,X). + */ + double sinpxx = sin(M_PI*x)/x; + double sinpx2 = sin(0.5*M_PI*x); + double t1 = sinpxx/tan(M_PI*b); + double t2 = 2.0*sinpx2*(sinpx2/x); + double trig = t1 - t2; + double value = dpoch1 * (1.0 + x*trig) + trig; + + if (err) { + double tmp_err = (fabs(dpoch1*x) + 1.0) * GSL_DBL_EPSILON * (fabs(t1) + fabs(t2)); + tmp_err += 2.0 * GSL_DBL_EPSILON * (fabs(incr) + 1.0) * fabs(value); + *err = tmp_err; + } + return value; +} + +/* functions I still need to define: */ + +ohana_gsl_sf_gammainv +ohana_gsl_sf_lngamma +ohana_gsl_sf_pow_int +ohana_gsl_sf_psi +ohana_gsl_sf_expm1 Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/legendre.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/legendre.c (revision 37533) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/legendre.c (revision 37533) @@ -0,0 +1,366 @@ +# include +# include + +/* some useful functions for legendre polynomials, spherical harmonics, and vector spherical harmonics */ + +/* Some of the functions below are derived based on code from the gnu scientific library (gsl) v1.11 + + The corresponding copyright from that code is given in doc/gsl.txt + + * From: + * specfunc/legendre_poly.c + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002 Gerard Jungman + * Author: G. Jungman + + * adaptations by Eugene Magnier + */ + +/* Calculate P_m^m(x) from the analytic result: + * P_m^m(x) = (-1)^m (2m-1)!! (1-x^2)^(m/2) , m > 0 + * = 1 , m = 0 + * from gsl v1.11 specfunc/legendre_poly.c: + */ +static double legendre_Pmm(int m, double x) { + + if (m == 0) return 1.0; + + double value = 1.0; + double factor = sqrt(1.0 - x) * sqrt(1.0 + x); + double coeff = 1.0; + + int i; + for (i = 1; i <= m; i++) { + value *= -coeff * factor; + coeff += 2.0; + } + return value; +} + +/** these functions also return the accumulated floating-point error if requested */ + +static double legendre_P0(double x, double *err) { + + double value = 1.0; + if (err) *err = 0.0; + + return value; +} + +static double legendre_P1(double x, double *err) { + + double value = x; + if (err) *err = 0.0; + + return value; +} + +static double legendre_P2(double x, double *err) { + + double value = 0.5*(3.0*x*x - 1.0); + if (err) *err = GSL_DBL_EPSILON * (fabs(3.0*x*x) + 1.0); + + return value; +} + +static double legendre_P3(double x, double *err) { + + double value = 0.5*x*(5.0*x*x - 3.0); + if (err) *err = GSL_DBL_EPSILON * (fabs(value) + 0.5 * fabs(x) * (fabs(5.0*x*x) + 3.0)); + + return value; +} + +double legendre_Pl (int l, double x, double *err) { + + double value; + + if (l < 0) return NAN; + if (x < -1.0) return NAN; + if (x > +1.0) return NAN; + + // pre-defined l-values: + if (l == 0) { + value = legendre_P0(x, err); + return value; + } + if (l == 1) { + value = legendre_P1(x, err); + return value; + } + if (l == 2) { + value = legendre_P2(x, err); + return value; + } + if (l == 3) { + value = legendre_P3(x, err); + return value; + } + + if (x == 1.0) { + value = 1.0; + if (err) *err = 0.0; + return value; + } + if (x == -1.0) { + result->val = (GSL_IS_ODD(l) ? -1.0 : 1.0); + if (err) *err = 0.0; + return value; + } + + if (l < 100000) { + /* upward recurrence: l P_l = (2l-1) z P_{l-1} - (l-1) P_{l-2} */ + + double p_LM2 = 1.0; /* P_0(x) */ + double p_LM1 = x; /* P_1(x) */ + double p_L = p_LM1; + + double e_LM2 = GSL_DBL_EPSILON; + double e_LM1 = fabs(x)*GSL_DBL_EPSILON; + double e_L = e_LM1; + + int L; + for (L = 2; L <= l; L++){ + p_L = (x*(2*L-1)*p_LM1 - (L-1)*p_LM2) / L; + p_LM2 = p_LM1; + p_LM1 = p_L; + + if (err) { + e_L = 0.5*(fabs(x)*(2*L-1.0) * e_LM1 + (L-1.0)*e_LM2)/L; + e_LM2 = e_LM1; + e_LM1 = e_L; + } + } + + value = p_L; + if (err) *err = e_L + l*fabs(p_L)*GSL_DBL_EPSILON; + return value; + } + + /* last resort : Asymptotic expansion. + * [Olver, p. 473] + */ +# if (0) + double u = l + 0.5; + double th = acos(x); + gsl_sf_result J0; + gsl_sf_result Jm1; + int stat_J0 = gsl_sf_bessel_J0_e(u*th, &J0); + int stat_Jm1 = gsl_sf_bessel_Jn_e(-1, u*th, &Jm1); + double pre; + double B00; + double c1; + + /* B00 = 1/8 (1 - th cot(th) / th^2 + * pre = sqrt(th/sin(th)) + */ + if(th < GSL_ROOT4_DBL_EPSILON) { + B00 = (1.0 + th*th/15.0)/24.0; + pre = 1.0 + th*th/12.0; + } + else { + double sin_th = sqrt(1.0 - x*x); + double cot_th = x / sin_th; + B00 = 1.0/8.0 * (1.0 - th * cot_th) / (th*th); + pre = sqrt(th/sin_th); + } + + c1 = th/u * B00; + + result->val = pre * (J0.val + c1 * Jm1.val); + result->err = pre * (J0.err + fabs(c1) * Jm1.err); + result->err += GSL_SQRT_DBL_EPSILON * fabs(result->val); + + return GSL_ERROR_SELECT_2(stat_J0, stat_Jm1); +# else + fprintf (stderr, "legendre polynomials are not defined for l >= 10000\n"); + return NAN; +# endif +} + +/* If l is large and m is large, then we have to worry + * about overflow. Calculate an approximate exponent which + * measures the normalization of this thing. + */ + +double legendre_Plm (int l, int m, double x, double *err) { + + // check on the domain + if (m < 0) return NAN; + if (l < m) return NAN; + if (x < -1.0) return NAN; + if (x > +1.0) return NAN; + + const double dif = l-m; + const double sum = l+m; + const double t_d = ( dif == 0.0 ? 0.0 : 0.5 * dif * (log(dif)-1.0) ); + const double t_s = ( dif == 0.0 ? 0.0 : 0.5 * sum * (log(sum)-1.0) ); + const double exp_check = 0.5 * log(2.0*l+1.0) + t_d - t_s; + + if (exp_check < GSL_LOG_DBL_MIN + 10.0) return NAN; + + /* Account for the error due to the + * representation of 1-x. + */ + const double err_amp = 1.0 / (GSL_DBL_EPSILON + fabs(1.0-fabs(x))); + + /* P_m^m(x) and P_{m+1}^m(x) */ + double p_mm = legendre_Pmm(m, x); + double p_mmp1 = x * (2*m + 1) * p_mm; + + double value; + + if (l == m){ + value = p_mm; + if (err) *err = err_amp * 2.0 * GSL_DBL_EPSILON * fabs(p_mm); + return value; + } + + if (l == m + 1) { + value = p_mmp1; + if (err) *err = err_amp * 2.0 * GSL_DBL_EPSILON * fabs(p_mmp1); + return value; + } + + /* upward recurrence: (l-m) P(l,m) = (2l-1) z P(l-1,m) - (l+m-1) P(l-2,m) + * start at P(m,m), P(m+1,m) + */ + + double p_LM2 = p_mm; + double p_LM1 = p_mmp1; + double p_L = 0.0; + + int L; + for (L = m + 2; L <= l; L++){ + p_L = (x*(2*L-1)*p_LM1 - (L+m-1)*p_LM2) / (L-m); + p_LM2 = p_LM1; + p_LM1 = p_L; + } + value = p_L; + + if (err) *err = err_amp * (0.5*(l-m) + 1.0) * GSL_DBL_EPSILON * fabs(p_L); + + return value; +} + +double legendre_Plm_sphere(int l, int m, double x, double *err) { + + // check on the domain + if (m < 0) return NAN; + if (l < m) return NAN; + if (x < -1.0) return NAN; + if (x > +1.0) return NAN; + + if (m == 0) { + + double Pvalue = legendre_Pl (l, x, err); + + double pre = sqrt((2.0*l + 1.0)/(4.0*M_PI)); + + double value = pre * Pvalue; + if (err) { + double tmp_err = pre * *err; + *err = tmp_err + 2.0 * GSL_DBL_EPSILON * fabs(value); + } + return value; + } + + if (x == 1.0 || x == -1.0) { + /* m > 0 here */ + double value = 0.0; + if (err) *err = 0.0; + return value; + } + + + /* m > 0 and |x| < 1 here */ + + /* Starting value for recursion. + * Y_m^m(x) = sqrt( (2m+1)/(4pi m) gamma(m+1/2)/gamma(m) ) (-1)^m (1-x^2)^(m/2) / pi^(1/4) + */ + + const double sgn = ( GSL_IS_ODD(m) ? -1.0 : 1.0); + const double y_mmp1_factor = x * sqrt(2.0*m + 3.0); + + double y_mmp1, y_mmp1_err; + + double *errptr = NULL; + double lncirc_err, lnpoch_err; + + // NOTE: only evaluate the error if requested + if (err) errptr = &lncirc_err + double lncirc = ohana_gsl_sf_log_1plusx(-x*x, errptr); + + if (err) errptr = &lnpoch_err + double lnpoch = ohana_gsl_sf_lnpoch(m, 0.5, errptr); /* Gamma(m+1/2)/Gamma(m) */ + + double lnpre_val = -0.25*M_LNPI + 0.5 * (lnpoch + m*lncirc); + double lnpre_err = NAN; + + if (err) { + lnpre_err = 0.25*M_LNPI*GSL_DBL_EPSILON + 0.5 * (lnpoch_err + fabs(m)*lncirc_err); + } + + /* Compute exp(ln_pre) with error term, avoiding call to gsl_sf_exp_err BJG */ + double ex_pre_val = exp(lnpre_val); + + double ex_pre_err = NAN; + if (err) { + ex_pre_err = 2.0*(sinh(lnpre_err) + GSL_DBL_EPSILON)*ex_pre_val; + } + + double sr = sqrt((2.0+1.0/m)/(4.0*M_PI)); + double y_mm = sgn * sr * ex_pre_val; + double y_mmp1 = y_mmp1_factor * y_mm; + + double y_mm_err = NAN; + double y_mmp1_err = NAN; + if (err) { + y_mm_err = 2.0 * GSL_DBL_EPSILON * fabs(y_mm) + sr * ex_pre_err; + y_mm_err *= 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-x)); + y_mmp1_err = fabs(y_mmp1_factor) * y_mm_err; + } + + if (l == m){ + double value = y_mm; + if (err) { + *err = y_mm_err; + *err += 2.0 * GSL_DBL_EPSILON * fabs(y_mm); + } + return value; + } + + if (l == m + 1) { + double value = y_mmp1; + if (err) { + *err = y_mmp1_err; + *err += 2.0 * GSL_DBL_EPSILON * fabs(y_mmp1); + } + return value; + } + + double y_L = 0.0; + double y_L_err = NAN; + + /* Compute Y_l^m, l > m+1, upward recursion on l. */ + int L; + for(L=m+2; L <= l; L++){ + const double rat1 = (double)(L-m)/(double)(L+m); + const double rat2 = (L-m-1.0)/(L+m-1.0); + const double factor1 = sqrt(rat1*(2.0*L+1.0)*(2.0*L-1.0)); + const double factor2 = sqrt(rat1*rat2*(2.0*L+1.0)/(2.0*L-3.0)); + y_L = (x*y_mmp1*factor1 - (L+m-1.0)*y_mm*factor2) / (L-m); + y_mm = y_mmp1; + y_mmp1 = y_L; + + if (err) { + y_L_err = 0.5*(fabs(x*factor1)*y_mmp1_err + fabs((L+m-1.0)*factor2)*y_mm_err) / fabs(L-m); + y_mm_err = y_mmp1_err; + y_mmp1_err = y_L_err; + } + } + + value = y_L; + if (err) *err = y_L_err + (0.5*(l-m) + 1.0) * GSL_DBL_EPSILON * fabs(y_L); + + return value; +}