; We consider IEEE double float division of positive finite reals. ; We denote such a number as (cons exp frac) where both exp and frac are Scheme integers. ; (cons e f) denotes the real (2^(e-1023-52))(2^52+f). No denormals just yet. ; This is required: http://cap-lore.com/code/Scheme/reposIntro (define (assert n b) (if (not b) (begin (write (list 'Bad n)) (newline)))) (define (hex x) (if (zero? x) "0" (let y ((x x) (j 0)) (if (zero? (modulo x 16)) (y (/ x 16) (+ j 4)) (string-append (number->string x 16) (if (> j 0) (string-append "∙2^" (number->string j 10)) "")))))) (define divide (let* ((ls arithmetic-shift) (p2 (lambda (n) (ls 1 n))) (Dos (fileVal "Do")) (DoV (Dos 'DoV)) (Do (Dos 'Do)) (2^52 (p2 52)) (cap (λ (x) (assert 'args (and (integer? x) (exact? x) (< x 2^52) (<= 0 x)))))) (lambda (M p) (let* ((2^M (p2 M)) (p (lambda (x) (if p (begin (write (x))(newline))))) (T (DoV 2^M (λ (j) (- (truncate (/ (p2 (+ (* 2 M) 1)) (+ 2^M j 1))) 2^M))))) ; (if (= M 7) (vector-set! T 64 (- 170 2^M))) ; CHEAT!!!! ; (if (= M 5) (vector-set! T 16 (- 42 2^M))) ; CHEAT!!!! (Do 2^M (lambda (j) (assert 'vecs (let ((t (vector-ref T j))) (and (integer? t) (< t 2^M) (> t -1)))))) (lambda (n d) (let* ((df (cdr d)) (nf (cdr n)) (ad (if (< nf df) 55 54)) (ed (+ 1023 54 (- (- (car n) (car d)) ad))) (i (ls df (- M 52))) (T_i (+ 2^M (vector-ref T i))) (T_id (* T_i (+ 2^52 df))) (T_in (ls (* T_i (+ 2^52 nf)) (- ad 1)))) (p (λ () (list M i (hex T_i) 'T_id (hex T_id) 'T_in (hex T_in)))) (cap df) (cap nf) (let iter ((r T_in)(q 0)(rb ad)) (p (λ () (list 'r (hex r) 'q (hex q) 'rb rb))) (assert 'invar (= (+ (* q T_id) r) T_in)) (if (>= rb 0) (let ((nq (ls (ls r (- (+ 53 rb))) (- rb M)))) (p (λ () (list 'nq (hex nq) 'prod (hex (* T_id nq))))) (iter (- r (* T_id nq)) (+ q nq) (- rb M))) (begin (if (= 0 1) (assert 'end (and (>= r 0) (< r T_id)))) (assert (cons 'q q) (= (ls q -53) 1)) (p (λ () (list 'quo (quotient r T_id) (remainder r T_id)))) (let s ((r r)(q q)) (if (>= r T_id) (begin (p (λ () 'Harumph))(s (- r T_id) (+ q 1))) (cons (cons ed (- (ls q -1) 2^52)) ; The one bit shift of q above compensates for keeping the rounding bit. (+ (ls (bitwise-and q 1) 1) (if (= r 0) 0 1)))))))))))))) ; scrounge round. (ignores "to even" clause) (define (gfdivr M p) (let ((dv (divide M p))) (lambda (n d) (let ((fq (dv n d))) (cons (caar fq) (+ (cdar fq) (if (> (cdr fq) 1) 1 0))))))) (define fdiv (divide 8 (= 0 0))) (define (ex a b) (write (list a b)) (newline) b) (define tx (let* ((q (λ (x) (* x x x x))) (2^-99 (expt 2 -99)) (2^4 (q 2.0)) (2^16 (q 2^4)) (2^64 (q 2^16)) (2^256 (q 2^64))) (lambda (a) (assert 'pos (positive? a)) (let e ((a a)(f 0)) ; (write (list 'tx a (/ (log a) (log 2)) f (* a (expt 2 f)))) (newline) (if (< a 1.0) (e (/ a 2^-99) (- f 99)) (if (>= a 2^256) (e (/ a 2^256)(+ f 256)) (if (>= a 2^64) (e (/ a 2^64)(+ f 64)) (if (>= a 2^16) (e (/ a 2^16)(+ f 16)) (if (>= a 2^4) (e (/ a 2^4)(+ f 4)) (if (>= a 2.0) (e (/ a 2.0)(+ f 1)) (cons (+ 1023 f) (- (truncate (inexact->exact (* a (expt 2 52)))) (expt 2 52))) )))))))))) (define (t a b) (let ((n (tx (* a b))) (d (tx a))) (write (list 't n d)) (newline) (list (fdiv n d) (tx b)))) (define (rfg s) (((fileVal "RC4") s) 'U)) (define one '(1023 . 0)) (define ses (cons 1023 (expt 2 51))) (fdiv one ses); => 5555555555555 (4/6) last q = (2^56)÷3 (fdiv ses one) (fdiv one one) (fdiv '(1025 . 0) (cons 1024 (expt 2 51))) ; 4/3 ---------- (let* ((fdiv (gfdivr 8 (= 0 1))) (rf (rfg "North")) (t (λ (a b) (let* ((iq (fdiv (tx a) (tx b)))(Q (tx (/ a b)))(t (equal? iq Q))) (if (not t) (begin (write (list a (tx a) b (tx b) iq Q)) (newline))) t)))) (and (t 2 3) (t 1 1) (t 3 2) (t 1 1/3) (t 1/3 1/3) (t 0.2575046593181226 0.7538716594077141) (let e ((n 100000)) (or (= n 0) (and (t (rf) (rf)) (e (- n 1))))))) (t 4/3 3) (t 55.4 78) (t 1 3) (fdiv (cons 1000 (expt 2 50)) (cons 1000 (expt 2 51))) ; => 35555555555555 (5/6) (fdiv '(1000 . 0) '(1000 . 0)) ; => fffffffffffff