; This version keeps low zeros of numerators thruout.
(define (P2 n) (expt 2 n))
(define l2 (let ((o (/ (log 2)))) (lambda (x) (* o (log x)))))
; (define (bn b n) (odd? (quotient b (P2 n))))

(define (assert m t) (if (not t) (begin (write (list "Fail" m))
  (newline))))
(define (Pb b n) (let* (
   (st (number->string b 2))(l (string-length st)))
  (string-append (substring st 0 (min l n))
  (if (< l n) "-" "...") (number->string l))))
(define (rp m v) (write (list m v (Pb v 13))) (newline) v)
(define (ph n) (number->string n 16))
(define (join m a b) (if (equal? a b) a (begin
   (write (list "Different" m (Pb a 54) (Pb b 54))) (newline) a)))
(define mx 0) ; mutable cel!!!

; div takes two integers, n, d.
; d is between 2^52 and 2^53.
; div returns the integer floor(2^53 * (n/d)).
; b is the number of quotient bits per cycle.
; tb is log of table size.
; tz is number of bits in table entry.
(define div (let ((b 9)(tb 7)(tz 9))
  (let ((ptb (P2 tb))(ptz (P2 tz))(pb (P2 b))
    (p52 (P2 52))(p53 (P2 53))(pbt (P2 (- 52 tb))))
  (let ((tab (make-vector ptb)))
  (let lp ((j (- ptb 1))) (if (< j 0) () (begin
    (vector-set! tab j (- (quotient (P2 (+ tb tz 1)) (+ ptb j)) 1))
    (lp (- j 1)))))
(write (list "t[76]" (vector-ref tab 76))) (newline)
  (lambda (n d) (assert 2 (and (integer? d) (>= d p52) (< d p53)))
(write (list "den=" (ph d))) (newline)
      (let* ((n (* n p53))
     (N (vector-ref tab (- (quotient d pbt) ptb)))
     (Nn (* (rp "N=" N) n))(Nd (* N d)))
(write (list "Zil" (Pb Nn 64) (Pb Nd 64))) (newline)
    (join "quo" (quotient Nn Nd)
     (let dz ((n (* ptz Nn))(divp (* ptz p53 ptb))(q 0))
     (if (= divp ptz) (begin
         (write (list "qw" (ph divp) (ph n) (ph Nd)))
         (newline)
         (quotient q (* ptz ptz)))
        (let ((td (quotient n (* p53 divp))))
        (assert 1 (>= td 0))
        (set! mx (max mx td))
(write (list "PLL" td mx (ph (* td Nd)))) (newline)
      (dz (- n (/ (* Nd td divp) ptz))
          (max (quotient divp pb) ptz)
          (+ (* td divp) q))))))))))))

; dv takes real positive Scheme reals between 2^52 and 2^53.
; dv returns 2^52 times their quotient.
(define dv (let* ((cv (lambda (x)
   (inexact->exact (floor (if (and (number? x) (real? x)
  (>= x (P2 52)) (< x (P2 53))) x
    (begin (write (list x "bam")) "foo")))))))
  (lambda (n d) (div (cv n) (cv d)))))

; dx takes real Scheme numbers between 1/2 and 1
; and returns the floating quotient.
(define (dx n d) (set! mx 0) (let ((y (join "end" 
  (+ 0. (* (P2 -53) (dv (* n (P2 53)) (* d (P2 53)))))
  (+ 0. (/ n d)))))
  (cons mx y)))

(dx 3/4 1/2)
(dx 0.6922436 0.8243664)
(dx 4/7 3/5)