Algorithm for simplifying decimal to fractions
I tried writing an algorithm to simplify a decimal to a fraction and realized it wasn't too simple.
Write 0.333333... as 1/3 for example.
Or 0.1666667, which is 1/6.
Surprisingly I looked online and all the codes I found where either too long, or wouldn't work in some cases. What was even more annoying was that they didn't work for recurring decimals. I was wondering however whether there would be a mathematician/programmer here who understands all the involved processes in simplifying a decimal to a fraction. Anyone?
Solution 1:
The algorithm that the other people have given you gets the answer by calculating the Continued Fraction of the number. This gives a fractional sequence which is guaranteed to converge very, very rapidly. However it is not guaranteed to give you the smallest fraction that is within a distance epsilon of a real number. To find that you have to walk the Stern-Brocot tree.
To do that you subtract off the floor to get the number in the range [0, 1), then your lower estimate is 0, and your upper estimate is 1. Now do a binary search until you are close enough. At each iteration if your lower is a/b and your upper is c/d your middle is (a+c)/(b+d). Test your middle against x, and either make the middle the upper, the lower, or return your final answer.
Here is some very non-idiomatic (and hence, hopefully, readable even if you don't know the language) Python that implements this algorithm.
def float_to_fraction (x, error=0.000001):
n = int(math.floor(x))
x -= n
if x < error:
return (n, 1)
elif 1 - error < x:
return (n+1, 1)
# The lower fraction is 0/1
lower_n = 0
lower_d = 1
# The upper fraction is 1/1
upper_n = 1
upper_d = 1
while True:
# The middle fraction is (lower_n + upper_n) / (lower_d + upper_d)
middle_n = lower_n + upper_n
middle_d = lower_d + upper_d
# If x + error < middle
if middle_d * (x + error) < middle_n:
# middle is our new upper
upper_n = middle_n
upper_d = middle_d
# Else If middle < x - error
elif middle_n < (x - error) * middle_d:
# middle is our new lower
lower_n = middle_n
lower_d = middle_d
# Else middle is our best fraction
else:
return (n * middle_d + middle_n, middle_d)
Solution 2:
(code improved Feb 2017 - scroll down to 'optimization'...)
(algorithm comparison table at the end of this answer)
I implemented btilly's answer in C# and...
- added support for negative numbers
- provide an
accuracyparameter to specify the max. relative error, not the max. absolute error;0.01would find a fraction within 1% of the value. - provide an optimization
-
Double.NaNandDouble.Infinityare not supported; you might want to handle those (example here).
public Fraction RealToFraction(double value, double accuracy)
{
if (accuracy <= 0.0 || accuracy >= 1.0)
{
throw new ArgumentOutOfRangeException("accuracy", "Must be > 0 and < 1.");
}
int sign = Math.Sign(value);
if (sign == -1)
{
value = Math.Abs(value);
}
// Accuracy is the maximum relative error; convert to absolute maxError
double maxError = sign == 0 ? accuracy : value * accuracy;
int n = (int) Math.Floor(value);
value -= n;
if (value < maxError)
{
return new Fraction(sign * n, 1);
}
if (1 - maxError < value)
{
return new Fraction(sign * (n + 1), 1);
}
// The lower fraction is 0/1
int lower_n = 0;
int lower_d = 1;
// The upper fraction is 1/1
int upper_n = 1;
int upper_d = 1;
while (true)
{
// The middle fraction is (lower_n + upper_n) / (lower_d + upper_d)
int middle_n = lower_n + upper_n;
int middle_d = lower_d + upper_d;
if (middle_d * (value + maxError) < middle_n)
{
// real + error < middle : middle is our new upper
upper_n = middle_n;
upper_d = middle_d;
}
else if (middle_n < (value - maxError) * middle_d)
{
// middle < real - error : middle is our new lower
lower_n = middle_n;
lower_d = middle_d;
}
else
{
// Middle is our best fraction
return new Fraction((n * middle_d + middle_n) * sign, middle_d);
}
}
}
The Fraction type is just a simple struct. Of course, use your own preferred type... (I like this one by Rick Davin.)
public struct Fraction
{
public Fraction(int n, int d)
{
N = n;
D = d;
}
public int N { get; private set; }
public int D { get; private set; }
}
Feb 2017 optimization
For certain values, like 0.01, 0.001, etc. the algorithm goes through hundreds or thousands of linear iterations. To fix this, I implemented a binary way of finding the final value -- thanks to btilly for this idea. Inside the if-statement substitute the following:
// real + error < middle : middle is our new upper
Seek(ref upper_n, ref upper_d, lower_n, lower_d, (un, ud) => (lower_d + ud) * (value + maxError) < (lower_n + un));
and
// middle < real - error : middle is our new lower
Seek(ref lower_n, ref lower_d, upper_n, upper_d, (ln, ld) => (ln + upper_n) < (value - maxError) * (ld + upper_d));
Here is the Seek method implementation:
/// <summary>
/// Binary seek for the value where f() becomes false.
/// </summary>
void Seek(ref int a, ref int b, int ainc, int binc, Func<int, int, bool> f)
{
a += ainc;
b += binc;
if (f(a, b))
{
int weight = 1;
do
{
weight *= 2;
a += ainc * weight;
b += binc * weight;
}
while (f(a, b));
do
{
weight /= 2;
int adec = ainc * weight;
int bdec = binc * weight;
if (!f(a - adec, b - bdec))
{
a -= adec;
b -= bdec;
}
}
while (weight > 1);
}
}
Algorithm comparison table
You may want to copy the table to your text editor for full screen viewing.
Accuracy: 1.0E-3 | Stern-Brocot OPTIMIZED | Eppstein | Richards
Input | Result Error Iterations Iterations | Result Error Iterations | Result Error Iterations
======================| =====================================================| =========================================| =========================================
0 | 0/1 (zero) 0 0 0 | 0/1 (zero) 0 0 | 0/1 (zero) 0 0
1 | 1/1 0 0 0 | 1001/1000 1.0E-3 1 | 1/1 0 0
3 | 3/1 0 0 0 | 1003/334 1.0E-3 1 | 3/1 0 0
-1 | -1/1 0 0 0 | -1001/1000 1.0E-3 1 | -1/1 0 0
-3 | -3/1 0 0 0 | -1003/334 1.0E-3 1 | -3/1 0 0
0.999999 | 1/1 1.0E-6 0 0 | 1000/1001 -1.0E-3 2 | 1/1 1.0E-6 0
-0.999999 | -1/1 1.0E-6 0 0 | -1000/1001 -1.0E-3 2 | -1/1 1.0E-6 0
1.000001 | 1/1 -1.0E-6 0 0 | 1001/1000 1.0E-3 1 | 1/1 -1.0E-6 0
-1.000001 | -1/1 -1.0E-6 0 0 | -1001/1000 1.0E-3 1 | -1/1 -1.0E-6 0
0.50 (1/2) | 1/2 0 1 1 | 999/1999 -5.0E-4 2 | 1/2 0 1
0.33... (1/3) | 1/3 0 2 2 | 999/2998 -3.3E-4 2 | 1/3 0 1
0.67... (2/3) | 2/3 0 2 2 | 999/1498 3.3E-4 3 | 2/3 0 2
0.25 (1/4) | 1/4 0 3 3 | 999/3997 -2.5E-4 2 | 1/4 0 1
0.11... (1/9) | 1/9 0 8 4 | 999/8992 -1.1E-4 2 | 1/9 0 1
0.09... (1/11) | 1/11 0 10 5 | 999/10990 -9.1E-5 2 | 1/11 0 1
0.62... (307/499) | 8/13 2.5E-4 5 5 | 913/1484 -2.2E-6 8 | 8/13 2.5E-4 5
0.14... (33/229) | 15/104 8.7E-4 20 9 | 974/6759 -4.5E-6 6 | 16/111 2.7E-4 3
0.05... (33/683) | 7/145 -8.4E-4 24 10 | 980/20283 1.5E-6 7 | 10/207 -1.5E-4 4
0.18... (100/541) | 17/92 -3.3E-4 11 10 | 939/5080 -2.0E-6 8 | 17/92 -3.3E-4 4
0.06... (33/541) | 5/82 -3.7E-4 19 8 | 995/16312 -1.9E-6 6 | 5/82 -3.7E-4 4
0.1 | 1/10 0 9 5 | 999/9991 -1.0E-4 2 | 1/10 0 1
0.2 | 1/5 0 4 3 | 999/4996 -2.0E-4 2 | 1/5 0 1
0.3 | 3/10 0 5 5 | 998/3327 -1.0E-4 4 | 3/10 0 3
0.4 | 2/5 0 3 3 | 999/2497 2.0E-4 3 | 2/5 0 2
0.5 | 1/2 0 1 1 | 999/1999 -5.0E-4 2 | 1/2 0 1
0.6 | 3/5 0 3 3 | 1000/1667 -2.0E-4 4 | 3/5 0 3
0.7 | 7/10 0 5 5 | 996/1423 -1.0E-4 4 | 7/10 0 3
0.8 | 4/5 0 4 3 | 997/1246 2.0E-4 3 | 4/5 0 2
0.9 | 9/10 0 9 5 | 998/1109 -1.0E-4 4 | 9/10 0 3
0.01 | 1/100 0 99 8 | 999/99901 -1.0E-5 2 | 1/100 0 1
0.001 | 1/1000 0 999 11 | 999/999001 -1.0E-6 2 | 1/1000 0 1
0.0001 | 1/9991 9.0E-4 9990 15 | 999/9990001 -1.0E-7 2 | 1/10000 0 1
1E-05 | 1/99901 9.9E-4 99900 18 | 1000/99999999 1.0E-8 3 | 1/99999 1.0E-5 1
0.33333333333 | 1/3 1.0E-11 2 2 | 1000/3001 -3.3E-4 2 | 1/3 1.0E-11 1
0.3 | 3/10 0 5 5 | 998/3327 -1.0E-4 4 | 3/10 0 3
0.33 | 30/91 -1.0E-3 32 8 | 991/3003 1.0E-5 3 | 33/100 0 2
0.333 | 167/502 -9.9E-4 169 11 | 1000/3003 1.0E-6 3 | 333/1000 0 2
0.7777 | 7/9 1.0E-4 5 4 | 997/1282 -1.1E-5 4 | 7/9 1.0E-4 3
0.101 | 10/99 1.0E-4 18 10 | 919/9099 1.1E-6 5 | 10/99 1.0E-4 3
0.10001 | 1/10 -1.0E-4 9 5 | 1/10 -1.0E-4 4 | 1/10 -1.0E-4 2
0.100000001 | 1/10 -1.0E-8 9 5 | 1000/9999 1.0E-4 3 | 1/10 -1.0E-8 2
0.001001 | 1/999 1.0E-6 998 11 | 1/999 1.0E-6 3 | 1/999 1.0E-6 1
0.0010000001 | 1/1000 -1.0E-7 999 11 | 1000/999999 9.0E-7 3 | 1/1000 -1.0E-7 2
0.11 | 10/91 -1.0E-3 18 9 | 1000/9091 -1.0E-5 4 | 10/91 -1.0E-3 2
0.1111 | 1/9 1.0E-4 8 4 | 1000/9001 -1.1E-5 2 | 1/9 1.0E-4 1
0.111111111111 | 1/9 1.0E-12 8 4 | 1000/9001 -1.1E-4 2 | 1/9 1.0E-12 1
1 | 1/1 0 0 0 | 1001/1000 1.0E-3 1 | 1/1 0 0
-1 | -1/1 0 0 0 | -1001/1000 1.0E-3 1 | -1/1 0 0
-0.5 | -1/2 0 1 1 | -999/1999 -5.0E-4 2 | -1/2 0 1
3.14 | 22/7 9.1E-4 6 4 | 964/307 2.1E-5 3 | 22/7 9.1E-4 1
3.1416 | 22/7 4.0E-4 6 4 | 732/233 9.8E-6 3 | 22/7 4.0E-4 1
3.14... (pi) | 22/7 4.0E-4 6 4 | 688/219 -1.3E-5 4 | 22/7 4.0E-4 1
0.14 | 7/50 0 13 7 | 995/7107 2.0E-5 3 | 7/50 0 2
0.1416 | 15/106 -6.4E-4 21 8 | 869/6137 9.2E-7 5 | 16/113 -5.0E-5 2
2.72... (e) | 68/25 6.3E-4 7 7 | 878/323 -5.7E-6 8 | 87/32 1.7E-4 5
0.141592653589793 | 15/106 -5.9E-4 21 8 | 991/6999 -7.0E-6 4 | 15/106 -5.9E-4 2
-1.33333333333333 | -4/3 2.5E-15 2 2 | -1001/751 -3.3E-4 2 | -4/3 2.5E-15 1
-1.3 | -13/10 0 5 5 | -992/763 1.0E-4 3 | -13/10 0 2
-1.33 | -97/73 -9.3E-4 26 8 | -935/703 1.1E-5 3 | -133/100 0 2
-1.333 | -4/3 2.5E-4 2 2 | -1001/751 -8.3E-5 2 | -4/3 2.5E-4 1
-1.33333337 | -4/3 -2.7E-8 2 2 | -999/749 3.3E-4 3 | -4/3 -2.7E-8 2
-1.7 | -17/10 0 5 5 | -991/583 -1.0E-4 4 | -17/10 0 3
-1.37 | -37/27 2.7E-4 7 7 | -996/727 1.0E-5 7 | -37/27 2.7E-4 5
-1.33337 | -4/3 -2.7E-5 2 2 | -999/749 3.1E-4 3 | -4/3 -2.7E-5 2
0.047619 | 1/21 1.0E-6 20 6 | 1000/21001 -4.7E-5 2 | 1/21 1.0E-6 1
12.125 | 97/8 0 7 4 | 982/81 -1.3E-4 2 | 97/8 0 1
5.5 | 11/2 0 1 1 | 995/181 -5.0E-4 2 | 11/2 0 1
0.1233333333333 | 9/73 -3.7E-4 16 8 | 971/7873 -3.4E-6 4 | 9/73 -3.7E-4 2
0.7454545454545 | 38/51 -4.8E-4 15 8 | 981/1316 -1.9E-5 6 | 38/51 -4.8E-4 4
0.01024801004 | 2/195 8.2E-4 98 9 | 488/47619 2.0E-8 13 | 2/195 8.2E-4 3
0.99011 | 91/92 -9.9E-4 91 8 | 801/809 1.3E-6 5 | 100/101 -1.1E-5 2
0.9901134545 | 91/92 -9.9E-4 91 8 | 601/607 1.9E-6 5 | 100/101 -1.5E-5 2
0.19999999 | 1/5 5.0E-8 4 3 | 1000/5001 -2.0E-4 2 | 1/5 5.0E-8 1
0.20000001 | 1/5 -5.0E-8 4 3 | 1000/4999 2.0E-4 3 | 1/5 -5.0E-8 2
5.0183168565E-05 | 1/19908 9.5E-4 19907 16 | 1000/19927001 -5.0E-8 2 | 1/19927 5.2E-12 1
3.909E-07 | 1/2555644 1.0E-3 2555643 23 | 1/1 2.6E6 (!) 1 | 1/2558199 1.1E-8 1
88900003.001 |88900003/1 -1.1E-11 0 0 |88900004/1 1.1E-8 1 |88900003/1 -1.1E-11 0
0.26... (5/19) | 5/19 0 7 6 | 996/3785 -5.3E-5 4 | 5/19 0 3
0.61... (37/61) | 17/28 9.7E-4 8 7 | 982/1619 -1.7E-5 8 | 17/28 9.7E-4 5
| | |
Accuracy: 1.0E-4 | Stern-Brocot OPTIMIZED | Eppstein | Richards
Input | Result Error Iterations Iterations | Result Error Iterations | Result Error Iterations
======================| =====================================================| =========================================| =========================================
0.62... (307/499) | 227/369 -8.8E-5 33 11 | 9816/15955 -2.0E-7 8 | 299/486 -6.7E-6 6
0.05... (33/683) | 23/476 6.4E-5 27 12 | 9989/206742 1.5E-7 7 | 23/476 6.4E-5 5
0.06... (33/541) | 28/459 6.6E-5 24 12 | 9971/163464 -1.9E-7 6 | 33/541 0 5
1E-05 | 1/99991 9.0E-5 99990 18 | 10000/999999999 1.0E-9 3 | 1/99999 1.0E-5 1
0.333 | 303/910 -9.9E-5 305 12 | 9991/30003 1.0E-7 3 | 333/1000 0 2
0.7777 | 556/715 -1.0E-4 84 12 | 7777/10000 0 8 | 1109/1426 -1.8E-7 4
3.14... (pi) | 289/92 -9.2E-5 19 8 | 9918/3157 -8.1E-7 4 | 333/106 -2.6E-5 2
2.72... (e) | 193/71 1.0E-5 10 9 | 9620/3539 6.3E-8 11 | 193/71 1.0E-5 7
0.7454545454545 | 41/55 6.1E-14 16 8 | 9960/13361 -1.8E-6 6 | 41/55 6.1E-14 5
0.01024801004 | 7/683 8.7E-5 101 12 | 9253/902907 -1.3E-10 16 | 7/683 8.7E-5 5
0.99011 | 100/101 -1.1E-5 100 8 | 901/910 -1.1E-7 6 | 100/101 -1.1E-5 2
0.9901134545 | 100/101 -1.5E-5 100 8 | 8813/8901 1.6E-8 7 | 100/101 -1.5E-5 2
0.26... (5/19) | 5/19 0 7 6 | 9996/37985 -5.3E-6 4 | 5/19 0 3
0.61... (37/61) | 37/61 0 10 8 | 9973/16442 -1.6E-6 8 | 37/61 0 7
Performance comparison
I performed detailed speed tests and plotted the results. Not looking at quality and only speed:
- The Stern-Brocot optimization slows it down by at most a factor 2, but the original Stern-Brocot can be hundreds or thousands times slower when it hits the unlucky values mentioned. That's still only a couple of microseconds though per call.
- Richards is consistently fast.
- Eppstein is around 3 times slower than the others.
Stern-Brocot and Richards compared:
- Both return nice fractions.
- Richards often results in a smaller error. It is also a bit faster.
- Stern-Brocot walks down the S-B tree. It finds the fraction of the lowest denominator that meets the required accuracy, then stops.
If you do not require the lowest denominator fraction, Richards is a good choice.
Solution 3:
I know you said you searched online, but if you missed the following paper it might be of some help. It includes a code example in Pascal.
Algorithm To Convert A Decimal To A Fraction*
Alternatively, as part of it's standard library, Ruby has code that deals with rational numbers. It can convert from floats to rationals and vice versa. I believe you can look through the code as well. The documentation is found here. I know you're not using Ruby, but it might help to look at the algorithms.
Additionally, you can call Ruby code from C# (or even write Ruby code inside a C# code file) if you use IronRuby, which runs on top of the .net framework.
*Updated to a new link as it appears the original URL is broken (http://homepage.smc.edu/kennedy_john/DEC2FRAC.pdf)