@@ -108,6 +108,92 @@ We get $x_{1} = 1$. Update the triple $(p_{1}, q_{1}, m_{1}) = ( \frac{1 \cdot 3
108108
109109##Implementation
110110``` cpp
111+ bool isSquare (long long n) {
112+ long long sqrtN = (long long)sqrt(n);
113+ return sqrtN * sqrtN == n;
114+ }
115+
116+ long long mod(long long a, long long b) {
117+ return (a % b + b) % b;
118+ }
119+
120+ long long modInv(long long a, long long b) {
121+ long long b0 = b, x0 = 0, x1 = 1;
122+ if (b == 1) return 1;
123+ while (a > 1) {
124+ long long q = a / b;
125+ long long temp = b;
126+ b = a % b;
127+ a = temp;
128+ temp = x0;
129+ x0 = x1 - q * x0;
130+ x1 = temp;
131+ }
132+ if (x1 < 0) x1 += b0;
133+ return x1;
134+ }
135+
136+
137+ // Chakravala method for solving Pell's equation
138+ pair<long long, long long> chakravala(int n) {
139+ // Check if n is a perfect square
140+ if (isSquare(n)) {
141+ throw invalid_argument("n is a perfect square. No solutions exist for Pell's equation.");
142+ }
143+
144+ // Initial values
145+ double sqrt_n = sqrt(n);
146+ long long a = (long long)floor(sqrt_n); // Initial a
147+ long long b = 1; // Initial b
148+ long long k = a * a - n; // Initial k
149+
150+ int steps = 0; // Step counter for iterations
151+
152+ // Repeat until k = 1
153+ while (k != 1) {
154+ long long absK = abs(k);
155+
156+ // Find m such that k | (a + bm), and minimize |m^2 - n|
157+ long long m;
158+ if (absK == 1) {
159+ m = (long long)round(sqrt(n)); // round to nearest integer
160+ } else {
161+ long long r = mod(-a, absK); // Negative of a mod(k) // (a + m*b)/|k|
162+ long long s = modInv(b, absK); // Modular inverse of b mod(k)
163+ m = mod(r * s, absK); // Compute m for (a + b*m) mod(k) = 0
164+
165+ // Approximate value of m
166+ // m = m + ((long long)floor((sqrt_n - m) / absK)) * absK;
167+
168+ // Adjust m to ensure m < sqrt(n) < m + k
169+ while (m > sqrt(n)) m -= absK;
170+ while (m + absK < sqrt_n) m += absK;
171+
172+ // Select closest value to n
173+ if (abs(m * m - n) > abs((m + absK) * (m + absK) - n)) {
174+ m = m + absK;
175+ }
176+ }
177+
178+ // Print the current triple
179+ cout << "[a = " << a << ", b = " << b << ", k = " << k << "]" << endl;
180+
181+ // Update a, b, k using the recurrence relations
182+ long long alpha = a;
183+ a = (m * a + n * b) / absK;
184+ b = (alpha + m * b) / absK;
185+ k = (m * m - n) / k;
186+
187+ // Increment step counter
188+ steps++;
189+ }
190+
191+ // Print final result
192+ cout << a << "^2 - " << n << " x " << b << "^2 = 1 in " << steps << " calculations." << endl;
193+
194+ // Return the solution as a pair (a, b)
195+ return {a, b};
196+ }
111197
112198```
113199