; compute sine
; Better compute sin(x/10(10^2))*10^(10^2) for integer x.
; http://cap-lore.com/MathPhys/PI.html
x - x^3/6 + x^5/120
x(1 - (1/6)x^2(1-(1/20)x^2(1 - ...)))
(define (M a b) (quotient (* a b) D))
(define (id x y) (quotient y x))
(define D (expt 10 100))
(define x (* 3 D))
(define x2 (M x x))
(M x (- D (id (* 2 3) (M x2
     (- D (id (* 4 5) (M x2 1)))))))
(define (s n) (if (> n 88) 0 (- D (id (* n (+ n 1)) (M x2 (s (+ n 2)))))))
(define (Sin x) (M x (s 2)))
(+ x (M x (s 2)))

(set! x (+ x (Sin x)))
(- 0.1 (/ 0.001 6))
(* x(- 1 (s 2)))
(sin 0.1)
-------
(define Sin (lambda (D) (let ((M (lambda (a b) (quotient (* a b) D))))
  (lambda (x) (let ((x2 (M x x))) (letrec ((s (lambda (n) 
      (if (> n 94) 0 (- D (quotient (M x2 (s (+ n 2))) (* n (+ n 1))))))))
    (M x (s 2))))))))
(define D (expt 10 100))
(define sine (Sin D))
(define (T x) (+ x (sine x)))
(T (T (T (T (* 3 D)))))
; => 31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170676 which is (10^100)pi-3.
-------
(define (Sin D t) (let ((M (lambda (a b) (quotient (* a b) D))))
  (lambda (x) (let ((x2 (M x x))) (letrec ((s (lambda (n) 
      (if (> n t) 0 (- D (quotient (M x2 (s (+ n 2))) (* n (+ n 1))))))))
    (M x (s 2)))))))
(define (pi n) (let more ((pi 31) (D 10) (t 4)) (let* ((D2 (* D D))
    (DN (* D D2))
    (np (* pi D2))
    (nt (* 3 t))
    (x (+ np ((Sin DN nt) np))))
;   (display (list x np DN)) (newline)
   (if (> (/ (log DN) (log 10)) n) x (more x DN nt)))))
(pi 1000) ; yields pi to 2005 places in about 1.35 seconds.
; digits "275900994"
------
(define (pi n) (let ((tn (expt 10 n)))
  (let more ((pi 6) (D 2) (t 3)) (let* ((D2 (* D D))
    (DN (* D D2))
    (np (* pi D2))
    (nt (* 3 t))
    (x (+ np ((Sin DN nt) np))))
;   (display (list DN np)) (newline)
   (if (> DN tn) (quotient (* x tn) DN)
        (more x DN nt))))))
;(pi 1000) => ... 66111959092164201989
------
(define as arithmetic-shift)
(define (Sin d t) (let ((M (lambda (a b) (as (* a b) (- d)))))
  (lambda (x) (let ((x2 (M x x))) (letrec ((s (lambda (n) 
      (if (> n t) 0 (- (as 1 d) (quotient (M x2 (s (+ n 2))) (* n (+ n 1))))))))
    (M x (s 2)))))))
(define (pi n) (let ((tn (expt 10 n)))
  (let more ((pi 6) (d 1) (t 2)) (let* ((d2 (+ d d))
    (dn (+ d d2))
    (np (as pi d2))
    (nt (* 3 t))
    (x (+ np ((Sin dn nt) np))))
;   (display (list dn np)) (newline)
   (if (> (as 1 dn) tn) (as (* x tn) (- dn))
        (more x dn nt))))))
;(pi 1000) => ... 66111959092164201989