R. William Gosper <rwg at osots.com> · Tue, 17 Jan 2006 13:37:36 -0800 (PST)

The full implementation of 2F1 numerics, sfloat, dfloat, and bfloat,
including all screw cases, is something of a thesis project, at
least how I did it.  The following message to Frank Olver contains
a sketch of the approach, plus my trick for doing the hardest part.
Let me know if you get to the point of attacking Problem 2.
Once completed, the whole shebang can be slightly tweaked to do
3F2[1], of which 2F1[z] is a limiting case.
.....
--------------
Dear Professor Olver,
Dick Askey told me you'd know who should get this.  I'll be happy to provide
clarifications, and even happier if I can provide more formal assistance at
this late stage of the project.
--Bill Gosper

    Methods for numerically difficult cases of 2F1(a,b;c|z)

(Note:  this message must be viewed with its original linebreaks using a
 fixed pitch font, e.g. Courier New, in a window at least 88 characters
 wide.)

Problem 1:  z near cis(+-pi/3).  The linear transformations

  z  <-  1 - z,  z  <-  1/z,   z  <-  z/(z - 1), . . .

(A&S 15.3.3, ..., 15.3.9) all leave z near cis(+-pi/3), where the
series converge slowly, or not at all.  If we were very lucky,
c = 2 a  or 2 b  mod 1, and we could use the quadratic transformation

    z  <-  z^2/(4 z - 4),

(A&s 15.3.15), but this is only applicable in a two-dimensional,
measure zero subset of the three-dimensional parameter space.

Here is how Macsyma (and Mathematica, at least as of a few years ago)
handle this.  Rejoice:  we don't need to get lucky--we can have our
z^2/(4 z - 4) convergence unconditionally if we switch from a series
to a first order recurrence on three variables.  Specifically, let

                      [d , e , f ] = [0, 1, 0],
                        0   0   0

              [                              (k + c - b - a) d  z        ]
              [                                               k          ]
              [      (k + a) (k + b) z (e  - --------------------)       ]
              [                          k          1 - z                ]
              [      ---------------------------------------------       ]
              [                             c       c + 1                ]
              [              4 (k + 1) (k + -) (k + -----)               ]
              [                             2         2                  ]
 [ d      ]   [                                                          ]
 [  k + 1 ]   [                           a b d  z                       ]
 [        ]   [                                k                         ]
 [ e      ] = [        (k + a) (k + b) z (-------- + (k + c) e )         ] .
 [  k + 1 ]   [                            1 - z              k          ]
 [        ]   [        -----------------------------------------         ]
 [ f      ]   [                             c       c + 1                ]
 [  k + 1 ]   [              4 (k + 1) (k + -) (k + -----)               ]
              [                             2         2                  ]
              [                                                          ]
              [      d  (k ((c - b - a) z + k (z - 2) - c) - a b z)      ]
              [       k                                                  ]
              [ f  - ---------------------------------------------- + e  ]
              [  k                        c                            k ]
              [                    2 (k + -) (1 - z)                     ]
              [                           2                              ]

Then d_k and e_k approach 0 like (z^2/(4z-4))^k, and the running sum f_k
approaches 2F1(a,b;c|z).  E.g., to numerically test Macsyma's valuation
of 2F1(1/2,1/6;1/3|cis(pi/3)):
                                                               %i %pi
                                        %i %pi                 ------
                                        ------      1/3          12
                             1  1  1      3       (2    + 1) %e
(d200)             hyper_2f1(-, -, -, %e      ) = -------------------
                             2  6  3                      3/4
                                                         3

which Macsyma deems
(c201) dfloat(%)

(d201) 0.25659521709148d0 %i + 0.95762638716458d0 = 

                                   0.25659521709148d0 %i + 0.95762638716458d0

we implement the foregoing recurrence,
(c202) olv(a,b,c,z):=block([d:0,e:1,f:0,n],for k from 0 unless
 f=last([d,e,f]:[(k+a)*(k+b)*z*(d*(k+c-b-a)*z/(z-1)+e)/(4*(k+1)*(k+c/2)*(k+(c+1)/2)),
  (k+a)*(k+b)*z*(e*(k+c)-a*b*d*z/(z-1))/(4*(k+1)*(k+c/2)*(k+(c+1)/2)),
  d*(k*((c-b-a)*z-c)-a*b*z+k^2*(z-2))/(2*(k+c/2)*(z-1))+f+e]) do n:k,
[d,e,f,n])$

and then [d_n,e_n,f_n,n], where n:=max(k), are
(c203) olv(1/2,1/6,1/3,cis(dfloat(%pi/3)))

(d203) [2.47871813091201d-19 %i + 2.2745179020215d-21, 

  - 2.56722341361481d-19 %i - 9.66582652669975d-19, 

   0.25659521709148d0 %i + 0.95762638716458d0, 25]

I.e., about 18 decimals in 25 steps.  (With only rational functions.)

Incidentally, (d200) above is the special case a=1/2 of

                              %i %pi
                              ------
                    a  2 a      3
(d194) hyper_2f1(a, -, ---, %e      ) = 
                    3   3

                               %i %pi a
                               --------
                 a   1   a/3      6                1                       1
 sqrt(%pi) Gamma(- + -) 4    %e         (--------------------- + ---------------------)
                 3   2                         1        a + 2          2        a + 1
                                         Gamma(-) Gamma(-----)   Gamma(-) Gamma(-----)
                                               3          3            3          3
 --------------------------------------------------------------------------------------,
                                          a + 1
                                          -----
                                            2
                                         3
a result that would be painful to derive from A&S 15.1.31 for lack of a contiguous
companion.  I hope DLMF provides contiguous pairs wherever possible to facilitate
both numeric and symbolic computation of 2F1s.  On that note, the accelerator
recurrence above produces a b z 2F1(a+1,b+1;c+1|z)/c when given the alternative
initialization
                      [d , e , f ] = [1, 0, 0].
                        0   0   0
Problem 2:
To evaluate 2F1[z] when |z| > 1 - epsilon, one needs the (rational)
linear transformations into a pair of 2F1[1-z] or 2F1[1/z]
(A&S 15.3.6, 15.3.7).  But when certain linear combinations of the
parameters are integers, the transformations give (F - F)/0, and
when you take a limit, you effectively differentiate the terms of the
hypergeometric series wrt the summation index.  (A&S 15.3.11, 15.3.12,
 15.3.14).  Since the general term is a product of four Gamma functions
and z^n, this effectively multiplies the terms by a linear combination of four
diGamma functions and log z.

I can provide a completely rational recurrence (omitted for "brevity"), isomorphic
to the acceleration above, which computes the requisite

   inf              n
   ====  (a)  (b)  x  (log(x) - psi (n + c) + psi (n + b) + psi (n + a) - psi (n + 1))
   \        n    n                 0             0             0             0        
    >    -----------------------------------------------------------------------------
   /                                        (c)  n!
   ====                                        n
   n = 0

(and companion), circumventing all four digammas and the log.

Please let me know if you find any of this of interest.
--Bill Gosper
PS, the Airy Functions section looks gorgeous.

Subject:	Re: [Maxima] numerical evaluation of 2F1
Date:	Tue, 17 Jan 2006 13:37:36 -0800 (PST)
From:	R. William Gosper <rwg at osots.com>
To:	fateman at cs.berkeley.edu