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