Krishna Myneni writes Re: Riemann zeta function in Forth

> It's probably time to implement the Riemann zeta function in Forth.

> http://physics.aps.org/synopsis-for/...ett.108.170601

See below!

-marcel
-- ---------------------

ANEW -zeta

-- Borwein's algorithm

#20 =: N_valid_digits
N_valid_digits #13 #10 */ =: N
CREATE Zc_ N 1+ FLOATS ALLOT

CREATE Zref_
+NaN F,
+NaN F,
1.6449340668482264364724151666460251892190e F,
1.2020569031595942853997381615114499907651e F,
1.0823232337111381915160036965411679027749e F,
1.0369277551433699263313654864570341680573e F,
1.0173430619844491397145179297909205279018e F,
1.0083492773819228268397975498497967595993e F,
1.0040773561979443393786852385086524652574e F,
1.0020083928260822144178527692324120604826e F,
1.0009945751278180853371459589003190170015e F,
1.0004941886041194645587022825264699364626e F,
1.0002460865533080482986379980477396709530e F,
1.0001227133475784891467518365263573957058e F,
1.0000612481350587048292585451051353337383e F,
1.0000305882363070204935517285106450625779e F,
1.0000152822594086518717325714876367220132e F,
1.0000076371976378997622736002935630292028e F,
1.0000038172932649998398564616446219397200e F,
1.0000019082127165539389256569577951013428e F,
1.0000009539620338727961131520386834493354e F,

: []Zref@ ( ix -- ) ( F: -- r ) Zref_ []FLOAT F@ ;
: []Zc@ ( ix -- ) ( F: -- r ) Zc_ []FLOAT F@ ;
: []Zc! ( ix -- ) ( F: r -- ) Zc_ []FLOAT F! ;

:NONAME ( -- )
1e 1e 1e FLOCALS| x y z |
N 1+
0 DO
x I []Zc!
N I - S>F N I + S>F F* 4e F* y F* FDUP TO y
I 2* 1+ S>F I 2* 2+ S>F F* z F* FDUP TO z
( y z ) F/ +TO x
LOOP ; EXECUTE

: ZETA ( F: x -- z )
FLOCAL x
0e
0 N 1-
?DO
N []Zc@ I []Zc@ F-
I 1+ S>F x F** F/ FSWAP F-
-1 +LOOP
( sum) N []Zc@ F/ x FNEGATE F2^x F2* 1e FSWAP F- F/ ;

: .ZETA ( -- )
N_valid_digits 1+
2 DO CR I 4 U.R
I S>F ZETA FDUP +E.
I []Zref@ FDUP +E.
FSWAP F-/ULP F>S ." (" 0 .R ." ulps)"
LOOP ;

:ABOUT CR ." Try: .ZETA -- Zeta(x) for x positive integer." ;

DOC
(*
FORTH> .zeta ( arg, computed, reference, error in ulps )
2 1.6449340668482264365e+0000 1.6449340668482264365e+0000 (0 ulps)
3 1.2020569031595942854e+0000 1.2020569031595942854e+0000 (0 ulps)
4 1.0823232337111381915e+0000 1.0823232337111381915e+0000 (0 ulps)
5 1.0369277551433699263e+0000 1.0369277551433699263e+0000 (0 ulps)
6 1.0173430619844491397e+0000 1.0173430619844491397e+0000 (0 ulps)
7 1.0083492773819228269e+0000 1.0083492773819228268e+0000 (-1 ulps)
8 1.0040773561979443395e+0000 1.0040773561979443394e+0000 (-1 ulps)
9 1.0020083928260822145e+0000 1.0020083928260822144e+0000 (-1 ulps)
10 1.0009945751278180854e+0000 1.0009945751278180853e+0000 (-1 ulps)
11 1.0004941886041194645e+0000 1.0004941886041194646e+0000 (1 ulps)
12 1.0002460865533080483e+0000 1.0002460865533080483e+0000 (0 ulps)
13 1.0001227133475784891e+0000 1.0001227133475784891e+0000 (0 ulps)
14 1.0000612481350587049e+0000 1.0000612481350587047e+0000 (-1 ulps)
15 1.0000305882363070205e+0000 1.0000305882363070205e+0000 (0 ulps)
16 1.0000152822594086518e+0000 1.0000152822594086519e+0000 (1 ulps)
17 1.0000076371976378998e+0000 1.0000076371976378998e+0000 (0 ulps)
18 1.0000038172932649999e+0000 1.0000038172932649998e+0000 (-1 ulps)
19 1.0000019082127165539e+0000 1.0000019082127165539e+0000 (0 ulps)
20 1.0000009539620338729e+0000 1.0000009539620338727e+0000 (-1 ulps)
*)
ENDDOC