#include #include #include #include const int size = 100; double *fraction; int pos = 0; void continued_fraction(int n) { double sqrtn, x, y, an, a0; memset(fraction, 0, size*sizeof(double)); pos = 0; sqrtn = sqrt((double)n); a0 = floor(sqrtn); if((int)sqrtn*sqrtn == n) // works only for rational square roots { fraction[pos] = sqrtn; return; } // // floor + ( sqrt(23) - x ) / y // nenner = (23 - x*x) // tmp = y * (sqrt(23) + x) / nenner) // an = floor(tmp) // x = an*nenner/y - x // x = a0; y = 1.0; while(1) { double denominator = (n - x*x); double tmp = y * (sqrtn + x) / denominator; an = floor(tmp); x = (an * denominator / y) - x; y = denominator / y; fraction[pos++] = an; if(an == 2*a0) break; } /* // not working because of double imprecision while(1) { an = floorl(a); r = 1/(a - an); a = r; fraction[++pos] = an; if(an == 2*a0) break; } */ } int main(void) { int i, count=0; fraction = (double*) malloc(size*sizeof(double)); for(i=2; i<=10000; i++) { continued_fraction(i); if(pos > 0 && pos&1 == 1) count++; } printf("%i\n", count); return 0; }