Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #11716 > unrolled thread

Riemann zeta function in Forth

Started byKrishna Myneni <krishna.myneni@ccreweb.org>
First post2012-04-28 07:48 -0700
Last post2012-05-05 02:55 -0700
Articles 18 — 7 participants

Back to article view | Back to comp.lang.forth

Contents

  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

#11716 — Riemann zeta function in Forth

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2012-04-28 07:48 -0700
SubjectRiemann zeta function in Forth
Message-ID<c81633a8-c351-4497-853d-62cb7d0c9f77@f17g2000yqj.googlegroups.com>
It's probably time to implement the Riemann zeta function in Forth.

http://physics.aps.org/synopsis-for/10.1103/PhysRevLett.108.170601

Krishna

[toc] | [next] | [standalone]

#11747

Frommhx@iae.nl
Date2012-04-29 09:18 -0700
Message-ID<23741374.234.1335716291651.JavaMail.geo-discussion-forums@vbez18>
In reply to#11716
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/10.1103/PhysRevLett.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

[toc] | [prev] | [next] | [standalone]

#11754

Frommhx@iae.nl
Date2012-04-29 11:49 -0700
Message-ID<10309024.425.1335725354962.JavaMail.geo-discussion-forums@vbmi19>
In reply to#11747
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/10.1103/PhysRevLett.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

[toc] | [prev] | [next] | [standalone]

#11763

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2012-04-29 17:08 -0700
Message-ID<73d31370-2292-4a41-b1c3-70563ddf8660@e42g2000yqa.googlegroups.com>
In reply to#11754
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/10.1103/PhysRevLett.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

[toc] | [prev] | [next] | [standalone]

#11767

Frommhx@iae.nl
Date2012-04-29 23:47 -0700
Message-ID<27338109.184.1335768456249.JavaMail.geo-discussion-forums@vbqq1>
In reply to#11763
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

[toc] | [prev] | [next] | [standalone]

#11768

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2012-04-30 04:17 -0700
Message-ID<fe90a237-66b8-4450-86d2-a26fb1628aed@2g2000yqk.googlegroups.com>
In reply to#11767
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_of_Riemann_Zeta_Function


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

Krishna

[toc] | [prev] | [next] | [standalone]

#11780

Frommhx@iae.nl
Date2012-04-30 07:25 -0700
Message-ID<23588120.127.1335795911306.JavaMail.geo-discussion-forums@vbfg3>
In reply to#11767
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

[toc] | [prev] | [next] | [standalone]

#11783

From"A. K." <akk@nospam.org>
Date2012-04-30 16:59 +0200
Message-ID<4f9ea897$0$9505$9b4e6d93@newsspool1.arcor-online.net>
In reply to#11780
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....

[toc] | [prev] | [next] | [standalone]

#11787

Frommhx@iae.nl
Date2012-04-30 09:33 -0700
Message-ID<1616278.277.1335803596666.JavaMail.geo-discussion-forums@vbuo17>
In reply to#11783
"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

[toc] | [prev] | [next] | [standalone]

#11793

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2012-04-30 15:21 -0700
Message-ID<b5d28a28-5f7b-4cd4-89e1-5a91faee2a2b@t23g2000yqd.googlegroups.com>
In reply to#11787
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 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.

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

[toc] | [prev] | [next] | [standalone]

#11799

FromAndrew Haley <andrew29@littlepinkcloud.invalid>
Date2012-05-01 02:53 -0500
Message-ID<OoCdnXWyf9QWCwLSnZ2dnUVZ_oGdnZ2d@supernews.com>
In reply to#11793
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.

[toc] | [prev] | [next] | [standalone]

#11804

From"Elizabeth D. Rather" <erather@forth.com>
Date2012-05-01 08:37 -1000
Message-ID<sJCdnf9RCeIfsD3SnZ2dnUVZ_hWdnZ2d@supernews.com>
In reply to#11799
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."
==================================================

[toc] | [prev] | [next] | [standalone]

#11824

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2012-05-02 07:06 -0700
Message-ID<36ee5d45-976c-45ee-8b61-8441ed43d7fa@h12g2000yqi.googlegroups.com>
In reply to#11767
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+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

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/examples/fsl/extras/zzeta.4th

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

Krishna

[toc] | [prev] | [next] | [standalone]

#11859

FromHans Bezemer <the.beez.speaks@gmail.com>
Date2012-05-03 17:43 +0200
Message-ID<4fa2a744$0$6881$e4fe514c@news2.news.xs4all.nl>
In reply to#11824
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/examples/fsl/extras/zzeta.4th
> 
Thanks, already ported to 4tH and committed to SVN.

Hans Bezemer

[toc] | [prev] | [next] | [standalone]

#11901

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2012-05-03 19:33 -0700
Message-ID<2c07c3de-6343-444b-9b04-d481d5ea7b86@t23g2000yqd.googlegroups.com>
In reply to#11859
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/examples/fsl/extras/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

[toc] | [prev] | [next] | [standalone]

#11907

FromThe Beez <the.beez.speaks@gmail.com>
Date2012-05-04 01:39 -0700
Message-ID<3b68e610-406f-4e0a-aeb4-8816861214cf@y15g2000yqh.googlegroups.com>
In reply to#11901
On 4 mei, 04:33, Krishna Myneni <krishna.myn...@ccreweb.org> wrote:
> 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.
Well, I even eliminated the FSL arrays. I do support them, but it
isn't strictly required to make the code run and a quick edit.

Hans

[toc] | [prev] | [next] | [standalone]

#11928

FromKrishna Myneni <krishna.myneni@ccreweb.org>
Date2012-05-04 15:26 -0700
Message-ID<f6cff72e-9713-41a7-b023-e359984736ed@ee2g2000vbb.googlegroups.com>
In reply to#11907
On May 4, 3:39 am, The Beez <the.beez.spe...@gmail.com> wrote:
> On 4 mei, 04:33, Krishna Myneni <krishna.myn...@ccreweb.org> wrote:> 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.
>
> Well, I even eliminated the FSL arrays. I do support them, but it
> isn't strictly required to make the code run and a quick edit.
>
> Hans

Yes, the FSL arrays are not required. I tend to use them for
consistency throughout my code, except in cases where efficiency is a
big issue.

Krishna

[toc] | [prev] | [next] | [standalone]

#11933

FromThe Beez <the.beez.speaks@gmail.com>
Date2012-05-05 02:55 -0700
Message-ID<510f339b-c20f-4f20-8b20-2195ecbb310c@ee2g2000vbb.googlegroups.com>
In reply to#11928
On 5 mei, 00:26, Krishna Myneni <krishna.myn...@ccreweb.org> wrote:
> Yes, the FSL arrays are not required. I tend to use them for
> consistency throughout my code, except in cases where efficiency is a
> big issue.
I know, I did the same thing with your Startrek program ;-)
Thanks again for that one too!

Hans

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.forth

csiph-web