# -*- org -*- #+CATEGORY: specfunc * Complex hypergeometric function 1F1 * Could probably return immediately for exact zeros in 3j,6j,9j functions. Easiest to implement for 3j. Note from Serge Winitzki : The package "matpack" (www.matpack.de) includes many special functions, also the 3j symbols. They refer to some quite complicated numerical methods using recursion relations to get the right answers for large momenta, and to 1975-1976 papers by Schulten and Gordon for the description of the algorithms. The papers can be downloaded for free at http://www.ks.uiuc.edu/Publications/Papers/ http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU76B http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU75A http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU75 * add Fresnel Integrals to specfunc. See TOMS 723 + 2 subsequent errata. * make mode variables consistent in specfunc -- some seem to be unnecessary from performance point of view since the speed difference is negligible. * From: "Alexander Babansky" To: "Brian Gough" Subject: Re: gsl-1.2 Date: Sun, 3 Nov 2002 14:15:15 -0500 Hi Brian, May I suggest you to add another function to gsl-1.2 ? It's a modified Ei(x) function: Em(x)=exp(-x)*Ei(x); As u might know, Ei(x) raises as e^x on the negative interval. Therefore, Ei(100) is very very large. But Ei(100)*exp(-100) = 0.010; Unfortunately, if u try x=800 u'll get overflow in Ei(800). but Ei(800)*exp(-800) should be around 0.0001; Modified function Em(x) is used in cos, sin integrals such as: int_0^\infinity dx sin(bx)/(x^2+z^2)=(1/2z)*(Em(bz)-Em(-bz)); int_0^\infinity dx x cos(bx)/(x^2+z^2)=(1/2)*(Em(bz)+Em(-bz)); One of possible ways to add it to the library is: Em(x) = - PV int_0^\infinity e^(-t)/(t+x) dt Sincerely, Alex DONE: Wed Nov 6 13:06:42 MST 2002 [GJ] ---------------------------------------------------------------------- The following should be finished before a 1.0 level release. * Implement the conicalP_sph_reg() functions. DONE: Fri Nov 6 23:33:53 MST 1998 [GJ] * Irregular (Q) Legendre functions, at least the integer order ones. More general cases can probably wait. DONE: Sat Nov 7 15:47:35 MST 1998 [GJ] * Make hyperg_1F1() work right. This is the last remaining source of test failures. The problem is with an unstable recursion in certain cases. Look for the recursion with the variable named "start_pair"; this is stupid hack to keep track of when the recursion result is going the wrong way for awhile by remembering the minimum value. An error estimate is amde from that. But it is just a hack. Somethign must be done abou that case. * Clean-up Coulomb wave functions. This does not mean completing a fully controlled low-energy evaluation, which is a larger project. DONE: Sun May 16 13:49:47 MDT 1999 [GJ] * Clean-up the Fermi-Dirac code. The full Fermi-Dirac functions can probably wait until a later release, but we should have at least the common j = integer and j = 1/2-integer cases for the 1.0 release. These are not too hard. DONE: Sat Nov 7 19:46:27 MST 1998 [GJ] * Go over the tests and make sure nothing is left out. * Sanitize all the error-checking, error-estimation, algorithm tuning, etc. * Fill out our scorecard, working from Lozier's "Software Needs in Special Functions" paper. * Final Seal of Approval This section has itself gone through several revisions (sigh), proving that the notion of done-ness is ill-defined. So it is worth stating the criteria for done-ness explicitly: o interfaces stabilized o error-estimation in place o all deprecated constructs removed o passes tests - airy.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - airy_der.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - airy_zero.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - atanint.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_I0.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_I1.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_In.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_Inu.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_J0.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_J1.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_Jn.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_Jnu.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_K0.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_K1.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_Kn.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_Knu.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_Y0.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_Y1.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_Yn.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_Ynu.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_amp_phase.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_i.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_j.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_k.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_olver.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_sequence.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_temme.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_y.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - bessel_zero.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - beta.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - chebyshev.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - clausen.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - coulomb.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - coulomb_bound.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - coupling.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - dawson.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - debye.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - dilog.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - elementary.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - ellint.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - elljac.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - erfc.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - exp.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - expint.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - expint3.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - fermi_dirac.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - gamma.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - gamma_inc.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - gegenbauer.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - hyperg.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - hyperg_0F1.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - hyperg_1F1.c - hyperg_2F0.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - hyperg_2F1.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - hyperg_U.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - laguerre.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - legendre_H3d.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - legendre_Qn.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - legendre_con.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - legendre_poly.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - log.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - poch.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - poly.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - pow_int.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - psi.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - result.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - shint.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - sinint.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - synchrotron.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - transport.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - trig.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: - zeta.c INTERFACES: ERRORESTIM: DEPRECATED: PASSTESTS: ---------------------------------------------------------------------- The following are important but probably will not see completion before a 1.0 level release. * Incomplete Fermi-Dirac functions. Other Fermi-Dirac functions, including the generic 1/2-integer case, which was not done. * Implement the low-energy regime for the Coulomb wave functions. This is fairly well understood in the recent literature but will require some detailed work. Specifically this means creating a drop-in replacement for coulomb_jwkb() which is controlled and extensible. * General Legendre functions (at least on the cut). This subsumes the toroidal functions, so we need not consider those separately. SLATEC code exists (originally due to Olver+Smith). * Characterize the algorithms. A significant fraction of the code is home-grown and it should be reviewed by other parties. ---------------------------------------------------------------------- The following are extra features which need not be implemented for a version 1.0 release. * Spheroidal wave functions. * Mathieu functions. * Weierstrass elliptic functions. ---------------------------------------------------------------------- Improve accuracy of ERF NNTP-Posting-Date: Thu, 11 Sep 2003 07:41:42 -0500 From: "George Marsaglia" Newsgroups: comp.lang.c References: Subject: Re: When (32-bit) double precision isn't precise enough Date: Thu, 11 Sep 2003 08:41:40 -0400 X-Priority: 3 X-MSMail-Priority: Normal X-Newsreader: Microsoft Outlook Express 6.00.2800.1158 X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2800.1165 Message-ID: Lines: 265 NNTP-Posting-Host: 68.35.247.101 X-Trace: sv3-4YY+jkhhdeQvGKAREa99vDBFHJoKVqVBdUTSuRxA71OwlgxX0uUFnKYs54FlnUs0Xb6BRngKigkd75d!tKin8l8rAQKylaP+4vzTI3AO33bivOw1lKDZUUtXe4lUMW1qn+goUp/Pfksstg== X-Complaints-To: abuse@comcast.net X-DMCA-Complaints-To: dmca@comcast.net X-Abuse-and-DMCA-Info: Please be sure to forward a copy of ALL headers X-Abuse-and-DMCA-Info: Otherwise we will be unable to process your complaint properly X-Postfilter: 1.1 Why most of those who deal with the normal integral in probability theory are still stuck with the historical baggage of the error function is a puzzle to me, as is the poor quality of the results one gets from standard library implementations of erf(). (One of the most common is based on ALGORITHM AS66, APPL. STATIST.(1973) Vol.22, .424 by HILL, which gives only 6-8 digit accuracy). Here is a listing of my method: /* Marsaglia Complementary Normal Distribution Function cPhi(x) = integral from x to infinity of exp(-.5*t^2)/sqrt(2*pi), x<15 15-digit accuracy for x<15, returns 0 for x>15. #include */ double cPhi(double x){ long double v[]={0.,.65567954241879847154L, .42136922928805447322L,.30459029871010329573L, .23665238291356067062L,.19280810471531576488L, .16237766089686746182L,.14010418345305024160L, .12313196325793229628L,.10978728257830829123L, .99028596471731921395e-1L,.90175675501064682280e-1L, .82766286501369177252e-1L,.76475761016248502993e-1L, .71069580538852107091e-1L,.66374235823250173591e-1L}; long double h,a,b,z,t,sum,pwr; int i,j; if(x>15.) return (0.); if(x<-15.) return (1.); j=fabs(x)+1.; z=j; h=fabs(x)-z; a=v[j]; b=z*a-1.; pwr=1.; sum=a+h*b; for(i=2;i<60;i+=2){ a=(a+z*b)/i; b=(b+z*a)/(i+1); pwr=pwr*h*h; t=sum; sum=sum+pwr*(a+h*b); if(sum==t) break; } sum=sum*exp(-.5*x*x-.91893853320467274178L); if(x<0.) sum=1.-sum; return ((double) sum); } */ end of listing */ The method is based on defining phi(x)=exp(-x^2)/sqrt(2pi) and R(x)=cPhi(x)/phi(x). The function R(x) is well-behaved and terms of its Taylor series are readily obtained by a two-term recursion. With an accurate representation of R(x) at ,say, x=0,1,2,...,15, a simple evaluation of the Taylor series at intermediate points provides up to 15 digits of accuracy. An article describing the method will be in the new version of my Diehard CDROM. A new version of the Diehard tests of randomness (but not yet the new DVDROM) is at http://www.csis.hku.hk/~diehard/ George Marsaglia