(     Title:  Complex gamma function with reliable floats
       File:  zgamma-rf.fs
    Started:  May 6, 2011
    Revised:  May 8, 2011  posted to c.l.f.
    Revised:  May 9, 2011  updated complex-rf names
    Revised:  May 11, 2011 adapted for forth-gmpfr 0.8.6

Adapted for forth-gmpfr reliable floats and complex-rf by
David N. Williams, from code published on comp.lang.forth,
Complex Gamma Function, by Marcel Hendrix, 2010-11-02, and
adapted for generic ANS Forth and the FSL by K. Myneni,
2010-11-03.

Reference:
http://groups.google.com/group/comp.lang.forth/msg/41acabc5dc66daf7
)
BASE @
DECIMAL

s" complex-rf.fs" included  \ loads gmpfr.fs

\ 53 1023 set-IEEE-754  \ binary64

\ Private:

\ Coefficients used by the GNU Scientific Library
7 CONSTANT g
0 value    reflect

create p
  rf"  0.99999999999980993e0" rf,  rf" 676.5203681218851e0"   rf,
  rf" -1259.1392167224028e0"  rf,  rf" 771.32342877765313e0"  rf,
  rf" -176.61502916214059e0"  rf,  rf" 12.507343278686905e0"  rf,
  rf" -0.13857109526572012e0" rf,  rf" 9.9843695780195716e-6" rf,
  rf"  1.5056327351493116e-7" rf,

zvariable temp

\ Public:

pi>rf 2 rfu* rf-sqrt rfconstant SQRT_2PI
rf" 0.0"	     rfconstant 0.0
rf" 0.5"             rfconstant 0.5
rf" 1" rf" 0"        zconstant  1+0i

: rfloat[] ( base ix -- addr ) cells + ;

\ Sum = p0 + p1/(1+z) + p2/(2+z) + ...
: zgam_series_sum (rf: z -- z zsum )
        p rf@ 0.0
        g 2 + 1 DO
          zover     I n>rf  0.0  z+ 1/z
          p I rfloat[] rf@  0.0  z* z+
        LOOP ;

: zgamma (rf: z1 -- z2 )
        false to reflect
        zdup real 0.5 rf< IF    \ reflect the argument about z=1
          1+0i zswap z-
          true to reflect
        THEN

        1+0i z- zgam_series_sum     \ rf: z zsum

        zover 0.5 0.0 z+  g n>rf 0.0 z+ zswap   \ rf: z t zsum
        zover znegate zexp  z*  SQRT_2PI z*f temp z!  \ rf: z t
        zover 0.5 0.0 z+    z^  temp z@  z*

        zswap
        reflect IF
          znegate
          pi>rf z*f zsin  z*
          pi>rf 0.0 zswap z/
        ELSE
          zdrop
        THEN
;

\ Reset-Search-Order
BASE !

1 [IF] \ test
14 set-output-rfprecision
cr rf"  0"   rf"  -2e-7" zgamma zs.
cr rf" -2"   rf"  -2e-7" zgamma zs.
cr rf" -2"   rf"  -2"    zgamma zs.
cr rf"  0.5" rf"   0"    zgamma zs.
cr rf"  1"   rf"   0"    zgamma zs.
cr rf"  2"   rf"   0"    zgamma zs.
cr rf"  5"   rf"   0"    zgamma zs.
cr
[THEN]

0 [IF]
First result is rzgamma, second is zgamma with gforth,
K. Myneni, 2010-11-3.

-0.57721566637525E0 + i0.49999999999998E7
-0.57721566526503   + i4999999.99999981

 0.46236671242882E0 + i0.24999999999998E7
 0.462366680243462  + i2499999.99999982

 0.10897678196206E-1 - i0.52829660788900E-2
 0.0108976781962063 - i0.00528296607889003

 0.17724538509055E1 + i0.00000000000000E0
 1.77245385090552 + i0.  ( sqrt[pi] )

 0.10000000000000E1 + i0.00000000000000E0
 1. + i0.  ( 0! )

 0.10000000000000E1 + i0.00000000000000E0
 1. + i0.  ( 1! )
 
 0.24000000000000E2 + i0.00000000000000E0
 24. + i0. ( 4! )
[THEN]