1/* 2 * pell - solve Pell's equation 3 * 4 * Copyright (C) 1999,2021 David I. Bell 5 * 6 * Calc is open software; you can redistribute it and/or modify it under 7 * the terms of the version 2.1 of the GNU Lesser General Public License 8 * as published by the Free Software Foundation. 9 * 10 * Calc is distributed in the hope that it will be useful, but WITHOUT 11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 12 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General 13 * Public License for more details. 14 * 15 * A copy of version 2.1 of the GNU Lesser General Public License is 16 * distributed with calc under the filename COPYING-LGPL. You should have 17 * received a copy with calc; if not, write to Free Software Foundation, Inc. 18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 19 * 20 * Under source code control: 1990/02/15 01:50:34 21 * File existed as early as: before 1990 22 * 23 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ 24 */ 25 26/* 27 * Solve Pell's equation; Returns the solution X to: X^2 - D * Y^2 = 1. 28 * Type the solution to Pell's equation for a particular D. 29 */ 30 31 32define pell(D) 33{ 34 local X, Y; 35 36 X = pellx(D); 37 if (isnull(X)) { 38 print "D=":D:" is square"; 39 return; 40 } 41 Y = isqrt((X^2 - 1) / D); 42 print X : "^2 - " : D : "*" : Y : "^2 = " : X^2 - D*Y^2; 43} 44 45 46/* 47 * Function to solve Pell's equation 48 * Returns the solution X to: 49 * X^2 - D * Y^2 = 1 50 */ 51define pellx(D) 52{ 53 local R, Rp, U, Up, V, Vp, A, T, Q1, Q2, n; 54 local mat ans[2,2]; 55 local mat tmp[2,2]; 56 57 R = isqrt(D); 58 Vp = D - R^2; 59 if (Vp == 0) 60 return; 61 Rp = R + R; 62 U = Rp; 63 Up = U; 64 V = 1; 65 A = 0; 66 n = 0; 67 ans[0,0] = 1; 68 ans[1,1] = 1; 69 tmp[0,1] = 1; 70 tmp[1,0] = 1; 71 do { 72 T = V; 73 V = A * (Up - U) + Vp; 74 Vp = T; 75 A = U // V; 76 Up = U; 77 U = Rp - U % V; 78 tmp[0,0] = A; 79 ans *= tmp; 80 n++; 81 } while (A != Rp); 82 Q2 = ans[[1]]; 83 Q1 = isqrt(Q2^2 * D + 1); 84 if (isodd(n)) { 85 T = Q1^2 + D * Q2^2; 86 Q2 = Q1 * Q2 * 2; 87 Q1 = T; 88 } 89 return Q1; 90} 91