diff --git a/projecteuler/066/main.cpp b/projecteuler/066/main.cpp
index f9e67c8..bec9b24 100644
--- a/projecteuler/066/main.cpp
+++ b/projecteuler/066/main.cpp
@@ -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;