Egyptian Fractions
June 4, 2013
The largest unit fraction that will fit a given ratio n ÷ d is ⌈ d / n ⌉. Given that, the needed function is simple:
(define (egypt n d)
(let loop ((n n) (d d) (xs (list)))
(if (= n 1) (reverse (cons d xs))
(let* ((x (ceiling (/ d n))) (y (- (/ n d) (/ x)))
(n (numerator y)) (d (denominator y)))
(loop n d (cons x xs))))))
We return the fraction as a list of denominators in ascending order. Here are some examples:
> (egypt 5 6)
(2 3)
> (egypt 7 15)
(3 8 120)
> (egypt 5 121)
(25 757 763309 873960180913 1527612795642093418846225)
You can run the program at http://programmingpraxis.codepad.org/eopIofz7. If you are interested, the study of Egyptian fractions is a fascinating topic; let Google be your guide.
Programming Praxis 
[…] today’s Programming Praxis exercise, our goal is to convert a given fraction to a sum of fractions with […]
My Haskell solution (see http://bonsaicode.wordpress.com/2013/06/04/programming-praxis-egyptian-fractions/ for a version with comments):
import Data.Ratio egypt :: Integer -> Integer -> [Integer] egypt 1 d = [d] egypt n d = e : egypt (numerator r) (denominator r) where (e,r) = (div (d+n-1) n, n%d - 1%e)here’s my general idea.
while numerator > 1:
subtract the largest factor that’s smaller than the numerator; distribute.
if no such factor:
multiply top and bottom by the smallest value that will put the numerator greater than the smallest factor bigger than it.
for example:
7/15 = (2+5)/15 = 1/3 +2/15 = 1/3 +4/30 = 1/3 +1/10 +1/30
25/37 = 50/(2*37) = (37+13)/(2*37) = 1/2 +13/(2*37) = 1/2 +1/37 +11/(2*37) = 1/2 +1/37 +44/(2*4*37) = 1/2 +1/37 +1/8 +7/(2*4*37)
= 1/2 +1/37 +1/8 +1/74 +3/(2*4*37) = 1/2 +1/37 +1/8 +1/74 +1/148 +1/296.
http://pastebin.com/y9KgqXs1
phil’s program loops indefinitely with inputs of 2/4 and 6/8 (perhaps others.) Maybe it needs to reduce fraction to lowest terms before proceeding.
okay thank you david. i tried to fix the program but there seems to be too many road blocks preventing this algorithm form working. currently
31/466 can’t be done.
Clojure solution for a pairing algorithm. Represent x/y as a list of x copies of the fraction 1/y. Reduce the list each step by looking for pairs with the same denominator. Replace each pair with 2/y (when y is even) or 2/(y+1) + 2/(y(y+1)) when y is odd. In each case, the 2 can be factored out, so the numerator will still be 1. The algorithm halts when there are no more equal denominators in the list.
(defn reduce-pairs [l] (cond (< (count l) 2) l (= (first l) (second l)) (let [den (denominator (first l))] (if (even? den) (cons (/ 2 den) (reduce-pairs (rest (rest l)))) (cons (/ 2 (inc den)) (cons (/ 2 (* den (inc den))) (reduce-pairs (rest (rest l))))))) :else (cons (first l) (reduce-pairs (rest l))))) (defn egypt [x y] (loop [l (repeat x (/ 1 y))] (let [l' (sort > (reduce-pairs l))] (if (= l l') l (recur l'))))) user=> (egypt 18 23) (1/2 1/6 1/12 1/35 1/276 1/2415)The Egyptians probably used the pair reduction method because they have a lot of ways in which they would reduce a fraction 2/n where n is odd. but they employed various heuristics to keep the denominators small when n is odd. So this version is also a pair reduction, but uses a lookup table for small odd denominators that looks up the fractions that were used on the Rhind Papyrus. For denominators greater than 101, the identity 1/n + 1/n = 2/(n+1) + 2/(n(n+1) is used.
(def rhind-papyrus [[1/2 1/6], [1/3 1/15], [1/4 1/28] ; 2/3, 2/5, 2/7 [1/6 1/18], [1/6 1/66], [1/8 1/52 1/104] ; 2/9, 2/11, 2/13 [1/10 1/30], [1/12 1/51 1/68], [1/12 1/76 1/114] ; 2/15, 2/17, 2/19 [1/14 1/42], [1/12 1/276], [1/15 1/75] ; 2/21, 2/23, 2/25 [1/18 1/54], [1/24 1/58 1/174 1/232] ; 2/27 2/29 [1/20 1/124 1/155], [1/22 1/66], [1/30 1/42] ; 2/31, 2/33, 2/35 [1/24 1/111 1/296], [1/26 1/78], [1/24 1/246 1/328] ; 2/37, 2/39, 2/41 [1/42 1/86 1/129 1/301], [1/30 1/90] ; 2/43, 2/45 [1/30 1/141 1/470], [1/28 1/196], [1/34 1/102] ; 2/47, 2/49, 2/51 [1/30 1/318 1/795], [1/30 1/330], [1/38 1/114] ; 2/53, 2/55, 2/57 [1/36 1/236 1/531], [1/40 1/244 1/488 1/610] ; 2/59, 2/61 [1/42 1/126], [1/39 1/195], [1/40 1/335 1/536] ; 2/63, 2/65, 2/67 [1/46 1/138], [1/40 1/568 1/710], [1/60 1/219 1/292 1/365] ; 2/69, 2/71 2/73 [1/50 1/150], [1/44 1/308], [1/60 1/237 1/316 1/790] ; 2/75, 2/77, 2/79 [1/54 1/162], [1/60 1/332 1/415 1/498] ; 2/81, 2/83 [1/39 1/195], [1/58 1/174], [1/60 1/356 1/534 1/890] ; 2/85, 2/87, 2/89 [1/70 1/130], [1/62 1/186], [1/60 1/380 1/570] ; 2/91 2/93 2/95 [1/56 1/679 1/776], [1/66 1/198], [1/101 1/202 1/303 1/606]] ; 2/97 2/99 2/101 ) (defn ahmes [n] (rhind-papyrus (dec (/ (dec n) 2)))) (defn reduce-pairs [l] (cond (< (count l) 2) l (= (first l) (second l)) (let [den (denominator (first l))] (cond (even? den) (cons (/ 2 den) (reduce-pairs (rest (rest l)))) (<= den 101) (reduce conj (ahmes den) (reduce-pairs (rest (rest l)))) :else (cons (/ 2 (inc den)) (cons (/ 2 (* den (inc den))) (reduce-pairs (rest (rest l))))))) :else (cons (first l) (reduce-pairs (rest l))))) (defn egypt [x y] (loop [l (repeat x (/ 1 y))] (let [l' (sort > (reduce-pairs l))] (if (= l l') l (recur l'))))) user=> (egypt 18 23) (1/2 1/6 1/12 1/46 1/138 1/276)There is a binary method which takes advantage of the fact that a non-repeating fraction in binary notation is already a sum of unit fractions (e.g. 5/16 = 0.0101 = 1/16 + 1/4, so the method splits a fraction x/y into a sum of inverse powers of two and a remainder portion that is also inverse powers of two multiplied by 1/y. See http://www.ics.uci.edu/~eppstein/numth/egypt/binary.html. Cook thing is that the number of terms is bounded by log(x) + log(y) and denominator is no larger than y^2.
(def powers_of_two (iterate (partial * 2) 1)) (defn bin_frac [x y] ; create Egyptian fraction from x/y; y MUST BE power of 2 (loop [x x, y y, result ()] (cond (= x 0) result (odd? x) (recur (quot x 2) (quot y 2) (cons (/ 1 y) result)) :else (recur (quot x 2) (quot y 2) result)))) (defn bin_divisor [x y] (first (drop-while (fn [p] (>= (mod (* x p) y) (* p 2))) powers_of_two))) (defn egypt [x y] (if (> x y) (cons (quot x y) (egypt (rem x y) y)) ;else (let [p (bin_divisor x y), r (rem (* x p) y), s (quot (* x p) y)] (concat (bin_frac s p) (map #(/ %1 y) (bin_frac r p)))))) user=> (egypt 5 121) (1/32 1/121 1/968 1/1936 1/3872) user=> (egypt 355 113) ; pi (3 1/8 1/113 1/226 1/452 1/904) user=> (egypt 18 23) (1/2 1/4 1/46 1/92)[…] This post presents C# solutions to a coin change problem as described in https://programmingpraxis.com/2013/06/04/egyptian-fractions. […]