Fraction approximation in C++

README.md

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)

main.cpp

    #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);
    	}
    }

main2.cpp

// 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;
	}
}

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.