Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]
Groups > comp.lang.forth > #9369
| From | mhx@iae.nl (Marcel Hendrix) |
|---|---|
| Subject | Re: ZGAMMA revisited |
| Newsgroups | comp.lang.forth |
| Message-ID | <01921399018435@frunobulax.edu> (permalink) |
| Date | 2012-02-04 22:11 +0200 |
| References | <97140304028435@frunobulax.edu> |
| Organization | SunSITE.dk - Supporting Open source |
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
Back to comp.lang.forth | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll 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