; 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