project-euler

https://projecteuler.net/
Log | Files | Refs | README

Euler_66.cpp (1640B)

      1 #include "Euler.h"
      2 
      3 bool valueOfDiophantine(cpp_int x, cpp_int y, cpp_int n)
      4 {
      5     return EulerUtility::power(x, 2) - (n * EulerUtility::power(y, 2)) == cpp_int(1);
      6 }
      7 
      8 std::vector<cpp_int> recurseFraction(std::vector<cpp_int> period, cpp_int n, std::vector<cpp_int> fraction)
      9 {
     10     if (n > period.size() - 1)
     11         return fraction;
     12 
     13     std::swap(fraction[0], fraction[1]);
     14 
     15     fraction[0] = (fraction[1] * period[period.size() - n.toInt() - 1]) + fraction[0];
     16 
     17     return recurseFraction(period, n + 1, fraction);
     18 }
     19 
     20 std::vector<cpp_int> recurseFraction(std::vector<cpp_int> period, cpp_int n)
     21 {
     22     std::vector<cpp_int> fraction;
     23 
     24     fraction.push_back(period[period.size() - 1]);
     25     fraction.push_back(1);
     26 
     27     return recurseFraction(period, 1, fraction);
     28 }
     29 
     30 cpp_int valueOfX(cpp_int n)
     31 {
     32     double n2 = std::sqrtl(n.toInt());
     33     cpp_int a = (int)n2, p = 0, q = 1;
     34 
     35     std::vector<cpp_int> period;
     36     std::vector<cpp_int> approx;
     37 
     38     period.push_back(a);
     39 
     40     do
     41     {
     42         p = a * q - p;
     43         q = ( n - p * p ) / q;
     44         a = (long)(( p.toLong() + n2 ) /q.toLong());
     45         
     46         period.push_back(a);
     47         approx = recurseFraction(period, 0);
     48 
     49     } while (valueOfDiophantine(approx[0], approx[1], n) != 1);
     50 
     51     return approx[0];
     52 }
     53 
     54 int Euler::Diophantine()
     55 {
     56     cpp_int currentMax = 0;
     57     int valueOfD = 0;
     58 
     59     for (int i = 2; i <= 1000; ++i)
     60     {
     61         if (!EulerUtility::isPerfectSquare(i))
     62         {
     63             cpp_int x = valueOfX(i);
     64 
     65             if (x > currentMax)
     66             {
     67                 currentMax = x;
     68                 valueOfD = i;
     69             }
     70         }
     71     }
     72 
     73     return valueOfD;
     74 }