egyptian_greedy.py

from fractions import Fraction as F, gcd
from os import environ as env

__all__ = ["greedy_solve"]


def _find_next(t, c):
    p = c
    zero = F(0, 1)
    while t - F(1, c) < zero:
        p = c
        c = 2 * c

    while c - p > 1:
        mid = (p + c) / 2
        if t - F(1, mid) >= zero:
            c = mid
        else:
            p = mid

    return c        


def greedy_solve(t):
    coll = []
    zero = F(0, 1)
    c = 1
    while True:
        if t - F(1, c) >= zero:
            coll.append(c)
            t -= F(1, c)
            if t == zero:
                break
        else:
            c = _find_next(t, c)

    return coll


TEST = "TEST" in env

if __name__ == "__main__":
    if not TEST:
        from sys import argv
        nom, denom = map(int, argv[1:3])
        g = gcd(nom, denom)
        nom, denom = nom / g, denom / g
        print greedy_solve(F(nom, denom))
    else:
        from random import randint as ri
        for i in xrange(30):
            nom, denom = sorted([ri(1, 1000), ri(1, 1000)])
            g = gcd(nom, denom)
            nom, denom = nom / g, denom / g
            res = greedy_solve(F(nom, denom))
            print "test %s: %s / %s = %s" % (i, nom, denom, " + ".join(map(lambda x: "1 / %s" % x, res)))
            assert F(nom, denom) == sum([F(1, el) for el in res])