Algorithm for simplifying decimal to fractions
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)
Solution 4
I found the same paper that Matt referenced, and I took a second and implemented it in Python. Maybe seeing the same idea in code will make it clearer. Granted, you requested an answer in C# and I'm giving it to you in Python, but it's a fairly trivial program, and I'm sure it would be easy to translate. The parameters are num (the decimal number you'd like to convert to a rational) and epsilon (the maximum allowed difference between num and the calculated rational). Some quick test runs find that it usually only takes two or three iterations to converge when epsilon is around 1e-4.
def dec2frac(num, epsilon, max_iter=20):
d = [0, 1] + ([0] * max_iter)
z = num
n = 1
t = 1
while num and t < max_iter and abs(n/d[t] - num) > epsilon:
t += 1
z = 1/(z - int(z))
d[t] = d[t-1] * int(z) + d[t-2]
# int(x + 0.5) is equivalent to rounding x.
n = int(num * d[t] + 0.5)
return n, d[t]
Edit: I just noticed your note about wanting them to work with recurring decimals. I don't know any languages that have syntax to support recurring decimals, so I'm not sure how one would go about handling them, but running 0.6666666 and 0.166666 through this method return the correct results (2/3 and 1/6, respectively).
Another Edit (I didn't think this would be so interesting!): If you want to know more about the theory behind this algorithm, Wikipedia has an excellent page on the Euclidian algorithm
Solution 5
The most populair solutions to this problem are Richards’ algorithm and the Stern-Brocot algorithm, implemented by btilly with speed optimalization by btilly and Jay Zed. Richards’ algorithm is the fastest, but does not guarantee to return the best fraction.
I have an solution to this problem which always gives the best fraction and is also faster than all of the algorithms above. Here is the algorithm in C# (explanation and speed test below).
This is a short algorithm without comments. An complete version is provided in the source code at the end.
public static Fraction DoubleToFractionSjaak(double value, double accuracy)
{
int sign = value < 0 ? -1 : 1;
value = value < 0 ? -value : value;
int integerpart = (int)value;
value -= integerpart;
double minimalvalue = value - accuracy;
if (minimalvalue < 0.0) return new Fraction(sign * integerpart, 1);
double maximumvalue = value + accuracy;
if (maximumvalue > 1.0) return new Fraction(sign * (integerpart + 1), 1);
int a = 0;
int b = 1;
int c = 1;
int d = (int)(1 / maximumvalue);
while (true)
{
int n = (int)((b * minimalvalue - a) / (c - d * minimalvalue));
if (n == 0) break;
a += n * c;
b += n * d;
n = (int)((c - d * maximumvalue) / (b * maximumvalue - a));
if (n == 0) break;
c += n * a;
d += n * b;
}
int denominator = b + d;
return new Fraction(sign * (integerpart * denominator + (a + c)), denominator);
}
Where Fraction is a simple class to store a fraction, like the following:
public class Fraction
{
public int Numerator { get; private set; }
public int Denominator { get; private set; }
public Fraction(int numerator, int denominator)
{
Numerator = numerator;
Denominator = denominator;
}
}
How it works
Like the other solutions mentioned, my solution is based on continued fraction. Other solutions like the one from Eppstein or solutions based on repeating decimals proved to be slower and/or give suboptimal results.
Continued fraction
Solutions based on continued fraction are mostly based on two algorithms, both described in an article by Ian Richards published here in 1981. He called them the “slow continued fraction algorithm” and the “fast continued fraction algorithm”. The first is known as the the Stern-Brocot algorithm while the latter is known as Richards’ algorithm.
My algorithm (short explanation)
To fully understand my algorithm, you need to have read the article by Ian Richards or at least understand what a Farey pair is. Furthermore, read the algorithm with comments at the end of this article.
The algorithm is using a Farey pair, containing a left and a right fraction. By repeatedly taking the mediant it is closing in on the target value. This is just like the slow algorithm but there are two major differences:
- Multiple iterations are performed at once as long as the mediant stay on one side of the target value.
- The left and right fraction cannot come closer to the target value than the given accuracy.
Alternately the right and left side of the target value are checked. If the algorithm cannot produce a result closer to the target value, the process ends. The resulting mediant is the optimal solution.
Speed test
I did some speed tests on my laptop with the following algorithms:
- Improved slow algorithm by Kay Zed and btilly
- John Kennedy’s implementation of the Fast algorithm, converted to C# by Kay Zed
- My implementation of the Fast algorithm (close to the original by Ian Richards)
- Jeremy Herrman’s implementation of the Fast algorithm
- My algorithm above
I omitted the original slow algorithm by btilly, because of its bad worst-case performance.
Test set
I choose a set of target values (very arbitrary) and calculated the fraction 100000 times with 5 different accuracies. Because possible some (future) algorithms couldn't handle improper fractions, only target values from 0.0 to 1.0 were tested. The accuracy was taken from the range from 2 to 6 decimal places (0.005 to 0.0000005). The following set was used:
0.999999, 0.000001, 0.25
0.33, 0.333, 0.3333, 0.33333, 0.333333, 0.333333333333,
0.666666666666, 0.777777777777, 0.090909090909, 0.263157894737,
0.606557377049, 0.745454545454, 0.000050183168565,
pi - 3, e - 2.0, sqrt(2) - 1
Results
I did 13 test runs. The result is in milliseconds needed for the whole data set.
Run 1 Run 2 Run 3 Run 4 Run 5 Run 6 Run 7 Run 8 Run 9 Run 10 Run 11 Run 12 Run 13
1. 9091 9222 9070 9111 9091 9108 9293 9118 9115 9113 9102 9143 9121
2. 7071 7125 7077 6987 7126 6985 7037 6964 7023 6980 7053 7050 6999
3. 6903 7059 7062 6891 6942 6880 6882 6918 6853 6918 6893 6993 6966
4. 7546 7554 7564 7504 7483 7529 7510 7512 7517 7719 7513 7520 7514
5. 6839 6951 6882 6836 6854 6880 6846 7017 6874 6867 6828 6848 6864
Conclusion (skipping the analysis)
Even without a statistical analysis, it's easy to see that my algorithm is faster than the other tested algorithms. The difference with the fastest variant of “fast algorithm” however is less than 1 percent. The Improved slow algorithm is 30%-35% slower than the fastest algorithm”.
On the other hand, even the slowest algorithm performs a calculation on average in less than a microsecond. So under normal circumstances speed is not really an issue. In my opinion the best algorithm is mainly a matter of taste, so choose any of the tested algorithms on other criteria.
- Does the algorithm gives the best result?
- Is the algorithm available in my favorite language?
- What is the code size of the algorithm?
- Is the algorithm readable, understandable?
Source code
The source code below contains all used algorithms. It includes:
- My original algorithm (with comments)
- A even faster version of my algorithm (but less readable)
- The original slow algorithm
- All tested algorithms
public class DoubleToFraction
{
// ===================================================
// Sjaak algorithm - original version
//
public static Fraction SjaakOriginal(double value, double accuracy)
{
// Split value in a sign, an integer part, a fractional part
int sign = value < 0 ? -1 : 1;
value = value < 0 ? -value : value;
int integerpart = (int)value;
value -= integerpart;
// check if the fractional part is near 0
double minimalvalue = value - accuracy;
if (minimalvalue < 0.0) return new Fraction(sign * integerpart, 1);
// check if the fractional part is near 1
double maximumvalue = value + accuracy;
if (maximumvalue > 1.0) return new Fraction(sign * (integerpart + 1), 1);
// The left fraction (a/b) is initially (0/1), the right fraction (c/d) is initially (1/1)
// Together they form a Farey pair.
// We will keep the left fraction below the minimumvalue and the right fraction above the maximumvalue
int a = 0;
int b = 1;
int c = 1;
int d = (int)(1 / maximumvalue);
// The first interation is performed above. Calculate maximum n where (n*a+c)/(n*b+d) >= maximumvalue
// This is the same as n <= 1/maximumvalue - 1, d will become n+1 = floor(1/maximumvalue)
// repeat forever (at least until we cannot close in anymore)
while (true)
{
// Close in from the left n times.
// Calculate maximum n where (a+n*c)/(b+n*d) <= minimalvalue
// This is the same as n <= (b * minimalvalue - a) / (c-d*minimalvalue)
int n = (int)((b * minimalvalue - a) / (c - d * minimalvalue));
// If we cannot close in from the left (and also not from the right anymore) the loop ends
if (n == 0) break;
// Update left fraction
a += n * c;
b += n * d;
// Close in from the right n times.
// Calculate maximum n where (n*a+c)/(n*b+d) >= maximumvalue
// This is the same as n <= (c - d * maximumvalue) / (b * maximumvalue - a)
n = (int)((c - d * maximumvalue) / (b * maximumvalue - a));
// If we cannot close in from the right (and also not from the left anymore) the loop ends
if (n == 0) break;
// Update right fraction
c += n * a;
d += n * b;
}
// We cannot close in anymore
// The best fraction will be the mediant of the left and right fraction = (a+c)/(b+d)
int denominator = b + d;
return new Fraction(sign * (integerpart * denominator + (a + c)), denominator);
}
// ===================================================
// Sjaak algorithm - faster version
//
public static Fraction SjaakFaster(double value, double accuracy)
{
int sign = value < 0 ? -1 : 1;
value = value < 0 ? -value : value;
int integerpart = (int)value;
value -= integerpart;
double minimalvalue = value - accuracy;
if (minimalvalue < 0.0) return new Fraction(sign * integerpart, 1);
double maximumvalue = value + accuracy;
if (maximumvalue > 1.0) return new Fraction(sign * (integerpart + 1), 1);
//int a = 0;
int b = 1;
//int c = 1;
int d = (int)(1 / maximumvalue);
double left_n = minimalvalue; // b * minimalvalue - a
double left_d = 1.0 - d * minimalvalue; // c - d * minimalvalue
double right_n = 1.0 - d * maximumvalue; // c - d * maximumvalue
double right_d = maximumvalue; // b * maximumvalue - a
while (true)
{
if (left_n < left_d) break;
int n = (int)(left_n / left_d);
//a += n * c;
b += n * d;
left_n -= n * left_d;
right_d -= n * right_n;
if (right_n < right_d) break;
n = (int)(right_n / right_d);
//c += n * a;
d += n * b;
left_d -= n * left_n;
right_n -= n * right_d;
}
int denominator = b + d;
int numerator = (int)(value * denominator + 0.5);
return new Fraction(sign * (integerpart * denominator + numerator), denominator);
}
// ===================================================
// Original Farley - Implemented by btilly
//
public static Fraction OriginalFarley(double value, double accuracy)
{
// Split value in a sign, an integer part, a fractional part
int sign = value < 0 ? -1 : 1;
value = value < 0 ? -value : value;
int integerpart = (int)value;
value -= integerpart;
// check if the fractional part is near 0
double minimalvalue = value - accuracy;
if (minimalvalue < 0.0) return new Fraction(sign * integerpart, 1);
// check if the fractional part is near 1
double maximumvalue = value + accuracy;
if (maximumvalue > 1.0) return new Fraction(sign * (integerpart + 1), 1);
// The lower fraction is 0/1
int lower_numerator = 0;
int lower_denominator = 1;
// The upper fraction is 1/1
int upper_numerator = 1;
int upper_denominator = 1;
while (true)
{
// The middle fraction is (lower_numerator + upper_numerator) / (lower_denominator + upper_denominator)
int middle_numerator = lower_numerator + upper_numerator;
int middle_denominator = lower_denominator + upper_denominator;
if (middle_denominator * maximumvalue < middle_numerator)
{
// real + error < middle : middle is our new upper
upper_numerator = middle_numerator;
upper_denominator = middle_denominator;
}
else if (middle_numerator < minimalvalue * middle_denominator)
{
// middle < real - error : middle is our new lower
lower_numerator = middle_numerator;
lower_denominator = middle_denominator;
}
else
{
return new Fraction(sign * (integerpart * middle_denominator + middle_numerator), middle_denominator);
}
}
}
// ===================================================
// Modified Farley - Implemented by btilly, Kay Zed
//
public static Fraction ModifiedFarley(double value, double accuracy)
{
// Split value in a sign, an integer part, a fractional part
int sign = value < 0 ? -1 : 1;
value = value < 0 ? -value : value;
int integerpart = (int)value;
value -= integerpart;
// check if the fractional part is near 0
double minimalvalue = value - accuracy;
if (minimalvalue < 0.0) return new Fraction(sign * integerpart, 1);
// check if the fractional part is near 1
double maximumvalue = value + accuracy;
if (maximumvalue > 1.0) return new Fraction(sign * (integerpart + 1), 1);
// The lower fraction is 0/1
int lower_numerator = 0;
int lower_denominator = 1;
// The upper fraction is 1/1
int upper_numerator = 1;
int upper_denominator = 1;
while (true)
{
// The middle fraction is (lower_numerator + upper_numerator) / (lower_denominator + upper_denominator)
int middle_numerator = lower_numerator + upper_numerator;
int middle_denominator = lower_denominator + upper_denominator;
if (middle_denominator * maximumvalue < middle_numerator)
{
// real + error < middle : middle is our new upper
ModifiedFarleySeek(ref upper_numerator, ref upper_denominator, lower_numerator, lower_denominator, (un, ud) => (lower_denominator + ud) * maximumvalue < (lower_numerator + un));
}
else if (middle_numerator < minimalvalue * middle_denominator)
{
// middle < real - error : middle is our new lower
ModifiedFarleySeek(ref lower_numerator, ref lower_denominator, upper_numerator, upper_denominator, (ln, ld) => (ln + upper_numerator) < minimalvalue * (ld + upper_denominator));
}
else
{
return new Fraction(sign * (integerpart * middle_denominator + middle_numerator), middle_denominator);
}
}
}
private static void ModifiedFarleySeek(ref int a, ref int b, int ainc, int binc, Func<int, int, bool> f)
{
// Binary seek for the value where f() becomes false
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);
}
}
// ===================================================
// Richards implementation by Jemery Hermann
//
public static Fraction RichardsJemeryHermann(double value, double accuracy, int maxIterations = 20)
{
// Split value in a sign, an integer part, a fractional part
int sign = value < 0 ? -1 : 1;
value = value < 0 ? -value : value;
int integerpart = (int)value;
value -= integerpart;
// check if the fractional part is near 0
double minimalvalue = value - accuracy;
if (minimalvalue < 0.0) return new Fraction(sign * integerpart, 1);
// check if the fractional part is near 1
double maximumvalue = value + accuracy;
if (maximumvalue > 1.0) return new Fraction(sign * (integerpart + 1), 1);
// Richards - Implemented by Jemery Hermann
double[] d = new double[maxIterations + 2];
d[1] = 1;
double z = value;
double n = 1;
int t = 1;
while (t < maxIterations && Math.Abs(n / d[t] - value) > accuracy)
{
t++;
z = 1 / (z - (int)z);
d[t] = d[t - 1] * (int)z + d[t - 2];
n = (int)(value * d[t] + 0.5);
}
return new Fraction(sign * (integerpart * (int)d[t] + (int)n), (int)d[t]);
}
// ===================================================
// Richards implementation by Kennedy
//
public static Fraction RichardsKennedy(double value, double accuracy)
{
// Split value in a sign, an integer part, a fractional part
int sign = value < 0 ? -1 : 1;
value = value < 0 ? -value : value;
int integerpart = (int)value;
value -= integerpart;
// check if the fractional part is near 0
double minimalvalue = value - accuracy;
if (minimalvalue < 0.0) return new Fraction(sign * integerpart, 1);
// check if the fractional part is near 1
double maximumvalue = value + accuracy;
if (maximumvalue > 1.0) return new Fraction(sign * (integerpart + 1), 1);
// Richards
double z = value;
int previousDenominator = 0;
int denominator = 1;
int numerator;
do
{
z = 1.0 / (z - (int)z);
int temp = denominator;
denominator = denominator * (int)z + previousDenominator;
previousDenominator = temp;
numerator = (int)(value * denominator + 0.5);
}
while (Math.Abs(value - (double)numerator / denominator) > accuracy && z != (int)z);
return new Fraction(sign * (integerpart * denominator + numerator), denominator);
}
// ===================================================
// Richards implementation by Sjaak
//
public static Fraction RichardsOriginal(double value, double accuracy)
{
// Split value in a sign, an integer part, a fractional part
int sign = value < 0 ? -1 : 1;
value = value < 0 ? -value : value;
int integerpart = (int)value;
value -= integerpart;
// check if the fractional part is near 0
double minimalvalue = value - accuracy;
if (minimalvalue < 0.0) return new Fraction(sign * integerpart, 1);
// check if the fractional part is near 1
double maximumvalue = value + accuracy;
if (maximumvalue > 1.0) return new Fraction(sign * (integerpart + 1), 1);
// Richards
double z = value;
int denominator0 = 0;
int denominator1 = 1;
int numerator0 = 1;
int numerator1 = 0;
int n = (int)z;
while (true)
{
z = 1.0 / (z - n);
n = (int)z;
int temp = denominator1;
denominator1 = denominator1 * n + denominator0;
denominator0 = temp;
temp = numerator1;
numerator1 = numerator1 * n + numerator0;
numerator0 = temp;
double d = (double)numerator1 / denominator1;
if (d > minimalvalue && d < maximumvalue) break;
}
return new Fraction(sign * (integerpart * denominator1 + numerator1), denominator1);
}
}
Chibueze Opata
Updated on July 08, 2022Comments
-
Chibueze Opata almost 2 years
I tried writing an algorithm to simplify a decimal to a fraction and realized it wasn't too simple.
Write
0.333333...as1/3for example.Or
0.1666667, which is1/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?