[eesmith]

I don't think it's that good. At the very least I implemented it and got different answers for the test case of 0.263157894737 = 5/19 depending on the choice of numerical representation.

float64 ended up at 87714233939/333314088968 while decimal and fraction ended up at 87719298244/333333333327.

Here is the code.

    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument("Q", nargs="?", default="0.263157894737")
    parser.add_argument("--method", choices = ("float", "decimal", "fraction"), default="float")
    args = parser.parse_args()

    if args.method == "float":
        print(f"Using float64 for {args.Q!r}")
        one = 1.0
        Q = float(args.Q)
        int2float = float
    elif args.method == "decimal":
        print(f"Using decimal for {args.Q!r}")
        import decimal
        one = decimal.Decimal(1)
        Q = decimal.Decimal(args.Q)
        int2float = decimal.Decimal
    elif args.method == "fraction":
        print(f"Using fractions for {args.Q!r}")
        import fractions
        one = fractions.Fraction(1)
        Q = fractions.Fraction(args.Q)
        int2float = fractions.Fraction

    def to_frac(X):
        Da = 0
        Db = 1
        Za = X
        while 1:
            Zb = one / (Za - int(Za))
            Dc = Db * int(Zb) + Da
            N = round(X * Dc)
            frac = N / int2float(Dc)
            print(f"  {N}/{Dc} = {frac} (diff: {X-frac})")
            if float(N) / float(Dc) == float(X):
                return (N, Dc)
            Da, Db = Db, Dc
            Za = Zb

    print("solution:", to_frac(Q))

Here's the output for 0.263157894737

  % python frac.py 0.263157894737 --method float
  Using float64 for '0.263157894737'
    1/3 = 0.3333333333333333 (diff: -0.0701754385963333)
    1/4 = 0.25 (diff: 0.01315789473700002)
    4/15 = 0.26666666666666666 (diff: -0.003508771929666643)
    5/19 = 0.2631578947368421 (diff: 1.5792922525292852e-13)
    87714233939/333314088968 = 0.263157894737 (diff: 0.0)
  solution: (87714233939, 333314088968)

  % python frac.py 0.263157894737 --method decimal
  Using decimal for '0.263157894737'
    1/3 = 0.3333333333333333333333333333 (diff: -0.0701754385963333333333333333)
    1/4 = 0.25 (diff: 0.013157894737)
    4/15 = 0.2666666666666666666666666667 (diff: -0.0035087719296666666666666667)
    5/19 = 0.2631578947368421052631578947 (diff: 1.578947368421053E-13)
    87719298244/333333333327 = 0.2631578947370000000000030000 (diff: -3.0000E-24)
  solution: (87719298244, 333333333327)

  % python frac.py 0.263157894737 --method fraction
  Using fractions for '0.263157894737'
    1/3 = 1/3 (diff: -210526315789/3000000000000)
    1/4 = 1/4 (diff: 13157894737/1000000000000)
    4/15 = 4/15 (diff: -10526315789/3000000000000)
    5/19 = 5/19 (diff: 3/19000000000000)
    87719298244/333333333327 = 87719298244/333333333327 (diff: -1/333333333327000000000000)
  solution: (87719298244, 333333333327)

For the pi = 3.14159265359 test case the solutions are all 226883371/72219220 . The sequence diverges at the point marked with "<---- here".

  Using float64
    22/7 = 3.142857142857143 (diff: -0.0012644892671427321)
    333/106 = 3.141509433962264 (diff: 8.321962773605307e-05)
    355/113 = 3.1415929203539825 (diff: -2.667639824593948e-07)
    103993/33102 = 3.1415926530119025 (diff: 5.780975698144175e-10)
    104348/33215 = 3.141592653921421 (diff: -3.3142111277584263e-10)
    208341/66317 = 3.1415926534674368 (diff: 1.2256329284809908e-10)
    312689/99532 = 3.1415926536189365 (diff: -2.893640882462023e-11)
    833719/265381 = 3.141592653581078 (diff: 8.922196315097608e-12)
    1146408/364913 = 3.141592653591404 (diff: -1.403765992336048e-12)
    5419351/1725033 = 3.1415926535898153 (diff: 1.8474111129762605e-13) <---- here
    6565759/2089946 = 3.141592653590093 (diff: -9.281464485866309e-14)
    11985110/3814979 = 3.141592653589967 (diff: 3.2862601528904634e-14)
    18550869/5904925 = 3.1415926535900116 (diff: -1.1546319456101628e-14)
    30535979/9719904 = 3.1415926535899943 (diff: 5.773159728050814e-15)
    49086848/15624829 = 3.141592653590001 (diff: -8.881784197001252e-16)
    226883371/72219220 = 3.14159265359 (diff: 0.0)
  solution: (226883371, 72219220)

For 0.10000000000000002 it gives 562949953421313/5629499534213129 (for float64) or 500000000000000/4999999999999999 (using fractions):

  Using float64 for '0.10000000000000002'
    1/9 = 0.1111111111111111 (diff: -0.011111111111111086)
    1/10 = 0.1 (diff: 1.3877787807814457e-17)
    562949953421313/5629499534213129 = 0.10000000000000002 (diff: 0.0)
  solution: (562949953421313, 5629499534213129)

  Using fractions for '0.10000000000000002'
    1/9 = 1/9 (diff: -4999999999999991/450000000000000000)
    1/10 = 1/10 (diff: 1/50000000000000000)
    500000000000000/4999999999999999 = 500000000000000/4999999999999999 (diff: -1/249999999999999950000000000000000)
  solution: (500000000000000, 4999999999999999)

FWIW, the exact solution is:

  >>> import fractions
  >>> fractions.Fraction("0.10000000000000002")
  Fraction(5000000000000001, 50000000000000000)

[dunham]

The pascal code cuts when it hits a certain accuracy (passed in as an argument).

[eesmith]

The main issue I see is that the algorithm - unlike the mediant version - does not generate maximally successive accurate approximations.

Yes, you can stop the algorithm at a certain accuracy, but that doesn't mean you can't get better for that given accuracy.

Consider the pi = 3.14159265359 case, and you want it precise to 1/1,000,000. The float64 algorithm gives 1146408/364913 or 5419351/1725033 because:

     ...
    1146408/364913 = 3.141592653591404 (diff: -1.403765992336048e-12)
      --- want a solution here ---
    5419351/1725033 = 3.1415926535898153 (diff: 1.8474111129762605e-13)
     ...

The mediant method, on the other hand, gives an intermediate solution:

  >>> import fractions
  >>> fractions.Fraction("3.14159265359").limit_denominator(1000000)
  Fraction(3126535, 995207)
  >>> float(fractions.Fraction("3.14159265359").limit_denominator(1000000))
  3.1415926535886505

That's a difference of -1.3495871087343403e-12 which is more accurate than 1146408/364913, and is not a solution found by the other algorithm.

Or, if you want a denominator of 364913 or 1725033 you can do that with mediants:

  >>> fractions.Fraction("3.14159265359").limit_denominator(364913)
  Fraction(1146408, 364913)
  >>> fractions.Fraction(3.14159265359).limit_denominator(1725033)
  Fraction(5419351, 1725033)

Another issue is the numerical range. Consider the input "1.5e-318". It causes overflow in the float, decimal, and fraction implementations I gave:

  % python frac.py 1.5e-318 --method float
  Using float64 for '1.5e-318'
  Traceback (most recent call last):
    File "frac.py", line 40, in <module>
      print("solution:", to_frac(Q))
    File "tmp.py", line 31, in to_frac
      Dc = Db * int(Zb) + Da
  OverflowError: cannot convert float infinity to integer


  % python frac.py 1.5e-318 --method fraction
  Using fractions for '1.5e-318'
   1/666666666... many 6s removed ...6666 = 1/666666666... many 6s removed ...6666
      (diff: -1/666 .. more 6s removed ... 6600.. even more zeros removed ...00)
  Traceback (most recent call last):
    File "frac.py", line 40, in <module>
      print("solution:", to_frac(Q))
    File "frac.py", line 35, in to_frac
      if float(N) / float(Dc) == float(X):
  OverflowError: int too large to convert to float

while with the mediant solution I don't need to worry about the input range beyond infinity and NaN:

  >>> from fractions import Fraction
  >>> Fraction(1.5e-318).limit_denominator(1000000)
  Fraction(0, 1)

(I used the float() calls to ensure the fraction and decimal methods stop at the limits of float64. If I remove them I still end up with numbers like 3/2E319 which would not be representable as a Turbo Pascal integer, while the Turbo Pascal mediant implementation would not have a problem.)

Finally, the mediant solution is easier to implement.