From aea7cccb4e77399516d5aa91af39037869d045ac Mon Sep 17 00:00:00 2001
From: Reiner Herrmann <reiner@reiner-h.de>
Date: Sun, 10 Oct 2010 23:41:55 +0200
Subject: projecteuler solution 064

---
 src/projecteuler/064.c | 84 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 84 insertions(+)
 create mode 100644 src/projecteuler/064.c

diff --git a/src/projecteuler/064.c b/src/projecteuler/064.c
new file mode 100644
index 0000000..0e60211
--- /dev/null
+++ b/src/projecteuler/064.c
@@ -0,0 +1,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;
+}
+
-- 
cgit v1.2.3