## 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.

Pages: 1 2

### 9 Responses to “Egyptian Fractions”

1. […] today’s Programming Praxis exercise, our goal is to convert a given fraction to a sum of fractions with […]

```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)
```
3. phil said

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

4. David said

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.

5. phil said

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.

6. David said

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)
```
7. David said

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)
```
8. David said

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)
```
9. […] This post presents C# solutions to a coin change problem as described inÂ http://programmingpraxis.com/2013/06/04/egyptian-fractions. […]