Solved 66, made some initialization error and had to gradually build up the chain instead of recomputing it everytime

This commit is contained in:
2021-11-13 01:52:05 +01:00
parent d2d19d3efa
commit 8f5217e793

View File

@@ -25,12 +25,13 @@ The value of x can be found by calculating the convergent of the continued fract
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 solution.
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;
@@ -45,6 +46,13 @@ 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);
@@ -83,18 +91,75 @@ vector<int> findCF(int n){
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();
int maxX = 0;
for(int d = 1; d <= 1000; ++d){
auto cf = findCF(d);
InfInt maxX = 0;
int maxD = 1;
for(int d = 1; d <= 1000; ++d){
if(sqrt(d) == sqrt((float)d)){
continue;
}
// Insert problem 065 code once I have that fixed for big integers (cpp needs a dedicated internal library for that)
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;