Pell's equation
You are encouraged to solve this task according to the task description, using any language you may know.
Pell's equation is a Diophantine equation of the form
- x2 - ny2 = 1
with integer solutions for x and y, where n is a given nonsquare 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.
C#
<lang csharp>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); } }
}</lang>
- 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
FreeBASIC
for n = 277 the result is wrong, I do not know if you can represent such large numbers in FreeBasic! <lang 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 </lang>
Go
<lang 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) }
}</lang>
- 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
Julia
<lang 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
</lang>
- 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
Perl
<lang 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);
}</lang>
- 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
Perl 6
<lang perl6>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 -> $n {
my ($x, $y) = pell($n); printf("x² - %3dy² = 1 for x = %-21s and y = %s\n", $n, $x, $y);
}</lang>
- 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
Phix
<lang Phix>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(), u = mpz_init(), a = mpz_init(), b = mpz_init()
while true 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 if mpz_cmp_si(t,1)==0 then return {a, b} end if end while
end function
sequence ns = {61, 109, 181, 277} 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) printf(1,"x^2 - %3d*y^2 = 1 for x = %27s and y = %25s\n", {n, xs, ys})
end for</lang>
- 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
REXX
<lang rexx>/*REXX program to solve Pell's equation for the smallest solution of positive integers. */ numeric digits 2200 /*ensure enought decimal digits for ans*/ parse arg $ /*obtain optinal arguments from the CL.*/ if $= | $="," then $= 61 109 181 277 /*Not specified? Then use the defaults*/ d=22 /*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.*/ say 'x^2 -'right(#,max(4,length(#))) "* y^2 == 1 when x="right(x, max(d,length(x))), ' and y='right(y, max(d,length(y))) end /*j*/
exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ 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*/
r= x + x parse value 1 0 with e1 e2 1 f2 f1 /*assign values for: E1, E2, F1, F2. */ z= 1 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</lang>
- output when using the default inputs:
x^2 - 61 * y^2 == 1 when x= 1766319049 and y= 226153980 x^2 - 109 * y^2 == 1 when x= 158070671986249 and y= 15140424455100 x^2 - 181 * y^2 == 1 when x= 2469645423824185801 and y= 183567298683461940 x^2 - 277 * y^2 == 1 when x= 159150073798980475849 and y= 9562401173878027020
Ruby
<lang 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)]} </lang>
- 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
Sidef
<lang ruby>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)
}</lang>
- 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
Visual Basic .NET
<lang vbnet>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</lang>
- 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
zkl
GNU Multiple Precision Arithmetic Library
<lang zkl>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); }
}</lang> <lang zkl>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));
}</lang>
- 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