1 files changed, 72 insertions, 0 deletions
diff --git a/chudnovsky.scm b/chudnovsky.scm
new file mode 100644
index 0000000..7f143ef
--- /dev/null
+++ b/chudnovsky.scm
@@ -0,0 +1,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)