Groups | Search | Server Info | Login | Register

Groups > comp.lang.forth > #23459

Re: FSL contribution available for review

Show all headers | View raw

C G Montgomery <cgm@physics.utoledo.edu> 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

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

Thread

FSL contribution available for review C G Montgomery <cgm@physics.utoledo.edu> - 2013-06-09 19:43 -0400
  Re: FSL contribution available for review mhx@iae.nl (Marcel Hendrix) - 2013-06-11 23:29 +0200
    Re: FSL contribution available for review C G Montgomery <cgm@physics.utoledo.edu> - 2013-06-11 19:38 -0400

csiph-web