Pell Equation
C++. This is the problem 66 from Project Euler.
Question: 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:
3^2 – 2×2^2 = 1
2^2 – 3×1^2 = 1
9^2 – 5×4^2 = 1
5^2 – 6×2^2 = 1
8^2 – 7×3^2 = 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. Link to the page.
Solution: We could call it “Pell Equation” as well. To be honest, it was pretty exhaustive to write an efficient algorithm for its solution. It was harder than I thought. I implemented “Chakravala Method”. Main problem was rounding error in C++ when the numbers getting really high (particularly after 16 digit). That’s why, I decided to use double data types for x, and y. Algorithm starts with setting up m value for a [x,y,k] trivial triple, then goes up to search x, y, and k values. When k reaches 1, that gives x and y values, and stops there. Also, there is a shortcut calculation for k = -1 and abs(k) = 2 values inside the algorithm. The answer is 661.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 | #include <iostream> #include <cmath> #include <vector> void CalculateCoefficientValue(const int _D); void pellEquation(const int _d, std::vector<double> & _xValues); int main() { const int D = 1000; CalculateCoefficientValue(D); return 0; } //////////////////////////////////////////////////////////////////////////////// void CalculateCoefficientValue(const int _D) { std::vector<double> xValues; //storing x values for each constant int coefficient = 1; double maxValue = 1; xValues.push_back(0); //indexing the first element to allign for (int d = 1; d <= _D; ++d) { double root = sqrt(d); //checking if d is square or not if (root == floor(root)) //if it is square, go to next integer { xValues.push_back(0); continue; } pellEquation(d, xValues); //searching for x, y if (xValues[d] >= maxValue) //indexing max value and coefficient d { maxValue = xValues[d]; coefficient = d; } } std::cout << "x domain will be created when the constant D goes from 1 to "<< _D << std::endl; std::cout << "The constant D is " << coefficient << " where the max x value lies in the x domain" << std::endl; std::cout << "x value is " << maxValue << " ,where D is " << coefficient << std::endl; } //////////////////////////////////////////////////////////////////////////////// void pellEquation(const int _d, std::vector<double> & _xValues) { int m = 1; int k = 1; double x = 0; double y = 0; double xref = 1; //trivial triple (m, 1, m^2 – D) to get a new triple (mx + Dy, x + my, k(m^2 – D)). //scaling down by k, using Bhaskara’s lemma, we get x^2 – Dy^2 = k while (k != 1 || y == 0) { //calculating m value to get [x, y, k] triple m = k * (m / k + 1) - m; int z = (m - sqrt(_d)) / k; m = m - z * k; x = (m * xref + _d * y) / abs(k); y = (m * y + xref) / abs(k); k = (pow(m, 2) - _d) / k; //Brahmagupta Lemma (shortcut) if (k == -1 || abs(k) == 2) { x = (pow(x, 2) + _d * pow(y, 2)) / abs(k); y = (2 * x * y) / abs(k); break; } xref = x; } _xValues.push_back(x); } |