; 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 (number? 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)) "")))))
   (if (pair? x) (cons (hex (car x)) (hex (cdr x))) x)))

(define divide (let* ((ls arithmetic-shift) (p2 (lambda (n) (ls 1 n))) (÷ quotient)
    (Dos (fileVal "Do")) (DoV (Dos 'DoV)) (Do (Dos 'Do)) (2^52 (p2 52)))
  (lambda (M p) (let* ((2^M (p2 M)) (num (p2 (+ (* 2 M) 2)))
      (p (lambda (x) (if p (begin (write (x))(newline)))))
      (mx (DoV (+ 1 (÷ 53 M)) (λ (j) 0))) (mxb 0)
      (T (DoV 2^M (λ (j) (- (÷ num (+ 2^M j 1)) (* 2 2^M))))))
    (Do 2^M (lambda (j) (assert 'vecs
      (let ((t (vector-ref T j))) (and (integer? t) (< t (* 2 2^M)) (> t -1))))))
  (cons (lambda (n d) (let* ((nf (cdr n)) (df (cdr d)) (ad (if (< nf df) 1 0))
      (ed (+ 1023 (car n) (- (car d)) (- ad)))
      (i (ls df (- M 52)))
      (N (if #f (+ (* 2 2^M) (vector-ref T i)) (let* ((z 2))
        (÷ (ls num z) (+ (ls (+ 2^52 df) (+ z M -52)) 1)))))

      (Nd (* N (+ 2^52 df))) (Nn (ls (* N (+ 2^52 nf)) ad)))
    (p (λ () (list M i (hex N) 'Nn (hex Nn) 'Nd (hex Nd))))
    (let ((cap (λ (x) (assert 'args (and (integer? x) (exact? x) (< x 2^52) (<= 0 x))))))
      (cap df) (cap nf))
    (let iter ((r Nn)(q 0)(rb 53))
      (p (λ () (list 'r (hex r) 'q (hex q) 'rb rb)))
      (assert 'Invar (= (+ (* q Nd) r) (ls Nn (- 53 rb))))
      (assert 'Invb (<= 0 r))
      (if (>= rb 0) (let ((nq (ls r -54)))
          (p (λ () (list 'nq (hex nq))))
          (let ((i (÷ rb M))) (vector-set! mx i (max nq (vector-ref mx i))))
          (iter (- (ls r M) (* Nd nq)) (+ (ls q M) nq) (- rb M)))
        (begin (set! mxb (max mxb (÷ r Nd))) (let s ((r r)(q q))
          (assert (cons 'invar r) (= (+ (* q Nd) r) (ls Nn (- 53 rb))))
          (if (>= r Nd) (begin (p (λ () 'Harumph)) (s (- r Nd) (+ q 1)))
                 (begin (p (λ () (list 'qr (hex q) (hex r) mxb)))
                   (cons (cons (+ 1023 (car n) (- (car d)) (- ad))
                         (- (ls q (- rb 1)) 2^52))
        ; The one bit shift of q finally discards the rounding bit.
                     (+ (ls (bitwise-and (ls q rb) 1) 1)
                        (if (and (= 0 (bitwise-and q (- (p2 (- rb)) 1))) (= r 0)) 0 1)
                      ))))))))))
            (cons mx (lambda () mxb)))))))

; scrounge round. (rounds to nearest. ignores "to even" clause, probably unneeded!)
(define (gfdivr M p) (let* ((dvx (divide M p)) (dv (car dvx)))
  (cons (lambda (n d) (let ((fq (dv n d)))
  (cons (caar fq) (+ (cdar fq) (if (> (cdr fq) 1) 1 0))))) (cdr dvx))))

(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 fdiv (car (divide 8 (= 0 0))))
(define div (lambda (a b) (let ((p (divide 8 #t))) (cons ((car p) a b)
  (cons (cadr p) ((cddr p)))))))
(define (t a b) (let ((n (tx (* a b))) (d (tx a)))
  (write (list 't n d)) (newline) (hex (cons (fdiv n d) (tx b)))))
(define one '(1023 . 0))
(define ses (cons 1023 (expt 2 51)))
(define onep (cons 1023 (expt 2 0)))
(define onem (cons 1022 (- (expt 2 52) 1)))
; ----------
(let* ((p (gfdivr 8 (= 0 1))) (fdiv (car p)) (rf (((fileVal "RC4") "East") 'U))
    (t (λ (a b) (let* ((iq (fdiv (tx a) (tx b)))(Q (tx (/ a b)))
        (t (and (equal? iq Q) (< (vector-ref (cadr p) 0) 3000))))
        (if (not t) (begin (write (list a (tx a) b (tx b) iq Q)) (newline))) t))))
  (cons (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))))))
   (cons (cadr p) ((cddr p)))))
; => (#t #7(564 529 538 532 545 542 511) . 2)

(let ((c (λ (x) (cons 1024 (string->number x 16)))))
  (div (c "dc3829aaf2f6d") (c "e085d0264e5b1"))) 

(div one one)
(div onep one)
(div one onep)
(div one ses)
(div ses one)
(div (tx 2) (tx 3))
((car (gfdivr 7 (= 0 0))) '(1021 . 135191923634558) '(1022 . 2286672621417021))
(div one ses); => 5555555555555 (4/6) last q = (2^56)÷3
(div ses one); => 8∙2^48
(div '(1025 . 0) (cons 1024 (expt 2 51))) ; 4/3
(t 4/3 3)
(t 55.4 78)

(div (cons 1000 (expt 2 50)) (cons 1000 (expt 2 51))) ; => 35555555555555 (5/6)
(div '(1000 . 0) '(1000 . 0)) ; => fffffffffffff
(div '(1022 . 2751306630312438) '(1019 . 2373797444546464))