Path: csiph.com!v102.xanadu-bbs.net!xanadu-bbs.net!news.glorb.com!dotsrc.org!filter.dotsrc.org!news.dotsrc.org!not-for-mail X-Newsreader: iForth 2.0 console (October 21, 2006) From: mhx@iae.nl (Marcel Hendrix) Subject: Re: FSL contribution available for review Newsgroups: comp.lang.forth Message-ID: <65061292978434@frunobulax.edu> Date: Tue, 11 Jun 2013 23:29:38 +0200 References: Lines: 64 Organization: SunSITE.dk - Supporting Open source NNTP-Posting-Host: 82.157.147.127 X-Trace: news.sunsite.dk DXC=O3J:=4LSCo1oDTKD0DTdP>YSB=nbEKnk;O5A8C=`f[d1hhA9SOoCJO8e7IB49][_92Y7S096R@MT<^aiVk_ writes: > From Hans Bezemer, a version of the error function which is intermediate, > in precision and speed, between the versions provided in Algorithm #62. > (The algorithm used led to some discussions in comp.lang.forth in October, > 2011.) The submitted code is definitely and unnecessarily weird. I suggest: create _L 7.618009172947146e1 f, -8.650532032941677e1 f, 2.401409824083091e1 f, -1.231739572450155e f, 0.1208650973866179e-2 f, -0.5395239384953e-5 f, : gammaln_new ( f1 -- f2) 1.000000000190015e fover 6 0 do 1e f+ _L i floats + f@ fover f/ frot f+ fswap loop fdrop fover fdup 5.5e f+ fdup fln frot 0.5e f+ f* f- fnegate fswap 2.5066282746310005e f* frot f/ fln f+ ; A test, comparing gammaln_new to the zgamma results: : test ( -- ) 1e 0e FLOCALS| sum arg | #30 0 DO CR ." \ " arg +e. 2 spaces arg 0e R,I->Z zgamma REAL FLN fdup arg gammaln_new f- fswap f/ +e. 1e +TO arg LOOP ; \ x gammaln_new \ 1.0000000000000000000e+0000 -2.6619999999999999994e+0002 \ 2.0000000000000000000e+0000 1.4566666666666666671e+0002 \ 3.0000000000000000000e+0000 1.5156837315413355634e-0016 \ 4.0000000000000000000e+0000 5.8634650642649096784e-0017 \ 5.0000000000000000000e+0000 3.2136599976015110008e-0017 \ 6.0000000000000000000e+0000 2.0019558716152115772e-0017 \ 7.0000000000000000000e+0000 -2.7457496325686340266e-0015 \ 8.0000000000000000000e+0000 -1.0006774080483571122e-0014 \ 9.0000000000000000000e+0000 -2.2116712234461730942e-0014 \ 1.0000000000000000000e+0001 -3.8521827653372058092e-0014 \ 1.1000000000000000000e+0001 -5.8262560404684105464e-0014 \ 1.2000000000000000000e+0001 -8.0298384254259491232e-0014 \ 1.3000000000000000000e+0001 -1.0367055194912764301e-0013 \ 1.4000000000000000000e+0001 -1.2757593891011606686e-0013 \ 1.5000000000000000000e+0001 -1.5138223371142210353e-0013 \ 1.6000000000000000000e+0001 -1.7461517233283822270e-0013 \ 1.7000000000000000000e+0001 -1.9693448721728147118e-0013 \ 1.8000000000000000000e+0001 -2.1811326046325099310e-0013 \ 1.9000000000000000000e+0001 -2.3800667844593437190e-0013 \ 2.0000000000000000000e+0001 -2.5653866524181842926e-0013 \ 2.1000000000000000000e+0001 -2.7367993864864199716e-0013 \ 2.2000000000000000000e+0001 -2.8944059304657091458e-0013 \ 2.3000000000000000000e+0001 -3.0385158276675865048e-0013 \ 2.4000000000000000000e+0001 -3.1696484957374142414e-0013 \ 2.5000000000000000000e+0001 -3.2884041033705138642e-0013 \ 2.6000000000000000000e+0001 -3.3954736559091735848e-0013 \ 2.7000000000000000000e+0001 -3.4915706165453221168e-0013 \ 2.8000000000000000000e+0001 -3.5774134954325300864e-0013 \ 2.9000000000000000000e+0001 -3.6537242328887360372e-0013 \ 3.0000000000000000000e+0001 -3.7211905235738305660e-0013 ok -marcel