From 1ef25a125370b40b183fa37eb2d73b8c60c9794b Mon Sep 17 00:00:00 2001 From: Thomas Guillermo Albers Raviola Date: Wed, 21 Jan 2026 16:15:50 +0100 Subject: Initial commit --- chudnovsky.scm | 72 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 chudnovsky.scm (limited to 'chudnovsky.scm') 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 . + +(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) -- cgit v1.2.3