Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #9299 > unrolled thread

ZGAMMA revisited

Started bymhx@iae.nl (Marcel Hendrix)
First post2012-01-30 01:05 +0200
Last post2012-02-05 08:01 -0800
Articles 5 — 2 participants

Back to article view | Back to comp.lang.forth

Contents

  ZGAMMA revisited mhx@iae.nl (Marcel Hendrix) - 2012-01-30 01:05 +0200
    Re: ZGAMMA revisited Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-01-30 05:21 -0800
      Re: ZGAMMA revisited mhx@iae.nl (Marcel Hendrix) - 2012-01-31 00:21 +0200
        Re: ZGAMMA revisited mhx@iae.nl (Marcel Hendrix) - 2012-02-04 22:11 +0200
          Re: ZGAMMA revisited Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-02-05 08:01 -0800

#9299 — ZGAMMA revisited

Frommhx@iae.nl (Marcel Hendrix)
Date2012-01-30 01:05 +0200
SubjectZGAMMA revisited
Message-ID<77980205028435@frunobulax.edu>
In 2011, there was an extended CLF thread on computing the ZGAMMA function.

Here's ZGAMMA once GMP, MPFR and MPC are made available to the Forth system.
The best set of precomputed constants (with known error bounds) I could find 
limit the ZGAMMA absolute accuracy to about 10^-32. Results are compared 
to Wolfram Alpha (Mathematica).

With respect to the syntax, note that the iForth MPC implementation prefixes
words with 'Z#' instead of 'Z'.

The series coefficients are found with Victor Toth's program (See Wikipedia
references). Implementing the latter code in Forth would mean adding vector 
manipulation to the complex number lexicon, which I found a bit too much for
now.

-marcel

-- ---------------------------
(*
 * LANGUAGE    : Forth
 * PROJECT     : iForth mpc binding
 * DESCRIPTION : ZGAMMA function
 * CATEGORY    : Multiple precision
 * AUTHOR      : Marcel Hendrix (Viktor Toth helped)
 * STARTED     : Saturday, January 21, 2012, 22:23
 * LAST CHANGE : Monday, January 23, 2012, 22:50, debugged with help of Viktor Toth
 * LAST CHANGE : Friday, January 27, 2012, 22:36, stack version
 *)

	NEEDS -mpc

	REVISION -zgamma "--- ZGAMMA              Version 2.00 ---"

	PRIVATES
DOC
(*
 Adapted for mpc reliable zfloats by Marcel Hendrix, from code 
 published on comp.lang.forth, Complex Gamma Function, 
 by Marcel Hendrix, 2010-11-02
 
 Coefficients computed with Paul Godfrey's algorithm.
 See Wikipedia article on the "Lanczos Approximation". 
 Viktor Toth's C++ implementation was used (after a debug session) to find suitable coefficients.
 Error should be less than 2.2e-34. At 19+17i the error is less than 1.5e-42.
*)
ENDDOC

CREATE p PRIVATE 
    HERE
	S" +2.506628274631000502415765284787870637797744129562922113629468995990881110091579" 	  Z#$,
	S" +28036067819.74631058890623306093743315094258174562411893178167594785029820343691" 	  Z#$,
	S" -216764945447.9755869479070651364242600873819229772036422932539103338394587590359" 	  Z#$,
	S" +761012054330.03945277892380396516383568687720603462649943464239776666289557561" 	  Z#$,
	S" -1605877031381.312566068303370014018139099488755887743760531057224792844245569506" 	  Z#$,
	S" +2272529090534.289746939568578464540304068026284134278423934821554819668011902319" 	  Z#$,
	S" -2278670107037.21881759388335364276809181501856030856262702587383408764591438163" 	  Z#$,
	S" +1667908002357.944532310778965036913434221785993961596904697164266287442696107954" 	  Z#$,
	S" -904943925880.4240647886338479311253638077507629604475254919583359513331257753362" 	  Z#$,
	S" +365912700027.1305257455310369346060251700082941010098710132183307816607289906837" 	  Z#$,
	S" -109997616204.9781865177477607977184817713664125903839197749902468856938625042953" 	  Z#$,
	S" +24348551356.55137104951171973668680413058186082622921754170066500983412568773559" 	  Z#$,
	S" -3901718994.80839765409884884266155117512356586927533266330933734889969290900669" 	  Z#$,
	S" +441185605.4397725284127049593852816175572412254568625243428474643320223860023753" 	  Z#$,
	S" -33945791.00736390050136485026584415793211077997267451669192064550871126457828082" 	  Z#$,
	S" +1689153.692762895209404107416511190726656623120943719309447720663543853259001793" 	  Z#$,
	S" -50600.94012073931556806265270806257866457444269544345034780336250443507642702667" 	  Z#$,
	S" +822.933843432667148223554485844040777863273975456492530425954140368874871160083"      Z#$,
	S" -6.221544351999599605772962778662317260603729847447341502208732474834838624590193"     Z#$,
	S" +0.01707716022400587822088240928659086747168574885167707124958428446802555083484932"   Z#$,
	S" -1.107293955647243700304586842608602996682084371750795405523313414922213678200919e-5"  Z#$,
	S" +7.200761897886692682089256765671462285946945437486892400210078640805620336060679e-10" Z#$,
     HERE SWAP - /MPC / =: szP  PRIVATE

S" 0.5"  F#IN F#CONSTANT F#0.5  PRIVATE
S" 0.5"  Z#IN Z#CONSTANT Z#0.5  PRIVATE
S" 23.0" Z#IN Z#CONSTANT Z#g    PRIVATE
F#PI F#LN  
F#0.0 R,I->Z# Z#CONSTANT Z#LNPI PRIVATE

-- sum = p0 + p1/(z+1) + p2/(z+2) + ...
: zgam_series ( z: z1 -- z2 )
	p Z#@ 
	1 szP 1- 
	       DO
		  Z#OVER I Z#N+ Z#1/Z
		  p I Z#[@] Z#* Z#+  
	 -1 +LOOP 
	 Z#NIP ; PRIVATE

: Z#GAMMA ( z -- gamma )
	FALSE LOCAL reflect 
	Z#DUP Z#REAL  F#0< IF ( reflect the argument about z=1 )
				  1+0i Z#SWAP Z#-
				  TRUE TO reflect
			ENDIF  

	1+0i Z#- Z#DUP zgam_series Z#LN ( -- z ln[R] ) 

	Z#OVER Z#0.5 Z#+  Z#g Z#+	( -- z ln[R] z+0.5+g ) \ "r"
	Z#DUP Z#LN 		  	( -- z ln[R] r ln[r] )
	3 Z#PICK Z#0.5 Z#+  Z#*	 	( -- z ln[R] r ln[r]*[z+0.5] )
	Z#SWAP Z#- Z#+			( -- z ln[R]+ln[r]*[z+0.5]-r ) \ "x"

        reflect IF 
		   Z#LNPI Z#SWAP Z#- 		( -- z ln[PI]-x )
		   Z#OVER F#PI Z#*F	     	( -- z ln[PI]-x z*pi )
		   Z#SIN Z#NEGATE Z#LN	Z#-	( -- z ln[PI]-x-ln[-sin[z*pi]] )
	     ENDIF 

	Z#EXP Z#NIP ;

NESTING @ 1 = 
 [IF] 
	CREATE ref
	S" ( -5.77215664901496561443468856202033878880433664190180e-1   4.99999999999980218880093441334274561554458356403277e+6  )" Z#$,
	S" (  4.613921675492047501466772884439451366959286255331092e-1  2.499999999999812676750202441499863808634588595751551e+6 )" Z#$,
	S" (  1.089767819620633575931431799910101541141356665833077e-2 -5.282966078890027084097283638947329682805866071495392e-3 )" Z#$,
	S" (  0.999999999999960437760186882668549123108916712806555     1.154431329802993122886937712404067757760867328380360e-7 )" Z#$,
	S" (  1.519040026700361374481609505450015036681862641859509e-1 -1.980488016185498197191013167096389454801612622462159e-2 )" Z#$,
	S" 1.772453850905516027298167483341145182795171595237654e+0" 								    Z#$,
	S"  1" 															    Z#$,
	S"  1" 															    Z#$,
	S" 24" 															    Z#$,

	CR #40 TO Z#DECIMALS

	: tt ( n -- ) ( Z: a -- ) 
		CR ." iForth      " Z#DUP Z#OUT
		CR ." Mathematica " ref [@]Z# Z#DUP Z#OUT 
		CR ." error       " Z#- Z#OUT ;

	CR S" (0   -2e-7)" Z#IN Z#GAMMA  0 tt 
	CR S" (-2  -2e-7)" Z#IN Z#GAMMA  1 tt
	CR S" (-2  -2   )" Z#IN Z#GAMMA  2 tt
	CR S" (1   -2e-7)" Z#IN Z#GAMMA  3 tt
	CR S" (1   -2   )" Z#IN Z#GAMMA  4 tt
	CR S" (0.5  0   )" Z#IN Z#GAMMA  5 tt
	CR S" (1    0   )" Z#IN Z#GAMMA  6 tt
	CR S" (2    0   )" Z#IN Z#GAMMA  7 tt
	CR S" (5    0   )" Z#IN Z#GAMMA  8 tt
[THEN]

DOC
(*
	iForth      -5.7721566490149656144346885620203387878506e-0001  4.9999999999998021888009344133427456154536e+0006 i
	Mathematica -5.7721566490149656144346885620203387888043e-0001  4.9999999999998021888009344133427456155445e+0006 i
	error        9.5367146522275018400089657754154308883458e-0038 -9.0914641332667993341422709305607169289612e-0032 i

	iForth       4.6139216754920475014667728844394513669592e-0001  2.4999999999998126767502024414998638086345e+0006 i
	Mathematica  4.6139216754920475014667728844394513669592e-0001  2.4999999999998126767502024414998638086345e+0006 i
	error        1.2565947847219021411854512457667087712921e-0043  7.7302023839869337108258467173788247688924e-0046 i

	iForth       1.0897678196206335759314317999101015411413e-0002 -5.2829660788900270840972836389473296828059e-0003 i
	Mathematica  1.0897678196206335759314317999101015411413e-0002 -5.2829660788900270840972836389473296828058e-0003 i
	error        1.9914747832878888598574857110895654796432e-0044 -1.0415876539349819141434225526834450092155e-0043 i

	iForth       9.9999999999996043776018688266854912310891e-0001  1.1544313298029931228869377124040677577602e-0007 i
	Mathematica  9.9999999999996043776018688266854912310891e-0001  1.1544313298029931228869377124040677577608e-0007 i
	error        2.0953967662390846752010419497533152636595e-0052 -6.3559675457835218371148308564662885002114e-0047 i

	iForth       1.5190400267003613744816095054500150366871e-0001 -1.9804880161854981971910131670963894548073e-0002 i
	Mathematica  1.5190400267003613744816095054500150366818e-0001 -1.9804880161854981971910131670963894548016e-0002 i
	error        5.3198244412282749886214173119895814375263e-0040 -5.7662922798013562829869479065228888979193e-0041 i

	iForth       1.7724538509055160272981674833411451827951e+0000  0.0000000000000000000000000000000000000000e-0001 i
	Mathematica  1.7724538509055160272981674833411451827951e+0000  0.0000000000000000000000000000000000000000e-0001 i
	error       -8.9801078728205359821693358142637695652489e-0053  0.0000000000000000000000000000000000000000e-0001 i

	iForth       1.0000000000000000000000000000000000000000e+0000  0.0000000000000000000000000000000000000000e-0001 i
	Mathematica  1.0000000000000000000000000000000000000000e+0000  0.0000000000000000000000000000000000000000e-0001 i
	error        1.0047885858726272361808591138953963396086e-0068  0.0000000000000000000000000000000000000000e-0001 i

	iForth       1.0000000000000000000000000000000000000000e+0000  0.0000000000000000000000000000000000000000e-0001 i
	Mathematica  1.0000000000000000000000000000000000000000e+0000  0.0000000000000000000000000000000000000000e-0001 i
	error        7.8875878690481839035512465908932195399229e-0068  0.0000000000000000000000000000000000000000e-0001 i

	iForth       2.4000000000000000000000000000000000000000e+0001  0.0000000000000000000000000000000000000000e-0001 i
	Mathematica  2.4000000000000000000000000000000000000000e+0001  0.0000000000000000000000000000000000000000e-0001 i
	error        6.6975457978667602462473517568207380320082e-0065  0.0000000000000000000000000000000000000000e-0001 i
*)
ENDDOC

:ABOUT 	CR ." ZGAMMA -- zgamma using arbitrary precision complex FP package based on GMP, MPFR and MPC" ;

NESTING @ 1 = [IF]	.ABOUT -zgamma CR   [THEN]
			DEPRIVE

                              (* End of Source *)

[toc] | [next] | [standalone]

#9308

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2012-01-30 05:21 -0800
Message-ID<bba646d9-97f7-4473-839a-6cf446a880f0@q8g2000yqa.googlegroups.com>
In reply to#9299
On Jan 29, 5:05 pm, m...@iae.nl (Marcel Hendrix) wrote:
> In 2011, there was an extended CLF thread on computing the ZGAMMA function.
>
> Here's ZGAMMA once GMP, MPFR and MPC are made available to the Forth system.

What's MPC? A multi-precision complex library. Is it part of the
package released by David Williams?

> The best set of precomputed constants (with known error bounds) I could find
> limit the ZGAMMA absolute accuracy to about 10^-32. Results are compared
> to Wolfram Alpha (Mathematica).
>
> With respect to the syntax, note that the iForth MPC implementation prefixes
> words with 'Z#' instead of 'Z'.
>
> The series coefficients are found with Victor Toth's program (See Wikipedia
> references). Implementing the latter code in Forth would mean adding vector
> manipulation to the complex number lexicon, which I found a bit too much for
> now.
>
> -marcel
>

Looks great. Have you checked its performance very near its poles? As
I recall, that's where the TOMS programs had serious issues.

Krishna

[toc] | [prev] | [next] | [standalone]

#9321

Frommhx@iae.nl (Marcel Hendrix)
Date2012-01-31 00:21 +0200
Message-ID<97140304028435@frunobulax.edu>
In reply to#9308
Krishna Myneni <krishna.myneni@ccreweb.org> writes Re: ZGAMMA revisited

> On Jan 29, 5:05=A0pm, mhx@iae.nl (Marcel Hendrix) wrote:
>> In 2011, there was an extended CLF thread on computing the ZGAMMA function.

>> Here's ZGAMMA once GMP, MPFR and MPC are made available to the Forth system.

> What's MPC? A multi-precision complex library. Is it part of the
> package released by David Williams?

No, but I ported it after I finished adapting David's package for Windows
(GMP and MPFR wrongly assume long ints are 64 bit wide). MPC is a nice companion 
package, but it does (did:-) not have complex gamma and error functions.
[..]
> Looks great. Have you checked its performance very near its poles? As
> I recall, that's where the TOMS programs had serious issues.

I have (now) tested it with the parameter values listed on 
ftp://ccreweb.org/software/fsl/extras/. I went and fetched 
71-decimal-place ref result values from Wolfram/Alpha. 

As you can see below, the maximum absolute error is 10^-31 around 
the origin, which is the maximum error Glendon Pugh designed them 
for (in 'An analysis of the Lanczos Gamma approximation' (PhD thesis, 2004)).

The ZWOFZ results are also listed, although I didn't show yet how to 
implement this word using MPC.

-marcel

-- ----------------------------------------------------

Reference values for the complex Gamma function
from http://wolframalpha.com [71 places]

z#gamma(  0.00e-0001 -1.99e-0007 i) error =  9.09e-0032
z#gamma( -2.00e+0000 -1.99e-0007 i) error =  1.25e-0043
z#gamma( -2.00e+0000 -2.00e+0000 i) error =  1.06e-0043
z#gamma(  1.00e+0000 -1.99e-0007 i) error =  6.35e-0047
z#gamma(  1.00e+0000 -2.00e+0000 i) error =  5.35e-0040
z#gamma(  5.00e-0001  0.00e-0001 i) error =  8.98e-0053
z#gamma(  1.00e+0000  0.00e-0001 i) error =  1.00e-0068
z#gamma(  2.00e+0000  0.00e-0001 i) error =  7.88e-0068
z#gamma(  5.00e+0000  0.00e-0001 i) error =  6.69e-0065
z#gamma( -2.50e+0000 -2.50e+0000 i) error =  1.43e-0044
z#gamma( -2.50e+0000 -2.00e+0000 i) error =  1.57e-0044
z#gamma( -2.50e+0000 -1.50e+0000 i) error =  1.69e-0044
z#gamma( -2.50e+0000 -1.00e+0000 i) error =  1.79e-0044
z#gamma( -2.50e+0000 -5.00e-0001 i) error =  1.84e-0044
z#gamma( -2.50e+0000  0.00e-0001 i) error =  1.86e-0044
z#gamma( -2.50e+0000  5.00e-0001 i) error =  1.84e-0044
z#gamma( -2.00e+0000 -3.00e+0000 i) error =  8.58e-0044
z#gamma( -2.00e+0000 -2.00e+0000 i) error =  1.06e-0043
z#gamma( -2.00e+0000 -1.00e+0000 i) error =  1.20e-0043
z#gamma( -2.00e+0000 -5.00e-0001 i) error =  1.24e-0043
z#gamma( -2.00e+0000 -1.00e-0001 i) error =  1.25e-0043
z#gamma( -2.00e+0000 -9.99e-0003 i) error =  1.25e-0043
z#gamma( -2.00e+0000 -1.00e-0004 i) error =  1.25e-0043
z#gamma( -2.00e+0000 -1.00e-0006 i) error =  1.25e-0043
z#gamma( -2.00e+0000 -1.99e-0007 i) error =  1.25e-0043
z#gamma( -2.00e+0000  4.00e+0000 i) error =  6.40e-0044
z#gamma( -1.50e+0000  0.00e-0001 i) error =  3.34e-0041
z#gamma( -5.00e-0001 -5.00e-0001 i) error =  4.29e-0041
z#gamma( -5.00e-0001  0.00e-0001 i) error =  5.55e-0041
z#gamma( -5.00e-0001  5.00e-0001 i) error =  4.29e-0041
z#gamma(  0.00e-0001 -1.50e+0000 i) error =  1.09e-0038
z#gamma(  0.00e-0001 -1.00e+0000 i) error =  1.74e-0038
z#gamma(  0.00e-0001 -5.00e-0001 i) error =  3.59e-0038
z#gamma(  0.00e-0001  5.00e-0001 i) error =  3.59e-0038
z#gamma(  0.00e-0001  1.00e+0000 i) error =  1.74e-0038
z#gamma(  5.00e-0001 -5.00e-0001 i) error =  2.35e-0039
z#gamma(  5.00e-0001  0.00e-0001 i) error =  2.32e-0039
z#gamma(  5.00e-0001  5.00e-0001 i) error =  2.35e-0039
z#gamma(  1.00e+0000  0.00e-0001 i) error =  1.00e-0068
z#gamma(  1.50e+0000  0.00e-0001 i) error =  1.55e-0041
z#gamma(  2.00e+0000  0.00e-0001 i) error =  7.88e-0068
z#gamma(  2.50e+0000  0.00e-0001 i) error =  6.16e-0041
z#gamma(  3.00e+0000 -1.00e+0000 i) error =  1.20e-0042
z#gamma(  3.00e+0000  0.00e-0001 i) error =  6.84e-0067
z#gamma(  4.00e+0000  0.00e-0001 i) error =  6.49e-0066
z#gamma(  5.00e+0000  0.00e-0001 i) error =  6.69e-0065

z#wofz( -2.00e+0000 -2.00e+0000 i) error =  1.98e-0056
z#wofz( -2.00e+0000 -1.50e+0000 i) error =  2.99e-0060
z#wofz( -2.00e+0000 -1.00e+0000 i) error =  2.75e-0057
z#wofz( -2.00e+0000 -5.00e-0001 i) error =  5.05e-0054
z#wofz( -2.00e+0000  0.00e-0001 i) error =  1.00e-0048
z#wofz( -2.00e+0000  5.00e-0001 i) error =  5.05e-0054
z#wofz( -2.00e+0000  1.00e+0000 i) error =  2.75e-0057
z#wofz( -2.00e+0000  1.50e+0000 i) error =  2.99e-0060
z#wofz( -2.00e+0000  2.00e+0000 i) error =  1.98e-0056
z#wofz( -1.50e+0000 -2.00e+0000 i) error =  3.83e-0058
z#wofz( -1.50e+0000 -1.50e+0000 i) error =  1.43e-0059
z#wofz( -1.50e+0000 -1.00e+0000 i) error =  3.12e-0058
z#wofz( -1.50e+0000 -5.00e-0001 i) error =  1.04e-0075
z#wofz( -1.50e+0000  0.00e-0001 i) error =  8.63e-0078
z#wofz( -1.50e+0000  5.00e-0001 i) error =  1.05e-0075
z#wofz( -1.50e+0000  1.00e+0000 i) error =  3.12e-0058
z#wofz( -1.50e+0000  1.50e+0000 i) error =  1.43e-0059
z#wofz( -1.50e+0000  2.00e+0000 i) error =  3.83e-0058
z#wofz( -1.00e+0000 -2.00e+0000 i) error =  2.25e-0057
z#wofz( -1.00e+0000 -1.50e+0000 i) error =  5.82e-0061
z#wofz( -1.00e+0000 -1.00e+0000 i) error =  3.86e-0077
z#wofz( -1.00e+0000 -5.00e-0001 i) error =  9.24e-0076
z#wofz( -1.00e+0000  0.00e-0001 i) error =  4.99e-0071
z#wofz( -1.00e+0000  5.00e-0001 i) error =  9.19e-0076
z#wofz( -1.00e+0000  1.00e+0000 i) error =  5.03e-0077
z#wofz( -1.00e+0000  1.50e+0000 i) error =  5.82e-0061
z#wofz( -1.00e+0000  2.00e+0000 i) error =  2.25e-0057
z#wofz( -5.00e-0001 -2.00e+0000 i) error =  6.34e-0057
z#wofz( -5.00e-0001 -1.50e+0000 i) error =  6.48e-0061
z#wofz( -5.00e-0001 -1.00e+0000 i) error =  0.00e-0001
z#wofz( -5.00e-0001 -5.00e-0001 i) error =  2.21e-0074
z#wofz( -5.00e-0001  0.00e-0001 i) error =  2.26e-0068
z#wofz( -5.00e-0001  5.00e-0001 i) error =  2.21e-0074
z#wofz( -5.00e-0001  1.00e+0000 i) error =  1.38e-0077
z#wofz( -5.00e-0001  1.50e+0000 i) error =  6.48e-0061
z#wofz( -5.00e-0001  2.00e+0000 i) error =  6.34e-0057
z#wofz(  0.00e-0001 -2.00e+0000 i) error =  1.51e-0059
z#wofz(  0.00e-0001 -1.50e+0000 i) error =  6.64e-0061
z#wofz(  0.00e-0001 -1.00e+0000 i) error =  0.00e-0001
z#wofz(  0.00e-0001 -5.00e-0001 i) error =  8.99e-0075
z#wofz(  0.00e-0001  0.00e-0001 i) error =  0.00e-0001
z#wofz(  0.00e-0001  5.00e-0001 i) error =  9.01e-0075
z#wofz(  0.00e-0001  1.00e+0000 i) error =  1.72e-0077
z#wofz(  0.00e-0001  1.50e+0000 i) error =  6.64e-0061
z#wofz(  0.00e-0001  2.00e+0000 i) error =  1.51e-0059
z#wofz(  5.00e-0001 -2.00e+0000 i) error =  6.34e-0057
z#wofz(  5.00e-0001 -1.50e+0000 i) error =  6.48e-0061
z#wofz(  5.00e-0001 -1.00e+0000 i) error =  0.00e-0001
z#wofz(  5.00e-0001 -5.00e-0001 i) error =  2.21e-0074
z#wofz(  5.00e-0001  0.00e-0001 i) error =  2.26e-0068
z#wofz(  5.00e-0001  5.00e-0001 i) error =  2.21e-0074
z#wofz(  5.00e-0001  1.00e+0000 i) error =  1.38e-0077
z#wofz(  5.00e-0001  1.50e+0000 i) error =  6.48e-0061
z#wofz(  5.00e-0001  2.00e+0000 i) error =  6.34e-0057
z#wofz(  1.00e+0000 -2.00e+0000 i) error =  2.25e-0057
z#wofz(  1.00e+0000 -1.50e+0000 i) error =  5.82e-0061
z#wofz(  1.00e+0000 -1.00e+0000 i) error =  3.86e-0077
z#wofz(  1.00e+0000 -5.00e-0001 i) error =  9.24e-0076
z#wofz(  1.00e+0000  0.00e-0001 i) error =  4.99e-0071
z#wofz(  1.00e+0000  5.00e-0001 i) error =  9.19e-0076
z#wofz(  1.00e+0000  1.00e+0000 i) error =  5.03e-0077
z#wofz(  1.00e+0000  1.50e+0000 i) error =  5.82e-0061
z#wofz(  1.00e+0000  2.00e+0000 i) error =  2.25e-0057
z#wofz(  1.50e+0000 -2.00e+0000 i) error =  3.83e-0058
z#wofz(  1.50e+0000 -1.50e+0000 i) error =  1.43e-0059
z#wofz(  1.50e+0000 -1.00e+0000 i) error =  3.12e-0058
z#wofz(  1.50e+0000 -5.00e-0001 i) error =  1.04e-0075
z#wofz(  1.50e+0000  0.00e-0001 i) error =  8.63e-0078
z#wofz(  1.50e+0000  5.00e-0001 i) error =  1.05e-0075
z#wofz(  1.50e+0000  1.00e+0000 i) error =  3.12e-0058
z#wofz(  1.50e+0000  1.50e+0000 i) error =  1.43e-0059
z#wofz(  1.50e+0000  2.00e+0000 i) error =  3.83e-0058
z#wofz(  2.00e+0000 -2.00e+0000 i) error =  1.98e-0056
z#wofz(  2.00e+0000 -1.50e+0000 i) error =  2.99e-0060
z#wofz(  2.00e+0000 -1.00e+0000 i) error =  2.75e-0057
z#wofz(  2.00e+0000 -5.00e-0001 i) error =  5.05e-0054
z#wofz(  2.00e+0000  0.00e-0001 i) error =  1.00e-0048
z#wofz(  2.00e+0000  5.00e-0001 i) error =  5.05e-0054
z#wofz(  2.00e+0000  1.00e+0000 i) error =  2.75e-0057
z#wofz(  2.00e+0000  1.50e+0000 i) error =  2.99e-0060
z#wofz(  2.00e+0000  2.00e+0000 i) error =  1.98e-0056

[toc] | [prev] | [next] | [standalone]

#9369

Frommhx@iae.nl (Marcel Hendrix)
Date2012-02-04 22:11 +0200
Message-ID<01921399018435@frunobulax.edu>
In reply to#9321
mhx@iae.nl (Marcel Hendrix) wrote Re: ZGAMMA revisited

> Krishna Myneni <krishna.myneni@ccreweb.org> writes Re: ZGAMMA revisited
[..]
>> Looks great. Have you checked its performance very near its poles? As
>> I recall, that's where the TOMS programs had serious issues.

> I have (now) tested it with the parameter values listed on 
> ftp://ccreweb.org/software/fsl/extras/. I went and fetched 
> 71-decimal-place ref result values from Wolfram/Alpha. 

> As you can see below, the maximum absolute error is 10^-31 around 
> the origin, which is the maximum error Glendon Pugh designed them 
> for (in 'An analysis of the Lanczos Gamma approximation' (PhD thesis, 2004)).

The accuracy has increased to 1e-41. I have added Godfrey's/Toth's 
P-generator, so it is now, in principal (*), possible to configure for almost 
arbitrary precision.

For brevity, I have deleted the routine to find the optimal g constant
(g should be somewhat smaller than n). G depends on n and only needs to
be computed once.

I find that 41 decimal digits is about the maximum what can be done with
reasonable computation speed. According to Godfrey, the matrix P is very
difficult to compute accurately (quad-precision for double precision 
results). Even with 2048 bit variables I don't get much more than 41 digits
of precision. However, the results are better than those of Victor Toth's 
program. Maybe one of the constants or literals has too low precision, but
I can't seem to find it.

-marcel

-- -----------------------------
(*
 * LANGUAGE    : Forth
 * PROJECT     : iForth mpc binding
 * DESCRIPTION : Z#GAMMA and Z#WOFZ (Faddeeva) functions, Z#ERF, Z#ERFC
 * CATEGORY    : Multiple precision
 * AUTHOR      : Marcel Hendrix
 * AUTHOR      : Krishna Myneni, ZWOFZ V1.0, 24 January 2011 
 * AUTHOR      : Godfrey, Paul (2001). "Lanczos Implementation of the Gamma Function"
 * AUTHOR      : Toth, Viktor (2005). "Programmable Calculators: The Lanczos Approximation"
 * STARTED     : Saturday, January 21, 2012, 22:23
 * LAST CHANGE : Friday, January 27, 2012, 22:36, stack version
 * LAST CHANGE : Saturday, February 04, 2012, 20:54; error < 1e-41, testing, documentation
 *)

	NEEDS -mpc

DOC
(*
 Adapted for mpc reliable zfloats by Marcel Hendrix, from code 
 published on comp.lang.forth, Complex Gamma Function, 
 by Marcel Hendrix, 2010-11-02
 
 Coefficients computed with Paul Godfrey's algorithm, 
   http://www.numericana.com/answer/info/godfrey.htm, 
 on the "Lanczos Approximation". 
 
 Error should be less than 4.67e-0041. 
 At 19+17i the error is less than 1.45e-41.

 Paul Godfrey's comments on c=B*B*C*F:
 -------------------------------------
 c is calculated from c=D*B*C*F where D is [1 -Sloane's sequence A002457],
 fact(2n+2)/(2fact(n)fact(n+1)), diagonal, and B is rows from the odd 
 cols of Abramowitz and Stegun table 24.1 (Binomial), unit Upper 
 triangular, and C is cols from the even cols of Abramowitz and Stegun 
 table 22.3 (Cheby), C(1)=1/2, Lower triangular, and F is 
 sqrt(2)/pi*gamma(z+0.5).*exp(z+g+0.5).*(z+g+0.5).^-(z+0.5).
 
 gamma(z+0.5) can be computed using factorials, 2^z, and sqrt(pi).
 Abramowitz and Stegun formulas 6.1.8 and 6.1.12.
*)
ENDDOC


#22 =: n  		        	
S" 0.5"  Z#IN Z#CONSTANT Z#0.5  
S" 21.8799999999999990052401699358597397804260253906250" 
	 Z#IN Z#VALUE    Z#g      
F#PI F#LN  
F#0.0 R,I->Z# Z#CONSTANT Z#LNPI 

n     Z#BLOCK P    
: P[r]@    ( r -- )		  P [@]Z# ;  
: P[r]!    ( r -- )		  P [!]Z# ;

n n * F#BLOCK B    
: B[r,c]@  ( r c -- ) SWAP n * +  B [@]F# ;  
: B[r,c]!  ( r c -- ) SWAP n * +  B [!]F# ;

n n * F#BLOCK C    
: C[r,c]@  ( r c -- ) SWAP n * +  C [@]F# ;  
: C[r,c]!  ( r c -- ) SWAP n * +  C [!]F# ;

n n * F#BLOCK D    
: D[r,c]@  ( r c -- ) SWAP n * +  D [@]F# ;  
: D[r,c]!  ( r c -- ) SWAP n * +  D [!]F# ;

n     F#BLOCK F    
: F[r]@    ( r -- )   	          F [@]F# ;  
: F[r]!    ( r -- )	          F [!]F# ;  

n     F#BLOCK tmp1 
: tmp1[r]@ ( r -- )   	       tmp1 [@]F# ;  
: tmp1[r]! ( r -- )            tmp1 [!]F# ;

n     F#BLOCK tmp2 
: tmp2[r]@ ( r -- )   	       tmp2 [@]F# ;	
: tmp2[r]! ( r -- )            tmp2 [!]F# ;

: COMBINE ( n k -- )
	DUP 0< IF  2DROP F#0.0 EXIT  ENDIF
	LOCALS| k n |
	F#1.0 n k - 1+ 1 ?DO  n 1+ I - F#N*  I F#N/  LOOP ;

: makeB ( -- )
	n 0 DO  F#1.0  0 I B[r,c]!  LOOP
	n 1 DO 	
		n 0 DO  
			J I + 1-  I J - COMBINE  
			J I + 1 AND IF  F#NEGATE  ENDIF
			J I B[r,c]!
		  LOOP
	  LOOP ;

: makeC ( -- )
	n 1 DO 
		I 1+ 0 DO  
			   F#0.0 
			   J 1+ 0 DO  
			   	      K 2*  I 2*      COMBINE 
				      I     I J + K - COMBINE F#*
				      F#+
			   	LOOP
			   J I - 1 AND IF  F#NEGATE  ENDIF
			   J I C[r,c]!  
		     LOOP
	  LOOP 
	F#0.5  0 0 C[r,c]! ;

: makeD ( -- )
	F#1.0           0 0 D[r,c]!
	F#1.0  F#NEGATE 1 1 D[r,c]!
	n 2 DO	
		I 1- I 1- D[r,c]@  I 2* 1- 2* F#N*  
		I 1- F#N/
		I I D[r,c]! 
	  LOOP ;

: makeF ( -- )
	n 0 DO 
		F#2.0  I F[r]!
		I 2* 1+  I 1+ ?DO  
				  J F[r]@ I F#N* 4 F#N/
				  J F[r]!
			     LOOP
		I N->F#  Z#g Z#REAL F#+  F#0.5 F#+
		I F[r]@	
			F#OVER F#EXP   F#*
			F#OVER I F#^N  F#/
			F#SWAP F#SQRT  F#/
		I F[r]! 
	  LOOP ;

: makeP  ( -- )	\ P = D*B*C*F
		makeB makeC makeD makeF
( C*F    )	n 0 DO 	F#0.0  n 0 DO  J I C[r,c]@  I    F[r]@ F#* F#+  LOOP  I tmp1[r]!  LOOP
( B*tmp1 )	n 0 DO 	F#0.0  n 0 DO  J I B[r,c]@  I tmp1[r]@ F#* F#+  LOOP  I tmp2[r]!  LOOP
( D*tmp2 )	n 0 DO 	F#0.0  n 0 DO  J I D[r,c]@  I tmp2[r]@ F#* F#+  LOOP  F#0.0 R,I->Z#  I P[r]!  LOOP ;

makeP  FORGET B

-- sum = p0 + p1/(z+1) + p2/(z+2) + ...
: zgam_series ( z: z1 -- z2 )
	0 P[r]@
	1 n 1- DO
		  Z#OVER I Z#N+ Z#1/Z
		  I P[r]@ Z#* Z#+  
	 -1 +LOOP 
	 Z#NIP ; 

: Z#GAMMA ( z -- gamma )
	Z#DUP Z#IMAG F#0= IF  Z#REAL F#GAMMA F#0.0 R,I->Z# EXIT  ENDIF
	FALSE LOCAL reflect 
	Z#DUP Z#REAL  F#0< IF ( reflect the argument about z=1 )
				  1+0i Z#SWAP Z#-
				  TRUE TO reflect
			ENDIF  

	1+0i Z#- Z#DUP zgam_series Z#LN ( -- z ln[R] ) 

	Z#OVER Z#0.5 Z#+  Z#g Z#+	( -- z ln[R] z+0.5+g ) \ "r"
	Z#DUP Z#LN 		  	( -- z ln[R] r ln[r] )
	3 Z#PICK Z#0.5 Z#+  Z#*	 	( -- z ln[R] r ln[r]*[z+0.5] )
	Z#SWAP Z#- Z#+			( -- z ln[R]+ln[r]*[z+0.5]-r ) \ "x"

        reflect IF 
		   Z#LNPI Z#SWAP Z#- 		( -- z ln[PI]-x )
		   Z#OVER F#PI Z#*F	     	( -- z ln[PI]-x z*pi )
		   Z#SIN Z#NEGATE Z#LN	Z#-	( -- z ln[PI]-x-ln[-sin[z*pi]] )
	     ENDIF 

	Z#EXP Z#NIP ;


1 [IF] 
	CREATE ref
	S" ( -0.5772156649015038444865707308395844477160077449056784923056111472009192607334784612103705009326        5.592405333333156476339539856304986168520135413649795772037789887519931076442722028285448224756e6        )" Z#$,
	S" ( -0.5772156649014965614434688562020338788804336641901804033533344132555560391778608380309006531595        4.999999999999802188800934413342745615544583564032775719864604802968974077221470016830354437666e6        )" Z#$,
	S" (  0.4613921675492047501466772884439451366959286255331092488308107054807886345899049001054408402412744     2.499999999999812676750202441499863808634588595751551773116557499745300126549085677827331229076830e6     )" Z#$,
	S" (  1.089767819620633575931431799910101541141356665833077e-2                                               -5.282966078890027084097283638947329682805866071495392e-3                                                 )" Z#$,
	S" (  0.999999999999960437760186882668549123108916712806555143972920960593794815444294003366070887533322259   1.154431329802993122886937712404067757760867328380360806706668826511112078355721676061801306319107296e-7 )" Z#$,
	S" (  0.1519040026700361374481609505450015036681862641859509057437762329750161509579256698043879906351668913 -0.01980488016185498197191013167096389454801612622462159210135851109453146547488858604011337290335010243  )" Z#$,
	S"  1.772453850905516027298167483341145182797549456122387128213807789852911284591032181374950656738544665416"  Z#$,
	S"  1" 	Z#$,
	S"  1" 	Z#$,
	S" 24" 	Z#$,

	CREATE ref2in
	S" (  19   17   )" Z#$,
	S" (  -2.5 -2.5 )" Z#$,	S" (  -2.5 -2	)" Z#$,	S" (  -2.5 -1.5	)" Z#$,	S" (  -2.5 -1	)" Z#$,	S" (  -2.5 -0.5	)" Z#$,	S" (  -2.5  0	)" Z#$,	S" (  -2.5  0.5 )" Z#$,	
	S" (  -2 -3	)" Z#$,	S" (  -2 -2	)" Z#$,	S" (  -2 -1	)" Z#$,	S" (  -2 -0.5	)" Z#$,	
	S" (  -2 -0.1	)" Z#$,	
	S" (  -2 -1e-2  )" Z#$,	
	S" (  -2 -1e-4  )" Z#$,	
	S" (  -2 -1e-6  )" Z#$,	
	S" (  -2 -2e-7  )" Z#$,	
	S" (  -2  4	)" Z#$,	
	S" (  -1.5  0	)" Z#$,	
	S" (  -0.5 -0.5	)" Z#$,	S" (  -0.5  0	)" Z#$,	S" (  -0.5  0.5	)" Z#$,	
	S" (   0 -1.5	)" Z#$,	S" (   0   -1	)" Z#$,	S" (   0   -0.5	)" Z#$,	S" (   0  0.5	)" Z#$,	S" (   0  1	)" Z#$,	
	S" ( 0.5 -0.5	)" Z#$,	S" (   0.5  0	)" Z#$,	S" (   0.5  0.5 )" Z#$, 
	S" (   1  0	)" Z#$,
	S" ( 1.5  0	)" Z#$,	
	S" (   2  0	)" Z#$,	S" (   2.5  0	)" Z#$,	
	S" (   3 -1	)" Z#$,	S" (   3    0	)" Z#$,	
	S" (   4  0	)" Z#$,	
	S" (   5  0	)" Z#$,

	CREATE ref2out
	S" (  1.66800622476031324355011890409209528835533764602489185010162635832930601822372682728912412120895e12 5.77782923985566521886203411555624804818352915811475661386507682039004814385018353984493480058434e12)" Z#$,
	S" (  0.0018596262288378007605140298772990657201833713854374686047506913935015066    0.0002712831060525439837242160962034060623576043192526247279453237439586030 )" Z#$,
	S" (  0.0045440424708998928448112357397385947420081062226643492249511950571780757    0.0047378548180488181976190390406144493612489679782855445314880866512142289 )" Z#$,
	S" (  0.003412139564239149028570842363649156500364328040716332665211214563015526     0.024053490434664735984426343338570327436132783210141801089253024645474196  )" Z#$,
	S" ( -0.041736625807893613744760138309780403748109093691845660691035034841702707     0.086369107369763484694186279347028210540939438970587338511624117856634249  )" Z#$,
	S" ( -0.33387520352243233740327727033956558807270634792337725625220954388092623      0.20645730796360841491828760756387298838346687701254812467624343946129235   )" Z#$,
	S" ( -0.945308720482941881225689324448610764158693043265273135047364154588219351     0   									 )" Z#$,
	S" ( -0.33387520352243233740327727033956558807270634792337725625220954388092623     -0.20645730796360841491828760756387298838346687701254812467624343946129235   )" Z#$,
	S" ( -0.00016317241827260774828275925258031819741311893448576591773591666572864989  -0.00112849591701795310547369969387328899536850056333539571044604506666388670)" Z#$,
	S" (  0.0108976781962063357593143179991010154114135666583307781282492546957629247   -0.0052829660788900270840972836389473296828058660714953920511962941564652637 )" Z#$,
	S" (  0.13390971760532574430161182219027076706630092506836253484120485881965141      0.09628651530237880980885565089138579075406019238256248836902837853950407   )" Z#$,
	S" (  0.32119401555445918511256007517484607980513263907752688365797023886197469      0.64091266903299752430095572023545107067767722234949743905269500723026429   )" Z#$,
	S" (  0.4542668315760802669659523134765663422166489510101463202539662319518575       4.9074302745784021858738556409335325139548120408090619790076212695993421    )" Z#$,
	S" (  0.461320126710157722258539005266352314016062534354685835412811483640234       49.990634940684132929533412542345511751978854511837161546043339886107730     )" Z#$,
	S" (  0.461392160344346133398958874035181401084136291681805580190                 4999.99990633837620450537169436781525862350137375010754841                     )" Z#$,
	S" (  0.461392167548513080945078063631450565143289051963880409811               499999.999999063383751013266656803608415650019482531266866                       )" Z#$,
	S" (  0.4613921675492047501466772884439451366959286255331092488308107054807886345899049001054408402412744     2.499999999999812676750202441499863808634588595751551773116557499745300126549085677827331229076830e6  )" Z#$,
	S" ( -0.000126872852422279549793478603992800073156747521802152331403311837522350325 -3.842307699536186907890816438311212389947979180635544628153796964917616e-6  )" Z#$,
	S" (  2.3632718012073547030642233111215269103967e+0000  0 )" Z#$,
	S" ( -1.5814778282557300107480456613162255350214059717569166462352202825058265       0.0548501708277647774074520857944242831967213736971639238624329503526707    )" Z#$,
	S" ( -3.5449077018110320545963349666822903655950e+0000  0 )" Z#$,
	S" ( -1.5814778282557300107480456613162255350214059717569166462352202825058265      -0.0548501708277647774074520857944242831967213736971639238624329503526707    )" Z#$,
	S" ( -0.03146902230831199010727996346048436773520532733941691082484206663955588      0.19142063360054488759000008529803738813356696165020554888022742167470925   )" Z#$,
	S" ( -0.15494982830181068512495513048388660519587965207932493026588027679886080      0.49801566811835604271369111746219809195296296758765009289264295499845830   )" Z#$,
	S" ( -0.3992794763291927125044534487971959663575337151385741121775945928369407       1.6033881941394344451955126231743084933936340977279108438291716209459245    )" Z#$,
	S" ( -0.3992794763291927125044534487971959663575337151385741121775945928369407      -1.6033881941394344451955126231743084933936340977279108438291716209459245    )" Z#$,
	S" ( -0.15494982830181068512495513048388660519587965207932493026588027679886080     -0.49801566811835604271369111746219809195296296758765009289264295499845830   )" Z#$,
	S" (  0.81816399954174739407774887355532490910906367272704028504882661642924863      0.76331382871398261667029678776090062591234229902987636118639366607657788   )" Z#$,
	S" (  1.7724538509055160272981674833411451827975e+0000  0 )" Z#$,
	S" (  0.81816399954174739407774887355532490910906367272704028504882661642924863     -0.76331382871398261667029678776090062591234229902987636118639366607657788   )" Z#$,
	S" (  1.0                        0.0 			  )" Z#$,
	S" (  8.8622692545275801364908374167057259139877e-0001  0 )" Z#$,
	S" (  1.0                        0.0 			  )" Z#$,
	S" (  1.3293403881791370204736256125058588870981e+0000  0 )" Z#$,
	S" (  0.9628651530237880980885565089138579075406019238256248836902837853950407      -1.3390971760532574430161182219027076706630092506836253484120485881965141    )" Z#$,
	S" (  2.0                        0.0  			  )" Z#$,
	S" (  6.0                        0.0 			  )" Z#$,
	S" (  24.0                       0.0 			  )" Z#$,

	: F#MAX ( f# a b -- c ) F#2DUP F#< IF F#NIP ELSE F#DROP ENDIF ;
	
	S" -Inf" F#IN F#CONSTANT F#-INF  
	S" +Inf" F#IN F#CONSTANT F#+INF  
	F#-INF        F#VALUE    ferror

	: tt ( n c-addr u -- )  
		Z#IN 
		Z#DUP CR ." z#gamma( " 2 Z#OUT.R
		Z#GAMMA ref [@]Z# Z#TUCK Z#- ." ) = (" Z#DUP #22 Z#OUT.R
		Z#ABS Z#ABS F#/ ferror F#OVER F#MAX TO ferror ." ) error = " 2 F#OUT.R ;

	: testz#gamma 
		F#-INF TO ferror
		CR ." *** Z#GAMMA Error Table ***"
		0 S" (0 -1.78813934326171875e-0007)" tt
		1 S" (0   -2e-7)" tt 
		2 S" (-2  -2e-7)" tt
		3 S" (-2  -2   )" tt
		4 S" (1   -2e-7)" tt
		5 S" (1   -2   )" tt
		6 S" (0.5  0   )" tt
		7 S" (1    0   )" tt
		8 S" (2    0   )" tt
		9 S" (5    0   )" tt
		#38 0 DO  ref2in  I Z#[@] CR ." z#gamma( " Z#DUP 2 Z#OUT.R ." ) = ("  Z#GAMMA 
			  ref2out I Z#[@] Z#TUCK Z#- Z#DUP #22 Z#OUT.R ." ) error = " 
			  Z#ABS Z#ABS F#/  ferror F#OVER F#MAX TO ferror 2 F#OUT.R 
		    LOOP 
		CR ." Maximum absolute error = " ferror 2 F#OUT.R ;

	CR testz#gamma
[THEN]

: ABOUT	CR ." Z#GAMMA -- Gamma using GMP, MPFR and MPC" ;

	ABOUT

[toc] | [prev] | [next] | [standalone]

#9415

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2012-02-05 08:01 -0800
Message-ID<5e8552aa-7438-4300-ba74-8982523c406e@g27g2000yqa.googlegroups.com>
In reply to#9369
On Feb 4, 2:11 pm, m...@iae.nl (Marcel Hendrix) wrote:
> m...@iae.nl (Marcel Hendrix) wrote Re: ZGAMMA revisited
...
> The accuracy has increased to 1e-41. I have added Godfrey's/Toth's
> P-generator, so it is now, in principal (*), possible to configure for almost
> arbitrary precision.
...

Can you post your mpc library bindings? They would be a useful
starting point for generating bindings for other Forth systems.

Krishna

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.forth

csiph-web