; 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 fdiv (let* ((M 8) (ls arithmetic-shift) (p2 (lambda (n) (ls 1 n)))
     (2^M (p2 M)) (Do (fileVal "Do")) (2^52 (p2 52)) (DoV (Do 'DoV))
     (T (DoV 2^M (λ (j) (- (truncate (/ (p2 (+ (* 2 M) 1)) (+ 2^M j 1))) 2^M))))
     (cap (λ (x) (assert 'args (and (integer? x) (exact? x) (< x 2^52) (<= 0 x))))))
;  (if (= M 7) (vector-set! T 64 (- 170 2^M))) ; CHEAT!!!!
 ; (if (= M 5) (vector-set! T 16 (- 42 2^M))) ; CHEAT!!!!
  ((Do '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))))
    (write (list M i (hex T_i) 'T_id (hex T_id) 'T_in (hex T_in))) (newline)
    (cap df) (cap nf)
    (let iter ((r T_in)(q 0)(rb ad))
      (write (list 'r (hex r) 'q (hex q) 'rb rb)) (newline)
      (assert 'invar (= (+ (* q T_id) r) T_in))
      (if (>= rb 0) (let ((nq (ls (ls r (- (+ 53 rb))) (- rb M))))
               (write (list 'nq (hex nq) 'prod (hex (* T_id nq))))(newline)
               (iter (- r (* T_id nq)) (+ q nq) (- rb M)))
            (begin (assert 'end (and (>= r 0) (< r T_id)))
                   (assert 'q (= (ls q -53) 1))
                   (write (list 'quo (quotient r T_id) (remainder  r T_id))) (newline)
                   (let s ((r r)(q q)) (if (>= r T_id) (begin 
                     (write 'Harumph)(newline)(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))))))))))))

(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 one '(1023 . 0))
(define ses (cons 1023 (expt 2 51)))
----------
(fdiv one ses); => 2aaaaaaaaaaaaa (4/6) last q = (2^56)÷3
(fdiv ses one)
(fdiv one one)
(fdiv '(1025 . 0) (cons 1024 (expt 2 51))) ; 4/3
(t 4/3 3)
(t 55.4 78)

(fdiv (cons 1000 (expt 2 50)) (cons 1000 (expt 2 51))) ; => aaaaaaaaaaaaa (5/6)
(fdiv '(1000 . 0) '(1000 . 0)) ; => one