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

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

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

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

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

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

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

## Phix

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

This now ignores the nonsquare part of the task spec, returning {1,0}.

`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*2mpz 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 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)    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
```

## Python

Translation of: D
`import math def fun(a, b, c):    t = a[0]    a[0] = b[0]    b[0] = b[0] * c + t def solvePell(n, a, b):    x = int(math.sqrt(n))    y = x    z = 1    r = x << 1    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[0] = f2[0]        b[0] = e2[0]        fun(b, a, x)        if a[0] * a[0] - n * b[0] * b[0] == 1:            return x = [0]y = [0]for n in [61, 109, 181, 277]:    solvePell(n, x, y)    print("x^2 - %3d * y^2 = 1 for x = %27d and y = %25d" % (n, x[0], y[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 =       159150073798980475849 and y =       9562401173878027020```

## REXX

`/*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= 22                                            /*used for aligning the output numbers.*/     do j=1  for words(\$);    #= word(\$, j)      /*process all the numbers in the list. */     parse value   pells(#)   with   x  y        /*extract the two values of  X  and  Y.*/     say 'x^2 -'right(#,max(4,length(#))) "* y^2 == 1  when x="right(x, max(d,length(x))),                                                      ' and y='right(y, max(d,length(y)))     end   /*j*/exit                                             /*stick a fork in it,  we're all done. *//*──────────────────────────────────────────────────────────────────────────────────────*/floor: procedure; parse arg x;  _= x % 1;          return  _   -    (x < 0)   *   (x \= _)/*──────────────────────────────────────────────────────────────────────────────────────*/iSqrt: procedure; parse arg x;  r= 0;     q= 1;           do  while q<=x;  q= q * 4;   end         do  while q>1; q= q%4; _= x-r-q; r= r%2; if _>=0  then do; x= _; r= r+q; end; end       return r                                  /*R:  is the integer square root of X. *//*──────────────────────────────────────────────────────────────────────────────────────*/pells: procedure; parse arg n; x= iSqrt(n);  y=x /*obtain arg;  obtain integer sqrt of N*/       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=            1766319049  and y=             226153980
x^2 - 109 * y^2 == 1  when x=       158070671986249  and y=        15140424455100
x^2 - 181 * y^2 == 1  when x=   2469645423824185801  and y=    183567298683461940
x^2 - 277 * y^2 == 1  when x= 159150073798980475849  and y=   9562401173878027020
```

## Ruby

Translation of: Sidef
`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

```

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

## zkl

Library: GMP
GNU Multiple Precision Arithmetic Library
Translation of: Perl6
`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
```