167 lines
4.1 KiB
C++
167 lines
4.1 KiB
C++
/*
|
||
Consider quadratic Diophantine equations of the form:
|
||
|
||
x2 – Dy2 = 1
|
||
|
||
For example, when D=13, the minimal solution in x is 6492 – 13×1802 = 1.
|
||
|
||
It can be assumed that there are no solutions in positive integers when D is square.
|
||
|
||
By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:
|
||
|
||
32 – 2×22 = 1
|
||
22 – 3×12 = 1
|
||
92 – 5×42 = 1
|
||
52 – 6×22 = 1
|
||
82 – 7×32 = 1
|
||
|
||
Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.
|
||
|
||
Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.
|
||
|
||
-----------------------------------------------------------------------------------------------------
|
||
|
||
The value of x can be found by calculating the convergent of the continued fraction of D,
|
||
right before it repeats.
|
||
|
||
So if y/x is the corresponding approximation of sqrt(D) for the value of a_n of the continued fraction right before it repeats,
|
||
we immediately know x. So in order to solve this we just have to combine the previous two problem solutions.
|
||
|
||
https://mathshistory.st-andrews.ac.uk/HistTopics/Pell/
|
||
*/
|
||
|
||
#include <bits/stdc++.h>
|
||
#include "../bigint.h"
|
||
|
||
using namespace std;
|
||
|
||
int gcd(int a, int b){
|
||
if(b == 0){
|
||
return a;
|
||
}
|
||
return gcd(b, a % b);
|
||
}
|
||
|
||
int gcd(int a, int b, int c){
|
||
return gcd(a, gcd(b, c));
|
||
}
|
||
|
||
InfInt gcd(InfInt a, InfInt b){
|
||
if(b == 0){
|
||
return a;
|
||
}
|
||
return gcd(b, a % b);
|
||
}
|
||
|
||
// Code yanked from problem 64
|
||
vector<int> findCF(int n){
|
||
float x = sqrt((float)n);
|
||
|
||
vector<int> cf;
|
||
|
||
if(x != sqrt(n)){
|
||
int a, b = 1, c = 1, d = 0;
|
||
int bn, cn, dn, g;
|
||
int b1, c1, d1;
|
||
|
||
for(int i = 0; ; ++i){
|
||
a = floor((floor(b * x) + d) / c);
|
||
cf.push_back(a);
|
||
|
||
bn = b*c;
|
||
cn = b*b*n - d*d - a*a*c*c + 2*a*c*d;
|
||
dn = a*c*c - c*d;
|
||
|
||
g = gcd(bn, cn, dn);
|
||
|
||
b = bn / g;
|
||
c = cn / g;
|
||
d = dn / g;
|
||
|
||
if(i == 0){
|
||
b1 = b;
|
||
c1 = c;
|
||
d1 = d;
|
||
} else if(b1 == b && c1 == c && d1 == d){
|
||
break;
|
||
}
|
||
}
|
||
}
|
||
|
||
return cf;
|
||
}
|
||
|
||
pair<InfInt, InfInt> fracAdd(const pair<InfInt, InfInt>& f1, const pair<InfInt, InfInt>& f2){
|
||
pair<InfInt, InfInt> result = {f1.first*f2.second + f2.first*f1.second, f1.second * f2.second};
|
||
InfInt g = gcd(result.first, result.second);
|
||
|
||
result.first /= g;
|
||
result.second /= g;
|
||
|
||
return result;
|
||
}
|
||
|
||
pair<InfInt, InfInt> fracInv(const pair<InfInt, InfInt>& f){
|
||
return pair<InfInt, InfInt>(f.second, f.first);
|
||
}
|
||
|
||
pair<InfInt, InfInt> cfApprox(const vector<int>& seq){
|
||
pair<InfInt, InfInt> result = {seq[0], 1};
|
||
pair<InfInt, InfInt> carry = {0, 1};
|
||
if(seq.size() > 1){
|
||
carry = {1, seq[1]};
|
||
for(size_t i = seq.size() - 2; i > 0; --i){
|
||
pair<InfInt, InfInt> newp = {seq[i], 1};
|
||
carry = fracInv(fracAdd(carry, newp));
|
||
}
|
||
}
|
||
|
||
result = fracAdd(result, carry);
|
||
|
||
return result;
|
||
}
|
||
|
||
int main(){
|
||
|
||
cout << "Hello this is Patrick" << endl;
|
||
auto start = chrono::high_resolution_clock::now();
|
||
|
||
InfInt maxX = 0;
|
||
int maxD = 1;
|
||
for(int d = 1; d <= 1000; ++d){
|
||
if(sqrt(d) == sqrt((float)d)){
|
||
continue;
|
||
}
|
||
|
||
vector<int> cf = findCF(d);
|
||
|
||
pair<InfInt, InfInt> init = {cf[0], 1};
|
||
pair<InfInt, InfInt> carry = {0, 1};
|
||
int n = 0;
|
||
InfInt x;
|
||
|
||
while(true){
|
||
auto result = fracAdd(init, carry);
|
||
|
||
if(result.first * result.first - (InfInt)d * result.second * result.second == 1){
|
||
x = result.first;
|
||
break;
|
||
}
|
||
|
||
carry = fracInv(fracAdd(carry, {cf[(n % (cf.size() - 1)) + 1], 1}));
|
||
n++;
|
||
}
|
||
|
||
maxX = max(x, maxX);
|
||
if(maxX == x){
|
||
maxD = d;
|
||
}
|
||
}
|
||
|
||
cout << "Maximum: " << maxD << endl;
|
||
|
||
auto duration = chrono::duration_cast<chrono::milliseconds>(chrono::high_resolution_clock::now() - start);
|
||
cout << (float)duration.count()/1000 << endl;
|
||
return 0;
|
||
}
|