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

Krishna

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

On Sunday, April 29, 2012 6:18:11 PM UTC+2, m...@iae.nl wrote:
> 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
> -- ---------------------

Of course, it works for non-integers and complex floats too:

: ZZETA ( F: z1 -- z2 )
ZLOCAL x
0+0i
0 N 1-
?DO
N []Zc@ I []Zc@ F- 0e R,I->Z
I 1+ S>F 0e R,I->Z x Z** Z/ ZSWAP Z-
-1 +LOOP
( sum) N []Zc@ Z/F
2+0i x ZNEGATE Z** Z2* 1+0i ZSWAP Z- Z/ ;

2.5e 1.25e zzeta 22 +ze.r (1.0766731004736795343e+0000 -2.5667157846102185318e-0001)

[..]

-marcel

On Apr 29, 1:49*pm, m...@iae.nl wrote:
> On Sunday, April 29, 2012 6:18:11 PM UTC+2, m...@iae.nl wrote:
> > 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
> > -- ---------------------
>
> Of course, it works for non-integers and complex floats too:
>
> : ZZETA ( F: z1 -- z2 )
> * * * * ZLOCAL x
> * * * * 0+0i
> * * * * 0 N 1-
> * * * * * * *?DO
> * * * * * * * * N []Zc@ *I []Zc@ F- *0e R,I->Z
> * * * * * * * * I 1+ S>F 0e R,I->Z x Z** *Z/ *ZSWAP Z-
> * * * * -1 +LOOP
> * * * * ( sum) N []Zc@ Z/F
> * * * * 2+0i x ZNEGATE Z** *Z2* *1+0i ZSWAP Z- *Z/ ;
>
> 2.5e 1.25e zzeta 22 +ze.r (1.0766731004736795343e+0000 -2.5667157846102185318e-0001)
>
> [..]
>
> -marcel

Nice. As with the complex gamma function, I expect trouble with
accuracy near the poles. I'll have to poke around for some reference
values.

Krishna

On Monday, April 30, 2012 2:08:27 AM UTC+2, Krishna Myneni wrote:
>> [..]
> Nice. As with the complex gamma function, I expect trouble with
> accuracy near the poles. I'll have to poke around for some reference
> values.

Here is the (complex) Riemann Zeta function for arbitrary precision, using iForth's MPC and MPFR extensions. As far as I know, the only pole is at s = 1+0i. I did not implement the results for s a negative integer (which involves Bernouilli numbers). Mathematica gives results for arbitrary s.

-marcel

-- -----------------
(*
* LANGUAGE : ANS Forth with extensions
* PROJECT : Forth Environments
* DESCRIPTION : The Riemann Zeta function
* CATEGORY : Special functions
* AUTHOR : Marcel Hendrix
* LAST CHANGE : April 29, 2012, Marcel Hendrix
*)



NEEDS -miscutil
NEEDS -mpfr
NEEDS -mpc

REVISION -zeta "--- Riemann Zeta fnc Version 1.00 ---"

PRIVATES

DOC
(*
In Canadian Mathematical Society, Conference Proceedings, 1991.
'An Efficient Algorithm for the Riemann Zeta Function,' P. Borwein

Borwein's algorithm 2:

k
.--- (n+i-1)! * 4^i
d_k = n * > ---------------
'--- (n-i)! * (2*i)!
i=0

n-1
-1 .---- (-1)^k (d_k - d_n)
Zeta(s) = --------------- * > ------------------
d_n*(1-2^(1-s)) '---- (k+1)^s
k=0

The algorithm only works for Re(s) > 1. A special case is when
Im(s) = 0.5, but it relates to efficiency, not to accuracy.
*)
ENDDOC

#60 =: N_valid_digits PRIVATE
N_valid_digits #14 #10 */ =: N PRIVATE
N 1+ F#ARRAY Zc_ PRIVATE

: []Zc@ ( ix -- ) ( F: -- r ) Zc_ [@]F# ; PRIVATE
: []Zc! ( ix -- ) ( F: r -- ) Zc_ [!]F# ; PRIVATE

F#0.0 F#VALUE x PRIVATE
F#0.0 F#VALUE y PRIVATE
F#0.0 F#VALUE z PRIVATE
0+0i Z#VALUE z#x PRIVATE

:NONAME ( -- )
F#1.0 TO x
F#1.0 TO y
F#1.0 TO z
N 1+
0 DO
x I []Zc!
N I - N->F# N I + N->F# F#* F#4.0 F#* y F#* F#DUP TO y
I 2* 1+ N->F# I 2* 2+ N->F# F#* z F#* F#DUP TO z
( y z ) F#/ +TO x
LOOP ; EXECUTE

-- Valid for x > 1
: ZETA ( F: x -- z )
TO x
x F#0= IF F#-0.5 EXIT ENDIF
x F#1.0 F#<= IF F#NaN EXIT ENDIF
F#0.0
0 N 1-
?DO
N []Zc@ I []Zc@ F#-
I 1+ N->F# x F#** F#/ F#SWAP F#-
-1 +LOOP
( sum) N []Zc@ F#/ x F#NEGATE F#2^x F#2* F#1.0 F#SWAP F#- F#/ ;

: ZZETA ( F: z1 -- z2 )
TO z#x
z#x 0+0i Z#= IF -0.5+0i EXIT ENDIF
z#x Z#REAL F#1.0 F#<= IF NaN+NaNi EXIT ENDIF
0+0i
0 N 1-
?DO
N []Zc@ I []Zc@ F#- F#0.0 R,I->Z#
I 1+ N->F# F#0.0 R,I->Z# z#x Z#** Z#/ Z#SWAP Z#-
-1 +LOOP
( sum) N []Zc@ Z#/F
2+0i z#x Z#NEGATE Z#** Z#2* 1+0i Z#SWAP Z#- Z#/ ;

NESTING @ 1 =
[IF]

#21 F#ARRAY Zref_
Zref_ #21 1 }}F#READ
-0.5
+Inf
1.6449340668482264364724151666460251892190
1.2020569031595942853997381615114499907651
1.0823232337111381915160036965411679027749
1.0369277551433699263313654864570341680573
1.0173430619844491397145179297909205279018
1.0083492773819228268397975498497967595993
1.0040773561979443393786852385086524652574
1.0020083928260822144178527692324120604826
1.0009945751278180853371459589003190170015
1.0004941886041194645587022825264699364626
1.0002460865533080482986379980477396709530
1.0001227133475784891467518365263573957058
1.0000612481350587048292585451051353337383
1.0000305882363070204935517285106450625779
1.0000152822594086518717325714876367220132
1.0000076371976378997622736002935630292028
1.0000038172932649998398564616446219397200
1.0000019082127165539389256569577951013428
1.0000009539620338727961131520386834493354

: []Zref@ ( ix -- ) ( F: -- r ) Zref_ [@]F# ;

: .ZETA ( -- )
#DECIMALS #40 TO #DECIMALS
#21 2 DO CR I 4 U.R
I N->F# F#0.0 R,I->Z# ZZETA Z#REAL F#.
I []Zref@ F#.
LOOP
TO #DECIMALS ;

DOC
(*
FORTH> .zeta ( More accurate than UBASIC )
2 1.6449340668482264364724151666460251892189e+0000 1.6449340668482264364724151666460251892190e+0000
3 1.2020569031595942853997381615114499907649e+0000 1.2020569031595942853997381615114499907651e+0000
4 1.0823232337111381915160036965411679027747e+0000 1.0823232337111381915160036965411679027748e+0000
5 1.0369277551433699263313654864570341680570e+0000 1.0369277551433699263313654864570341680573e+0000
6 1.0173430619844491397145179297909205279018e+0000 1.0173430619844491397145179297909205279017e+0000
7 1.0083492773819228268397975498497967595998e+0000 1.0083492773819228268397975498497967595993e+0000
8 1.0040773561979443393786852385086524652589e+0000 1.0040773561979443393786852385086524652574e+0000
9 1.0020083928260822144178527692324120604856e+0000 1.0020083928260822144178527692324120604825e+0000
10 1.0009945751278180853371459589003190170060e+0000 1.0009945751278180853371459589003190170015e+0000
11 1.0004941886041194645587022825264699364686e+0000 1.0004941886041194645587022825264699364626e+0000
12 1.0002460865533080482986379980477396709604e+0000 1.0002460865533080482986379980477396709529e+0000
13 1.0001227133475784891467518365263573957142e+0000 1.0001227133475784891467518365263573957057e+0000
14 1.0000612481350587048292585451051353337474e+0000 1.0000612481350587048292585451051353337382e+0000
15 1.0000305882363070204935517285106450625876e+0000 1.0000305882363070204935517285106450625779e+0000
16 1.0000152822594086518717325714876367220232e+0000 1.0000152822594086518717325714876367220132e+0000
17 1.0000076371976378997622736002935630292130e+0000 1.0000076371976378997622736002935630292028e+0000
18 1.0000038172932649998398564616446219397304e+0000 1.0000038172932649998398564616446219397200e+0000
19 1.0000019082127165539389256569577951013532e+0000 1.0000019082127165539389256569577951013427e+0000
20 1.0000009539620338727961131520386834493459e+0000 1.0000009539620338727961131520386834493354e+0000 ok
*)
ENDDOC

:ABOUT CR ." Try: .ZETA -- Test real Zeta(x) for x > 1." ;

[ELSE]

:ABOUT CR ." Try: ( F#: x1 -- x2 ) ZETA -- Zeta(x) for x > 1."
CR ." ( Z#: z1 -- z2 ) ZZETA -- Zeta(z1) " ;

[THEN]

NESTING @ 1 = [IF] .ABOUT -zeta CR [THEN]
DEPRIVE

(* End of Source *)


-marcel

On Apr 30, 1:47*am, m...@iae.nl wrote:
> On Monday, April 30, 2012 2:08:27 AM UTC+2, Krishna Myneni wrote:
> >> [..]
> > Nice. As with the complex gamma function, I expect trouble with
> > accuracy near the poles. I'll have to poke around for some reference
> > values.
>
> Here is the (complex) Riemann Zeta function for arbitrary precision, using iForth's MPC and MPFR extensions. As far as I know, the only pole is at s= 1+0i. I did not implement the results for s a negative integer (which involves Bernouilli numbers). ..

Looks like that's correct,

http://www.proofwiki.org/wiki/Poles_..._Zeta_Function


I guess it must have local minima and maxima (peaks and valleys) in
the complex plane.

Krishna

On Monday, April 30, 2012 8:47:36 AM UTC+2, m...@iae.nl wrote:
[..]
Please remove the lines
z#x 0+0i Z#= IF -0.5+0i EXIT ENDIF
z#x Z#REAL F#1.0 F#<= IF NaN+NaNi EXIT ENDIF
from ZZETA (and ZETA) -- they are unnecessary. The continuation
works everywhere but the accuracy decreases outside of the
recommended domain. I have now implemented P.Borwein's Algorithm 3,
but it is very much slower as it needs a few hundred terms for 60
significant digits.

ZZETA matches Mathematica for the negative real axis and some
selected poles with Re(s) < 1.

-marcel

On 30.04.2012 16:25, mhx@iae.nl wrote:
> On Monday, April 30, 2012 8:47:36 AM UTC+2, m...@iae.nl wrote:
> [..]
> Please remove the lines
> z#x 0+0i Z#= IF -0.5+0i EXIT ENDIF
> z#x Z#REAL F#1.0 F#<= IF NaN+NaNi EXIT ENDIF
> from ZZETA (and ZETA) -- they are unnecessary. The continuation
> works everywhere but the accuracy decreases outside of the
> recommended domain. I have now implemented P.Borwein's Algorithm 3,
> but it is very much slower as it needs a few hundred terms for 60
> significant digits.
>
> ZZETA matches Mathematica for the negative real axis and some
> selected poles with Re(s)< 1.
>
> -marcel
>

Under the somewhat simplifying assumption that Forth is good for small
devices, but not altoo well suited for mathematics (there are better
competitors in that realm), what could one practically achieve with
numbercrunching the Zetafunction in Forth??? Crypto-stuff???

Of course, I understand that it can be a lot of fun....

"A. K." <akk@..rg> writes Re: Riemann zeta function in Forth

> On 30.04.2012 16:25, mh...nl wrote:
[..]
>> ZZETA matches Mathematica for the negative real axis and some
>> selected poles with Re(s)< 1.

> Under the somewhat simplifying assumption that Forth is good for small
> devices, but not altoo well suited for mathematics (there are better
> competitors in that realm),

Well, I do not agree that Forth is good only for small devices.
For instance, I have Forth software that hooks both into SPICE
simulator engines (LTSpice and NGSPICE) and Simulink, and allows
digital state-feedback control of switched mode power supplies
that are being designed at the device-level. I think I am on
the point to have gained sufficient knowledge to need neither
SPICE nor Simulink to do what I want (I can do it all in Forth).

On the other hand, I am installing Code Composer Studio right now
to work with the Piccolo controlSTICK DSP (40 Euros, can do FP and
fixed-point). I have no Forth software for that device (yet :-)

For each task imaginable there exists a better tool than I am
currently using. However, I'd like to know more or less exactly HOW
such tools work, as I observe they break rather easily (The TI tools
took 10 minutes diskaccess, 1 GB memory and 20 minutes internet contact
to install), leaving the user with no other option than try something
else. I have found that Forth is quite good at helping me find out how
things work. This is becoming more and more important now less and less
people seem to have such (saleable) knowledge anymore.

> what could one practically achieve with
> numbercrunching the Zetafunction in Forth??? Crypto-stuff???

I don't know (isn't it nice to gain some understanding of the zeta
function while using Forth?). For me it is a way to test the MPFR and MPC interfaces to iForth. Indeed, I improved the functionality and fixed a
few bugs during the Riemann Zeta coding. And I learned quite a few things
about approximating functions and complex numbers that would be quite an
unpleasant and prohibitively time-expensive study in other settings.

> Of course, I understand that it can be a lot of fun....

The driving factor for most things I do.

-marcel

On Apr 30, 11:33*am, m...@iae.nl wrote:
> "A. K." <akk@..rg> writes Re: Riemann zeta function in Forth
>
>
>
> > On 30.04.2012 16:25, mh...nl wrote:
> [..]
> >> ZZETA matches Mathematica for the negative real axis and some
> >> selected poles with Re(s)< *1.
> > Under the somewhat simplifying assumption that Forth is good for small
> > devices, but not altoo well suited for mathematics (there are better
> > competitors in that realm),
>
> Well, I do not agree that Forth is good only for small devices.

Ditto. I didn't know anyone actually claimed this anymore, with modern
Forth systems. AFAIK Gforth doesn't advertise that one should use it
for small apps on small systems.

>

> For each task imaginable there exists a better tool than I am
> currently using. However, I'd like to know more or less exactly HOW
> such tools work, as I observe they break rather easily ...

My observations also. Not only do they break easily, they are also too
cumbersome to use, and take too much time to learn to use effectively,
even if they can calculate more efficiently. And, who can afford the
time to learn exactly the most efficient tool for a given task, unless
it offers an absolutely compelling benefit. Forth is very handy for
exploratory computing, much better than Fortran or BASIC for such a
purpose.

> took 10 minutes diskaccess, 1 GB memory and 20 minutes internet contact
> to install), leaving the user with no other option than try something
> else. I have found that Forth is quite good at helping me find out how
> things work. This is becoming more and more important now less and less
> people seem to have such (saleable) knowledge anymore.

Yes. More and more people know how to use the same canned software,
and don't know better than to question nonsensical results when they
occur.

>
> > * * * * * * * * * * * * * * what could one practically achieve with
> > numbercrunching the Zetafunction in Forth??? Crypto-stuff???
>
> I don't know (isn't it nice to gain some understanding of the zeta
> function while using Forth?). For me it is a way to test the MPFR and MPCinterfaces to iForth. Indeed, I improved the functionality and fixed a
> few bugs during the Riemann Zeta coding. And I learned quite a few things
> about approximating functions and complex numbers that would be quite an
> unpleasant and prohibitively time-expensive study in other settings.

I may want to do a freezing transition in glass simulation someday,
look at the density of prime numbers along the number line, or perform
some other calculation where the Riemann zeta function comes into
play... Having a ready Forth module allows me to integrate it with my
existing code. The existing library code for scientific computation in
Forth is adequate to handle a significant portion of my needs. People
still use Fortran for calculations because of the large set of
available libraries, so providing libraries in different languages
just makes sense.
>
> > Of course, I understand that it can be a lot of fun....
>
> The driving factor for most things I do.
>

Yes, and it can be practical too in terms of time and money.

Krishna

Krishna Myneni <krishna.myneni@ccreweb.org> wrote:
> On Apr 30, 11:33?am, m...@iae.nl wrote:
>> "A. K." <akk@..rg> writes Re: Riemann zeta function in Forth
>>
>> > On 30.04.2012 16:25, mh...nl wrote:
>> [..]
>> >> ZZETA matches Mathematica for the negative real axis and some
>> >> selected poles with Re(s)< ?1.
>> > Under the somewhat simplifying assumption that Forth is good for small
>> > devices, but not altoo well suited for mathematics (there are better
>> > competitors in that realm),
>>
>> Well, I do not agree that Forth is good only for small devices.
>
> Ditto. I didn't know anyone actually claimed this anymore, with
> modern Forth systems. AFAIK Gforth doesn't advertise that one should
> use it for small apps on small systems.

I'm not really sure what constitutes a small device: even some tiny
systems have quite powerful processors these days. I think Forth
still has some virtue, even on these systems.

Andrew.

On 4/30/12 9:53 PM, Andrew Haley wrote:
> Krishna Myneni<krishna.myneni@ccreweb.org> wrote:
>> On Apr 30, 11:33?am, m...@iae.nl wrote:
>>> "A. K."<akk@..rg> writes Re: Riemann zeta function in Forth
>>>
>>>> On 30.04.2012 16:25, mh...nl wrote:
>>> [..]
>>>>> ZZETA matches Mathematica for the negative real axis and some
>>>>> selected poles with Re(s)< ?1.
>>>> Under the somewhat simplifying assumption that Forth is good for small
>>>> devices, but not altoo well suited for mathematics (there are better
>>>> competitors in that realm),
>>>
>>> Well, I do not agree that Forth is good only for small devices.
>>
>> Ditto. I didn't know anyone actually claimed this anymore, with
>> modern Forth systems. AFAIK Gforth doesn't advertise that one should
>> use it for small apps on small systems.
>
> I'm not really sure what constitutes a small device: even some tiny
> systems have quite powerful processors these days. I think Forth
> still has some virtue, even on these systems.

Forth is good for many things on systems of all sizes.

I take great exception to the notion that powerful hardware means there
is no need for efficient programming tools or systems that facilitate
the programer's task. The history of technology is that whatever the
size of the platform, users' desires to use it will grow. There will
always be a place for technology that can use platforms more
efficiently, and if it can get the work done quicker in many cases, so
much the better.

Cheers,
Elizabeth

--
==================================================
Elizabeth D. Rather (US & Canada) 800-55-FORTH
FORTH Inc. +1 310.999.6784
5959 West Century Blvd. Suite 700
Los Angeles, CA 90045
http://www.forth.com

"Forth-based products and Services for real-time
applications since 1973."
==================================================

On Apr 30, 1:47*am, m...@iae.nl wrote:
> On Monday, April 30, 2012 2:08:27 AM UTC+2, Krishna Myneni wrote:
> >> [..]
> > Nice. As with the complex gamma function, I expect trouble with
> > accuracy near the poles. I'll have to poke around for some reference
> > values.
>
> Here is the (complex) Riemann Zeta function for arbitrary precision, using iForth's MPC and MPFR extensions. As far as I know, the only pole is at s= 1+0i. I did not implement the results for s a negative integer (which involves Bernouilli numbers). Mathematica gives results for arbitrary s.
>
> -marcel
>
> -- -----------------
> (*
> ** LANGUAGE * *: ANS Forth with extensions
> ** PROJECT * * : Forth Environments
> ** DESCRIPTION : The Riemann Zeta function
> ** CATEGORY * *: Special functions
> ** AUTHOR * * *: Marcel Hendrix
> ** LAST CHANGE : April 29, 2012, Marcel Hendrix
> **)
>
> * * * * NEEDS -miscutil
> * * * * NEEDS -mpfr
> * * * * NEEDS -mpc
>
> * * * * REVISION -zeta "--- Riemann Zeta fnc * *Version 1.00 ---"
>
> * * * * PRIVATES
>
> DOC
> (*
> * In Canadian Mathematical Society, Conference Proceedings, 1991.
> * 'An Efficient Algorithm for the Riemann Zeta Function,' P. Borwein
>
> * Borwein's algorithm 2:
>
> * * * * * * * k
> * * * * * * .--- (n+i-1)! * 4^i
> * d_k = n * *> * ---------------
> * * * * * * '--- (n-i)! * (2*i)!
> * * * * * * *i=0
>
> * * * * * * * * * * * * * * * * n-1
> * * * * * * * * *-1 * * * * * .---- (-1)^k (d_k - d_n)
> * Zeta(s) = --------------- * *> * *------------------
> * * * * * * d_n*(1-2^(1-s)) * '---- * * *(k+1)^s
> * * * * * * * * * * * * * * * * k=0
>
> * The algorithm only works for Re(s) > 1. A special case is when
> * Im(s) = 0.5, but it relates to efficiency, not to accuracy.
> *)
> ENDDOC
>
> #60 =: N_valid_digits * * * * * PRIVATE
> N_valid_digits #14 #10 */ =: N *PRIVATE
> N 1+ F#ARRAY Zc_ * * * * * * * *PRIVATE
>
> : []Zc@ * ( ix -- ) ( F: -- r ) Zc_ [@]F# ; PRIVATE
> : []Zc! * ( ix -- ) ( F: r -- ) Zc_ [!]F# ; PRIVATE
>
> F#0.0 F#VALUE x * PRIVATE
> F#0.0 F#VALUE y * PRIVATE
> F#0.0 F#VALUE z * PRIVATE
> 0+0i *Z#VALUE z#x PRIVATE
>
> :NONAME ( -- )
> * * * * F#1.0 TO x
> * * * * F#1.0 TO y
> * * * * F#1.0 TO z
> * * * * N 1+
> * * * * 0 DO
> * * * * * * * * x I []Zc!
> * * * * * * * * N I - * N->F# *N I + * N->F# F#* *F#4.0 F#* *y F#* *F#DUP TO y
> * * * * * * * * I 2* 1+ N->F# *I 2* 2+ N->F# F#* * * * * * * z F#* *F#DUP TO z
> * * * * * * * * ( y z ) F#/ +TO x
> * * * * LOOP ; EXECUTE
>
> -- Valid for x > 1
> : ZETA ( F: x -- z )
> * * * * TO x
> * * * * x F#0= * * * IF *F#-0.5 EXIT *ENDIF
> * * * * x F#1.0 F#<= IF *F#NaN *EXIT *ENDIF
> * * * * F#0.0
> * * * * 0 N 1-
> * * * * * * *?DO
> * * * * * * * * N []Zc@ *I []Zc@ F#-
> * * * * * * * * I 1+ N->F# x F#** *F#/ *F#SWAP F#-
> * * * * -1 +LOOP
> * * * * ( sum) N []Zc@ F#/ * x F#NEGATE F#2^x F#2* *F#1.0 F#SWAP F#- *F#/ ;
>
> : ZZETA ( F: z1 -- z2 )
> * * * * TO z#x
> * * * * z#x 0+0i * * * * Z#= *IF *-0.5+0i *EXIT *ENDIF
> * * * * z#x Z#REAL F#1.0 F#<= IF *NaN+NaNi EXIT *ENDIF
> * * * * 0+0i
> * * * * 0 N 1-
> * * * * * * *?DO
> * * * * * * * * N []Zc@ *I []Zc@ F#- *F#0.0 R,I->Z#
> * * * * * * * * I 1+ N->F# *F#0.0 R,I->Z# z#x Z#** *Z#/ *Z#SWAP Z#-
> * * * * -1 +LOOP
> * * * * ( sum) N []Zc@ Z#/F
> * * * * 2+0i z#x Z#NEGATE Z#** *Z#2* *1+0i Z#SWAP Z#- *Z#/ ;
>
> NESTING @ 1 =
> * [IF]
>
> * * * * #21 F#ARRAY Zref_
> * * * * Zref_ #21 1 }}F#READ
> * * * * * *-0.5
> * * * * * *+Inf
> * * * * * *1.6449340668482264364724151666460251892190
> * * * * * *1.2020569031595942853997381615114499907651
> * * * * * *1.0823232337111381915160036965411679027749
> * * * * * *1.0369277551433699263313654864570341680573
> * * * * * *1.0173430619844491397145179297909205279018
> * * * * * *1.0083492773819228268397975498497967595993
> * * * * * *1.0040773561979443393786852385086524652574
> * * * * * *1.0020083928260822144178527692324120604826
> * * * * * *1.0009945751278180853371459589003190170015
> * * * * * *1.0004941886041194645587022825264699364626
> * * * * * *1.0002460865533080482986379980477396709530
> * * * * * *1.0001227133475784891467518365263573957058
> * * * * * *1.0000612481350587048292585451051353337383
> * * * * * *1.0000305882363070204935517285106450625779
> * * * * * *1.0000152822594086518717325714876367220132
> * * * * * *1.0000076371976378997622736002935630292028
> * * * * * *1.0000038172932649998398564616446219397200
> * * * * * *1.0000019082127165539389256569577951013428
> * * * * * *1.0000009539620338727961131520386834493354
>
> * * * * : []Zref@ ( ix -- ) ( F: -- r ) Zref_ [@]F# ;
>
> * * * * : .ZETA ( -- )
> * * * * * * * * #DECIMALS #40 TO #DECIMALS
> * * * * * * * * #21 2 DO CR I 4 U.R
> * * * * * * * * * * * * * * I N->F# F#0.0 R,I->Z# ZZETA Z#REAL F#.
> * * * * * * * * * * * * * * I []Zref@ * ** * * * * * * * * * *F#.
> * * * * * * * * * * LOOP
> * * * * * * * * TO #DECIMALS ;
>
> * * * * DOC
> * * * * (*
> * * * * FORTH> .zeta ( More accurate than UBASIC )
> * * * * * *2 1.6449340668482264364724151666460251892189e+00001. 6449340668482264364724151666460251892190e+0000
> * * * * * *3 1.2020569031595942853997381615114499907649e+00001. 2020569031595942853997381615114499907651e+0000
> * * * * * *4 1.0823232337111381915160036965411679027747e+00001. 0823232337111381915160036965411679027748e+0000
> * * * * * *5 1.0369277551433699263313654864570341680570e+00001. 0369277551433699263313654864570341680573e+0000
> * * * * * *6 1.0173430619844491397145179297909205279018e+00001. 0173430619844491397145179297909205279017e+0000
> * * * * * *7 1.0083492773819228268397975498497967595998e+00001. 0083492773819228268397975498497967595993e+0000
> * * * * * *8 1.0040773561979443393786852385086524652589e+00001. 0040773561979443393786852385086524652574e+0000
> * * * * * *9 1.0020083928260822144178527692324120604856e+00001. 0020083928260822144178527692324120604825e+0000
> * * * * * 10 1.0009945751278180853371459589003190170060e+0000 1..0009945751278180853371459589003190170015e+0000
> * * * * * 11 1.0004941886041194645587022825264699364686e+0000 1..0004941886041194645587022825264699364626e+0000
> * * * * * 12 1.0002460865533080482986379980477396709604e+0000 1..0002460865533080482986379980477396709529e+0000
> * * * * * 13 1.0001227133475784891467518365263573957142e+0000 1..0001227133475784891467518365263573957057e+0000
> * * * * * 14 1.0000612481350587048292585451051353337474e+0000 1..0000612481350587048292585451051353337382e+0000
> * * * * * 15 1.0000305882363070204935517285106450625876e+0000 1..0000305882363070204935517285106450625779e+0000
> * * * * * 16 1.0000152822594086518717325714876367220232e+0000 1..0000152822594086518717325714876367220132e+0000
> * * * * * 17 1.0000076371976378997622736002935630292130e+0000 1..0000076371976378997622736002935630292028e+0000
> * * * * * 18 1.0000038172932649998398564616446219397304e+0000 1..0000038172932649998398564616446219397200e+0000
> * * * * * 19 1.0000019082127165539389256569577951013532e+0000 1..0000019082127165539389256569577951013427e+0000
> * * * * * 20 1.0000009539620338727961131520386834493459e+0000 1..0000009539620338727961131520386834493354e+0000 ok
> * * * * *)
> * * * * ENDDOC
>
> * * * * :ABOUT *CR ." Try: .ZETA * -- Test real Zeta(x) for x> 1." ;
>
> [ELSE]
>
> :ABOUT *CR ." Try: ( F#: x1 -- x2 ) ZETA * -- Zeta(x) for x > 1."
> * * * * CR ." * * *( Z#: z1 -- z2 ) ZZETA *-- Zeta(z1) " ;
>
> [THEN]
>
> NESTING @ 1 = [IF] * * *.ABOUT -zeta CR * * * * [THEN]
> * * * * * * * * * * * * DEPRIVE
>
> * * * * * * * * * * * * * * * (* End of Source *)
>
> -marcel

I've ported the above code to a double precision version, compatible
with the Forth Scientific Library (FSL). The double precision version
is accurate to machine precision, about 2e-16. The ported version may
be found at,

ftp://ccreweb.org/software/kforth/ex...tras/zzeta.4th

We need some reference values for ZZETA which are off the real axis.

Krishna

Krishna Myneni wrote:
> I've ported the above code to a double precision version, compatible
> with the Forth Scientific Library (FSL). The double precision version
> is accurate to machine precision, about 2e-16. The ported version may
> be found at,
>
> ftp://ccreweb.org/software/kforth/ex...tras/zzeta.4th
>
Thanks, already ported to 4tH and committed to SVN.

Hans Bezemer

On May 3, 10:43*am, Hans Bezemer <the.beez.spe...@gmail.com> wrote:
> Krishna Myneni wrote:
> > I've ported the above code to a double precision version, compatible
> > with the Forth Scientific Library (FSL). The double precision version
> > is accurate to machine precision, about 2e-16. The ported version may
> > be found at,
>
> >ftp://ccreweb.org/software/kforth/ex...tras/zzeta.4th
>
> Thanks, already ported to 4tH and committed to SVN.
>
> Hans Bezemer

You're welcome, ... and thanks to Marcel for developing and posting
this code. It was relatively easy to make it work under plain ANS
Forth (without the MPC extensions). I think the code should load
without issues on ANS systems, except for maybe some filename fix up
in the include statements.

Krishna