| * Could probably return immediately for exact zeros in 3j,6j,9j |
| functions. Easiest to implement for 3j. |
| |
| Note from Serge Winitzki <serge@cosmos.phy.tufts.edu>: |
| |
| 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" <babansky@mail.ru> |
| To: "Brian Gough" <bjg@network-theory.co.uk> |
| 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" <geo@stat.fsu.edu> |
| Newsgroups: comp.lang.c |
| References: <t4J7b.18514$98.4310@nwrddc03.gnilink.net> |
| 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: <wq2dnQBikNwb8P2iU-KYvg@comcast.com> |
| 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 <math.h> |
| */ |
| |
| 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 |