Topic: -={ Challenge }=- Calculate Pi

[highflyingbird]

about Calculating Pi, you can see here

Here is my routine :

Code: [Select]

;;; Calculate Pi                                                       
;;; n---precision,e.g :800                                             
(defun CalPi (n / b c d e f g h r s x digits)
  (setq digits 10000)                            ;long a=10000,a group every four numbers
  (setq c (/ (1+ n) (/ (log 2) (log 10))))  ;the times of iteration
  (setq c (fix c))                     
  (setq e 0)
  (setq r "")                                ;return value
  (setq h (/ digits 5))                         
  (repeat c                                     
    (setq f (cons h f))                          ;f[b++]=a/5
  )
  (repeat (1+ (/ n 4))                          ;repeat 1+ n/4 times
    (setq d 0)                                  ;d=0
    (setq g (+ c c))                            ;g=c*2
    (setq b c) ;b=c
    (setq x nil)
    (while (> b 0)
      (setq d (* d b)) ;d*=b
      (setq b (1- b)) ;--b
      (setq d (+ d (* (car f) digits))) ;d+=f[b]*a
      (setq f (cdr f))
      (setq g (1- g))
      (setq x (cons (rem d g) x)) ;f[b]=d%--g
      (setq d (/ d g)) ;d/=g--
      (setq g (1- g))
    )
    (setq f (reverse x))
    (repeat 13
      (setq f (cdr f))
    )
    (setq s (+ e (/ d digits)))
    (setq s (substr (itoa (+ digits s)) 2)) ;printf("%.4d",e+d/a)
    (setq r (strcat r s))
    (setq e (rem d digits))
    (setq c (- c 13)) ;c-=14,I changed it a little
  )
  (setq r (substr r 2 n))          ;
  (setq r (strcat "3." r))        ;
)
you can test it like:

Code: [Select]

;;;highflybird    2008.4.20 Haikou, lst Edition,2011.5.11 2nd Edition   
(defun C:Pi (/ n Result t1 t2 time file path)
  (initget 7)    ;No zero,No null,No negative
  (setq n (getint "\nPlease Enter The Precision<Less than 16000>:")) ;Input
  (if (> n 16000)                               
    (alert "Your Input is exceed of 16000!")
    (progn 
      (setq t1     (getvar "TDUSRTIMER"))        ;Timer begins
      (setq Result (CalPi n))              ;Calculate Pi
      (setq t2     (getvar "TDUSRTIMER"))        ;End timer.
      (setq time   (* 86400 (- t2 t1)))
      (setq path   (getvar "DWGPREFIX"))
      (if (and (setq file (getfiled "Save to File" path "txt" 1))
       (setq file (Open file "W"))
  )
(progn
  (princ Result file)
  (close file)
)
(princ Result)
      )
      (princ "\nIt takes <second>:")                           
      (princ time)                              ;print time
    )
  )   
  (princ)                                        ;quit program
)
base on this:

Code: [Select]

long a=10000,b,c=2800,d,e,f[2801],g;
main()
{
   for(;b-c;) f[b++]=a/5;
   for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
   for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}

[ElpanovEvgeniy]

my old versions:

use Leibniz formula for Pi

Code: [Select]

(defun Leibniz_pi (/ time)
 ;; ElpanovEvgeniy
 ;; Calculation of Pi
 ;; Last edit: 25.04.2006
 ;; (Leibniz_pi)
 ;; >> Pi = 3.1415926  spent time 72.056 sec.
 ;; use Leibniz formula for Pi
 ;| (+ (- (/ 4 1.) (/ 4 3.))
       (- (/ 4 5.) (/ 4 7.))
       (- (/ 4 9.) (/ 4 11.))
       (- (/ 4 13.) (/ 4 15.))
       (- (/ 4 17.) (/ 4 19.))
       (- (/ 4 21.) (/ 4 23.))
       ...
    )|;
 (defun rec-pi_1 (i0)
  (if (< i i0)
   (+ (- (/ 4. i) (/ 4. (+ 2. i))) (rec-pi_1 (progn (setq i (+ 4. i)) i0)))
   0.
  ) ;_  if
 ) ;_  defun
 (defun rec-pi (i)
  (if (< i 99999990)
   (+ (rec-pi_1 (+ i 10000.)) (rec-pi i))
   0.
  ) ;_  if
 )
 (setq time (car (_VL-TIMES)))
 (princ (strcat "Pi = "
                (rtos (rec-pi 1) 2 7)
                "  spent time "
                (rtos (/ (- (car (_VL-TIMES)) time) 1000.) 2 3)
                " sec."
        )
 )
 (princ)
)
use Monte Carlo method for Pi

Code: [Select]

(defun monte-carlo-pi (n / i)
 ;; ElpanovEvgeniy
 ;; Calculation of Pi
 ;; Last edit: 10.02.2008
 ;; use Monte Carlo method for Pi
 ;|

(repeat 5 (princ (strcat "\n " (rtos (monte-carlo-pi 1000000.) 2 5)))(princ))

>>
 3.14258
 3.14138
 3.14159
 3.14105
 3.14331
 |;
 (defun randnum (/ MODULUS RANDOM)
                ;|
http://intervision.hjem.wanadoo.dk/lisps/randnum.lsp
 Randnum.lsp
 Returns a random number.
 Written by Paul Furman, 1996.
 Based on algorithm by Doug Cooper, 1982.
 (defun randnum (/ modulus multiplier increment random)
  (if (not seed)
   (setq seed (getvar "DATE"))
  ) ;_  if
  (setq modulus    65536
        multiplier 25173
        increment  13849
        seed       (rem (+ (* multiplier seed) increment) modulus)
        random     (/ seed modulus)
  ) ;_  setq
 ) ;_  defun
(setq modulus    4294967296.0
      multiplier 1664525
      increment  1
)
Edit ElpanovEvgeniy
|;
  (/ (setq seed (rem (+ (* 25173.
                           (if seed
                            seed
                            (setq seed (getvar "DATE"))
                           ) ;_  if
                        ) ;_  *
                        13849.
                     ) ;_  +
                     65536
                ) ;_  rem
     ) ;_  setq
     65536
  ) ;_  /
 ) ;_  defun
 (setq i 0.
       y (sqrt 0.5)
       x (- 1. y)
 ) ;_  setq
 (repeat (fix n)
  (if (< (sqrt (+ (expt (+ y (* (randnum) x)) 2.) (expt (* (randnum) y) 2.))) 1.)
   (setq i (1+ i))
  ) ;_  if
 ) ;_  repeat
 (+ (* x y (/ i n) 8.) 2.)
)

[highflyingbird]

ElpanovEvgeniy,Very good introduce for  Monte Carlo method.

Leibniz formula  is too slow.
My code is base on ：
(Sloane's A054387 and A054388). Using Euler's convergence improvement transformation gives
Pi/2= 1+1/3(1+2/5(1+3/7(1+4/9(...(1+k/(2k+1)))))   
      

Code: [Select]

(defun floatPi (k / x i)
  (setq x 2.)
  (setq i (1+ k))
  (repeat k
    (setq i (1- i))
    (setq x (+ 2 (* x (/ i (+ i i 1.)))))
  )
)

[ElpanovEvgeniy]

These programs are not written for the competition - as an educational tool, they have completely fulfilled their tasks! 

[ribarm]

This is an old topic, but I feel I had to revive it...

Can someone confirm me error in Leibniz formula...

I have strange incongruity on 8th decimal - perhaps this is why Evgeniy typed (rtos (...) 2 7), but I've tried with (rtos (...) 2 16) and did get results in witch all decimals are correct except 8th one... Really odd... I am starting to believe that CAD is making miscalculation...

With Euler it's expected that others decimals that follow from 8th are bad as CAD can't calculate so long array of numbers witch are going to infinity without cutting calculation after major array passed...

So can someone confirm my results :

Code - Auto/Visual Lisp: [Select]

1. (defun Leibniz_pi (/ time)
2.  ;; ElpanovEvgeniy
3.  ;; Calculation of Pi
4.  ;; Last edit: 25.04.2006
5.  ;; (Leibniz_pi)
6.  ;; >> Pi = 3.1415926  spent time 72.056 sec.
7.  ;; use Leibniz formula for Pi
8.  ;| (+ (- (/ 4 1.) (/ 4 3.))
9.        (- (/ 4 5.) (/ 4 7.))
10.        (- (/ 4 9.) (/ 4 11.))
11.        (- (/ 4 13.) (/ 4 15.))
12.        (- (/ 4 17.) (/ 4 19.))
13.        (- (/ 4 21.) (/ 4 23.))
14.        ...
15.     )|;
16.  (defun rec-pi_1 (i0)
17.   (if (< i i0)
18.    (+ (- (/ 4. i) (/ 4. (+ 2. i))) (rec-pi_1 (progn (setq i (+ 4. i)) i0)))
19.    0.
20.   ) ;_  if
21.  ) ;_  defun
22.  (defun rec-pi (i)
23.   (if (< i 99999990)
24.    (+ (rec-pi_1 (+ i 10000.)) (rec-pi i))
25.    0.
26.   ) ;_  if
27.  )
28.  (setq time (car (_VL-TIMES)))
29.  (princ (strcat "Pi = "
30.                 (rtos (rec-pi 1) 2 16)
31.                 "  spent time "
32.                 (rtos (/ (- (car (_VL-TIMES)) time) 1000.) 2 3)
33.                 " sec."
34.         )
35.  )
36.  (princ)
37. )
38.  
39. ;; Command: (Leibniz_pi)
40. ;; Pi = 3.141592633589793  spent time 39.874 sec.
41.  
42. ;; Command: (rtos pi 2 16)
43. ;; "3.141592653589793"
44.  
45. ;; !!! ERROR IN Leibniz on 8th decimal !!! ;;
46. ;; !!! Unexpected incongruity !!! ;;
47.  
48. ;;; Euler ;;;
49. ;;; PI=(sqrt (* 6.0 (+ 1.0 (/ 1.0 (expt 2.0 2)) (/ 1.0 (expt 3.0 2))... (/ 1.0 (expt r 2))))) ;;;
50.  
51. (defun Euler_pi (/ time)
52.   (defun sum_1 (i0)
53.     (if (< i i0)
54.       (+ (/ 1.0 (expt (float i) 2)) (sum_1 (progn (setq i (1+ i)) i0)))
55.       0.
56.     )
57.   )
58.   (defun sum (i)
59.     (if (< i 89999990)
60.       (+ (sum_1 (+ i 10000.)) (sum i))
61.       0.
62.     )
63.   )
64.   (setq time (car (_VL-TIMES)))
65.   (princ (strcat "Pi = "
66.                  (rtos (sqrt (* 6.0 (sum 1))) 2 16)
67.                  "  spent time "
68.                  (rtos (/ (- (car (_VL-TIMES)) time) 1000.) 2 3)
69.                  " sec."
70.          )
71.   )
72.   (princ)
73. )
74.  
75. ;; Command: (Euler_pi)
76. ;; Pi = 3.141592642979463  spent time 146.110 sec.
77.  
78. ;; Command: (rtos pi 2 16)
79. ;; "3.141592653589793"
80.  
81. ;; !!! ERROR IN Euler from 8th decimal !!! ;;
82. ;; !!! Expected incongruity !!! ;;
83.  
84. (defun floatPi (k / x i)
85.   (setq x 2.)
86.   (setq i (1+ k))
87.   (repeat k
88.     (setq i (1- i))
89.     (setq x (+ 2 (* x (/ i (+ i i 1.)))))
90.   )
91. )
92.  
93. ;; Command: (rtos (floatPi 999999) 2 16)
94. ;; "3.141592653589793"
95.  
96. ;; Command: (rtos pi 2 16)
97. ;; "3.141592653589793"
98.  
99. ;; Here everything is OK ;;
100.  

M.R.

[ribarm]

Weird, now I changed this :

Code: [Select]

(defun rec-pi (i)
 (if (< i 99999990)
  (+ (rec-pi_1 (+ i 10000.)) (rec-pi i))
  0.
 ) ;_  if
)

To this :

Code: [Select]

(defun rec-pi (i)
 (if (< i 9999999)
  (+ (rec-pi_1 (+ i 10000.)) (rec-pi i))
  0.
 ) ;_  if
)

And now I get wrong 7th digit, but also rest is OK...
Really weird...

Command: (Leibniz_pi)
Pi = 3.141592453589793  spent time 15.187 sec.

Command: (rtos pi 2 16)
"3.141592653589793"

Note - this is my netbook's time - it should be much faster if on my previous testing PC.

[WILL HATCH]

check this out:
http://www.lispworks.com/documentation/lcl50/aug/aug-170.html

If an implementation of Common Lisp uses only one internal floating-point format, that format is designated as single-float by convention

If the implementation here is a single precision floating point value then you're going to have some limitations from there, the biggest mantissa which can be represented in a single is 1.99999988079071 and the smallest increment is 0.00000011920928955078125 which would lead to the problems you've noticed. But I do believe AutoLisp uses a double precision float which would not explain this...

also check out the section on unusual behavior:
http://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80

I did some playing and ran into a similar problem, the following will not converge any further. the sacrifice for speed...

Code - Auto/Visual Lisp: [Select]

1. ;; Will Hatch
2. ;; Calculate pi using Nilakantha method
3. ;; Converges after 4532 iterations
4. ;; Accurate to 12 decimal places
5. ;; 3.141592653589525
6. ;; 2013-08-30
7. (defun Nilakantha_pi (its / i pi time)
8.   (setq time (car (_VL-TIMES)))
9.   (setq i 2)
10.   (setq pi 3.0)
11.   (repeat its
12.     (setq pi (+ pi (-  
13.       (/ 4. (* (+ i 0.) (+ i 1.) (+ i 2.)))
14.       (/ 4. (* (+ i 2.) (+ i 3.) (+ i 4.)))
15.     )))
16.     (setq i (+ i 4))
17.   )
18.  (princ (strcat "Pi = "
19.                 (rtos pi 2 16)
20.                 "  spent time "
21.                 (rtos (/ (- (car (_VL-TIMES)) time) 1000.) 2 3)
22.                 " sec."
23.         )
24.  )
25.  (princ)
26. )