;; 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)