064.c


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>


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