Pell's equation: Difference between revisions
Langurmonkey (talk | contribs) (→{{header|langur}}: format number string with commas) |
Langurmonkey (talk | contribs) |
||
(38 intermediate revisions by 21 users not shown) | |||
Line 3: | Line 3: | ||
'''Pell's equation''' (also called the '''Pell–Fermat''' equation) is a [https://en.wikipedia.org/wiki/Diophantine_equation <u>Diophantine equation</u>] of the form: |
'''Pell's equation''' (also called the '''Pell–Fermat''' equation) is a [https://en.wikipedia.org/wiki/Diophantine_equation <u>Diophantine equation</u>] of the form: |
||
:::::: <big> <b> x<sup>2</sup> - ny<sup>2</sup> = 1 </b> </big> |
:::::: <big> <b> x<sup>2</sup> - ny<sup>2</sup> = 1 </b> </big> |
||
with integer solutions for '''x''' and '''y''', where '''n''' is a given non-square positive integer. |
with integer solutions for '''x''' and '''y''', where '''n''' is a given non-square positive integer. |
||
Line 15: | Line 15: | ||
:* Wikipedia entry: [https://en.wikipedia.org/wiki/Pell%27s_equation <u>Pell's equation</u>]. |
:* Wikipedia entry: [https://en.wikipedia.org/wiki/Pell%27s_equation <u>Pell's equation</u>]. |
||
<br><br> |
<br><br> |
||
=={{header|11l}}== |
|||
{{trans|Python}} |
|||
<syntaxhighlight lang="11l">F solvePell(n) |
|||
V x = Int(sqrt(n)) |
|||
V (y, z, r) = (x, 1, x << 1) |
|||
BigInt e1 = 1 |
|||
BigInt e2 = 0 |
|||
BigInt f1 = 0 |
|||
BigInt f2 = 1 |
|||
L |
|||
y = r * z - y |
|||
z = (n - y * y) I/ z |
|||
r = (x + y) I/ z |
|||
(e1, e2) = (e2, e1 + e2 * r) |
|||
(f1, f2) = (f2, f1 + f2 * r) |
|||
V (a, b) = (f2 * x + e2, f2) |
|||
I a * a - n * b * b == 1 |
|||
R (a, b) |
|||
L(n) [61, 109, 181, 277] |
|||
V (x, y) = solvePell(n) |
|||
print(‘x^2 - #3 * y^2 = 1 for x = #27 and y = #25’.format(n, x, y))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 |
|||
x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 |
|||
x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 |
|||
x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020 |
|||
</pre> |
|||
=={{header|Ada}}== |
|||
{{Trans|C}}{{works with|Ada 2022}} |
|||
<syntaxhighlight lang="ada">with Ada.Text_Io; |
|||
with Ada.Numerics.Elementary_Functions; |
|||
with Ada.Numerics.Big_Numbers.Big_Integers; |
|||
procedure Pells_Equation is |
|||
use Ada.Numerics.Big_Numbers.Big_Integers; |
|||
type Pair is |
|||
record |
|||
V1, V2 : Big_Integer; |
|||
end record; |
|||
procedure Solve_Pell (N : Natural; X, Y : out Big_Integer) is |
|||
use Ada.Numerics.Elementary_Functions; |
|||
Big_N : constant Big_Integer := To_Big_Integer (N); |
|||
XX : constant Big_Integer := To_Big_Integer (Natural (Float'Floor (Sqrt (Float (N))))); |
|||
begin |
|||
if XX**2 = Big_N then |
|||
X := 1; Y := 0; |
|||
return; |
|||
end if; |
|||
declare |
|||
YY : Big_Integer := XX; |
|||
Z : Big_Integer := 1; |
|||
R : Big_Integer := 2 * XX; |
|||
E : Pair := Pair'(V1 => 1, V2 => 0); |
|||
F : Pair := Pair'(V1 => 0, V2 => 1); |
|||
begin |
|||
loop |
|||
YY := R * Z - YY; |
|||
Z := (Big_N - YY**2) / Z; |
|||
R := (XX + YY) / Z; |
|||
E := Pair'(V1 => E.V2, V2 => R * E.V2 + E.V1); |
|||
F := Pair'(V1 => F.V2, V2 => R * F.V2 + F.V1); |
|||
X := E.V2 + XX * F.V2; |
|||
Y := F.V2; |
|||
exit when X**2 - Big_N * Y**2 = 1; |
|||
end loop; |
|||
end; |
|||
end Solve_Pell; |
|||
procedure Test (N : Natural) is |
|||
package Natural_Io is new Ada.Text_Io.Integer_Io (Natural); |
|||
use Ada.Text_Io, Natural_Io; |
|||
X, Y : Big_Integer; |
|||
begin |
|||
Solve_Pell (N, X => X, Y => Y); |
|||
Put ("X**2 - "); |
|||
Put (N, Width => 3); |
|||
Put (" * Y**2 = 1 for X = "); |
|||
Put (To_String (X, Width => 22)); |
|||
Put (" and Y = "); |
|||
Put (To_String (Y, Width => 20)); |
|||
New_Line; |
|||
end Test; |
|||
begin |
|||
Test (61); |
|||
Test (109); |
|||
Test (181); |
|||
Test (277); |
|||
end Pells_Equation;</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
X**2 - 61 * Y**2 = 1 for X = 1766319049 and Y = 226153980 |
|||
X**2 - 109 * Y**2 = 1 for X = 158070671986249 and Y = 15140424455100 |
|||
X**2 - 181 * Y**2 = 1 for X = 2469645423824185801 and Y = 183567298683461940 |
|||
X**2 - 277 * Y**2 = 1 for X = 159150073798980475849 and Y = 9562401173878027020 |
|||
</pre> |
|||
=={{header|ALGOL 68}}== |
=={{header|ALGOL 68}}== |
||
{{Trans|Sidef}} Also tests for a trival solution only (if n is a perfect square only 1, 0 is solution). |
{{Trans|Sidef}} Also tests for a trival solution only (if n is a perfect square only 1, 0 is solution). |
||
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}} |
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}} |
||
< |
<syntaxhighlight lang="algol68">BEGIN |
||
# find solutions to Pell's eqauation: x^2 - ny^2 = 1 for integer x, y, n # |
# find solutions to Pell's eqauation: x^2 - ny^2 = 1 for integer x, y, n # |
||
MODE BIGINT = LONG LONG INT; |
MODE BIGINT = LONG LONG INT; |
||
Line 57: | Line 166: | ||
print( ("x^2 - ", whole( n, -3 ), " * y^2 = 1 for x = ", whole( v1 OF r, -21), " and y = ", whole( v2 OF r, -21 ), newline ) ) |
print( ("x^2 - ", whole( n, -3 ), " * y^2 = 1 for x = ", whole( v1 OF r, -21), " and y = ", whole( v2 OF r, -21 ), newline ) ) |
||
OD |
OD |
||
END</ |
END</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 65: | Line 174: | ||
x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020 |
x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020 |
||
</pre> |
</pre> |
||
=={{header|Arturo}}== |
|||
{{trans|Python}} |
|||
<syntaxhighlight lang="rebol">solvePell: function [n][ |
|||
x: to :integer sqrt n |
|||
[y, z, r]: @[x, 1, shl x 1] |
|||
[e1, e2]: [1, 0] |
|||
[f1, f2]: [0, 1] |
|||
while [true][ |
|||
y: (r * z) - y |
|||
z: (n - y * y) / z |
|||
r: (x + y) / z |
|||
[e1, e2]: @[e2, e1 + e2 * r] |
|||
[f1, f2]: @[f2, f1 + f2 * r] |
|||
[a, b]: @[e2 + f2 * x, f2] |
|||
if 1 = (a*a) - n*b*b -> |
|||
return @[a, b] |
|||
] |
|||
] |
|||
loop [61 109 181 277] 'n [ |
|||
[x, y]: solvePell n |
|||
print ["x² -" n "* y² = 1 for (x,y) =" x "," y] |
|||
]</syntaxhighlight> |
|||
{{out}} |
|||
<pre>x² - 61 * y² = 1 for (x,y) = 1766319049 , 226153980 |
|||
x² - 109 * y² = 1 for (x,y) = 158070671986249 , 15140424455100 |
|||
x² - 181 * y² = 1 for (x,y) = 2469645423824185801 , 183567298683461940 |
|||
x² - 277 * y² = 1 for (x,y) = 159150073798980475849 , 9562401173878027020</pre> |
|||
=={{header|C}}== |
=={{header|C}}== |
||
{{trans|ALGOL 68}} |
{{trans|ALGOL 68}} |
||
For n = 277, the x value is not correct because 64-bits is not enough to represent the value. |
For n = 277, the x value is not correct because 64-bits is not enough to represent the value. |
||
< |
<syntaxhighlight lang="c">#include <math.h> |
||
#include <stdbool.h> |
#include <stdbool.h> |
||
#include <stdint.h> |
#include <stdint.h> |
||
Line 120: | Line 264: | ||
void test(int n) { |
void test(int n) { |
||
struct Pair r = solvePell(n); |
struct Pair r = solvePell(n); |
||
printf("x^2 - %3d * y^2 = 1 for x = %21llu |
printf("x^2 - %3d * y^2 = 1 for x = %21llu and y = %21llu\n", n, r.v1, r.v2); |
||
} |
} |
||
Line 130: | Line 274: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 |
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 |
||
x^2 - 109 * y^2 = 1 for x = 158070671986249 |
x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 |
||
x^2 - 181 * y^2 = 1 for x = 2469645423824185801 |
x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 |
||
x^2 - 277 * y^2 = 1 for x = 11576121209304062921 |
x^2 - 277 * y^2 = 1 for x = 11576121209304062921 and y = 9562401173878027020</pre> |
||
=={{header|C++}}== |
=={{header|C++}}== |
||
{{trans|C}} |
{{trans|C}} |
||
As with the C solution, the output for n = 277 is not correct because 64-bits is not enough to represent the value. |
As with the C solution, the output for n = 277 is not correct because 64-bits is not enough to represent the value. |
||
< |
<syntaxhighlight lang="cpp">#include <iomanip> |
||
#include <iostream> |
#include <iostream> |
||
#include <tuple> |
#include <tuple> |
||
Line 189: | Line 333: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 |
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 |
||
Line 198: | Line 342: | ||
=={{header|C sharp|C#}}== |
=={{header|C sharp|C#}}== |
||
{{trans|Sidef}} |
{{trans|Sidef}} |
||
< |
<syntaxhighlight lang="csharp">using System; |
||
using System.Numerics; |
using System.Numerics; |
||
Line 228: | Line 372: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 |
<pre>x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 |
||
Line 237: | Line 381: | ||
=={{header|D}}== |
=={{header|D}}== |
||
{{trans|C#}} |
{{trans|C#}} |
||
< |
<syntaxhighlight lang="d">import std.bigint; |
||
import std.math; |
import std.math; |
||
import std.stdio; |
import std.stdio; |
||
Line 277: | Line 421: | ||
writefln("x^2 - %3d * y^2 = 1 for x = %27d and y = %25d", n, x, y); |
writefln("x^2 - %3d * y^2 = 1 for x = %27d and y = %25d", n, x, y); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 |
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 |
||
Line 283: | Line 427: | ||
x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 |
x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 |
||
x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020</pre> |
x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020</pre> |
||
=={{header|Delphi}}== |
|||
{{libheader| System.SysUtils}} |
|||
{{libheader| Velthuis.BigIntegers}} |
|||
{{Trans|Go}} |
|||
<syntaxhighlight lang="delphi"> |
|||
program Pells_equation; |
|||
{$APPTYPE CONSOLE} |
|||
uses |
|||
System.SysUtils, |
|||
Velthuis.BigIntegers; |
|||
type |
|||
TPellResult = record |
|||
x, y: BigInteger; |
|||
end; |
|||
function SolvePell(nn: UInt64): TPellResult; |
|||
var |
|||
n, x, y, z, r, e1, e2, f1, t, u, a, b: BigInteger; |
|||
begin |
|||
n := nn; |
|||
x := nn; |
|||
x := BigInteger.Sqrt(x); |
|||
y := BigInteger(x); |
|||
z := BigInteger.One; |
|||
r := x shl 1; |
|||
e1 := BigInteger.One; |
|||
e2 := BigInteger.Zero; |
|||
f1 := BigInteger.Zero; |
|||
b := BigInteger.One; |
|||
while True do |
|||
begin |
|||
y := (r * z) - y; |
|||
z := (n - (y * y)) div z; |
|||
r := (x + y) div z; |
|||
u := BigInteger(e1); |
|||
e1 := BigInteger(e2); |
|||
e2 := (r * e2) + u; |
|||
u := BigInteger(f1); |
|||
f1 := BigInteger(b); |
|||
b := r * b + u; |
|||
a := e2 + x * b; |
|||
t := (a * a) - (n * b * b); |
|||
if t = 1 then |
|||
begin |
|||
with Result do |
|||
begin |
|||
x := BigInteger(a); |
|||
y := BigInteger(b); |
|||
end; |
|||
Break; |
|||
end; |
|||
end; |
|||
end; |
|||
const |
|||
ns: TArray<UInt64> = [61, 109, 181, 277]; |
|||
fmt = 'x^2 - %3d*y^2 = 1 for x = %-21s and y = %s'; |
|||
begin |
|||
for var n in ns do |
|||
with SolvePell(n) do |
|||
writeln(format(fmt, [n, x.ToString, y.ToString])); |
|||
{$IFNDEF UNIX} readln; {$ENDIF} |
|||
end.</syntaxhighlight> |
|||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
{{trans|Sidef}} |
{{trans|Sidef}} |
||
< |
<syntaxhighlight lang="factor">USING: formatting kernel locals math math.functions sequences ; |
||
:: solve-pell ( n -- a b ) |
:: solve-pell ( n -- a b ) |
||
Line 317: | Line 536: | ||
dup solve-pell |
dup solve-pell |
||
"x^2 - %3d*y^2 = 1 for x = %-21d and y = %d\n" printf |
"x^2 - %3d*y^2 = 1 for x = %-21d and y = %d\n" printf |
||
] each</ |
] each</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 329: | Line 548: | ||
{{trans|Visual Basic .NET}} |
{{trans|Visual Basic .NET}} |
||
'''for n = 277 the result is wrong, I do not know if you can represent such large numbers in FreeBasic!''' |
'''for n = 277 the result is wrong, I do not know if you can represent such large numbers in FreeBasic!''' |
||
< |
<syntaxhighlight lang="freebasic"> |
||
Sub Fun(Byref a As LongInt, Byref b As LongInt, c As Integer) |
Sub Fun(Byref a As LongInt, Byref b As LongInt, c As Integer) |
||
Dim As LongInt t |
Dim As LongInt t |
||
Line 354: | Line 573: | ||
Print Using "x^2 - ### * y^2 = 1 for x = ##################### and y = #####################"; n(i); x; y |
Print Using "x^2 - ### * y^2 = 1 for x = ##################### and y = #####################"; n(i); x; y |
||
Next i |
Next i |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 366: | Line 585: | ||
=={{header|Go}}== |
=={{header|Go}}== |
||
{{trans|Sidef}} |
{{trans|Sidef}} |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 427: | Line 646: | ||
fmt.Printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", n, x, y) |
fmt.Printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", n, x, y) |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 435: | Line 654: | ||
x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 |
x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 |
||
x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020 |
x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020 |
||
</pre> |
|||
=={{header|Haskell}}== |
|||
{{trans|Julia}} |
|||
<syntaxhighlight lang="haskell">pell :: Integer -> (Integer, Integer) |
|||
pell n = go (x, 1, x * 2, 1, 0, 0, 1) |
|||
where |
|||
x = floor $ sqrt $ fromIntegral n |
|||
go (y, z, r, e1, e2, f1, f2) = |
|||
let y' = r * z - y |
|||
z' = (n - y' * y') `div` z |
|||
r' = (x + y') `div` z' |
|||
(e1', e2') = (e2, e2 * r' + e1) |
|||
(f1', f2') = (f2, f2 * r' + f1) |
|||
(a, b) = (f2', e2') |
|||
(b', a') = (a, a * x + b) |
|||
in if a' * a' - n * b' * b' == 1 |
|||
then (a', b') |
|||
else go (y', z', r', e1', e2', f1', f2')</syntaxhighlight> |
|||
<pre>λ> mapM_ print $ pell <$> [61,109,181,277] |
|||
(1766319049,226153980) |
|||
(158070671986249,15140424455100) |
|||
(2469645423824185801,183567298683461940) |
|||
(159150073798980475849,9562401173878027020)</pre> |
|||
=={{header|J}}== |
|||
<syntaxhighlight lang="j">NB. sqrt representation for continued fraction |
|||
sqrt_cf =: 3 : 0 |
|||
rep=. '' [ 'm d'=. 0 1 [ a =. a0=. <. %: y |
|||
while. a ~: +: a0 do. |
|||
rep=. rep , a=. <. (a0+m) % d=. d %~ y - *: m=. m -~ a*d |
|||
end. a0;rep |
|||
) |
|||
NB. find x,y such that x^2 - n*y^2 = 1 using continued fractions |
|||
pell =: 3 : 0 |
|||
n =. 1 [ 'a0 as' =. x: &.> sqrt_cf y |
|||
while. 1 do. cs =. 2 x: (+%)/\ a0, n$as NB. convergents |
|||
if. # sols =. I. 1 = (*: cs) +/ . * 1 , -y do. cs {~ {. sols return. end. |
|||
n =. +: n |
|||
end. |
|||
) |
|||
</syntaxhighlight> |
|||
Check that task is actually solved |
|||
<syntaxhighlight lang="j">verify =: 3 : 0 |
|||
assert. 1 = (*: xy) +/ . * 1 _61 [ echo 61 ; xy =. pell 61 |
|||
assert. 1 = (*: xy) +/ . * 1 _109 [ echo 109 ; xy =. pell 109 |
|||
assert. 1 = (*: xy) +/ . * 1 _181 [ echo 181 ; xy =. pell 181 |
|||
assert. 1 = (*: xy) +/ . * 1 _277 [ echo 277 ; xy =. pell 277 |
|||
) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> verify '' |
|||
┌──┬────────────────────┐ |
|||
│61│1766319049 226153980│ |
|||
└──┴────────────────────┘ |
|||
┌───┬──────────────────────────────┐ |
|||
│109│158070671986249 15140424455100│ |
|||
└───┴──────────────────────────────┘ |
|||
┌───┬──────────────────────────────────────┐ |
|||
│181│2469645423824185801 183567298683461940│ |
|||
└───┴──────────────────────────────────────┘ |
|||
┌───┬─────────────────────────────────────────┐ |
|||
│277│159150073798980475849 9562401173878027020│ |
|||
└───┴─────────────────────────────────────────┘ |
|||
</pre> |
</pre> |
||
=={{header|Java}}== |
=={{header|Java}}== |
||
< |
<syntaxhighlight lang="java"> |
||
import java.math.BigInteger; |
import java.math.BigInteger; |
||
import java.text.NumberFormat; |
import java.text.NumberFormat; |
||
Line 524: | Line 811: | ||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
Line 547: | Line 834: | ||
x = 2,565,007,112,872,132,129,669,406,439,503,954,211,359,492,684,749,762,901,360,167,370,740,763,715,001,557,789,090,674,216,330,243,703,833,040,774,221,628,256,858,633,287,876,949,448,689,668,281,446,637,464,359,482,677,366,420,261,407,112,316,649,010,675,881,349,744,201 |
x = 2,565,007,112,872,132,129,669,406,439,503,954,211,359,492,684,749,762,901,360,167,370,740,763,715,001,557,789,090,674,216,330,243,703,833,040,774,221,628,256,858,633,287,876,949,448,689,668,281,446,637,464,359,482,677,366,420,261,407,112,316,649,010,675,881,349,744,201 |
||
y = 27,126,610,172,119,035,540,864,542,981,075,550,089,190,381,938,849,116,323,732,855,930,990,771,728,447,597,698,969,628,164,719,475,714,805,646,913,222,890,277,024,408,337,458,564,351,161,990,641,948,210,581,361,708,373,955,113,191,451,102,494,265,278,824,127,994,180 |
y = 27,126,610,172,119,035,540,864,542,981,075,550,089,190,381,938,849,116,323,732,855,930,990,771,728,447,597,698,969,628,164,719,475,714,805,646,913,222,890,277,024,408,337,458,564,351,161,990,641,948,210,581,361,708,373,955,113,191,451,102,494,265,278,824,127,994,180 |
||
</pre> |
|||
=={{header|jq}}== |
|||
{{trans|Wren}} |
|||
'''Works with gojq, the Go implementation of jq''' |
|||
'''Preliminaries''' |
|||
<syntaxhighlight lang="jq"># If $j is 0, then an error condition is raised; |
|||
# otherwise, assuming infinite-precision integer arithmetic, |
|||
# if the input and $j are integers, then the result will be an integer. |
|||
def idivide($i; $j): |
|||
($i % $j) as $mod |
|||
| ($i - $mod) / $j ; |
|||
def idivide($j): |
|||
idivide(.; $j); |
|||
# input should be a non-negative integer for accuracy |
|||
# but may be any non-negative finite number |
|||
def isqrt: |
|||
def irt: |
|||
. as $x |
|||
| 1 | until(. > $x; . * 4) as $q |
|||
| {$q, $x, r: 0} |
|||
| until( .q <= 1; |
|||
.q |= idivide(4) |
|||
| .t = .x - .r - .q |
|||
| .r |= idivide(2) |
|||
| if .t >= 0 |
|||
then .x = .t |
|||
| .r += .q |
|||
else . |
|||
end) |
|||
| .r ; |
|||
if type == "number" and (isinfinite|not) and (isnan|not) and . >= 0 |
|||
then irt |
|||
else "isqrt requires a non-negative integer for accuracy" | error |
|||
end ;</syntaxhighlight> |
|||
'''The Task''' |
|||
<syntaxhighlight lang="jq">def solvePell: |
|||
. as $n |
|||
| ($n|isqrt) as $x |
|||
| { $x, |
|||
y : $x, |
|||
z : 1, |
|||
r : ($x * 2), |
|||
v1 : 1, |
|||
v2 : 0, |
|||
f1 : 0, |
|||
f2 : 1 } |
|||
| until(.emit; |
|||
.y = .r*.z - .y |
|||
| .z = idivide($n - .y*.y; .z) |
|||
| .r = idivide(.x + .y; .z) |
|||
| .v1 as $t |
|||
| .v1 = .v2 |
|||
| .v2 = .r*.v2 + $t |
|||
| .f1 as $t |
|||
| .f1 = .f2 |
|||
| .f2 = .r*.f2 + $t |
|||
| (.v2 + .x*.f2) as $a |
|||
| .f2 as $b |
|||
| if ($a*$a - $n*$b*$b == 1) then .emit = [$a, $b] else . end |
|||
).emit ; |
|||
(61, 109, 181, 277) |
|||
| solvePell as $res |
|||
| "x² - \(.)y² = 1 for x = \($res[0]) and y = \($res[1])"</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
x² - 61y² = 1 for x = 1766319049 and y = 226153980 |
|||
x² - 109y² = 1 for x = 158070671986249 and y = 15140424455100 |
|||
x² - 181y² = 1 for x = 2469645423824185801 and y = 183567298683461940 |
|||
x² - 277y² = 1 for x = 159150073798980475849 and y = 9562401173878027020 |
|||
</pre> |
</pre> |
||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
{{trans|C#}} |
{{trans|C#}} |
||
< |
<syntaxhighlight lang="julia">function pell(n) |
||
x = BigInt(floor(sqrt(n))) |
x = BigInt(floor(sqrt(n))) |
||
y, z, r = x, BigInt(1), x << 1 |
y, z, r = x, BigInt(1), x << 1 |
||
Line 573: | Line 933: | ||
println("x\u00b2 - $target", "y\u00b2 = 1 for x = $x and y = $y") |
println("x\u00b2 - $target", "y\u00b2 = 1 for x = $x and y = $y") |
||
end |
end |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
x² - 61y² = 1 for x = 1766319049 and y = 226153980 |
x² - 61y² = 1 for x = 1766319049 and y = 226153980 |
||
Line 583: | Line 943: | ||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
{{trans|C#}} |
{{trans|C#}} |
||
< |
<syntaxhighlight lang="scala">import java.math.BigInteger |
||
import kotlin.math.sqrt |
import kotlin.math.sqrt |
||
Line 652: | Line 1,012: | ||
println("x^2 - %3d * y^2 = 1 for x = %,27d and y = %,25d".format(it, x.value, y.value)) |
println("x^2 - %3d * y^2 = 1 for x = %,27d and y = %,25d".format(it, x.value, y.value)) |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 |
<pre>x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 |
||
Line 658: | Line 1,018: | ||
x^2 - 181 * y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 |
x^2 - 181 * y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 |
||
x^2 - 277 * y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020</pre> |
x^2 - 277 * y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020</pre> |
||
=={{header|Lambdatalk}}== |
|||
Computing big numbers requires the lib_BN library. |
|||
<syntaxhighlight lang="Scheme"> |
|||
{def pell |
|||
{lambda {:n} |
|||
{let { {:n :n} |
|||
{:x {BN.intPart {BN.sqrt :n}}} // x=int(sqrt(n)) |
|||
} {pell.r :n :x :x 1 {* 2 :x} 1 0 0 1} |
|||
}}} |
|||
-> pell |
|||
{def pell.r |
|||
{lambda {:n :x :y :z :r :e1 :e2 :f1 :f2} |
|||
{let { {:n :n} {:x :x} {:z :z} {:r :r} // no closure -> |
|||
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2} // must reassign :( |
|||
{:y {BN.- {BN.* :r :z} :y}} // y=rz-y |
|||
} {let { {:n :n} {:x :x} {:y :y} {:r :r} |
|||
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2} |
|||
{:z {BN.intPart |
|||
{BN./ {BN.- :n {BN.* :y :y}} :z}}} // z=(n-y*y)//z |
|||
} {let { {:n :n} {:x :x} {:y :y} {:z :z} |
|||
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2} |
|||
{:r {BN.intPart {BN./ {BN.+ :x :y} :z}}} // r= (x+y)//z |
|||
} {let { {:n :n} {:x :x} {:y :y} {:z :z} {:r :r} |
|||
{:e1 :e2} // e1=e2 |
|||
{:e2 {BN.+ {BN.* :r :e2} :e1}} // e2=r*e2+e1 |
|||
{:f1 :f2} // f1=f2 |
|||
{:f2 {BN.+ {BN.* :r :f2} :f1}} // f2=r*f2+f1 |
|||
} {let { {:n :n} {:x :x} {:y :y} {:z :z} {:r :r} |
|||
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2} |
|||
{:a {BN.+ :e2 {BN.* :x :f2}}} // a=e2+x*f2 |
|||
{:b :f2} // b=f2 |
|||
} {if {= {BN.compare {BN.- {BN.* :a :a} |
|||
{BN.* :n {BN.* :b :b}}} |
|||
1} |
|||
0} // a*a-n*b*b == 1 |
|||
then {div}x{sup 2} - n*y{sup 2} = 1 for n=:n, x=:a, y=:b |
|||
else {pell.r :n :x :y :z :r :e1 :e2 :f1 :f2} // do it again |
|||
}}}}}}}} |
|||
-> pell.r |
|||
{S.map pell 61 109 181 277} |
|||
-> |
|||
x^2 - n*y^2 = 1 for n=61, x=1766319049, y=226153980 |
|||
x^2 - n*y^2 = 1 for n=109, x=158070671986249, y=15140424455100 |
|||
x^2 - n*y^2 = 1 for n=181, x=2469645423824185801, y=183567298683461940 |
|||
x^2 - n*y^2 = 1 for n=277, x=159150073798980475849, y=9562401173878027020 |
|||
</syntaxhighlight> |
|||
=={{header|langur}}== |
=={{header|langur}}== |
||
{{trans|D}} |
{{trans|D}} |
||
<syntaxhighlight lang="langur">val .fun = fn(.a, .b, .c) { [.b, .b * .c + .a] } |
|||
{{works with|langur|0.10}} |
|||
Prior to 0.10, multi-variable declaration/assignment would use parentheses around variable names and values. |
|||
<lang langur>val .fun = f [.b, .b x .c + .a] |
|||
val .solvePell = |
val .solvePell = fn(.n) { |
||
val .x = |
val .x = trunc .n ^/ 2 |
||
var .y, .z, .r = .x, 1, .x |
var .y, .z, .r = .x, 1, .x * 2 |
||
var .e1, .e2, .f1, .f2 = 1, 0, 0, 1 |
var .e1, .e2, .f1, .f2 = 1, 0, 0, 1 |
||
for { |
for { |
||
.y = .r |
.y = .r * .z - .y |
||
.z = (.n - .y |
.z = (.n - .y * .y) \ .z |
||
.r = (.x + .y) \ .z |
.r = (.x + .y) \ .z |
||
.e1, .e2 = .fun(.e1, .e2, .r) |
.e1, .e2 = .fun(.e1, .e2, .r) |
||
.f1, .f2 = .fun(.f1, .f2, .r) |
.f1, .f2 = .fun(.f1, .f2, .r) |
||
val .b, .a = .fun(.e2, .f2, .x) |
val .b, .a = .fun(.e2, .f2, .x) |
||
if .a^2 - .n |
if .a^2 - .n * .b^2 == 1: return [.a, .b] |
||
} |
} |
||
} |
} |
||
val .C = |
val .C = fn(.x) { |
||
# format number string with commas |
# format number string with commas |
||
var .neg, .s = |
var .neg, .s = "", string .x |
||
if .s[1] == '-' { |
if .s[1] == '-' { |
||
.s = |
.neg, .s = "-", rest .s |
||
.neg = "-" |
|||
} |
} |
||
.neg ~ join ",", split -3, .s |
.neg ~ join ",", split -3, .s |
||
Line 694: | Line 1,099: | ||
for .n in [61, 109, 181, 277, 8941] { |
for .n in [61, 109, 181, 277, 8941] { |
||
val .x, .y = .solvePell(.n) |
val .x, .y = .solvePell(.n) |
||
writeln $"x² - \.n;y² = 1 for:\n\tx = \.x: |
writeln $"x² - \.n;y² = 1 for:\n\tx = \.x:fn C;\n\ty = \.y:fn C;\n" |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
Line 718: | Line 1,123: | ||
x = 2,565,007,112,872,132,129,669,406,439,503,954,211,359,492,684,749,762,901,360,167,370,740,763,715,001,557,789,090,674,216,330,243,703,833,040,774,221,628,256,858,633,287,876,949,448,689,668,281,446,637,464,359,482,677,366,420,261,407,112,316,649,010,675,881,349,744,201 |
x = 2,565,007,112,872,132,129,669,406,439,503,954,211,359,492,684,749,762,901,360,167,370,740,763,715,001,557,789,090,674,216,330,243,703,833,040,774,221,628,256,858,633,287,876,949,448,689,668,281,446,637,464,359,482,677,366,420,261,407,112,316,649,010,675,881,349,744,201 |
||
y = 27,126,610,172,119,035,540,864,542,981,075,550,089,190,381,938,849,116,323,732,855,930,990,771,728,447,597,698,969,628,164,719,475,714,805,646,913,222,890,277,024,408,337,458,564,351,161,990,641,948,210,581,361,708,373,955,113,191,451,102,494,265,278,824,127,994,180 |
y = 27,126,610,172,119,035,540,864,542,981,075,550,089,190,381,938,849,116,323,732,855,930,990,771,728,447,597,698,969,628,164,719,475,714,805,646,913,222,890,277,024,408,337,458,564,351,161,990,641,948,210,581,361,708,373,955,113,191,451,102,494,265,278,824,127,994,180 |
||
</pre> |
|||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
|||
<syntaxhighlight lang="mathematica">FindInstance[x^2 - 61 y^2 == 1, {x, y}, PositiveIntegers] |
|||
FindInstance[x^2 - 109 y^2 == 1, {x, y}, PositiveIntegers] |
|||
FindInstance[x^2 - 181 y^2 == 1, {x, y}, PositiveIntegers] |
|||
FindInstance[x^2 - 277 y^2 == 1, {x, y}, PositiveIntegers]</syntaxhighlight> |
|||
{{out}} |
|||
<pre>{{x -> 1766319049, y -> 226153980}} |
|||
{{x -> 158070671986249, y -> 15140424455100}} |
|||
{{x -> 2469645423824185801, y -> 183567298683461940}} |
|||
{{x -> 159150073798980475849, y -> 9562401173878027020}}</pre> |
|||
=={{header|Nim}}== |
|||
{{trans|Python}} |
|||
{{libheader|bignum}} |
|||
<syntaxhighlight lang="nim">import math, strformat |
|||
import bignum |
|||
func solvePell(n: int): (Int, Int) = |
|||
let x = newInt(sqrt(n.toFloat).int) |
|||
var (y, z, r) = (x, newInt(1), x shl 1) |
|||
var (e1, e2) = (newInt(1), newInt(0)) |
|||
var (f1, f2) = (newInt(0), newInt(1)) |
|||
while true: |
|||
y = r * z - y |
|||
z = (n - y * y) div z |
|||
r = (x + y) div z |
|||
(e1, e2) = (e2, e1 + e2 * r) |
|||
(f1, f2) = (f2, f1 + f2 * r) |
|||
let (a, b) = (f2 * x + e2, f2) |
|||
if a * a - n * b * b == 1: |
|||
return (a, b) |
|||
for n in [61, 109, 181, 277]: |
|||
let (x, y) = solvePell(n) |
|||
echo &"x² - {n:3} * y² = 1 for (x, y) = ({x:>21}, {y:>19})"</syntaxhighlight> |
|||
{{out}} |
|||
<pre>x² - 61 * y² = 1 for (x, y) = ( 1766319049, 226153980) |
|||
x² - 109 * y² = 1 for (x, y) = ( 158070671986249, 15140424455100) |
|||
x² - 181 * y² = 1 for (x, y) = ( 2469645423824185801, 183567298683461940) |
|||
x² - 277 * y² = 1 for (x, y) = (159150073798980475849, 9562401173878027020)</pre> |
|||
=={{header|Pascal}}== |
|||
{{libheader|IntXLib4Pascal}} |
|||
A console application in Free Pascal, created with the Lazarus IDE. |
|||
Pascal has no built-in support for arbitrarily large integers. The program below requires integers larger than 64 bits, and therefore uses an external library. The code could easily be adapted to solve the negative Pell equation x^2 - n*y^2 = -1, or show that no solution exists. |
|||
<syntaxhighlight lang="pascal"> |
|||
program Pell_console; |
|||
uses SysUtils, |
|||
uIntX; // uIntX is a unit in the library IntXLib4Pascal. |
|||
// uIntX.TIntX is an arbitrarily large integer. |
|||
// For the given n: if there are non-trivial solutions of x^2 - n*y^2 = 1 |
|||
// in non-negative integers (x,y), return the smallest. |
|||
// Else return the trivial solution (x,y) = (1,0). |
|||
procedure SolvePell( n : integer; out x, y : uIntX.TIntX); |
|||
var |
|||
m, a, c, d : integer; |
|||
p, q, p_next, q_next, p_prev, q_prev : uIntX.TIntX; |
|||
evenNrSteps : boolean; |
|||
begin |
|||
if (n >= 0) then m := Trunc( Sqrt( 1.0*n + 0.5)) // or use Rosetta Code Isqrt |
|||
else m := 0; |
|||
if n <= m*m then begin // if n is not a positive non-square |
|||
x := 1; y := 0; exit; // return a trivial solution |
|||
end; |
|||
c := m; d := 1; |
|||
p := 1; q := 0; |
|||
p_prev := 0; q_prev := 1; |
|||
a := m; |
|||
evenNrSteps := true; |
|||
repeat |
|||
// Get the next convergent p/q in the continued fraction for sqrt(n) |
|||
p_next := a*p + p_prev; |
|||
q_next := a*q + q_prev; |
|||
p_prev := p; p := p_next; |
|||
q_prev := q; q := q_next; |
|||
// Get the next term a in the continued fraction for sqrt(n) |
|||
Assert((n - c*c) mod d = 0); // optional sanity check |
|||
d := (n - c*c) div d; |
|||
a := (m + c) div d; |
|||
c := a*d - c; |
|||
evenNrSteps := not evenNrSteps; |
|||
until (c = m) and (d = 1); |
|||
{ |
|||
If the first return to (c,d) = (m,1) occurs after an even number of steps, |
|||
then p^2 - n*q^2 = 1, and there is no solution to x^2 - n*y^2 = -1. |
|||
Else p^2 - n*q^2 = -1, and to get a solution to x^2 - n*y^2 = 1 we can |
|||
either continue until we return to (c,d) = (m,1) for the second time, |
|||
or use the short cut below. |
|||
} |
|||
if evenNrSteps then begin |
|||
x := p; y := q; |
|||
end |
|||
else begin |
|||
x := 2*p*p + 1; y := 2*p*q |
|||
end; |
|||
end; |
|||
// For the given n: show the Pell solution on the console. |
|||
procedure ShowPellSolution( n : integer); |
|||
var |
|||
x, y : uIntX.TIntX; |
|||
lineOut : string; |
|||
begin |
|||
SolvePell( n, x, y); |
|||
lineOut := SysUtils.Format( 'n = %d --> (', [n]); |
|||
lineOut := lineOut + x.ToString + ', ' + y.ToString + ')'; |
|||
WriteLn( lineOut); |
|||
end; |
|||
// Main routine |
|||
begin |
|||
ShowPellSolution( 61); |
|||
ShowPellSolution( 109); |
|||
ShowPellSolution( 181); |
|||
ShowPellSolution( 277); |
|||
end. |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
n = 61 --> (1766319049, 226153980) |
|||
n = 109 --> (158070671986249, 15140424455100) |
|||
n = 181 --> (2469645423824185801, 183567298683461940) |
|||
n = 277 --> (159150073798980475849, 9562401173878027020) |
|||
</pre> |
</pre> |
||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
< |
<syntaxhighlight lang="perl">sub solve_pell { |
||
my ($n) = @_; |
my ($n) = @_; |
||
Line 755: | Line 1,291: | ||
my ($x, $y) = solve_pell($n); |
my ($x, $y) = solve_pell($n); |
||
printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", $n, $x, $y); |
printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", $n, $x, $y); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 767: | Line 1,303: | ||
{{trans|C#}} |
{{trans|C#}} |
||
{{trans|Go}} |
{{trans|Go}} |
||
{{libheader|mpfr}} |
{{libheader|Phix/mpfr}} |
||
This now ignores the nonsquare part of the task spec, returning {1,0}. |
This now ignores the nonsquare part of the task spec, returning {1,0}. |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>include mpfr.e |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
procedure fun(mpz a,b,t, integer c) |
|||
<span style="color: #008080;">include</span> <span style="color: #7060A8;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> |
|||
-- {a,b} = {b,c*b+a} (and t gets trashed) |
|||
mpz_set(t,a) |
|||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">fun</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpz</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">)</span> |
|||
mpz_set(a,b) |
|||
<span style="color: #000080;font-style:italic;">-- {a,b} = {b,c*b+a} (and t gets trashed)</span> |
|||
mpz_mul_si(b,b,c) |
|||
<span style="color: #7060A8;">mpz_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> |
|||
mpz_add(b,b,t) |
|||
<span style="color: #7060A8;">mpz_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> |
|||
end procedure |
|||
<span style="color: #7060A8;">mpz_mul_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpz_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
|||
function SolvePell(integer n) |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
|||
integer x = floor(sqrt(n)), y = x, z = 1, r = x*2 |
|||
mpz e1 = mpz_init(1), e2 = mpz_init(), |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">SolvePell</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
f1 = mpz_init(), f2 = mpz_init(1), |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)),</span> <span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">2</span> |
|||
t = mpz_init(0), u = mpz_init(), |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">e1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">e2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span> |
|||
a = mpz_init(1), b = mpz_init(0) |
|||
<span style="color: #000000;">f1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span> <span style="color: #000000;">f2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span> |
|||
if x*x!=n then |
|||
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">u</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span> |
|||
while mpz_cmp_si(t,1)!=0 do |
|||
<span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">b</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> |
|||
y = r*z - y |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">n</span> <span style="color: #008080;">then</span> |
|||
z = floor((n-y*y)/z) |
|||
<span style="color: #008080;">while</span> <span style="color: #7060A8;">mpz_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">0</span> <span style="color: #008080;">do</span> |
|||
r = floor((x+y)/z) |
|||
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">y</span> |
|||
fun(e1,e2,t,r) -- {e1,e2} = {e2,r*e2+e1} |
|||
<span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">((</span><span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">y</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span> |
|||
fun(f1,f2,t,r) -- {f1,f2} = {f2,r*r2+f1} |
|||
<span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">((</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span> |
|||
mpz_set(a,f2) |
|||
<span style="color: #000000;">fun</span><span style="color: #0000FF;">(</span><span style="color: #000000;">e1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">e2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- {e1,e2} = {e2,r*e2+e1}</span> |
|||
mpz_set(b,e2) |
|||
<span style="color: #000000;">fun</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- {f1,f2} = {f2,r*r2+f1}</span> |
|||
fun(b,a,t,x) -- {b,a} = {f2,x*f2+e2} |
|||
<span style="color: #7060A8;">mpz_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f2</span><span style="color: #0000FF;">)</span> |
|||
mpz_mul(t,a,a) |
|||
<span style="color: #7060A8;">mpz_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">e2</span><span style="color: #0000FF;">)</span> |
|||
mpz_mul_si(u,b,n) |
|||
<span style="color: #000000;">fun</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- {b,a} = {f2,x*f2+e2}</span> |
|||
mpz_mul(u,u,b) |
|||
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> |
|||
mpz_sub(t,t,u) -- t = a^2-n*b^2 |
|||
<span style="color: #7060A8;">mpz_mul_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">u</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
end while |
|||
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">u</span><span style="color: #0000FF;">,</span><span style="color: #000000;">u</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> |
|||
end if |
|||
<span style="color: #7060A8;">mpz_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">u</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- t = a^2-n*b^2</span> |
|||
return {a, b} |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
end function |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">}</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">split_into_chunks</span><span style="color: #0000FF;">(</span><span style="color: #004080;">string</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">one</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rest</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">one</span><span style="color: #0000FF;">]}</span> |
|||
<span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">one</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$]</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">while</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">k</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">min</span><span style="color: #0000FF;">(</span><span style="color: #000000;">l</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rest</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">k</span><span style="color: #0000FF;">])</span> |
|||
<span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$]</span> |
|||
<span style="color: #000000;">l</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">k</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #7060A8;">join</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">&</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #008000;">' '</span><span style="color: #0000FF;">,</span><span style="color: #000000;">29</span><span style="color: #0000FF;">))&</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">&</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #008000;">' '</span><span style="color: #0000FF;">,</span><span style="color: #000000;">17</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">ns</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">61</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">109</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">181</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">277</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8941</span><span style="color: #0000FF;">}</span> |
|||
sequence ns = {4, 61, 109, 181, 277, 8941} |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ns</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
for i=1 to length(ns) do |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ns</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
integer n = ns[i] |
|||
<span style="color: #004080;">mpz</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">SolvePell</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
mpz {x, y} = SolvePell(n) |
|||
<span style="color: #004080;">string</span> <span style="color: #000000;">xs</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">comma_fill</span><span style="color: #0000FF;">:=</span><span style="color: #004600;">true</span><span style="color: #0000FF;">),</span> |
|||
string xs = mpz_get_str(x,comma_fill:=true), |
|||
<span style="color: #000000;">ys</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">,</span><span style="color: #000000;">comma_fill</span><span style="color: #0000FF;">:=</span><span style="color: #004600;">true</span><span style="color: #0000FF;">)</span> |
|||
ys = mpz_get_str(y,comma_fill:=true) |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xs</span><span style="color: #0000FF;">)></span><span style="color: #000000;">97</span> <span style="color: #008080;">then</span> |
|||
printf(1,"x^2 - %3d*y^2 = 1 for x = %27s and y = %25s\n", {n, xs, ys}) |
|||
<span style="color: #000000;">xs</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">split_into_chunks</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xs</span><span style="color: #0000FF;">,</span><span style="color: #000000;">98</span><span style="color: #0000FF;">,</span><span style="color: #000000;">96</span><span style="color: #0000FF;">)</span> |
|||
end for</lang> |
|||
<span style="color: #000000;">ys</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">split_into_chunks</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ys</span><span style="color: #0000FF;">,</span><span style="color: #000000;">99</span><span style="color: #0000FF;">,</span><span style="color: #000000;">96</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"x^2 - %3d*y^2 = 1 for x = %27s and y = %25s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">xs</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ys</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 829: | Line 1,386: | ||
</pre> |
</pre> |
||
=={{header|Prolog}}== |
|||
pell(A, X, Y) finds all solutions X, Y s.t. X^2 - A*Y^2 = 1. Therefore the once() predicate can be used to only select the first one. |
|||
<syntaxhighlight lang="prolog"> |
|||
% Find the square root as a continued fraction |
|||
cf_sqrt(N, Sz, [A0, Frac]) :- |
|||
A0 is floor(sqrt(N)), |
|||
(A0*A0 =:= N -> |
|||
Sz = 0, Frac = [] |
|||
; |
|||
cf_sqrt(N, A0, A0, 0, 1, 0, [], Sz, Frac)). |
|||
cf_sqrt(N, A, A0, M0, D0, Sz0, L, Sz, R) :- |
|||
M1 is D0*A0 - M0, |
|||
D1 is (N - M1*M1) div D0, |
|||
A1 is (A + M1) div D1, |
|||
(A1 =:= 2*A -> |
|||
succ(Sz0, Sz), revtl([A1|L], R, R) |
|||
; |
|||
succ(Sz0, Sz1), cf_sqrt(N, A, A1, M1, D1, Sz1, [A1|L], Sz, R)). |
|||
revtl([], Z, Z). |
|||
revtl([A|As], Bs, Z) :- revtl(As, [A|Bs], Z). |
|||
% evaluate an infinite continued fraction as a lazy list of convergents. |
|||
% |
|||
convergents([A0, As], Lz) :- |
|||
lazy_list(next_convergent, eval_state(1, 0, A0, 1, As), Lz). |
|||
next_convergent(eval_state(P0, Q0, P1, Q1, [Term|Ts]), eval_state(P1, Q1, P2, Q2, Ts), R) :- |
|||
P2 is Term*P1 + P0, |
|||
Q2 is Term*Q1 + Q0, |
|||
R is P1 rdiv Q1. |
|||
% solve Pell's equation |
|||
% |
|||
pell(N, X, Y) :- |
|||
cf_sqrt(N, _, D), convergents(D, Rs), |
|||
once((member(R, Rs), ratio(R, P, Q), P*P - N*Q*Q =:= 1)), |
|||
pell_seq(N, P, Q, X, Y). |
|||
ratio(N, N, 1) :- integer(N). |
|||
ratio(P rdiv Q, P, Q). |
|||
pell_seq(_, X, Y, X, Y). |
|||
pell_seq(N, X0, Y0, X2, Y2) :- |
|||
pell_seq(N, X0, Y0, X1, Y1), |
|||
X2 is X0*X1 + N*Y0*Y1, |
|||
Y2 is X0*Y1 + Y0*X1. |
|||
</syntaxhighlight> |
|||
{{Out}} |
|||
<pre> |
|||
% show how we can keep generating solutions for x^2 - 3y^2 = 1 |
|||
?- pell(3,X,Y). |
|||
X = 2, |
|||
Y = 1 ; |
|||
X = 7, |
|||
Y = 4 ; |
|||
X = 26, |
|||
Y = 15 ; |
|||
X = 97, |
|||
Y = 56 ; |
|||
X = 362, |
|||
Y = 209 ; |
|||
X = 1351, |
|||
Y = 780 ; |
|||
X = 5042, |
|||
Y = 2911 . |
|||
% solve the task |
|||
?- forall((member(A, [61, 109, 181, 277]), once(pell(A, X, Y))), (write(X**2-A*Y**2=1), nl)). |
|||
1766319049**2-61*226153980**2=1 |
|||
158070671986249**2-109*15140424455100**2=1 |
|||
2469645423824185801**2-181*183567298683461940**2=1 |
|||
159150073798980475849**2-277*9562401173878027020**2=1 |
|||
true. |
|||
</pre> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
||
{{trans|D}} |
{{trans|D}} |
||
< |
<syntaxhighlight lang="python">import math |
||
def |
def solvePell(n): |
||
t = a[0] |
|||
a[0] = b[0] |
|||
b[0] = b[0] * c + t |
|||
def solvePell(n, a, b): |
|||
x = int(math.sqrt(n)) |
x = int(math.sqrt(n)) |
||
y = x |
y, z, r = x, 1, x << 1 |
||
e1, e2 = 1, 0 |
|||
f1, f2 = 0, 1 |
|||
e1 = [1] |
|||
e2 = [0] |
|||
f1 = [0] |
|||
f2 = [1] |
|||
while True: |
while True: |
||
y = r * z - y |
y = r * z - y |
||
z = |
z = (n - y * y) // z |
||
r = (x + y) // z |
r = (x + y) // z |
||
fun(e1, e2, r) |
|||
fun(f1, f2, r) |
|||
a[0] = f2[0] |
|||
b[0] = e2[0] |
|||
fun(b, a, x) |
|||
if a[0] * a[0] - n * b[0] * b[0] == 1: |
|||
return |
|||
e1, e2 = e2, e1 + e2 * r |
|||
x = [0] |
|||
f1, f2 = f2, f1 + f2 * r |
|||
y = [0] |
|||
a, b = f2 * x + e2, f2 |
|||
if a * a - n * b * b == 1: |
|||
return a, b |
|||
for n in [61, 109, 181, 277]: |
for n in [61, 109, 181, 277]: |
||
x, y = solvePell(n) |
|||
print("x^2 - %3d * y^2 = 1 for x = %27d and y = %25d" % (n, x |
print("x^2 - %3d * y^2 = 1 for x = %27d and y = %25d" % (n, x, y))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 |
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 |
||
Line 875: | Line 1,500: | ||
{{trans|Perl}} |
{{trans|Perl}} |
||
<lang |
<syntaxhighlight lang="raku" line>use Lingua::EN::Numbers; |
||
sub pell (Int $n) { |
sub pell (Int $n) { |
||
Line 907: | Line 1,532: | ||
my ($x, $y) = pell($n); |
my ($x, $y) = pell($n); |
||
printf "x² - %sy² = 1 for:\n\tx = %s\n\ty = %s\n\n", $n, |($x, $y)»., |
printf "x² - %sy² = 1 for:\n\tx = %s\n\ty = %s\n\n", $n, |($x, $y)»., |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>x² - 61y² = 1 for: |
<pre>x² - 61y² = 1 for: |
||
Line 930: | Line 1,555: | ||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
A little extra code was added to align and commatize the gihugeic numbers for readability. |
|||
<lang rexx>/*REXX program to solve Pell's equation for the smallest solution of positive integers. */ |
|||
<syntaxhighlight lang="rexx">/*REXX program to solve Pell's equation for the smallest solution of positive integers. */ |
|||
numeric digits 2200 /*ensure enough decimal digs for answer*/ |
numeric digits 2200 /*ensure enough decimal digs for answer*/ |
||
parse arg $ /*obtain optional arguments from the CL*/ |
parse arg $ /*obtain optional arguments from the CL*/ |
||
if $=='' | $=="," then $= 61 109 181 277 /*Not specified? Then use the defaults*/ |
if $=='' | $=="," then $= 61 109 181 277 /*Not specified? Then use the defaults*/ |
||
d= |
d= 28 /*used for aligning the output numbers.*/ |
||
do j=1 for words($); #= word($, j) |
do j=1 for words($); #= word($, j) /*process all the numbers in the list. */ |
||
parse value pells(#) with x y |
parse value pells(#) with x y /*extract the two values of X and Y.*/ |
||
cx= comma(x); Lcx= length(cx); cy= comma(y); Lcy= length(cy) |
|||
say 'x^2 -'right(#, max(4, length(#))) "* y^2 == 1" , |
|||
' when x='right(cx, max(d, Lcx)) " and y="right(cy, max(d, Lcy)) |
|||
end /*j*/ |
|||
end /*j*/ |
|||
exit /*stick a fork in it, we're all done. */ |
|||
exit 0 /*stick a fork in it, we're all done. */ |
|||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
comma: parse arg ?; do jc=length(?)-3 to 1 by -3; ?= insert(',', ?, jc); end; return ? |
|||
floor: procedure; parse arg x; _= x % 1; return _ - (x < 0) * (x \= _) |
floor: procedure; parse arg x; _= x % 1; return _ - (x < 0) * (x \= _) |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
Line 952: | Line 1,580: | ||
z= 1; r= x + x |
z= 1; r= x + x |
||
do until ( (e2 + x*f2)**2 - n*f2*f2) == 1 |
do until ( (e2 + x*f2)**2 - n*f2*f2) == 1 |
||
y= r*z |
y= r*z - y; z= floor( (n - y*y) / z) |
||
r= floor( (x + y ) / z) |
|||
parse value e2 r*e2 + e1 with e1 e2 |
|||
parse value |
parse value f2 r*f2 + f1 with f1 f2 |
||
parse value f2 r*f2 + f1 with f1 f2 |
|||
end /*until*/ |
end /*until*/ |
||
return e2 + x * f2 f2</ |
return e2 + x * f2 f2</syntaxhighlight> |
||
{{out|output|text= when using the default inputs:}} |
{{out|output|text= when using the default inputs:}} |
||
<pre> |
<pre> |
||
x^2 - 61 * y^2 == 1 when x= |
x^2 - 61 * y^2 == 1 when x= 1,766,319,049 and y= 226,153,980 |
||
x^2 - 109 * y^2 == 1 when x= |
x^2 - 109 * y^2 == 1 when x= 158,070,671,986,249 and y= 15,140,424,455,100 |
||
x^2 - 181 * y^2 == 1 when x= |
x^2 - 181 * y^2 == 1 when x= 2,469,645,423,824,185,801 and y= 183,567,298,683,461,940 |
||
x^2 - 277 * y^2 == 1 when x= |
x^2 - 277 * y^2 == 1 when x= 159,150,073,798,980,475,849 and y= 9,562,401,173,878,027,020 |
||
</pre> |
</pre> |
||
=={{header|Ruby}}== |
=={{header|Ruby}}== |
||
{{trans|Sidef}} |
{{trans|Sidef}} |
||
< |
<syntaxhighlight lang="ruby">def solve_pell(n) |
||
x = Integer.sqrt(n) |
x = Integer.sqrt(n) |
||
y = x |
y = x |
||
Line 989: | Line 1,616: | ||
[61, 109, 181, 277].each {|n| puts "x*x - %3s*y*y = 1 for x = %-21s and y = %s" % [n, *solve_pell(n)]} |
[61, 109, 181, 277].each {|n| puts "x*x - %3s*y*y = 1 for x = %-21s and y = %s" % [n, *solve_pell(n)]} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 998: | Line 1,625: | ||
</pre> |
</pre> |
||
=={{header|Rust}}== |
|||
<syntaxhighlight lang="rust"> |
|||
use num_bigint::{ToBigInt, BigInt}; |
|||
use num_traits::{Zero, One}; |
|||
//use std::mem::replace in the loop if you want this to be more efficient |
|||
fn main() { |
|||
test(61u64); |
|||
test(109u64); |
|||
test(181u64); |
|||
test(277u64); |
|||
} |
|||
struct Pair { |
|||
v1: BigInt, |
|||
v2: BigInt, |
|||
} |
|||
impl Pair { |
|||
pub fn make_pair(a: &BigInt, b: &BigInt) -> Pair { |
|||
Pair { |
|||
v1: a.clone(), |
|||
v2: b.clone(), |
|||
} |
|||
} |
|||
} |
|||
fn solve_pell(n: u64) -> Pair{ |
|||
let x: BigInt = ((n as f64).sqrt()).to_bigint().unwrap(); |
|||
if x.clone() * x.clone() == n.to_bigint().unwrap() { |
|||
Pair::make_pair(&One::one(), &Zero::zero()) |
|||
} else { |
|||
let mut y: BigInt = x.clone(); |
|||
let mut z: BigInt = One::one(); |
|||
let mut r: BigInt = ( &z + &z) * x.clone(); |
|||
let mut e: Pair = Pair::make_pair(&One::one(), &Zero::zero()); |
|||
let mut f: Pair = Pair::make_pair(&Zero::zero() ,&One::one()); |
|||
let mut a: BigInt = Zero::zero(); |
|||
let mut b: BigInt = Zero::zero(); |
|||
while &a * &a - n * &b * &b != One::one() { |
|||
//println!("{} {} {}", y, z, r); |
|||
y = &r * &z - &y; |
|||
z = (n - &y * &y) / &z; |
|||
r = (&x + &y) / &z; |
|||
e = Pair::make_pair(&e.v2, &(&r * &e.v2 + &e.v1)); |
|||
f = Pair::make_pair(&f.v2, &(&r * &f.v2 + &f.v1)); |
|||
a = &e.v2 + &x * &f.v2; |
|||
b = f.v2.clone(); |
|||
} |
|||
let pa = &a; |
|||
let pb = &b; |
|||
Pair::make_pair(&pa.clone(), &pb.clone()) |
|||
} |
|||
} |
|||
fn test(n: u64) { |
|||
let r: Pair = solve_pell(n); |
|||
println!("x^2 - {} * y^2 = 1 for x = {} and y = {}", n, r.v1, r.v2); |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 |
|||
x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 |
|||
x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 |
|||
x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020 |
|||
</pre> |
|||
=={{header|Scala}}== |
|||
<syntaxhighlight lang="scala">def pellFermat(n: Int): (BigInt,BigInt) = { |
|||
import scala.math.{sqrt, floor} |
|||
val x = BigInt(floor(sqrt(n)).toInt) |
|||
var i = 0 |
|||
// Use the Continued Fractions method |
|||
def converge(y:BigInt, z:BigInt, r:BigInt, e1:BigInt, e2:BigInt, f1:BigInt, f2:BigInt ) : (BigInt,BigInt) = { |
|||
val a = f2 * x + e2 |
|||
val b = f2 |
|||
if (a * a - n * b * b == 1) { |
|||
return (a, b) |
|||
} |
|||
val yh = r * z - y |
|||
val zh = (n - yh * yh) / z |
|||
val rh = (x + yh) / zh |
|||
converge(yh,zh,rh,e2,e1 + e2 * rh,f2,f1 + f2 * rh) |
|||
} |
|||
converge(x,BigInt("1"),x << 1,BigInt("1"),BigInt("0"),BigInt("0"),BigInt("1")) |
|||
} |
|||
val nums = List(61,109,181,277) |
|||
val solutions = nums.map{pellFermat(_)} |
|||
{ |
|||
println("For Pell's Equation, x\u00b2 - ny\u00b2 = 1\n") |
|||
(nums zip solutions).foreach{ case (n, (x,y)) => println(s"n = $n, x = $x, y = $y")} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>For Pell's Equation, x² - ny² = 1 |
|||
n = 61, x = 1766319049, y = 226153980 |
|||
n = 109, x = 158070671986249, y = 15140424455100 |
|||
n = 181, x = 2469645423824185801, y = 183567298683461940 |
|||
n = 277, x = 159150073798980475849, y = 9562401173878027020</pre> |
|||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
< |
<syntaxhighlight lang="ruby">func solve_pell(n) { |
||
var x = n.isqrt |
var x = n.isqrt |
||
Line 1,031: | Line 1,770: | ||
var (x, y) = solve_pell(n) |
var (x, y) = solve_pell(n) |
||
printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", n, x, y) |
printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", n, x, y) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,045: | Line 1,784: | ||
{{libheader|AttaSwift's BigInt}} |
{{libheader|AttaSwift's BigInt}} |
||
< |
<syntaxhighlight lang="swift">func solvePell<T: BinaryInteger>(n: T, _ a: inout T, _ b: inout T) { |
||
func swap(_ a: inout T, _ b: inout T, mul by: T) { |
func swap(_ a: inout T, _ b: inout T, mul by: T) { |
||
(a, b) = (b, b * by + a) |
(a, b) = (b, b * by + a) |
||
Line 1,084: | Line 1,823: | ||
print("x\u{00b2} - \(n)y\u{00b2} = 1 for x = \(x) and y = \(y)") |
print("x\u{00b2} - \(n)y\u{00b2} = 1 for x = \(x) and y = \(y)") |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,095: | Line 1,834: | ||
=={{header|Visual Basic .NET}}== |
=={{header|Visual Basic .NET}}== |
||
{{trans|Sidef}} |
{{trans|Sidef}} |
||
< |
<syntaxhighlight lang="vbnet">Imports System.Numerics |
||
Module Module1 |
Module Module1 |
||
Line 1,119: | Line 1,858: | ||
Next |
Next |
||
End Sub |
End Sub |
||
End Module</ |
End Module</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 |
<pre>x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 |
||
Line 1,125: | Line 1,864: | ||
x^2 - 181 * y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 |
x^2 - 181 * y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 |
||
x^2 - 277 * y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020</pre> |
x^2 - 277 * y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020</pre> |
||
=={{header|Wren}}== |
|||
{{trans|Sidef}} |
|||
{{libheader|Wren-big}} |
|||
{{libheader|Wren-fmt}} |
|||
<syntaxhighlight lang="wren">import "./big" for BigInt |
|||
import "./fmt" for Fmt |
|||
var solvePell = Fn.new { |n| |
|||
n = BigInt.new(n) |
|||
var x = n.isqrt |
|||
var y = x.copy() |
|||
var z = BigInt.one |
|||
var r = x * 2 |
|||
var e1 = BigInt.one |
|||
var e2 = BigInt.zero |
|||
var f1 = BigInt.zero |
|||
var f2 = BigInt.one |
|||
while (true) { |
|||
y = r*z - y |
|||
z = (n - y*y) / z |
|||
r = (x + y) / z |
|||
var t = e1.copy() |
|||
e1 = e2.copy() |
|||
e2 = r*e2 + t |
|||
t = f1.copy() |
|||
f1 = f2.copy() |
|||
f2 = r*f2 + t |
|||
var a = e2 + x*f2 |
|||
var b = f2.copy() |
|||
if (a*a - n*b*b == BigInt.one) return [a, b] |
|||
} |
|||
} |
|||
for (n in [61, 109, 181, 277]) { |
|||
var res = solvePell.call(n) |
|||
Fmt.print("x² - $3dy² = 1 for x = $-21i and y = $i", n, res[0], res[1]) |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
x² - 61y² = 1 for x = 1766319049 and y = 226153980 |
|||
x² - 109y² = 1 for x = 158070671986249 and y = 15140424455100 |
|||
x² - 181y² = 1 for x = 2469645423824185801 and y = 183567298683461940 |
|||
x² - 277y² = 1 for x = 159150073798980475849 and y = 9562401173878027020 |
|||
</pre> |
|||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
{{libheader|GMP}} GNU Multiple Precision Arithmetic Library |
{{libheader|GMP}} GNU Multiple Precision Arithmetic Library |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="zkl">var [const] BI=Import("zklBigNum"); // libGMP |
||
fcn solve_pell(n){ |
fcn solve_pell(n){ |
||
Line 1,145: | Line 1,930: | ||
if (A*A - B*B*n == 1) return(A,B); |
if (A*A - B*B*n == 1) return(A,B); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
< |
<syntaxhighlight lang="zkl">foreach n in (T(61, 109, 181, 277)){ |
||
x,y:=solve_pell(n); |
x,y:=solve_pell(n); |
||
println("x^2 - %3d*y^2 = 1 for x = %-21d and y = %d".fmt(n,x,y)); |
println("x^2 - %3d*y^2 = 1 for x = %-21d and y = %d".fmt(n,x,y)); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
Revision as of 19:27, 1 May 2024
You are encouraged to solve this task according to the task description, using any language you may know.
Pell's equation (also called the Pell–Fermat equation) is a Diophantine equation of the form:
- x2 - ny2 = 1
with integer solutions for x and y, where n is a given non-square positive integer.
- Task requirements
-
- find the smallest solution in positive integers to Pell's equation for n = {61, 109, 181, 277}.
- See also
-
- Wikipedia entry: Pell's equation.
11l
F solvePell(n)
V x = Int(sqrt(n))
V (y, z, r) = (x, 1, x << 1)
BigInt e1 = 1
BigInt e2 = 0
BigInt f1 = 0
BigInt f2 = 1
L
y = r * z - y
z = (n - y * y) I/ z
r = (x + y) I/ z
(e1, e2) = (e2, e1 + e2 * r)
(f1, f2) = (f2, f1 + f2 * r)
V (a, b) = (f2 * x + e2, f2)
I a * a - n * b * b == 1
R (a, b)
L(n) [61, 109, 181, 277]
V (x, y) = solvePell(n)
print(‘x^2 - #3 * y^2 = 1 for x = #27 and y = #25’.format(n, x, y))
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Ada
with Ada.Text_Io;
with Ada.Numerics.Elementary_Functions;
with Ada.Numerics.Big_Numbers.Big_Integers;
procedure Pells_Equation is
use Ada.Numerics.Big_Numbers.Big_Integers;
type Pair is
record
V1, V2 : Big_Integer;
end record;
procedure Solve_Pell (N : Natural; X, Y : out Big_Integer) is
use Ada.Numerics.Elementary_Functions;
Big_N : constant Big_Integer := To_Big_Integer (N);
XX : constant Big_Integer := To_Big_Integer (Natural (Float'Floor (Sqrt (Float (N)))));
begin
if XX**2 = Big_N then
X := 1; Y := 0;
return;
end if;
declare
YY : Big_Integer := XX;
Z : Big_Integer := 1;
R : Big_Integer := 2 * XX;
E : Pair := Pair'(V1 => 1, V2 => 0);
F : Pair := Pair'(V1 => 0, V2 => 1);
begin
loop
YY := R * Z - YY;
Z := (Big_N - YY**2) / Z;
R := (XX + YY) / Z;
E := Pair'(V1 => E.V2, V2 => R * E.V2 + E.V1);
F := Pair'(V1 => F.V2, V2 => R * F.V2 + F.V1);
X := E.V2 + XX * F.V2;
Y := F.V2;
exit when X**2 - Big_N * Y**2 = 1;
end loop;
end;
end Solve_Pell;
procedure Test (N : Natural) is
package Natural_Io is new Ada.Text_Io.Integer_Io (Natural);
use Ada.Text_Io, Natural_Io;
X, Y : Big_Integer;
begin
Solve_Pell (N, X => X, Y => Y);
Put ("X**2 - ");
Put (N, Width => 3);
Put (" * Y**2 = 1 for X = ");
Put (To_String (X, Width => 22));
Put (" and Y = ");
Put (To_String (Y, Width => 20));
New_Line;
end Test;
begin
Test (61);
Test (109);
Test (181);
Test (277);
end Pells_Equation;
- Output:
X**2 - 61 * Y**2 = 1 for X = 1766319049 and Y = 226153980 X**2 - 109 * Y**2 = 1 for X = 158070671986249 and Y = 15140424455100 X**2 - 181 * Y**2 = 1 for X = 2469645423824185801 and Y = 183567298683461940 X**2 - 277 * Y**2 = 1 for X = 159150073798980475849 and Y = 9562401173878027020
ALGOL 68
Also tests for a trival solution only (if n is a perfect square only 1, 0 is solution).
BEGIN
# find solutions to Pell's eqauation: x^2 - ny^2 = 1 for integer x, y, n #
MODE BIGINT = LONG LONG INT;
MODE BIGPAIR = STRUCT( BIGINT v1, v2 );
PROC solve pell = ( INT n )BIGPAIR:
IF INT x = ENTIER( sqrt( n ) );
x * x = n
THEN
# n is a erfect square - no solution otheg than 1,0 #
BIGPAIR( 1, 0 )
ELSE
# there are non-trivial solutions #
INT y := x;
INT z := 1;
INT r := 2*x;
BIGPAIR e := BIGPAIR( 1, 0 );
BIGPAIR f := BIGPAIR( 0, 1 );
BIGINT a := 0;
BIGINT b := 0;
WHILE
y := (r*z - y);
z := ENTIER ((n - y*y) / z);
r := ENTIER ((x + y) / z);
e := BIGPAIR( v2 OF e, r * v2 OF e + v1 OF e );
f := BIGPAIR( v2 OF f, r * v2 OF f + v1 OF f );
a := (v2 OF e + x*v2 OF f);
b := v2 OF f;
a*a - n*b*b /= 1
DO SKIP OD;
BIGPAIR( a, b )
FI # solve pell # ;
# task test cases #
[]INT nv = (61, 109, 181, 277);
FOR i FROM LWB nv TO UPB nv DO
INT n = nv[ i ];
BIGPAIR r = solve pell(n);
print( ("x^2 - ", whole( n, -3 ), " * y^2 = 1 for x = ", whole( v1 OF r, -21), " and y = ", whole( v2 OF r, -21 ), newline ) )
OD
END
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Arturo
solvePell: function [n][
x: to :integer sqrt n
[y, z, r]: @[x, 1, shl x 1]
[e1, e2]: [1, 0]
[f1, f2]: [0, 1]
while [true][
y: (r * z) - y
z: (n - y * y) / z
r: (x + y) / z
[e1, e2]: @[e2, e1 + e2 * r]
[f1, f2]: @[f2, f1 + f2 * r]
[a, b]: @[e2 + f2 * x, f2]
if 1 = (a*a) - n*b*b ->
return @[a, b]
]
]
loop [61 109 181 277] 'n [
[x, y]: solvePell n
print ["x² -" n "* y² = 1 for (x,y) =" x "," y]
]
- Output:
x² - 61 * y² = 1 for (x,y) = 1766319049 , 226153980 x² - 109 * y² = 1 for (x,y) = 158070671986249 , 15140424455100 x² - 181 * y² = 1 for (x,y) = 2469645423824185801 , 183567298683461940 x² - 277 * y² = 1 for (x,y) = 159150073798980475849 , 9562401173878027020
C
For n = 277, the x value is not correct because 64-bits is not enough to represent the value.
#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
struct Pair {
uint64_t v1, v2;
};
struct Pair makePair(uint64_t a, uint64_t b) {
struct Pair r;
r.v1 = a;
r.v2 = b;
return r;
}
struct Pair solvePell(int n) {
int x = (int) sqrt(n);
if (x * x == n) {
// n is a perfect square - no solution other than 1,0
return makePair(1, 0);
} else {
// there are non-trivial solutions
int y = x;
int z = 1;
int r = 2 * x;
struct Pair e = makePair(1, 0);
struct Pair f = makePair(0, 1);
uint64_t a = 0;
uint64_t b = 0;
while (true) {
y = r * z - y;
z = (n - y * y) / z;
r = (x + y) / z;
e = makePair(e.v2, r * e.v2 + e.v1);
f = makePair(f.v2, r * f.v2 + f.v1);
a = e.v2 + x * f.v2;
b = f.v2;
if (a * a - n * b * b == 1) {
break;
}
}
return makePair(a, b);
}
}
void test(int n) {
struct Pair r = solvePell(n);
printf("x^2 - %3d * y^2 = 1 for x = %21llu and y = %21llu\n", n, r.v1, r.v2);
}
int main() {
test(61);
test(109);
test(181);
test(277);
return 0;
}
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 11576121209304062921 and y = 9562401173878027020
C++
As with the C solution, the output for n = 277 is not correct because 64-bits is not enough to represent the value.
#include <iomanip>
#include <iostream>
#include <tuple>
std::tuple<uint64_t, uint64_t> solvePell(int n) {
int x = (int)sqrt(n);
if (x * x == n) {
// n is a perfect square - no solution other than 1,0
return std::make_pair(1, 0);
}
// there are non-trivial solutions
int y = x;
int z = 1;
int r = 2 * x;
std::tuple<uint64_t, uint64_t> e = std::make_pair(1, 0);
std::tuple<uint64_t, uint64_t> f = std::make_pair(0, 1);
uint64_t a = 0;
uint64_t b = 0;
while (true) {
y = r * z - y;
z = (n - y * y) / z;
r = (x + y) / z;
e = std::make_pair(std::get<1>(e), r * std::get<1>(e) + std::get<0>(e));
f = std::make_pair(std::get<1>(f), r * std::get<1>(f) + std::get<0>(f));
a = std::get<1>(e) + x * std::get<1>(f);
b = std::get<1>(f);
if (a * a - n * b * b == 1) {
break;
}
}
return std::make_pair(a, b);
}
void test(int n) {
auto r = solvePell(n);
std::cout << "x^2 - " << std::setw(3) << n << " * y^2 = 1 for x = " << std::setw(21) << std::get<0>(r) << " and y = " << std::setw(21) << std::get<1>(r) << '\n';
}
int main() {
test(61);
test(109);
test(181);
test(277);
return 0;
}
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 11576121209304062921 and y = 9562401173878027020
C#
using System;
using System.Numerics;
static class Program
{
static void Fun(ref BigInteger a, ref BigInteger b, int c)
{
BigInteger t = a; a = b; b = b * c + t;
}
static void SolvePell(int n, ref BigInteger a, ref BigInteger b)
{
int x = (int)Math.Sqrt(n), y = x, z = 1, r = x << 1;
BigInteger e1 = 1, e2 = 0, f1 = 0, f2 = 1;
while (true)
{
y = r * z - y; z = (n - y * y) / z; r = (x + y) / z;
Fun(ref e1, ref e2, r); Fun(ref f1, ref f2, r); a = f2; b = e2; Fun(ref b, ref a, x);
if (a * a - n * b * b == 1) return;
}
}
static void Main()
{
BigInteger x, y; foreach (int n in new[] { 61, 109, 181, 277 })
{
SolvePell(n, ref x, ref y);
Console.WriteLine("x^2 - {0,3} * y^2 = 1 for x = {1,27:n0} and y = {2,25:n0}", n, x, y);
}
}
}
- Output:
x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 x^2 - 109 * y^2 = 1 for x = 158,070,671,986,249 and y = 15,140,424,455,100 x^2 - 181 * y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 x^2 - 277 * y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020
D
import std.bigint;
import std.math;
import std.stdio;
void fun(ref BigInt a, ref BigInt b, int c) {
auto t = a;
a = b;
b = b * c + t;
}
void solvePell(int n, ref BigInt a, ref BigInt b) {
int x = cast(int) sqrt(cast(real) n);
int y = x;
int z = 1;
int r = x << 1;
BigInt e1 = 1;
BigInt e2 = 0;
BigInt f1 = 0;
BigInt f2 = 1;
while (true) {
y = r * z - y;
z = (n - y * y) / z;
r = (x + y) / z;
fun(e1, e2, r);
fun(f1, f2, r);
a = f2;
b = e2;
fun(b, a, x);
if (a * a - n * b * b == 1) {
return;
}
}
}
void main() {
BigInt x, y;
foreach(n; [61, 109, 181, 277]) {
solvePell(n, x, y);
writefln("x^2 - %3d * y^2 = 1 for x = %27d and y = %25d", n, x, y);
}
}
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Delphi
program Pells_equation;
{$APPTYPE CONSOLE}
uses
System.SysUtils,
Velthuis.BigIntegers;
type
TPellResult = record
x, y: BigInteger;
end;
function SolvePell(nn: UInt64): TPellResult;
var
n, x, y, z, r, e1, e2, f1, t, u, a, b: BigInteger;
begin
n := nn;
x := nn;
x := BigInteger.Sqrt(x);
y := BigInteger(x);
z := BigInteger.One;
r := x shl 1;
e1 := BigInteger.One;
e2 := BigInteger.Zero;
f1 := BigInteger.Zero;
b := BigInteger.One;
while True do
begin
y := (r * z) - y;
z := (n - (y * y)) div z;
r := (x + y) div z;
u := BigInteger(e1);
e1 := BigInteger(e2);
e2 := (r * e2) + u;
u := BigInteger(f1);
f1 := BigInteger(b);
b := r * b + u;
a := e2 + x * b;
t := (a * a) - (n * b * b);
if t = 1 then
begin
with Result do
begin
x := BigInteger(a);
y := BigInteger(b);
end;
Break;
end;
end;
end;
const
ns: TArray<UInt64> = [61, 109, 181, 277];
fmt = 'x^2 - %3d*y^2 = 1 for x = %-21s and y = %s';
begin
for var n in ns do
with SolvePell(n) do
writeln(format(fmt, [n, x.ToString, y.ToString]));
{$IFNDEF UNIX} readln; {$ENDIF}
end.
Factor
USING: formatting kernel locals math math.functions sequences ;
:: solve-pell ( n -- a b )
n sqrt >integer :> x!
x :> y!
1 :> z!
2 x * :> r!
1 0 :> ( e1! e2! )
0 1 :> ( f1! f2! )
0 0 :> ( a! b! )
[ a sq b sq n * - 1 = ] [
r z * y - y!
n y sq - z / floor z!
x y + z / floor r!
e2 r e2 * e1 + e2! e1!
f2 r f2 * f1 + f2! f1!
e2 x f2 * + a!
f2 b!
] until
a b ;
{ 61 109 181 277 } [
dup solve-pell
"x^2 - %3d*y^2 = 1 for x = %-21d and y = %d\n" printf
] each
- Output:
x^2 - 61*y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109*y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
FreeBASIC
for n = 277 the result is wrong, I do not know if you can represent such large numbers in FreeBasic!
Sub Fun(Byref a As LongInt, Byref b As LongInt, c As Integer)
Dim As LongInt t
t = a : a = b : b = b * c + t
End Sub
Sub SolvePell(n As Integer, Byref a As LongInt, Byref b As LongInt)
Dim As Integer z, r
Dim As LongInt x, y, e1, e2, f1, f2
x = Sqr(n) : y = x : z = 1 : r = 2 * x
e1 = 1 : e2 = 0 : f1 = 0 : f2 = 1
While True
y = r * z - y : z = (n - y * y) / z : r = (x + y) / z
Fun(e1, e2, r) : Fun(f1, f2, r) : a = f2 : b = e2 : Fun(b, a, x)
If a * a - n * b * b = 1 Then Exit Sub
Wend
End Sub
Dim As Integer i
Dim As LongInt x, y
Dim As Integer n(0 To 3) = {61, 109, 181, 277}
For i = 0 To 3 ''n In {61, 109, 181, 277}
SolvePell(n(i), x, y)
Print Using "x^2 - ### * y^2 = 1 for x = ##################### and y = #####################"; n(i); x; y
Next i
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = -6870622864405488695 and y = -8884342899831524596
Go
package main
import (
"fmt"
"math/big"
)
var big1 = new(big.Int).SetUint64(1)
func solvePell(nn uint64) (*big.Int, *big.Int) {
n := new(big.Int).SetUint64(nn)
x := new(big.Int).Set(n)
x.Sqrt(x)
y := new(big.Int).Set(x)
z := new(big.Int).SetUint64(1)
r := new(big.Int).Lsh(x, 1)
e1 := new(big.Int).SetUint64(1)
e2 := new(big.Int)
f1 := new(big.Int)
f2 := new(big.Int).SetUint64(1)
t := new(big.Int)
u := new(big.Int)
a := new(big.Int)
b := new(big.Int)
for {
t.Mul(r, z)
y.Sub(t, y)
t.Mul(y, y)
t.Sub(n, t)
z.Quo(t, z)
t.Add(x, y)
r.Quo(t, z)
u.Set(e1)
e1.Set(e2)
t.Mul(r, e2)
e2.Add(t, u)
u.Set(f1)
f1.Set(f2)
t.Mul(r, f2)
f2.Add(t, u)
t.Mul(x, f2)
a.Add(e2, t)
b.Set(f2)
t.Mul(a, a)
u.Mul(n, b)
u.Mul(u, b)
t.Sub(t, u)
if t.Cmp(big1) == 0 {
return a, b
}
}
}
func main() {
ns := []uint64{61, 109, 181, 277}
for _, n := range ns {
x, y := solvePell(n)
fmt.Printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", n, x, y)
}
}
- Output:
x^2 - 61*y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109*y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Haskell
pell :: Integer -> (Integer, Integer)
pell n = go (x, 1, x * 2, 1, 0, 0, 1)
where
x = floor $ sqrt $ fromIntegral n
go (y, z, r, e1, e2, f1, f2) =
let y' = r * z - y
z' = (n - y' * y') `div` z
r' = (x + y') `div` z'
(e1', e2') = (e2, e2 * r' + e1)
(f1', f2') = (f2, f2 * r' + f1)
(a, b) = (f2', e2')
(b', a') = (a, a * x + b)
in if a' * a' - n * b' * b' == 1
then (a', b')
else go (y', z', r', e1', e2', f1', f2')
λ> mapM_ print $ pell <$> [61,109,181,277] (1766319049,226153980) (158070671986249,15140424455100) (2469645423824185801,183567298683461940) (159150073798980475849,9562401173878027020)
J
NB. sqrt representation for continued fraction
sqrt_cf =: 3 : 0
rep=. '' [ 'm d'=. 0 1 [ a =. a0=. <. %: y
while. a ~: +: a0 do.
rep=. rep , a=. <. (a0+m) % d=. d %~ y - *: m=. m -~ a*d
end. a0;rep
)
NB. find x,y such that x^2 - n*y^2 = 1 using continued fractions
pell =: 3 : 0
n =. 1 [ 'a0 as' =. x: &.> sqrt_cf y
while. 1 do. cs =. 2 x: (+%)/\ a0, n$as NB. convergents
if. # sols =. I. 1 = (*: cs) +/ . * 1 , -y do. cs {~ {. sols return. end.
n =. +: n
end.
)
Check that task is actually solved
verify =: 3 : 0
assert. 1 = (*: xy) +/ . * 1 _61 [ echo 61 ; xy =. pell 61
assert. 1 = (*: xy) +/ . * 1 _109 [ echo 109 ; xy =. pell 109
assert. 1 = (*: xy) +/ . * 1 _181 [ echo 181 ; xy =. pell 181
assert. 1 = (*: xy) +/ . * 1 _277 [ echo 277 ; xy =. pell 277
)
- Output:
verify '' ┌──┬────────────────────┐ │61│1766319049 226153980│ └──┴────────────────────┘ ┌───┬──────────────────────────────┐ │109│158070671986249 15140424455100│ └───┴──────────────────────────────┘ ┌───┬──────────────────────────────────────┐ │181│2469645423824185801 183567298683461940│ └───┴──────────────────────────────────────┘ ┌───┬─────────────────────────────────────────┐ │277│159150073798980475849 9562401173878027020│ └───┴─────────────────────────────────────────┘
Java
import java.math.BigInteger;
import java.text.NumberFormat;
import java.util.ArrayList;
import java.util.List;
public class PellsEquation {
public static void main(String[] args) {
NumberFormat format = NumberFormat.getInstance();
for ( int n : new int[] {61, 109, 181, 277, 8941} ) {
BigInteger[] pell = pellsEquation(n);
System.out.printf("x^2 - %3d * y^2 = 1 for:%n x = %s%n y = %s%n%n", n, format.format(pell[0]), format.format(pell[1]));
}
}
private static final BigInteger[] pellsEquation(int n) {
int a0 = (int) Math.sqrt(n);
if ( a0*a0 == n ) {
throw new IllegalArgumentException("ERROR 102: Invalid n = " + n);
}
List<Integer> continuedFrac = continuedFraction(n);
int count = 0;
BigInteger ajm2 = BigInteger.ONE;
BigInteger ajm1 = new BigInteger(a0 + "");
BigInteger bjm2 = BigInteger.ZERO;
BigInteger bjm1 = BigInteger.ONE;
boolean stop = (continuedFrac.size() % 2 == 1);
if ( continuedFrac.size() == 2 ) {
stop = true;
}
while ( true ) {
count++;
BigInteger bn = new BigInteger(continuedFrac.get(count) + "");
BigInteger aj = bn.multiply(ajm1).add(ajm2);
BigInteger bj = bn.multiply(bjm1).add(bjm2);
if ( stop && (count == continuedFrac.size()-2 || continuedFrac.size() == 2) ) {
return new BigInteger[] {aj, bj};
}
else if (continuedFrac.size() % 2 == 0 && count == continuedFrac.size()-2 ) {
stop = true;
}
if ( count == continuedFrac.size()-1 ) {
count = 0;
}
ajm2 = ajm1;
ajm1 = aj;
bjm2 = bjm1;
bjm1 = bj;
}
}
private static final List<Integer> continuedFraction(int n) {
List<Integer> answer = new ArrayList<Integer>();
int a0 = (int) Math.sqrt(n);
answer.add(a0);
int a = -a0;
int aStart = a;
int b = 1;
int bStart = b;
while ( true ) {
//count++;
int[] values = iterateFrac(n, a, b);
answer.add(values[0]);
a = values[1];
b = values[2];
if (a == aStart && b == bStart) break;
}
return answer;
}
// array[0] = new part of cont frac
// array[1] = new a
// array[2] = new b
private static final int[] iterateFrac(int n, int a, int b) {
int x = (int) Math.floor((b * Math.sqrt(n) - b * a)/(n - a * a));
int[] answer = new int[3];
answer[0] = x;
answer[1] = -(b * a + x *(n - a * a)) / b;
answer[2] = (n - a * a) / b;
return answer;
}
}
- Output:
x^2 - 61 * y^2 = 1 for: x = 1,766,319,049 y = 226,153,980 x^2 - 109 * y^2 = 1 for: x = 158,070,671,986,249 y = 15,140,424,455,100 x^2 - 181 * y^2 = 1 for: x = 2,469,645,423,824,185,801 y = 183,567,298,683,461,940 x^2 - 277 * y^2 = 1 for: x = 159,150,073,798,980,475,849 y = 9,562,401,173,878,027,020 x^2 - 8941 * y^2 = 1 for: x = 2,565,007,112,872,132,129,669,406,439,503,954,211,359,492,684,749,762,901,360,167,370,740,763,715,001,557,789,090,674,216,330,243,703,833,040,774,221,628,256,858,633,287,876,949,448,689,668,281,446,637,464,359,482,677,366,420,261,407,112,316,649,010,675,881,349,744,201 y = 27,126,610,172,119,035,540,864,542,981,075,550,089,190,381,938,849,116,323,732,855,930,990,771,728,447,597,698,969,628,164,719,475,714,805,646,913,222,890,277,024,408,337,458,564,351,161,990,641,948,210,581,361,708,373,955,113,191,451,102,494,265,278,824,127,994,180
jq
Works with gojq, the Go implementation of jq
Preliminaries
# If $j is 0, then an error condition is raised;
# otherwise, assuming infinite-precision integer arithmetic,
# if the input and $j are integers, then the result will be an integer.
def idivide($i; $j):
($i % $j) as $mod
| ($i - $mod) / $j ;
def idivide($j):
idivide(.; $j);
# input should be a non-negative integer for accuracy
# but may be any non-negative finite number
def isqrt:
def irt:
. as $x
| 1 | until(. > $x; . * 4) as $q
| {$q, $x, r: 0}
| until( .q <= 1;
.q |= idivide(4)
| .t = .x - .r - .q
| .r |= idivide(2)
| if .t >= 0
then .x = .t
| .r += .q
else .
end)
| .r ;
if type == "number" and (isinfinite|not) and (isnan|not) and . >= 0
then irt
else "isqrt requires a non-negative integer for accuracy" | error
end ;
The Task
def solvePell:
. as $n
| ($n|isqrt) as $x
| { $x,
y : $x,
z : 1,
r : ($x * 2),
v1 : 1,
v2 : 0,
f1 : 0,
f2 : 1 }
| until(.emit;
.y = .r*.z - .y
| .z = idivide($n - .y*.y; .z)
| .r = idivide(.x + .y; .z)
| .v1 as $t
| .v1 = .v2
| .v2 = .r*.v2 + $t
| .f1 as $t
| .f1 = .f2
| .f2 = .r*.f2 + $t
| (.v2 + .x*.f2) as $a
| .f2 as $b
| if ($a*$a - $n*$b*$b == 1) then .emit = [$a, $b] else . end
).emit ;
(61, 109, 181, 277)
| solvePell as $res
| "x² - \(.)y² = 1 for x = \($res[0]) and y = \($res[1])"
- Output:
x² - 61y² = 1 for x = 1766319049 and y = 226153980 x² - 109y² = 1 for x = 158070671986249 and y = 15140424455100 x² - 181y² = 1 for x = 2469645423824185801 and y = 183567298683461940 x² - 277y² = 1 for x = 159150073798980475849 and y = 9562401173878027020
Julia
function pell(n)
x = BigInt(floor(sqrt(n)))
y, z, r = x, BigInt(1), x << 1
e1, e2, f1, f2 = BigInt(1), BigInt(0), BigInt(0), BigInt(1)
while true
y = r * z - y
z = div(n - y * y, z)
r = div(x + y, z)
e1, e2 = e2, e2 * r + e1
f1, f2 = f2, f2 * r + f1
a, b = f2, e2
b, a = a, a * x + b
if a * a - n * b * b == 1
return a, b
end
end
end
for target in BigInt[61, 109, 181, 277]
x, y = pell(target)
println("x\u00b2 - $target", "y\u00b2 = 1 for x = $x and y = $y")
end
- Output:
x² - 61y² = 1 for x = 1766319049 and y = 226153980 x² - 109y² = 1 for x = 158070671986249 and y = 15140424455100 x² - 181y² = 1 for x = 2469645423824185801 and y = 183567298683461940 x² - 277y² = 1 for x = 159150073798980475849 and y = 9562401173878027020
Kotlin
import java.math.BigInteger
import kotlin.math.sqrt
class BIRef(var value: BigInteger) {
operator fun minus(b: BIRef): BIRef {
return BIRef(value - b.value)
}
operator fun times(b: BIRef): BIRef {
return BIRef(value * b.value)
}
override fun equals(other: Any?): Boolean {
if (this === other) return true
if (javaClass != other?.javaClass) return false
other as BIRef
if (value != other.value) return false
return true
}
override fun hashCode(): Int {
return value.hashCode()
}
override fun toString(): String {
return value.toString()
}
}
fun f(a: BIRef, b: BIRef, c: Int) {
val t = a.value
a.value = b.value
b.value = b.value * BigInteger.valueOf(c.toLong()) + t
}
fun solvePell(n: Int, a: BIRef, b: BIRef) {
val x = sqrt(n.toDouble()).toInt()
var y = x
var z = 1
var r = x shl 1
val e1 = BIRef(BigInteger.ONE)
val e2 = BIRef(BigInteger.ZERO)
val f1 = BIRef(BigInteger.ZERO)
val f2 = BIRef(BigInteger.ONE)
while (true) {
y = r * z - y
z = (n - y * y) / z
r = (x + y) / z
f(e1, e2, r)
f(f1, f2, r)
a.value = f2.value
b.value = e2.value
f(b, a, x)
if (a * a - BIRef(n.toBigInteger()) * b * b == BIRef(BigInteger.ONE)) {
return
}
}
}
fun main() {
val x = BIRef(BigInteger.ZERO)
val y = BIRef(BigInteger.ZERO)
intArrayOf(61, 109, 181, 277).forEach {
solvePell(it, x, y)
println("x^2 - %3d * y^2 = 1 for x = %,27d and y = %,25d".format(it, x.value, y.value))
}
}
- Output:
x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 x^2 - 109 * y^2 = 1 for x = 158,070,671,986,249 and y = 15,140,424,455,100 x^2 - 181 * y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 x^2 - 277 * y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020
Lambdatalk
Computing big numbers requires the lib_BN library.
{def pell
{lambda {:n}
{let { {:n :n}
{:x {BN.intPart {BN.sqrt :n}}} // x=int(sqrt(n))
} {pell.r :n :x :x 1 {* 2 :x} 1 0 0 1}
}}}
-> pell
{def pell.r
{lambda {:n :x :y :z :r :e1 :e2 :f1 :f2}
{let { {:n :n} {:x :x} {:z :z} {:r :r} // no closure ->
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2} // must reassign :(
{:y {BN.- {BN.* :r :z} :y}} // y=rz-y
} {let { {:n :n} {:x :x} {:y :y} {:r :r}
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}
{:z {BN.intPart
{BN./ {BN.- :n {BN.* :y :y}} :z}}} // z=(n-y*y)//z
} {let { {:n :n} {:x :x} {:y :y} {:z :z}
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}
{:r {BN.intPart {BN./ {BN.+ :x :y} :z}}} // r= (x+y)//z
} {let { {:n :n} {:x :x} {:y :y} {:z :z} {:r :r}
{:e1 :e2} // e1=e2
{:e2 {BN.+ {BN.* :r :e2} :e1}} // e2=r*e2+e1
{:f1 :f2} // f1=f2
{:f2 {BN.+ {BN.* :r :f2} :f1}} // f2=r*f2+f1
} {let { {:n :n} {:x :x} {:y :y} {:z :z} {:r :r}
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}
{:a {BN.+ :e2 {BN.* :x :f2}}} // a=e2+x*f2
{:b :f2} // b=f2
} {if {= {BN.compare {BN.- {BN.* :a :a}
{BN.* :n {BN.* :b :b}}}
1}
0} // a*a-n*b*b == 1
then {div}x{sup 2} - n*y{sup 2} = 1 for n=:n, x=:a, y=:b
else {pell.r :n :x :y :z :r :e1 :e2 :f1 :f2} // do it again
}}}}}}}}
-> pell.r
{S.map pell 61 109 181 277}
->
x^2 - n*y^2 = 1 for n=61, x=1766319049, y=226153980
x^2 - n*y^2 = 1 for n=109, x=158070671986249, y=15140424455100
x^2 - n*y^2 = 1 for n=181, x=2469645423824185801, y=183567298683461940
x^2 - n*y^2 = 1 for n=277, x=159150073798980475849, y=9562401173878027020
langur
val .fun = fn(.a, .b, .c) { [.b, .b * .c + .a] }
val .solvePell = fn(.n) {
val .x = trunc .n ^/ 2
var .y, .z, .r = .x, 1, .x * 2
var .e1, .e2, .f1, .f2 = 1, 0, 0, 1
for {
.y = .r * .z - .y
.z = (.n - .y * .y) \ .z
.r = (.x + .y) \ .z
.e1, .e2 = .fun(.e1, .e2, .r)
.f1, .f2 = .fun(.f1, .f2, .r)
val .b, .a = .fun(.e2, .f2, .x)
if .a^2 - .n * .b^2 == 1: return [.a, .b]
}
}
val .C = fn(.x) {
# format number string with commas
var .neg, .s = "", string .x
if .s[1] == '-' {
.neg, .s = "-", rest .s
}
.neg ~ join ",", split -3, .s
}
for .n in [61, 109, 181, 277, 8941] {
val .x, .y = .solvePell(.n)
writeln $"x² - \.n;y² = 1 for:\n\tx = \.x:fn C;\n\ty = \.y:fn C;\n"
}
- Output:
x² - 61y² = 1 for: x = 1,766,319,049 y = 226,153,980 x² - 109y² = 1 for: x = 158,070,671,986,249 y = 15,140,424,455,100 x² - 181y² = 1 for: x = 2,469,645,423,824,185,801 y = 183,567,298,683,461,940 x² - 277y² = 1 for: x = 159,150,073,798,980,475,849 y = 9,562,401,173,878,027,020 x² - 8941y² = 1 for: x = 2,565,007,112,872,132,129,669,406,439,503,954,211,359,492,684,749,762,901,360,167,370,740,763,715,001,557,789,090,674,216,330,243,703,833,040,774,221,628,256,858,633,287,876,949,448,689,668,281,446,637,464,359,482,677,366,420,261,407,112,316,649,010,675,881,349,744,201 y = 27,126,610,172,119,035,540,864,542,981,075,550,089,190,381,938,849,116,323,732,855,930,990,771,728,447,597,698,969,628,164,719,475,714,805,646,913,222,890,277,024,408,337,458,564,351,161,990,641,948,210,581,361,708,373,955,113,191,451,102,494,265,278,824,127,994,180
Mathematica/Wolfram Language
FindInstance[x^2 - 61 y^2 == 1, {x, y}, PositiveIntegers]
FindInstance[x^2 - 109 y^2 == 1, {x, y}, PositiveIntegers]
FindInstance[x^2 - 181 y^2 == 1, {x, y}, PositiveIntegers]
FindInstance[x^2 - 277 y^2 == 1, {x, y}, PositiveIntegers]
- Output:
{{x -> 1766319049, y -> 226153980}} {{x -> 158070671986249, y -> 15140424455100}} {{x -> 2469645423824185801, y -> 183567298683461940}} {{x -> 159150073798980475849, y -> 9562401173878027020}}
Nim
import math, strformat
import bignum
func solvePell(n: int): (Int, Int) =
let x = newInt(sqrt(n.toFloat).int)
var (y, z, r) = (x, newInt(1), x shl 1)
var (e1, e2) = (newInt(1), newInt(0))
var (f1, f2) = (newInt(0), newInt(1))
while true:
y = r * z - y
z = (n - y * y) div z
r = (x + y) div z
(e1, e2) = (e2, e1 + e2 * r)
(f1, f2) = (f2, f1 + f2 * r)
let (a, b) = (f2 * x + e2, f2)
if a * a - n * b * b == 1:
return (a, b)
for n in [61, 109, 181, 277]:
let (x, y) = solvePell(n)
echo &"x² - {n:3} * y² = 1 for (x, y) = ({x:>21}, {y:>19})"
- Output:
x² - 61 * y² = 1 for (x, y) = ( 1766319049, 226153980) x² - 109 * y² = 1 for (x, y) = ( 158070671986249, 15140424455100) x² - 181 * y² = 1 for (x, y) = ( 2469645423824185801, 183567298683461940) x² - 277 * y² = 1 for (x, y) = (159150073798980475849, 9562401173878027020)
Pascal
A console application in Free Pascal, created with the Lazarus IDE.
Pascal has no built-in support for arbitrarily large integers. The program below requires integers larger than 64 bits, and therefore uses an external library. The code could easily be adapted to solve the negative Pell equation x^2 - n*y^2 = -1, or show that no solution exists.
program Pell_console;
uses SysUtils,
uIntX; // uIntX is a unit in the library IntXLib4Pascal.
// uIntX.TIntX is an arbitrarily large integer.
// For the given n: if there are non-trivial solutions of x^2 - n*y^2 = 1
// in non-negative integers (x,y), return the smallest.
// Else return the trivial solution (x,y) = (1,0).
procedure SolvePell( n : integer; out x, y : uIntX.TIntX);
var
m, a, c, d : integer;
p, q, p_next, q_next, p_prev, q_prev : uIntX.TIntX;
evenNrSteps : boolean;
begin
if (n >= 0) then m := Trunc( Sqrt( 1.0*n + 0.5)) // or use Rosetta Code Isqrt
else m := 0;
if n <= m*m then begin // if n is not a positive non-square
x := 1; y := 0; exit; // return a trivial solution
end;
c := m; d := 1;
p := 1; q := 0;
p_prev := 0; q_prev := 1;
a := m;
evenNrSteps := true;
repeat
// Get the next convergent p/q in the continued fraction for sqrt(n)
p_next := a*p + p_prev;
q_next := a*q + q_prev;
p_prev := p; p := p_next;
q_prev := q; q := q_next;
// Get the next term a in the continued fraction for sqrt(n)
Assert((n - c*c) mod d = 0); // optional sanity check
d := (n - c*c) div d;
a := (m + c) div d;
c := a*d - c;
evenNrSteps := not evenNrSteps;
until (c = m) and (d = 1);
{
If the first return to (c,d) = (m,1) occurs after an even number of steps,
then p^2 - n*q^2 = 1, and there is no solution to x^2 - n*y^2 = -1.
Else p^2 - n*q^2 = -1, and to get a solution to x^2 - n*y^2 = 1 we can
either continue until we return to (c,d) = (m,1) for the second time,
or use the short cut below.
}
if evenNrSteps then begin
x := p; y := q;
end
else begin
x := 2*p*p + 1; y := 2*p*q
end;
end;
// For the given n: show the Pell solution on the console.
procedure ShowPellSolution( n : integer);
var
x, y : uIntX.TIntX;
lineOut : string;
begin
SolvePell( n, x, y);
lineOut := SysUtils.Format( 'n = %d --> (', [n]);
lineOut := lineOut + x.ToString + ', ' + y.ToString + ')';
WriteLn( lineOut);
end;
// Main routine
begin
ShowPellSolution( 61);
ShowPellSolution( 109);
ShowPellSolution( 181);
ShowPellSolution( 277);
end.
- Output:
n = 61 --> (1766319049, 226153980) n = 109 --> (158070671986249, 15140424455100) n = 181 --> (2469645423824185801, 183567298683461940) n = 277 --> (159150073798980475849, 9562401173878027020)
Perl
sub solve_pell {
my ($n) = @_;
use bigint try => 'GMP';
my $x = int(sqrt($n));
my $y = $x;
my $z = 1;
my $r = 2 * $x;
my ($e1, $e2) = (1, 0);
my ($f1, $f2) = (0, 1);
for (; ;) {
$y = $r * $z - $y;
$z = int(($n - $y * $y) / $z);
$r = int(($x + $y) / $z);
($e1, $e2) = ($e2, $r * $e2 + $e1);
($f1, $f2) = ($f2, $r * $f2 + $f1);
my $A = $e2 + $x * $f2;
my $B = $f2;
if ($A**2 - $n * $B**2 == 1) {
return ($A, $B);
}
}
}
foreach my $n (61, 109, 181, 277) {
my ($x, $y) = solve_pell($n);
printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", $n, $x, $y);
}
- Output:
x^2 - 61*y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109*y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Phix
This now ignores the nonsquare part of the task spec, returning {1,0}.
with javascript_semantics include mpfr.e procedure fun(mpz a,b,t, integer c) -- {a,b} = {b,c*b+a} (and t gets trashed) mpz_set(t,a) mpz_set(a,b) mpz_mul_si(b,b,c) mpz_add(b,b,t) end procedure function SolvePell(integer n) integer x = floor(sqrt(n)), y = x, z = 1, r = x*2 mpz e1 = mpz_init(1), e2 = mpz_init(), f1 = mpz_init(), f2 = mpz_init(1), t = mpz_init(0), u = mpz_init(), a = mpz_init(1), b = mpz_init(0) if x*x!=n then while mpz_cmp_si(t,1)!=0 do y = r*z - y z = floor((n-y*y)/z) r = floor((x+y)/z) fun(e1,e2,t,r) -- {e1,e2} = {e2,r*e2+e1} fun(f1,f2,t,r) -- {f1,f2} = {f2,r*r2+f1} mpz_set(a,f2) mpz_set(b,e2) fun(b,a,t,x) -- {b,a} = {f2,x*f2+e2} mpz_mul(t,a,a) mpz_mul_si(u,b,n) mpz_mul(u,u,b) mpz_sub(t,t,u) -- t = a^2-n*b^2 end while end if return {a, b} end function function split_into_chunks(string x, integer one, rest) sequence res = {x[1..one]} x = x[one+1..$] integer l = length(x) while l do integer k = min(l,rest) res = append(res,x[1..k]) x = x[k+1..$] l -= k end while return join(res,"\n"&repeat(' ',29))&"\n"&repeat(' ',17) end function sequence ns = {4, 61, 109, 181, 277, 8941} for i=1 to length(ns) do integer n = ns[i] mpz {x, y} = SolvePell(n) string xs = mpz_get_str(x,comma_fill:=true), ys = mpz_get_str(y,comma_fill:=true) if length(xs)>97 then xs = split_into_chunks(xs,98,96) ys = split_into_chunks(ys,99,96) end if printf(1,"x^2 - %3d*y^2 = 1 for x = %27s and y = %25s\n", {n, xs, ys}) end for
- Output:
x^2 - 4*y^2 = 1 for x = 1 and y = 0 x^2 - 61*y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 x^2 - 109*y^2 = 1 for x = 158,070,671,986,249 and y = 15,140,424,455,100 x^2 - 181*y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 x^2 - 277*y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020 x^2 - 8941*y^2 = 1 for x = 2,565,007,112,872,132,129,669,406,439,503,954,211,359,492,684,749,762, 901,360,167,370,740,763,715,001,557,789,090,674,216,330,243,703,833, 040,774,221,628,256,858,633,287,876,949,448,689,668,281,446,637,464, 359,482,677,366,420,261,407,112,316,649,010,675,881,349,744,201 and y = 27,126,610,172,119,035,540,864,542,981,075,550,089,190,381,938,849,116, 323,732,855,930,990,771,728,447,597,698,969,628,164,719,475,714,805, 646,913,222,890,277,024,408,337,458,564,351,161,990,641,948,210,581, 361,708,373,955,113,191,451,102,494,265,278,824,127,994,180
Prolog
pell(A, X, Y) finds all solutions X, Y s.t. X^2 - A*Y^2 = 1. Therefore the once() predicate can be used to only select the first one.
% Find the square root as a continued fraction
cf_sqrt(N, Sz, [A0, Frac]) :-
A0 is floor(sqrt(N)),
(A0*A0 =:= N ->
Sz = 0, Frac = []
;
cf_sqrt(N, A0, A0, 0, 1, 0, [], Sz, Frac)).
cf_sqrt(N, A, A0, M0, D0, Sz0, L, Sz, R) :-
M1 is D0*A0 - M0,
D1 is (N - M1*M1) div D0,
A1 is (A + M1) div D1,
(A1 =:= 2*A ->
succ(Sz0, Sz), revtl([A1|L], R, R)
;
succ(Sz0, Sz1), cf_sqrt(N, A, A1, M1, D1, Sz1, [A1|L], Sz, R)).
revtl([], Z, Z).
revtl([A|As], Bs, Z) :- revtl(As, [A|Bs], Z).
% evaluate an infinite continued fraction as a lazy list of convergents.
%
convergents([A0, As], Lz) :-
lazy_list(next_convergent, eval_state(1, 0, A0, 1, As), Lz).
next_convergent(eval_state(P0, Q0, P1, Q1, [Term|Ts]), eval_state(P1, Q1, P2, Q2, Ts), R) :-
P2 is Term*P1 + P0,
Q2 is Term*Q1 + Q0,
R is P1 rdiv Q1.
% solve Pell's equation
%
pell(N, X, Y) :-
cf_sqrt(N, _, D), convergents(D, Rs),
once((member(R, Rs), ratio(R, P, Q), P*P - N*Q*Q =:= 1)),
pell_seq(N, P, Q, X, Y).
ratio(N, N, 1) :- integer(N).
ratio(P rdiv Q, P, Q).
pell_seq(_, X, Y, X, Y).
pell_seq(N, X0, Y0, X2, Y2) :-
pell_seq(N, X0, Y0, X1, Y1),
X2 is X0*X1 + N*Y0*Y1,
Y2 is X0*Y1 + Y0*X1.
- Output:
% show how we can keep generating solutions for x^2 - 3y^2 = 1 ?- pell(3,X,Y). X = 2, Y = 1 ; X = 7, Y = 4 ; X = 26, Y = 15 ; X = 97, Y = 56 ; X = 362, Y = 209 ; X = 1351, Y = 780 ; X = 5042, Y = 2911 . % solve the task ?- forall((member(A, [61, 109, 181, 277]), once(pell(A, X, Y))), (write(X**2-A*Y**2=1), nl)). 1766319049**2-61*226153980**2=1 158070671986249**2-109*15140424455100**2=1 2469645423824185801**2-181*183567298683461940**2=1 159150073798980475849**2-277*9562401173878027020**2=1 true.
Python
import math
def solvePell(n):
x = int(math.sqrt(n))
y, z, r = x, 1, x << 1
e1, e2 = 1, 0
f1, f2 = 0, 1
while True:
y = r * z - y
z = (n - y * y) // z
r = (x + y) // z
e1, e2 = e2, e1 + e2 * r
f1, f2 = f2, f1 + f2 * r
a, b = f2 * x + e2, f2
if a * a - n * b * b == 1:
return a, b
for n in [61, 109, 181, 277]:
x, y = solvePell(n)
print("x^2 - %3d * y^2 = 1 for x = %27d and y = %25d" % (n, x, y))
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Raku
(formerly Perl 6)
use Lingua::EN::Numbers;
sub pell (Int $n) {
my $y = my $x = Int(sqrt $n);
my $z = 1;
my $r = 2 * $x;
my ($e1, $e2) = (1, 0);
my ($f1, $f2) = (0, 1);
loop {
$y = $r * $z - $y;
$z = Int(($n - $y²) / $z);
$r = Int(($x + $y) / $z);
($e1, $e2) = ($e2, $r * $e2 + $e1);
($f1, $f2) = ($f2, $r * $f2 + $f1);
my $A = $e2 + $x * $f2;
my $B = $f2;
if ($A² - $n * $B² == 1) {
return ($A, $B);
}
}
}
for 61, 109, 181, 277, 8941 -> $n {
next if $n.sqrt.narrow ~~ Int;
my ($x, $y) = pell($n);
printf "x² - %sy² = 1 for:\n\tx = %s\n\ty = %s\n\n", $n, |($x, $y)».,
}
- Output:
x² - 61y² = 1 for: x = 1,766,319,049 y = 226,153,980 x² - 109y² = 1 for: x = 158,070,671,986,249 y = 15,140,424,455,100 x² - 181y² = 1 for: x = 2,469,645,423,824,185,801 y = 183,567,298,683,461,940 x² - 277y² = 1 for: x = 159,150,073,798,980,475,849 y = 9,562,401,173,878,027,020 x² - 8941y² = 1 for: x = 2,565,007,112,872,132,129,669,406,439,503,954,211,359,492,684,749,762,901,360,167,370,740,763,715,001,557,789,090,674,216,330,243,703,833,040,774,221,628,256,858,633,287,876,949,448,689,668,281,446,637,464,359,482,677,366,420,261,407,112,316,649,010,675,881,349,744,201 y = 27,126,610,172,119,035,540,864,542,981,075,550,089,190,381,938,849,116,323,732,855,930,990,771,728,447,597,698,969,628,164,719,475,714,805,646,913,222,890,277,024,408,337,458,564,351,161,990,641,948,210,581,361,708,373,955,113,191,451,102,494,265,278,824,127,994,180
REXX
A little extra code was added to align and commatize the gihugeic numbers for readability.
/*REXX program to solve Pell's equation for the smallest solution of positive integers. */
numeric digits 2200 /*ensure enough decimal digs for answer*/
parse arg $ /*obtain optional arguments from the CL*/
if $=='' | $=="," then $= 61 109 181 277 /*Not specified? Then use the defaults*/
d= 28 /*used for aligning the output numbers.*/
do j=1 for words($); #= word($, j) /*process all the numbers in the list. */
parse value pells(#) with x y /*extract the two values of X and Y.*/
cx= comma(x); Lcx= length(cx); cy= comma(y); Lcy= length(cy)
say 'x^2 -'right(#, max(4, length(#))) "* y^2 == 1" ,
' when x='right(cx, max(d, Lcx)) " and y="right(cy, max(d, Lcy))
end /*j*/
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
comma: parse arg ?; do jc=length(?)-3 to 1 by -3; ?= insert(',', ?, jc); end; return ?
floor: procedure; parse arg x; _= x % 1; return _ - (x < 0) * (x \= _)
/*──────────────────────────────────────────────────────────────────────────────────────*/
iSqrt: procedure; parse arg x; r= 0; q= 1; do while q<=x; q= q * 4; end
do while q>1; q= q%4; _= x-r-q; r= r%2; if _>=0 then do; x= _; r= r+q; end; end
return r /*R: is the integer square root of X. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
pells: procedure; parse arg n; x= iSqrt(n); y=x /*obtain arg; obtain integer sqrt of N*/
parse value 1 0 with e1 e2 1 f2 f1 /*assign values for: E1, E2, and F2, F1*/
z= 1; r= x + x
do until ( (e2 + x*f2)**2 - n*f2*f2) == 1
y= r*z - y; z= floor( (n - y*y) / z)
r= floor( (x + y ) / z)
parse value e2 r*e2 + e1 with e1 e2
parse value f2 r*f2 + f1 with f1 f2
end /*until*/
return e2 + x * f2 f2
- output when using the default inputs:
x^2 - 61 * y^2 == 1 when x= 1,766,319,049 and y= 226,153,980 x^2 - 109 * y^2 == 1 when x= 158,070,671,986,249 and y= 15,140,424,455,100 x^2 - 181 * y^2 == 1 when x= 2,469,645,423,824,185,801 and y= 183,567,298,683,461,940 x^2 - 277 * y^2 == 1 when x= 159,150,073,798,980,475,849 and y= 9,562,401,173,878,027,020
Ruby
def solve_pell(n)
x = Integer.sqrt(n)
y = x
z = 1
r = 2*x
e1, e2 = 1, 0
f1, f2 = 0, 1
loop do
y = r*z - y
z = (n - y*y) / z
r = (x + y) / z
e1, e2 = e2, r*e2 + e1
f1, f2 = f2, r*f2 + f1
a, b = e2 + x*f2, f2
break a, b if a*a - n*b*b == 1
end
end
[61, 109, 181, 277].each {|n| puts "x*x - %3s*y*y = 1 for x = %-21s and y = %s" % [n, *solve_pell(n)]}
- Output:
x*x - 61*y*y = 1 for x = 1766319049 and y = 226153980 x*x - 109*y*y = 1 for x = 158070671986249 and y = 15140424455100 x*x - 181*y*y = 1 for x = 2469645423824185801 and y = 183567298683461940 x*x - 277*y*y = 1 for x = 159150073798980475849 and y = 9562401173878027020
Rust
use num_bigint::{ToBigInt, BigInt};
use num_traits::{Zero, One};
//use std::mem::replace in the loop if you want this to be more efficient
fn main() {
test(61u64);
test(109u64);
test(181u64);
test(277u64);
}
struct Pair {
v1: BigInt,
v2: BigInt,
}
impl Pair {
pub fn make_pair(a: &BigInt, b: &BigInt) -> Pair {
Pair {
v1: a.clone(),
v2: b.clone(),
}
}
}
fn solve_pell(n: u64) -> Pair{
let x: BigInt = ((n as f64).sqrt()).to_bigint().unwrap();
if x.clone() * x.clone() == n.to_bigint().unwrap() {
Pair::make_pair(&One::one(), &Zero::zero())
} else {
let mut y: BigInt = x.clone();
let mut z: BigInt = One::one();
let mut r: BigInt = ( &z + &z) * x.clone();
let mut e: Pair = Pair::make_pair(&One::one(), &Zero::zero());
let mut f: Pair = Pair::make_pair(&Zero::zero() ,&One::one());
let mut a: BigInt = Zero::zero();
let mut b: BigInt = Zero::zero();
while &a * &a - n * &b * &b != One::one() {
//println!("{} {} {}", y, z, r);
y = &r * &z - &y;
z = (n - &y * &y) / &z;
r = (&x + &y) / &z;
e = Pair::make_pair(&e.v2, &(&r * &e.v2 + &e.v1));
f = Pair::make_pair(&f.v2, &(&r * &f.v2 + &f.v1));
a = &e.v2 + &x * &f.v2;
b = f.v2.clone();
}
let pa = &a;
let pb = &b;
Pair::make_pair(&pa.clone(), &pb.clone())
}
}
fn test(n: u64) {
let r: Pair = solve_pell(n);
println!("x^2 - {} * y^2 = 1 for x = {} and y = {}", n, r.v1, r.v2);
}
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Scala
def pellFermat(n: Int): (BigInt,BigInt) = {
import scala.math.{sqrt, floor}
val x = BigInt(floor(sqrt(n)).toInt)
var i = 0
// Use the Continued Fractions method
def converge(y:BigInt, z:BigInt, r:BigInt, e1:BigInt, e2:BigInt, f1:BigInt, f2:BigInt ) : (BigInt,BigInt) = {
val a = f2 * x + e2
val b = f2
if (a * a - n * b * b == 1) {
return (a, b)
}
val yh = r * z - y
val zh = (n - yh * yh) / z
val rh = (x + yh) / zh
converge(yh,zh,rh,e2,e1 + e2 * rh,f2,f1 + f2 * rh)
}
converge(x,BigInt("1"),x << 1,BigInt("1"),BigInt("0"),BigInt("0"),BigInt("1"))
}
val nums = List(61,109,181,277)
val solutions = nums.map{pellFermat(_)}
{
println("For Pell's Equation, x\u00b2 - ny\u00b2 = 1\n")
(nums zip solutions).foreach{ case (n, (x,y)) => println(s"n = $n, x = $x, y = $y")}
}
- Output:
For Pell's Equation, x² - ny² = 1 n = 61, x = 1766319049, y = 226153980 n = 109, x = 158070671986249, y = 15140424455100 n = 181, x = 2469645423824185801, y = 183567298683461940 n = 277, x = 159150073798980475849, y = 9562401173878027020
Sidef
func solve_pell(n) {
var x = n.isqrt
var y = x
var z = 1
var r = 2*x
var (e1, e2) = (1, 0)
var (f1, f2) = (0, 1)
loop {
y = (r*z - y)
z = floor((n - y*y) / z)
r = floor((x + y) / z)
(e1, e2) = (e2, r*e2 + e1)
(f1, f2) = (f2, r*f2 + f1)
var A = (e2 + x*f2)
var B = f2
if (A**2 - n*B**2 == 1) {
return (A, B)
}
}
}
for n in [61, 109, 181, 277] {
var (x, y) = solve_pell(n)
printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", n, x, y)
}
- Output:
x^2 - 61*y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109*y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Swift
func solvePell<T: BinaryInteger>(n: T, _ a: inout T, _ b: inout T) {
func swap(_ a: inout T, _ b: inout T, mul by: T) {
(a, b) = (b, b * by + a)
}
let x = T(Double(n).squareRoot())
var y = x
var z = T(1)
var r = x << 1
var e1 = T(1)
var e2 = T(0)
var f1 = T(0)
var f2 = T(1)
while true {
y = r * z - y
z = (n - y * y) / z
r = (x + y) / z
swap(&e1, &e2, mul: r)
swap(&f1, &f2, mul: r)
(a, b) = (f2, e2)
swap(&b, &a, mul: x)
if a * a - n * b * b == 1 {
return
}
}
}
var x = BigInt(0)
var y = BigInt(0)
for n in [61, 109, 181, 277] {
solvePell(n: BigInt(n), &x, &y)
print("x\u{00b2} - \(n)y\u{00b2} = 1 for x = \(x) and y = \(y)")
}
- Output:
x² - 61y² = 1 for x = 1766319049 and y = 226153980 x² - 109y² = 1 for x = 158070671986249 and y = 15140424455100 x² - 181y² = 1 for x = 2469645423824185801 and y = 183567298683461940 x² - 277y² = 1 for x = 159150073798980475849 and y = 9562401173878027020
Visual Basic .NET
Imports System.Numerics
Module Module1
Sub Fun(ByRef a As BigInteger, ByRef b As BigInteger, c As Integer)
Dim t As BigInteger = a : a = b : b = b * c + t
End Sub
Sub SolvePell(n As Integer, ByRef a As BigInteger, ByRef b As BigInteger)
Dim x As Integer = Math.Sqrt(n), y As Integer = x, z As Integer = 1, r As Integer = x << 1,
e1 As BigInteger = 1, e2 As BigInteger = 0, f1 As BigInteger = 0, f2 As BigInteger = 1
While True
y = r * z - y : z = (n - y * y) / z : r = (x + y) / z
Fun(e1, e2, r) : Fun(f1, f2, r) : a = f2 : b = e2 : Fun(b, a, x)
If a * a - n * b * b = 1 Then Exit Sub
End While
End Sub
Sub Main()
Dim x As BigInteger, y As BigInteger
For Each n As Integer In {61, 109, 181, 277}
SolvePell(n, x, y)
Console.WriteLine("x^2 - {0,3} * y^2 = 1 for x = {1,27:n0} and y = {2,25:n0}", n, x, y)
Next
End Sub
End Module
- Output:
x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 x^2 - 109 * y^2 = 1 for x = 158,070,671,986,249 and y = 15,140,424,455,100 x^2 - 181 * y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 x^2 - 277 * y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020
Wren
import "./big" for BigInt
import "./fmt" for Fmt
var solvePell = Fn.new { |n|
n = BigInt.new(n)
var x = n.isqrt
var y = x.copy()
var z = BigInt.one
var r = x * 2
var e1 = BigInt.one
var e2 = BigInt.zero
var f1 = BigInt.zero
var f2 = BigInt.one
while (true) {
y = r*z - y
z = (n - y*y) / z
r = (x + y) / z
var t = e1.copy()
e1 = e2.copy()
e2 = r*e2 + t
t = f1.copy()
f1 = f2.copy()
f2 = r*f2 + t
var a = e2 + x*f2
var b = f2.copy()
if (a*a - n*b*b == BigInt.one) return [a, b]
}
}
for (n in [61, 109, 181, 277]) {
var res = solvePell.call(n)
Fmt.print("x² - $3dy² = 1 for x = $-21i and y = $i", n, res[0], res[1])
}
- Output:
x² - 61y² = 1 for x = 1766319049 and y = 226153980 x² - 109y² = 1 for x = 158070671986249 and y = 15140424455100 x² - 181y² = 1 for x = 2469645423824185801 and y = 183567298683461940 x² - 277y² = 1 for x = 159150073798980475849 and y = 9562401173878027020
zkl
GNU Multiple Precision Arithmetic Library
var [const] BI=Import("zklBigNum"); // libGMP
fcn solve_pell(n){
x,y,z,r := BI(n).root(2), x.copy(), BI(1), x*2;
e1,e2, f1,f2 := BI(1), BI(0), BI(0), BI(1);
reg t; // a,b = c,d is a=c; b=d
do(30_000){ // throttle this in case of screw up
y,z,r = (r*z - y), (n - y*y)/z, (x + y)/z;
t,e2,e1 = e2, r*e2 + e1, t;
t,f2,f1 = f2, r*f2 + f1, t;
A,B := e2 + x*f2, f2;
if (A*A - B*B*n == 1) return(A,B);
}
}
foreach n in (T(61, 109, 181, 277)){
x,y:=solve_pell(n);
println("x^2 - %3d*y^2 = 1 for x = %-21d and y = %d".fmt(n,x,y));
}
- Output:
x^2 - 61*y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109*y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
- Programming Tasks
- Mathematics
- 11l
- Ada
- ALGOL 68
- Arturo
- C
- C++
- C sharp
- D
- Delphi
- System.SysUtils
- Velthuis.BigIntegers
- Factor
- FreeBASIC
- Go
- Haskell
- J
- Java
- Jq
- Julia
- Kotlin
- Lambdatalk
- Langur
- Mathematica
- Wolfram Language
- Nim
- Bignum
- Pascal
- IntXLib4Pascal
- Perl
- Phix
- Phix/mpfr
- Prolog
- Python
- Raku
- REXX
- Ruby
- Rust
- Scala
- Sidef
- Swift
- AttaSwift's BigInt
- Visual Basic .NET
- Wren
- Wren-big
- Wren-fmt
- Zkl
- GMP