projects / unix-history / blame
| Commit | Line | Data |
|---|---|---|
| 198b5bf9 C |
1 | /* |
| 2 | * Copyright (c) 1993 David I. Bell | |
| 3 | * Permission is granted to use, distribute, or modify this source, | |
| 4 | * provided that this copyright notice remains intact. | |
| 5 | * | |
| 6 | * Solve Pell's equation; Returns the solution X to: X^2 - D * Y^2 = 1. | |
| 7 | * Type the solution to pells equation for a particular D. | |
| 8 | */ | |
| 9 | ||
| 10 | define pell(D) | |
| 11 | { | |
| 12 | local X, Y; | |
| 13 | ||
| 14 | X = pellx(D); | |
| 15 | if (isnull(X)) { | |
| 16 | print "D=":D:" is square"; | |
| 17 | return; | |
| 18 | } | |
| 19 | Y = isqrt((X^2 - 1) / D); | |
| 20 | print X : "^2 - " : D : "*" : Y : "^2 = " : X^2 - D*Y^2; | |
| 21 | } | |
| 22 | ||
| 23 | ||
| 24 | /* | |
| 25 | * Function to solve Pell's equation | |
| 26 | * Returns the solution X to: | |
| 27 | * X^2 - D * Y^2 = 1 | |
| 28 | */ | |
| 29 | define pellx(D) | |
| 30 | { | |
| 31 | local R, Rp, U, Up, V, Vp, A, T, Q1, Q2, n; | |
| 32 | local mat ans[2,2]; | |
| 33 | local mat tmp[2,2]; | |
| 34 | ||
| 35 | R = isqrt(D); | |
| 36 | Vp = D - R^2; | |
| 37 | if (Vp == 0) | |
| 38 | return; | |
| 39 | Rp = R + R; | |
| 40 | U = Rp; | |
| 41 | Up = U; | |
| 42 | V = 1; | |
| 43 | A = 0; | |
| 44 | n = 0; | |
| 45 | ans[0,0] = 1; | |
| 46 | ans[1,1] = 1; | |
| 47 | tmp[0,1] = 1; | |
| 48 | tmp[1,0] = 1; | |
| 49 | do { | |
| 50 | T = V; | |
| 51 | V = A * (Up - U) + Vp; | |
| 52 | Vp = T; | |
| 53 | A = U // V; | |
| 54 | Up = U; | |
| 55 | U = Rp - U % V; | |
| 56 | tmp[0,0] = A; | |
| 57 | ans *= tmp; | |
| 58 | n++; | |
| 59 | } while (A != Rp); | |
| 60 | Q2 = ans[[1]]; | |
| 61 | Q1 = isqrt(Q2^2 * D + 1); | |
| 62 | if (isodd(n)) { | |
| 63 | T = Q1^2 + D * Q2^2; | |
| 64 | Q2 = Q1 * Q2 * 2; | |
| 65 | Q1 = T; | |
| 66 | } | |
| 67 | return Q1; | |
| 68 | } | |
| 69 | ||
| 70 | global lib_debug; | |
| 71 | if (lib_debug >= 0) { | |
| 72 | print "pell(D) defined"; | |
| 73 | print "pellx(D) defined"; | |
| 74 | } |