Pell's equation

Pell's equation is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Pell's equation is a Diophantine equation of the form

x² - ny² = 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.

Go

Translation of: Sidef

<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

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

Works with: Rakudo version 2018.12

Translation of: Perl

<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

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