Version 31 of pi

Updated 2012-11-26 15:38:27 by pooryorick

See Also

Discussion

Abstractly, pi is the ratio between the circumference and the diameter of a circle. It can be proven that the number is irrational, so you'll never get an exact representation of its value in Tcl.

The best way to use it in Tcl is to create a global variable (called pi of course) which contains a reasonably accurate approximation and then use that variable throughout your program.

DKF: If you don't want to use a global var, a proc can also be used:

proc Pi {} {return 3.1415926535897931} ;# RS

atan provides a handy way to ask Tcl for the value of pi:

% expr {atan(1) * 4}
3.1415926535897931

MGS Actually, using acos() is (slightly) more efficient:

% set tcl_precision 17
17
% expr {acos(-1)}
3.1415926535897931

Does anyone have any data on which method is preferable from a numerical point of view?

IDG Both contain the assumption that the transcendental functions are accurate to the last ulp. In many math libraries this is not so. I think you are safer with a string representation:

set pi 3.1415926535897931

AM: This is certainly recommendable: that way you can be sure to get a value which is platform-independent or at least highly predictable.

GWM: Evaluating the algorithm below or looking up on the web (direct link: http://www.joyofpi.com/pi.html ) shows the correct value of pi (to the precision given above) ends in the digit 2 (not 1). Pi = 3.141592653589793238462643383. The value given is however good enough to calculate the circumference of a sphere the mean radius of the earth to the size of an atom.


GS: (030927) Here is a small program ables to compute 2400 digits of pi:

# pi-2400.tcl 
# 2400 digits of pi with a spigot algorithm
    
set e 0
for {set b 0} {$b <= 8400} {incr b} {set f($b) 2000}
for {set c 8400} {$c > 0} {incr c -14} {
    set d 0
    for {set b $c} {$b > 0} {incr b -1} {
        set g [expr 2*$b -1]
        set d [expr ($d*$b) + ($f($b)*10000)]
        set f($b) [expr round([expr fmod($d,$g)])]
        set d [expr $d/$g]
    }
    puts -nonewline [format "%.4i" [expr $e+($d/10000)]]
    flush stdout
    set e [expr round([expr fmod($d,10000)])]
}  

It uses a spigot algorithm. More details in A spigot algorithm for the digits of pi, Stanley Rabinowitz and Stan Wagon, American Mathematical Monthly, March 1995, pp195-203.

MGS [2003/09/27]: Here's a more efficient version:

set e 0
for {set b 0} {$b <= 8400} {incr b} {set f($b) 2000}
for {set c 8400} {$c > 0} {incr c -14} {
    set d 0
    for {set b $c} {$b > 0} {incr b -1} {
        set g [expr {2 * $b - 1}]
        set d [expr {($d*$b) + ($f($b)*10000)}]
        set f($b) [expr {round(fmod($d,$g))}]
        set d [expr {$d / $g}]
    }
    puts -nonewline [format "%.4i" [expr {$e+($d/10000)}]]
    flush stdout
    set e [expr {round(fmod($d,10000))}]
}

AM: The efficiency is gained by putting the expressions in braces - then [expr] can parse it once and store the result.

KBK: The following code is several times faster still. The improvements that gain further performance are:

  • encapsulating code in a proc, to take advantage of the local variable table.
  • using integer arithmetic rather than floating point throughout.
  • using a list rather than a hash table to store the f array.

I'm sure that someone can squeeze another factor of 2 or so out of it. :^)

proc pi3 {} {
    set n 2400
    set e 0
    set f {}
    for { set b 0 } { $b <= 8400 } { incr b } {
        lappend f 2000
    }
    for { set c 8400 } { $c > 0 } { incr c -14 } {
        set d 0
        set g [expr { $c * 2 }]
        set b $c
        while 1 {
            incr d [expr { [lindex $f $b] * 10000 }]
            lset f $b [expr {$d % [incr g -1]}]
            set d [expr { $d / $g }]
            incr g -1
            if { [incr b -1] == 0 } break
            set d [expr { $d * $b }]
        }
        puts -nonewline [format %04d [expr { $e + $d / 10000 }]]
        flush stdout
        set e [expr { $d % 10000 }]
    }
    puts {}
}

puts [time pi3]

It might be amusing to implement the Bailey-Borwein-Plouffe algorithm, a digit-extraction algorithm for exactly calculating the digits of pi (http://mathworld.wolfram.com/BBP-TypeFormula.html , http://www.cecm.sfu.ca/~pborwein/ ). Fortran source code is available at the latter site under "My Papers on Pi" (direct link: http://www.cecm.sfu.ca/~pborwein/PISTUFF/Apistuff.html ). Note that the algorithm provides hexadecimal digits, so conversion to decimal would have to be done explicitly if desired.

Albert Davidson Chou HotFusionMan at Yahoo.com

Lars H: Note however that conversion hexadecimal to decimal requires knowing all digits more significant than the one you want to convert, so knowing some hexadecimal digits around e.g magnitude 10^-10000 (i.e., around the 8305th hexadecimal) is not much help for computing the 10000th decimal of pi. Though the algorithm is still a very nice piece of mathematics, this renders it not very useful as a party trick.


KPV: see also Buffon's Needle for one of the earliest Monte Carlo simulation to compute the value of pi.


TV: For the handy, dandy, quick-result requiring, scientifically oriented etc, the tcl/tk front-ended mathematics packageMaxima has a arbitrary precision possibility built in, which quickly gave me this result, which I didn't verify, but it seems handy (and quick):

 (C2)  block([FPPREC:30],bfloat(%pi));
 (D2)                                                3.14159265358979323846264338328B0

comments on the accuracy ?


Sarnold: With math::bigfloat, you can ask (it is in tcllib):

package require math::bigfloat
math::bigfloat::tostr [math::bigfloat::pi 100]
=> 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

The formuma used is :

Pi/4 = 44.atan(1/57) + 7.atan(1/239) - 12.atan(1/682) + 24.atan(1/12943)

with the development of :

atan(1/N) = 1/N - 1/N^3 + 1/N^5 - 1/N^7 ...

This is a formula of Strömer (http://fr.wikipedia.org/wiki/Pi ). Of the few formulaes I investigated, it was the fastest asymptotically. Any how, in tcllib there is a math::constants package which you might want to use and to improve.


Jaf: I kind of like this one (we owe it to Ramanujan I think):

  • Take the three first odd numbers 1,3,5
  • Write each of them twice 113355
  • Put a division sign in the middle: 113/355
  • the inverse 355/113 is pi - 3e-7

Check:

 expr acos(-1)
 => 3.141592653589793
 expr 355/113.0
 => 3.1415929203539825
 expr 355/113.0-acos(-1)
 => 2.667641894049666e-7

"pi" is also the short name of a user on this wiki, pi31415.