Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]
Groups > comp.lang.forth > #11767
| From | mhx@iae.nl |
|---|---|
| Newsgroups | comp.lang.forth |
| Subject | Re: Riemann zeta function in Forth |
| Date | 2012-04-29 23:47 -0700 |
| Organization | http://groups.google.com |
| Message-ID | <27338109.184.1335768456249.JavaMail.geo-discussion-forums@vbqq1> (permalink) |
| References | <c81633a8-c351-4497-853d-62cb7d0c9f77@f17g2000yqj.googlegroups.com> <23741374.234.1335716291651.JavaMail.geo-discussion-forums@vbez18> <10309024.425.1335725354962.JavaMail.geo-discussion-forums@vbmi19> <73d31370-2292-4a41-b1c3-70563ddf8660@e42g2000yqa.googlegroups.com> |
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
Back to comp.lang.forth | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread
Riemann zeta function in Forth Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-04-28 07:48 -0700
Re: Riemann zeta function in Forth mhx@iae.nl - 2012-04-29 09:18 -0700
Re: Riemann zeta function in Forth mhx@iae.nl - 2012-04-29 11:49 -0700
Re: Riemann zeta function in Forth Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-04-29 17:08 -0700
Re: Riemann zeta function in Forth mhx@iae.nl - 2012-04-29 23:47 -0700
Re: Riemann zeta function in Forth Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-04-30 04:17 -0700
Re: Riemann zeta function in Forth mhx@iae.nl - 2012-04-30 07:25 -0700
Re: Riemann zeta function in Forth "A. K." <akk@nospam.org> - 2012-04-30 16:59 +0200
Re: Riemann zeta function in Forth mhx@iae.nl - 2012-04-30 09:33 -0700
Re: Riemann zeta function in Forth Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-04-30 15:21 -0700
Re: Riemann zeta function in Forth Andrew Haley <andrew29@littlepinkcloud.invalid> - 2012-05-01 02:53 -0500
Re: Riemann zeta function in Forth "Elizabeth D. Rather" <erather@forth.com> - 2012-05-01 08:37 -1000
Re: Riemann zeta function in Forth Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-05-02 07:06 -0700
Re: Riemann zeta function in Forth Hans Bezemer <the.beez.speaks@gmail.com> - 2012-05-03 17:43 +0200
Re: Riemann zeta function in Forth Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-05-03 19:33 -0700
Re: Riemann zeta function in Forth The Beez <the.beez.speaks@gmail.com> - 2012-05-04 01:39 -0700
Re: Riemann zeta function in Forth Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-05-04 15:26 -0700
Re: Riemann zeta function in Forth The Beez <the.beez.speaks@gmail.com> - 2012-05-05 02:55 -0700
csiph-web