I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

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

x2   -   ny2   =   1

with integer solutions for   x   and   y,   where   n   is a given non-square positive integer.

•   find the smallest solution in positive integers to Pell's equation for   n = {61, 109, 181, 277}.

## 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
```

## 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 ) )    ODEND`
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

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 + tEnd 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    WendEnd Sub Dim As Integer iDim As LongInt x, yDim 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; yNext 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
```

${\displaystyle 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
```

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 fractionsqrt_cf =: 3 : 0rep=. '' [ 'm d'=. 0 1 [ a =. a0=. <. %: ywhile. a ~: +: a0 do.  rep=. rep , a=. <. (a0+m) % d=. d %~ y - *: m=. m -~ a*dend. a0;rep) NB. find x,y such that x^2 - n*y^2 = 1 using continued fractions pell =: 3 : 0n =. 1 [ 'a0 as' =. x: &.> sqrt_cf ywhile. 1 do. cs =. 2 x: (+%)/\ a0, n\$as NB. convergents  if. # sols =. I. 1 = (*: cs) +/ . * 1 , -y do. cs {~ {. sols return. end.  n =. +: nend.) `

Check that task is actually solved

`verify =: 3 : 0assert. 1 = (*: xy) +/ . * 1 _61  [ echo 61  ; xy =. pell 61assert. 1 = (*: xy) +/ . * 1 _109 [ echo 109 ; xy =. pell 109assert. 1 = (*: xy) +/ . * 1 _181 [ echo 181 ; xy =. pell 181assert. 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
```

## 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    endend 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.BigIntegerimport 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```

## langur

Translation of: D
Works with: langur version 0.10

Prior to 0.10, multi-variable declaration/assignment would use parentheses around variable names and values.

`val .fun = f [.b, .b x .c + .a] val .solvePell = f(.n) {    val .x = truncate .n ^/ 2    var .y, .z, .r = .x, 1, .x x 2    var .e1, .e2, .f1, .f2 = 1, 0, 0, 1     for {        .y = .r x .z - .y        .z = (.n - .y x .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 x .b^2 == 1: return [.a, .b]    }} val .C = f(.x) {    # format number string with commas    var .neg, .s = ZLS, toString .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:.C;\n\ty = \.y:.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, strformatimport 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)```

## 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)
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 .

?- 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)».&comma;}`
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  endend [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
```

## 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
`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 SubEnd 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 BigIntimport "/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
```