ZGAMMA revisited

From	mhx@iae.nl (Marcel Hendrix)
Subject	ZGAMMA revisited
Newsgroups	comp.lang.forth
Message-ID	<77980205028435@frunobulax.edu> (permalink)
Date	2012-01-30 01:05 +0200
Organization	SunSITE.dk - Supporting Open source

Show all headers | View raw

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

Back to comp.lang.forth | Previous | Next — Next in thread | Find similar | Unroll thread

Thread

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

csiph-web