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

(define (digits-of-pi n)
  ;; k=0 -> 3.2222222222 ...
  ;; k=1 -> 0.0212212212 ...
  ;; start vector with 3.243443443 ...
  (define digits (let ((v (make-vector n 4)))
		   (vector-set! v 0 2)
		   (do ((i 2 (+ i 3)))
		       ((>= i (vector-length v)) v)
		     (vector-set! v i 3))))
  (define (partial-term k)
    (let ([8k (* 8 k)])
      (- (/ 4 (+ 8k 1)) (/ 2 (+ 8k 4)) (/ 1 (+ 8k 5)) (/ 1 (+ 8k 6)))))
  (define (do-carry! v)
    (do ([i (- (vector-length v) 1) (- i 1)])
	((= i 0))
      (when (>= (vector-ref v i) 16)
	(vector-set! v (- i 1) (+ (vector-ref v (- i 1))
				  (fxdiv (vector-ref v i) 16)))
	(vector-set! v i (fxmod (vector-ref v i) 16)))))
  ;; k=2 -> 0.000acef6dc ...
  ;; k=3 -> 0.000055038c ...
  ;; k=4 -> 0.00000326fe ...
  (do ([k 2 (+ k 1)])
      ((= k n))
    (let* ([a (partial-term k)]
	   [den (denominator a)])
      (do ([i (+ k 1) (+ i 1)]
	   [num (fx* 16 (numerator a)) (fxmod (fx* 16 num) den)])
	  ((>= i n))
	(let ([v (vector-ref digits i)])
	  (vector-set! digits i (+ v (fxdiv (fx* 16 num) den)))))))
  (do-carry! digits)
  digits)

(define (display-digits-of-pi n)
  (let ([digits (digits-of-pi (+ n 4))])
    (display "3.")
    (do ([i 0 (+ i 1)])
	((= i n))
      (display (vector-ref digits i)))
    (newline)))

(display-digits-of-pi 100)

(define (digits-of-pi* n)
  (define digits (make-vector n 0))
  (define (term k)
    (let ([8k (* 8 k)])
      (/ (- (/ 4 (+ 8k 1)) (/ 2 (+ 8k 4)) (/ 1 (+ 8k 5)) (/ 1 (+ 8k 6)))
	 (expt 16 k))))
  ;; sum n+1 terms
  (define sum (let loop ([k 0] [s 0])
		(if (> k n) s (loop (+ k 1) (+ s (term k))))))
  (let ([den (denominator sum)])
    (do ([i 0 (+ i 1)]
	 [num (numerator sum) (* 16 (mod num den))])
	((>= i n) digits)
      (vector-set! digits i (div num den)))))