Skip to content

Instantly share code, notes, and snippets.

@mikeando
Last active February 4, 2023 01:02
Show Gist options
  • Save mikeando/7073d62385a34a61a6f7 to your computer and use it in GitHub Desktop.
Save mikeando/7073d62385a34a61a6f7 to your computer and use it in GitHub Desktop.
Fraction approximation in C++

Here's a version based on the answer at http://stackoverflow.com/questions/12098461/how-can-i-detect-if-a-float-has-a-repeating-decimal-expansion-in-c/12101996#12101996

It's hard coded to look for an approximation to 1.2734, and 10 iterations - its output is :

    0 1/1 = 1.000000  (eps = -0.273400)
    1 4/3 = 1.333333  (eps = 0.059933)
    2 5/4 = 1.250000  (eps = -0.023400)
    3 9/7 = 1.285714  (eps = 0.012314)
    4 14/11 = 1.272727  (eps = -0.000673)
    5 163/128 = 1.273438  (eps = 0.000038)
    6 177/139 = 1.273381  (eps = -0.000019)
    7 340/267 = 1.273408  (eps = 0.000008)
    8 517/406 = 1.273399  (eps = -0.000001)
    9 2925/2297 = 1.273400  (eps = 0.000000)

Or changing the value to 1.17171

0 1/1 = 1.000000  (eps = -0.171710)
1 6/5 = 1.200000  (eps = 0.028290)
2 7/6 = 1.166667  (eps = -0.005043)
3 34/29 = 1.172414  (eps = 0.000704)
4 41/35 = 1.171429  (eps = -0.000281)
5 116/99 = 1.171717  (eps = 0.000007)
6 1549/1322 = 1.171710  (eps = -0.000000)
7 1665/1421 = 1.171710  (eps = 0.000000)
8 18199/15532 = 1.171710  (eps = 0.000000)
9 19864/16953 = 1.171710  (eps = 0.000000)

main2.cpp contains a comparison of continued fraction method of main.cpp with an alternate algorithm based on searching the Stern-Brocot tree, adapted to C++ from http://stackoverflow.com/a/5128558/221955 .

The Stern Brocot method tends to find more approximations, and smaller denominators than the Continued Fraction approach but is computationally more expensive - but both are very fast. Here's the output at various tolerance levels.

tol=0.500000 SB: 1/1 = 1.000000 (err=-0.171710)    CF: 1/1 = 1.000000 (err=-0.171710)
tol=0.250000 SB: 1/1 = 1.000000 (err=-0.171710)    CF: 1/1 = 1.000000 (err=-0.171710)
tol=0.125000 SB: 5/4 = 1.250000 (err=0.078290)    CF: 6/5 = 1.200000 (err=0.028290)
tol=0.062500 SB: 6/5 = 1.200000 (err=0.028290)    CF: 6/5 = 1.200000 (err=0.028290)
tol=0.031250 SB: 6/5 = 1.200000 (err=0.028290)    CF: 6/5 = 1.200000 (err=0.028290)
tol=0.015625 SB: 7/6 = 1.166667 (err=-0.005043)    CF: 7/6 = 1.166667 (err=-0.005043)
tol=0.007812 SB: 7/6 = 1.166667 (err=-0.005043)    CF: 7/6 = 1.166667 (err=-0.005043)
tol=0.003906 SB: 27/23 = 1.173913 (err=0.002203)    CF: 34/29 = 1.172414 (err=0.000704)
tol=0.001953 SB: 34/29 = 1.172414 (err=0.000704)    CF: 34/29 = 1.172414 (err=0.000704)
tol=0.000977 SB: 34/29 = 1.172414 (err=0.000704)    CF: 34/29 = 1.172414 (err=0.000704)
tol=0.000488 SB: 41/35 = 1.171429 (err=-0.000281)    CF: 41/35 = 1.171429 (err=-0.000281)
tol=0.000244 SB: 75/64 = 1.171875 (err=0.000165)    CF: 116/99 = 1.171717 (err=0.000007)
tol=0.000122 SB: 116/99 = 1.171717 (err=0.000007)    CF: 116/99 = 1.171717 (err=0.000007)
tol=0.000061 SB: 116/99 = 1.171717 (err=0.000007)    CF: 116/99 = 1.171717 (err=0.000007)
tol=0.000031 SB: 116/99 = 1.171717 (err=0.000007)    CF: 116/99 = 1.171717 (err=0.000007)
tol=0.000015 SB: 116/99 = 1.171717 (err=0.000007)    CF: 116/99 = 1.171717 (err=0.000007)
tol=0.000008 SB: 116/99 = 1.171717 (err=0.000007)    CF: 116/99 = 1.171717 (err=0.000007)
tol=0.000004 SB: 1085/926 = 1.171706 (err=-0.000004)    CF: 1549/1322 = 1.171710 (err=-0.000000)
tol=0.000002 SB: 1317/1124 = 1.171708 (err=-0.000002)    CF: 1549/1322 = 1.171710 (err=-0.000000)
tol=0.000001 SB: 1549/1322 = 1.171710 (err=-0.000000)    CF: 1549/1322 = 1.171710 (err=-0.000000)
#include <stdio.h>
#include <math.h>
int main() {
float target = 1.2734;
int Aprev[2] = {1, 0};
int Bprev[2] = {0, 1};
float x = target;
for(int i=0;i<10;++i) {
int n = floor(x);
x-=n;
x=1/x;
int A = Aprev[0] + n * Aprev[1];
Aprev[0] = Aprev[1];
Aprev[1] = A;
int B = Bprev[0] + n * Bprev[1];
Bprev[0] = Bprev[1];
Bprev[1] = B;
double approx = (double)B/(double)A;
printf("%d %d/%d = %f (eps = %f)\n",i, B, A, approx, approx-target);
}
}
// Implements continued fraction and Stern Brocot tree search
#include <stdio.h>
#include <math.h>
struct fraction {
int n;
int d;
fraction(int n, int d) {
this->n = n;
this->d = d;
}
float asFloat() {
return float((double)n/(double)d);
}
};
fraction continued_fraction_approximation(float f, float tol) {
int Aprev[2] = {1, 0};
int Bprev[2] = {0, 1};
float x = f;
while(true) {
int n = floor(x);
x-=n;
x=1/x;
int A = Aprev[0] + n * Aprev[1];
Aprev[0] = Aprev[1];
Aprev[1] = A;
int B = Bprev[0] + n * Bprev[1];
Bprev[0] = Bprev[1];
Bprev[1] = B;
double approx = (double)B/(double)A;
if( fabs(approx - f) < tol)
return fraction(B,A);
}
}
fraction stern_brocot_search(float f, float tol) {
int base = floor(f);
f-=base;
if(f<tol)
return fraction(base,1);
if(1-tol < f)
return fraction(base+1,1);
fraction lower(0,1);
fraction upper(1,1);
while(true) {
fraction middle(lower.n+upper.n, lower.d+upper.d);
if( middle.d*(f+tol) < middle.n ) {
upper = middle;
} else if (middle.n < middle.d*(f-tol) ) {
lower = middle;
} else {
return fraction(middle.d*base + middle.n, middle.d);
}
}
}
int main() {
float target = 1.17171;
float tol=0.5;
for(int i=0; i<20; ++i) {
fraction sb = stern_brocot_search(target, tol);
fraction cf = continued_fraction_approximation(target,tol);
printf("tol=%f SB: %d/%d = %f (err=%f) CF: %d/%d = %f (err=%f)\n", tol, sb.n, sb.d, sb.asFloat(), sb.asFloat()-target, cf.n, cf.d, cf.asFloat(), cf.asFloat()-target );
tol/=2;
}
}
@tomalakgeretkal
Copy link

Unfortunately this goes horribly wrong when the whole part of the input is out of the range of int. Since floating-point values can hold massive numbers (albeit imprecisely) just whacking in an int64_t isn't a catch-all either.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment