; 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 (let ((DoV ((fileVal "Do") 'DoV))) (lambda (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)))
   (if (vector? x) (DoV (vector-length x) (λ (j) (hex (vector-ref x j))))))))))

(define divide (let* ((ls arithmetic-shift) (p2 (lambda (n) (ls 1 n))) (÷ quotient)
       (2^52 (p2 52)))
  (lambda (M) (let* ((2^M (p2 M)) (num (p2 (+ (* 2 M) 2))))
  (lambda (n d) (let* ((nf (cdr n)) (df (cdr d)) (ad (if (< nf df) 1 0))
      (N (let* ((z 1)) (÷ (ls num z) (+ (ls (+ 2^52 df) (+ z M -52)) 1))))
      (Nd (* N (+ 2^52 df))) (Nn (ls (* N (+ 2^52 nf)) ad)))
    (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))
      (assert 'Invar (= (+ (* q Nd) r) (ls Nn (- 53 rb))))
      (assert 'Invb (<= 0 r))
      (if (>= rb 0) (let ((nq (ls r -54)))
          (iter (- (ls r M) (* Nd nq)) (+ (ls q M) nq) (- rb M)))
        (let s ((r r)(q q))
          (assert (cons 'invar r) (= (+ (* q Nd) r) (ls Nn (- 53 rb))))
          (if (>= r Nd) (s (- r Nd) (+ q 1))
                 (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)
                      ))))))))))))

; scrounge round. (rounds to nearest. ignores "to even" clause, probably unneeded!)
(define (gfdivr M) (let* ((dv (divide M)))
  (lambda (n d) (let ((fq (dv n d)))
  (cons (caar fq) (+ (cdar fq) (if (> (cdr fq) 1) 1 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 fdiv (divide 8))
(define div fdiv)
(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* ((fdiv (gfdivr 8)) (rf (((fileVal "RC4") "full") 'U))
    (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 

(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) ; => ((("4∙2^8" . "8∙2^48") . "1") "4∙2^8" . "8∙2^48")
(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))