* 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