# -*- 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 <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