/*
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;
}