Pell's equation

Pell's equation
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:

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

Translation of: Python

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

Translation of: C

Works with: Ada 2022

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

Translation of: Sidef

Also tests for a trival solution only (if n is a perfect square only 1, 0 is solution).

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32

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

Translation of: Python

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

Translation of: ALGOL 68

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++

Translation of: 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#

Translation of: Sidef

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

Translation of: C#

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

Library: System.SysUtils

Library: Velthuis.BigIntegers

Translation of: Go

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

Translation of: Sidef

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

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!

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

$Insertformulahere$

Go

Translation of: Sidef

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

Translation of: Julia

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

Translation of: Wren

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

Translation of: C#

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

Translation of: C#

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

Translation of: D

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

Translation of: Python

Library: bignum

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

Library: 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.

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

Translation of: C#

Translation of: Go

Library: Phix/mpfr

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

Translation of: D

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)

Works with: Rakudo version 2018.12

Translation of: Perl

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

Translation of: Sidef

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

Translation of: Kotlin

Library: AttaSwift's BigInt

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

Translation of: Sidef

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

Translation of: Sidef

Library: Wren-big

Library: Wren-fmt

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

Library: GMP

GNU Multiple Precision Arithmetic Library

Translation of: Raku

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