1 files changed, 112 insertions, 0 deletions
diff --git a/066.c b/066.c
new file mode 100644
index 0000000..3c1014e
--- /dev/null
+++ b/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;
+}
+