# File: EgyptianFractions.py
# Author: Keith Schwarz (htiek@cs.stanford.edu)
#
# An implementation of the greedy algorithm for decomposing a fraction into an
# Egyptian fraction (a sum of distinct unit fractions). Egyptian fractions
# are a representation of fractions that dates back at least 3500 years (the
# Rhind Mathematical Papyrus contains a table of fractions written out this
# way). A number is expressed as an Egyptian fraction if it is written as a
# sum of unit fractions (fractions whose numerators are all zero) such that no
# fraction is duplicated. For example, 1/2 is already an Egyptian fraction,
# but 2/3 is not. However, 2/3 = 1/2 + 1/6 is an Egyptian fraction, though
# 2/3 = 1/3 + 1/3 is not because 1/3 is duplicated.
#
# Any rational number has at least one representation as an Egyptian fraction,
# and one simple algorithm for finding an Egyptian fraction representation of
# a rational number is given in Fibonacci's Liber Abaci (the same text that
# contains the eponymous Fibonacci sequence and the introduction of Hindu-
# Arabic numerals to Europe). This algorithm is a greedy algorithm that works
# as follows: if the rational number already has a numerator of 1, then we are
# done. Otherwise, subtract out the largest possible unit fraction from the
# number and repeat this process. For example, to compute an Egyptian fraction
# representation of 42/137, we would write
#
# 42/137 = 1/4 + 31/548
# = 1/4 + 1/18 + 5/4932
# = 1/4 + 1/18 + 1/987 + 1/1622628
#
# Notice that the denominator may grow very large; here, starting with 42/137,
# the last denominator is over a million!
#
# In order to prove correctness of this algorithm, we need to show that it
# eventually terminates and that it is correct. We restrict ourselves to
# rational numbers in the range (0, 1) and then prove, by induction on the
# numerator a, that for any rational number a/b in this range, the algorithm
# terminates with a legal Egyptian fraction.
#
# As our base case, if the numerator is 1, then our fraction is already an
# Egyptian fraction and we are done. For the inductive step, assume that for
# any numerator a' such that 1 <= a' <= a, that the claim holds. We need to
# show that for any rational number a / b with a < b that the greedy algorithm
# terminates and is correct. To do this, let 1 / n be the largest unit
# fraction smaller than a / b. This means that 1 / n < a / b < 1 / (n - 1).
# Now consider the difference a / b - 1 / n. This is given by
#
# a 1 an - b
# --- - --- = ------
# b n bn
#
# Now, since we have that 1 / n < a / b, this means that b < an, so 0 < an - b.
# Also, since a / b < 1 / (n - 1), we have that a(n - 1) < b, so an - a < b, so
# an - b < a. Thus an - b is a positive natural number smaller than a, so by
# our inductive hypothesis we can write (an - b) / bn as the sum of distict
# unit fractions. This sum, plus 1/n, is equal to a / b.
#
# To conclude that the algorithm is correct, we need to show that the sum of
# 1 / n and the unit fractions from (an - b) / bn do not contain any duplicate
# unit fractions. Our inductive hypothesis says that none of the fractions
# from the (an - b)/bn terms can be duplicates, so the only possible duplicates
# we could have must be the fraction 1 / n. This would only be possible if we
# had that 1 / n + 1 / n = 2 / n < a / b. Now, recall that a / b < 1 / (n - 1)
# because n is the denominator of the largest unit fraction no larger than a/b.
# This means tha 2 / n < 1 / (n - 1), so 2(n - 1) < n, so 2n - 1 < n, so n < 1.
# But this is impossible, because this would mean that the denominator is 0.
# We have reached a contradiction, so our assumption was wrong and no fraction
# is duplicated. Thus every rational number a / b in the range (0, 1) has an
# Egyptian fraction representation that can be found using the greedy
# algorithm.
#
# All that remains to get an efficient implementation of this algorithm is a
# way to find, given an arbitrary rational number a / b in the range (0, 1),
# the largest unit fraction smaller than a / b. Write b = ka + r using
# the division algorithm; we know that r != 0 because otherwise a / b = 1 / k,
# which is a unit fraction. This means that b = ka + r for some r in the range
# 1 < r < a. Consequently, we have that
#
# a a 1
# --- = ------ = -------
# b ka + r k + r/a
#
# Now, since 1 < r < a, 1/a < r/a < 1. This means that a / b = 1 / (k + r/a)
# is smaller than 1 / k, but it is larger than 1 / (k + 1). Consequently, the
# largest unit fraction smaller than a / b is 1 / (k + 1) =
# 1 / (floor(b / a) + 1). This gives the following algorithm:
#
# If a = 1, return a / b. (It's already an Egyptian fraction)
# Otherwise:
# Let k = floor(b / a) + 1.
# Append 1 / k to the list of unit fractions.
# Recursively compute a fraction for a / b - 1 / k.
#
#
# This implementation unrolls this recursion into a simple iterative solution,
# which is used in this code.
#
# The analysis of this algorithm is a bit tricky. We know from our correctness
# proof that the numerator of the fraction must decrease on each iteration, so
# the maximum possible work is O(a f(a, b)), where f(a, b) is the amount of
# work required on each iteration. Notice that if we begin the iteration with
# the fraction a / b, on the next iteration we have the fraction
#
# a 1 ak - b
# --- - ----- = --------
# b k bk
#
# Notice that this means the denominator is growing by a factor of k, where k
# is the unit fraction subtracted out. But this is (roughly) b / a, so if we
# don't simplify this fraction we get that the denominator has gone from b to
# (roughly) b^2 / a. In other words, if b is large relative to a, this may
# have twice as many bits in it as we began with.
#
# On each iteration of the algorithm, we compute one divide (to compute the
# denominator of the unit fraction to subtract out) and one rational number
# subtraction. We assume that dividing b by a takes O((lg b + lg a)^2) time
# (a reasonable assumption). Since we know that a < b (because we restrict
# ourselves to fractions in the range (0, 1), this means that the division step
# of dividing b by a takes O(lg^2 b) time. But we know that on each iteration
# the number of bits in b may double, so on iteration i this time requirement
# is given by O((2^i lg b))^2) = O(4^i lg^2 b). Summing this up over all O(a)
# iterations, we get that the runtime would be O(4^a lg^2 b), which is
# doubly-exponential in the number of bits of a. Since the work done to do the
# division is a constant fraction of the total work done, this would give a
# net runtime of O(4^a lg^2 b).
from fractions import Fraction
def greedyEgyptianFraction(rational):
# Sanity check: the rational number should be in the range (0, 1)
if rational <= 0 or rational >= 1:
raise Exception("Rational number out of range" , rational)
# Create a list to store the Egyptian fraction representation.
result = []
# Now, iteratively subtract out the largest unit fraction that may be
# subtracted out until we arrive at a unit fraction.
while True:
# If the rational number has numerator 1, we're done.
if rational.numerator == 1:
result.append(rational)
return result
# Otherwise, find the largest unit fraction less than the current
# rational number. This is given by the ceiling of the denominator
# divided by the numerator
unitFraction = Fraction(1, rational.denominator / rational.numerator + 1)
result.append(unitFraction)
# Subtract out this unit fraction.
rational = rational - unitFraction