From 410a4689002b0af693790a25121fb11e9f1e5ad5 Mon Sep 17 00:00:00 2001
From: Reiner Herrmann <reiner@reiner-h.de>
Date: Wed, 13 Oct 2010 03:44:37 +0200
Subject: projecteuler solution 066

---
 src/projecteuler/065.c |   2 +-
 src/projecteuler/066.c | 112 +++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 113 insertions(+), 1 deletion(-)
 create mode 100644 src/projecteuler/066.c

(limited to 'src')

diff --git a/src/projecteuler/065.c b/src/projecteuler/065.c
index 3e93085..94c17fd 100644
--- a/src/projecteuler/065.c
+++ b/src/projecteuler/065.c
@@ -23,7 +23,7 @@ int main(void)
 	mpz_init_set_ui(denominator, factors[limit-2]);
 	for(pos=limit-3; pos>=0; pos--)
 	{
-		mpz_addmul_ui(numerator, denominator, factors[pos]); // numerator *= denominator * factors[pos]
+		mpz_addmul_ui(numerator, denominator, factors[pos]); // numerator += denominator * factors[pos]
 		mpz_swap(numerator, denominator);
 	}
 
diff --git a/src/projecteuler/066.c b/src/projecteuler/066.c
new file mode 100644
index 0000000..3c1014e
--- /dev/null
+++ b/src/projecteuler/066.c
@@ -0,0 +1,112 @@
+#include <gmp.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+// numerator/denominator pairs of continued fractions are solutions to pell's equations
+
+/*
+ * taken from problem 65
+ *
+ * numerator, denominator: resulting num / denom are stored there
+ * limit: depth of continued fraction
+ * start: decimal part (pre-period)
+ * factors, factor_len: period of continued fraction with factor_len elements
+ *
+ */
+void compute(mpz_t numerator, mpz_t denominator, unsigned int limit,
+		unsigned long start, const unsigned int* factors, unsigned int factor_len)
+{
+	int i;
+	mpz_set_ui(numerator, 1);
+	mpz_set_ui(denominator, factors[limit % factor_len]);
+	for(i=limit-1; i>=0; i--)
+	{
+		mpz_addmul_ui(numerator, denominator, factors[i % factor_len]);
+		mpz_swap(numerator, denominator);
+	}
+	mpz_addmul_ui(numerator, denominator, start);
+}
+
+// taken from problem 64
+int continued_fraction(unsigned int* fraction, unsigned int n)
+{
+	double sqrtn, x, y;
+	unsigned int an, a0, pos=0;
+
+	sqrtn = sqrt((double)n);
+	a0 = (int)floor(sqrtn);
+
+	fraction[pos++] = sqrtn;
+
+	if((int)sqrtn*sqrtn == n) // works only for rational square roots
+		return 1;
+
+	x = (double)a0;
+	y = 1.0;
+
+	while(1)
+	{
+		double denominator = (n - x*x);
+		double tmp = y * (sqrtn + x) / denominator;
+		an = (int)floor(tmp);
+		x = (an * denominator / y) - x;
+		y = denominator / y;
+
+		fraction[pos++] = an;
+
+		if(an == 2*a0)
+			break;
+	}
+
+	return pos;
+}
+
+const int max_period_size = 100;
+
+int main(void)
+{
+	int d, i, maxd=0;
+	mpz_t maxx, numerator, denominator, tmp1, tmp2;
+	unsigned int* fraction = (unsigned int*) malloc(max_period_size*sizeof(unsigned int));
+
+	mpz_init_set_ui(maxx, 0);
+	mpz_init(numerator);
+	mpz_init(denominator);
+	mpz_init(tmp1);
+	mpz_init(tmp2);
+
+	for(d=2; d<=1000; d++)
+	{
+
+		memset(fraction, 0, max_period_size*sizeof(unsigned int));
+		int frac_len = continued_fraction(fraction, d);
+
+		if(frac_len == 1) // no solution for squares
+			continue;
+
+		for(i=0;;i++)
+		{
+			compute(numerator, denominator, i, fraction[0], &(fraction[1]), frac_len-1);
+
+			// numerator^2 - denominator^2 * d
+			mpz_mul(tmp1, numerator, numerator);
+			mpz_mul(tmp2, denominator, denominator);
+			mpz_submul_ui(tmp1, tmp2, d);
+			if(mpz_cmp_ui(tmp1, 1) == 0)
+				break;
+		}
+
+		if(mpz_cmp(numerator, maxx) > 0)
+		{
+			mpz_set(maxx, numerator);
+			maxd = d;
+		}
+	}
+
+	printf("%d\n", maxd);
+
+	return 0;
+}
+
-- 
cgit v1.2.3