'''
Python routine to find rational approximations to real numbers.  If
you have the mpmath module (get it from http://code.google.com/p/mpmath/),
it will be used and you can work to arbitrary precisions.

Copyright (C) 2009 Don Peterson

---------------------------------------------------------------------------

This program is free software; you can redistribute it
and/or modify it under the terms of the GNU General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be
useful, but WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.  See the GNU General Public License for more
details.

You should have received a copy of the GNU General Public
License along with this program; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330,
Boston, MA 02111-1307  USA

---------------------------------------------------------------------------

Translated from C code at:
    NIST Guide to Available Math Software.
    Source for module FRAC from package C.
    Retrieved from NETLIB on Tue Nov 11 18:08:26 1997.
 
/***************************************************************************
Given a number, v, this function outputs two integers, d and n, such that
 
    v = n / d
 
to accuracy
 
    epsilon = | (v - n/d) / v | <= error
 
input:  v     = decimal number you want replaced by fraction.
        error = accuracy to which the fraction should
                represent the number v.
 
output: n     = numerator of representing fraction.
        d     = denominator of representing fraction.
 
return value:  -1.0 if (v < MIN || v > MAX || error < 0.0)
             | (v - n/d) / v | otherwise.
 
Note:  This program only works for positive numbers, v.
 
reference:  Jerome Spanier and Keith B. Oldham, "An Atlas
    of Functions," Springer-Verlag, 1987, pp. 665-7.
***************************************************************************/
'''

from __future__ import division
import sys

try:
    from mpmath import *
    Dec = mpf
    mp.dps = 20
    def about(x):
        return nstr(x, 3)
except:
    from math import *
    Dec = float
    def about(x):
        return "%.2e" % x

def frac(x, error, max_iterations=500):
        sign = 1
        if x < 0:
            sign = -1
            x = abs(x)
        if x == 0:
            return 0, 1, 0
        d, D, n, r, epsilon, iterations = 1, 1, int(x), Dec(1), Dec(0), 0
        N = n + 1
        one = True
        while True:
            iterations += 1
            if iterations > max_iterations:
                raise Exception("Too many iterations")
            if not one:
                if r > 1.0:
                    N += n*int(r)
                    D += d*int(r)
                    n += N
                    d += D
                else:
                    r = 1.0/r
            else:
                one = False
            r = 0.0
            if x*d != n:
                r = (N - x*D)/(x*d - n)
                if r <= 1.0:
                    t = N
                    N = n
                    n = t
                    t = D
                    D = d
                    d = t
            epsilon = abs(1.0 - n/(x*d))
            if epsilon <= error:
                return sign*n, d, epsilon
            if r != 0.0:
                continue
            m = 10.0
            while m*epsilon < 1.0:
                m *= 10.0
            epsilon = (1.0/m)*(int(0.5 + m*epsilon))
            if r == 0.0:
                raise Exception("Error in frac() routine")

if __name__ == "__main__":
    if len(sys.argv) < 2:
        sys.stderr.write('''Usage:  %s decimal_number [n]
  Prints out a rational approximation for the decimal number that is 
  within 10^(-n) of the number.  If n is not given, it defaults to 8.

  The python math library is in scope, so decimal_number can also be
  expressions like 'pi - sqrt(3)/cos(2)'.  Example:  the golden ratio
  '(1 + sqrt(5))/2' is approximated by:
      10946/6765 error = 9.7719085516476242217e-9
''' % sys.argv[0])
        exit(1)
    try:
        x = Dec(sys.argv[1])
    except:
        try:
            x = Dec(eval(sys.argv[1]))
        except:
            sys.stderr.write("Don't know how to evaluate '%s'" % sys.argv[1])
            exit(1)
    N = 8
    if len(sys.argv) > 2:
        N = int(sys.argv[2])
    n, d, eps = frac(x, 1/Dec("10")**N)
    print "%d/%d" % (n, d), "error =", about(abs(x - Dec(n)/d))