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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
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;
}
|
