Last Hallowe'en (minus epsilon) I boasted
A trickier problem, for which the mediant algorithm seems better suited: What is the nearest fraction to x with numerator <= a and denominator <= b? e.g., (farint pi 999) 355 113 (farint pi 999 112) 333 106
but notice I said "for hand computation" For some reason (personal stupidity?), the iterated mediant is ugly to code. The following was maybe my sixth try: (defun farint (x &optional (nhi 400) (dhi nhi) (ln 0) (ld 1) (hn 1) (hd 0)) (if (> (+ ln hn) (* (+ ld hd) x)) (let* ((m (min (if (= 0 ln) nhi (floor(- nhi hn) ln)) (floor (- dhi hd) ld))) (d (- (* x ld) ln)) (k (if (= 0 d) m (ceiling (- hn (* x hd)) d)))) (if (< k m) (farint x nhi dhi (setf hn (+ (* k ln) hn)) (setf hd (+ (* k ld) hd)) (- hn ln) (- hd ld)) (let* ((n (+ (* m ln) hn)) (d (+ (* m ld) hd))) (if (< (* 2 d ld x) (+ (* ld n) (* ln d))) (values ln ld) (values n d))))) [...] After seven months and literally millions of calls (Danny Hillis is doing another clock), I cringe to announce that (farint 207727/146410 250) gives 227/160 = 1.41875, while the correct answer is 166/117 ~ 1.4188034+, and 207727/146410 ~ 1.41880336+. After another few tries, ;0 (defun farint (x &optional (nhi 400) (dhi nhi) (ln 0) (ld 1) (hn 1) (hd 0)) (if (> (+ ln hn) (* (+ ld hd) x)) ;1 (let* ((m (min (if (= 0 ln) nhi (floor(- nhi hn) ln)) (floor (- dhi hd) ld))) ;2 (d (- (* x ld) ln)) ;3 (k (if (= 0 d) m (ceiling (- hn (* x hd)) d)))) ;4 (if (< k m) (farint x nhi dhi (setf hn (+ (* k ln) hn)) ;5 (setf hd (+ (* k ld) hd)) (- hn ln) (- hd ld)) (let* ((n (+ (* m ln) hn)) (d (+ (* m ld) hd))) ;6 (if (< (* 2 d ld x) (+ (* ld n) (* ln d))) (values ln ld) ;7 (let* ((hn (- n ln)) (hd (- d ld))) ;8 (if (> (* 2 d hd x) (+ (* hd n) (* hn d))) (values hn hd) (values n d))))))) ;9 (let* ((m (min (floor (- nhi ln) hn) ;10 (if (= 0 hd) dhi (floor (- dhi ld) hd)))) (d (- hn (* x hd))) (k (if (= 0 d) m (ceiling (- (* x ld) ln) d)))) (if (< k m) (farint x nhi dhi (- (setf ln (+ (* k hn) ln)) hn) (- (setf ld (+ (* k hd) ld)) hd) ln ld) (let* ((n (+ (* m hn) ln)) (d (+ (* m hd) ld))) (if (> (* 2 d hd x) (+ (* hd n) (* hn d))) (values hn hd) (let* ((ln (- n hn)) (ld (- d hd))) (if (< (* 2 d ld x) (+ (* ld n) (* ln d))) (values ln ld) (values n d))))))))) Comments (= prayerful handwaving): ;0 Find closest n/d to given x with n <= nhi, d <= dhi. ; Narrow a bracketing interval, initially [0/1,1/0], via ; recursive bisection, using MEDIANT for the mean. ;1 If "midpoint" is also high, shade upper limit downward ; by combining with a multiple, k, of lower. ;2 m = largest k w/o exceeding nhi and dhi. ;3 d = denom of k estimate. ;4 k just sufficient to overshade to an underestimate. ;5 if k small enough, recurse with tightened brackets, k-1 yielding ; the new overestimate. ;6 if k too big, we're limited to m. ;7 But we also have the option of k = infinity, i.e. shade all the way ; to lower limit, which is sometimes better. ;8 And sometimes m-1 is better, because the CEILING overshot! ;9 But usually not. ;10 If "midpoint" was low, shade lower brackets upward instead ... Despite the ugliness, this is pretty fast. --rwg