Pell's equation

From Rosetta Code
Revision as of 14:08, 15 June 2019 by Thundergnat (talk | contribs) (→‎{{header|Perl 6}}: commas biiig ints)
Task
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

C#

Translation of: Sidef

<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

Translation of: Visual Basic .NET

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>

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

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


Julia

Translation of: C#

<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

Works with: Rakudo version 2018.12
Translation of: Perl

<lang perl6>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)».,

}</lang>

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

Phix

Translation of: C#
Translation of: Go
Library: mpfr

<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

Translation of: Sidef

<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

Translation of: Sidef

<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

Library: GMP

GNU Multiple Precision Arithmetic Library

Translation of: Perl6

<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
Retrieved from "https://rosettacode.org/wiki/Pell%27s_equation?oldid=287829"