[ create a new paste ] login | about

Link: http://codepad.org/XKfVY5Ua    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Jun 14:
; the digits of pi, again

(define (isqrt n)
  (if (not (and (positive? n) (integer? n)))
      (error 'isqrt "must be positive integer")
      (let loop ((x n))
        (let ((y (quotient (+ x (quotient n x)) 2)))
          (if (< y x) (loop y) x)))))

(define ch-A 13591409)
(define ch-B 545140134)
(define ch-C 640320)
(define ch-C^3 (expt 640320 3))
(define ch-D 12)

(define (ch-split a b)
  (display a) (display " ") (display b) (newline)
  (if (= 1 (- b a))
      (let ((g (* (- (* 6 b) 5) (- (* 2 b) 1) (- (* 6 b) 1))))
        (list g (quotient (* ch-C^3 (expt b 3)) 24)
                (* (expt -1 b) g (+ (* b ch-B) ch-A))))
      (let* ((mid (quotient (+ a b) 2))
             (gpq1 (ch-split a mid)) (gpq2 (ch-split mid b))
             (g1 (car gpq1)) (p1 (cadr gpq1)) (q1 (caddr gpq1))
             (g2 (car gpq2)) (p2 (cadr gpq2)) (q2 (caddr gpq2)))
        (list (* g1 g2) (* p1 p2) (+ (* q1 p2) (* q2 g1))))))

(define (pi digits)
  (let* ((num-terms (inexact->exact (floor (+ 2 (/ digits 14.181647462)))))
         (sqrt-C (isqrt (* ch-C (expt 100 digits)))))
    (let* ((gpq (ch-split 0 num-terms))
           (g (car gpq)) (p (cadr gpq)) (q (caddr gpq)))
      (quotient (* p ch-C sqrt-C) (* ch-D (+ q (* p ch-A)))))))

(display (pi 1000))


Output:
0 72
0 36
0 18
0 9
0 4
0 2
0 1
1 2
2 4
2 3
3 4
4 9
4 6
4 5
5 6
6 9
6 7
7 9
7 8
8 9
9 18
9 13
9 11
9 10
10 11
11 13
11 12
12 13
13 18
13 15
13 14
14 15
15 18
15 16
16 18
16 17
17 18
18 36
18 27
18 22
18 20
18 19
19 20
20 22
20 21
21 22
22 27
22 24
22 23
23 24
24 27
24 25
25 27
25 26
26 27
27 36
27 31
27 29
27 28
28 29
29 31
29 30
30 31
31 36
31 33
31 32
32 33
33 36
33 34
34 36
34 35
35 36
36 72
36 54
36 45
36 40
36 38
36 37
37 38
38 40
38 39
39 40
40 45
40 42
40 41
41 42
42 45
42 43
43 45
43 44
44 45
45 54
45 49
45 47
45 46
46 47
47 49
47 48
48 49
49 54
49 51
49 50
50 51
51 54
51 52
52 54
52 53
53 54
54 72
54 63
54 58
54 56
54 55
55 56
56 58
56 57
57 58
58 63
58 60
58 59
59 60
60 63
60 61
61 63
61 62
62 63
63 72
63 67
63 65
63 64
64 65
65 67
65 66
66 67
67 72
67 69
67 68
68 69
69 72
69 70
70 72
70 71
71 72
31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989


Create a new paste based on this one


Comments: