blob: 7f143ef88fc06255a11899d40224dfb03816e32c (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
|
;; Copyright (C) 2026 Thomas Guillermo Albers Raviola
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program. If not, see <https://www.gnu.org/licenses/>.
(import (rnrs)
(srfi srfi-1)
(srfi srfi-8)
(srfi srfi-11))
(define (binary-splitting n base-case)
(define (iteration a b)
(cond [(= a (- b 1))
(base-case b)]
[(> b a)
(let ([m (fxdiv (fx+ a b) 2)])
(let-values ([(p1 q1 r1) (iteration a m)]
[(p2 q2 r2) (iteration m b)])
(values (+ (* p1 q2) (* p2 r1))
(* q1 q2)
(* r1 r2))))]
[else
(error 'iteration "a >= b")]))
(iteration 0 n))
(define (sum n)
(define (base-case b)
(let ([r (* (- (* 2 b) 1) (- (* 6 b) 1) (- (* 6 b) 5))])
(values (* (expt -1 b) (+ (* 545140134 b) 13591409) r)
(* 10939058860032000 b b b)
r)))
(receive (p q r)
(binary-splitting n base-case)
(+ 13591409 (/ p q))))
(define (factor n)
;; compute sqrt(1+x) using taylor series and binary splitting
(define x 5/10000)
(define (base-case b)
(values (* (expt -1 (+ b 1)) 2 b x)
(* 4 b b)
(* 2 b (- (* 2 b) 1) x)))
(receive (p q r)
(binary-splitting n base-case)
(* 42688000 (+ 1 (/ p q)))))
(define (compute-pi a b)
(/ (factor a) (sum b)))
;; TODO: how to determine a and b from k?
(define (display-digits-of-pi a b k)
(let ([pi (compute-pi a b)])
(let loop ([i 0]
[n (numerator pi)])
(receive (d m) (div-and-mod n (denominator pi))
(display d)
(when (= i 0)
(display "."))
(unless (= i k)
(loop (+ i 1) (* 10 m)))))))
(display-digits-of-pi 10 10 10)
|
