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.
- 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
<lang 11l>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))</lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
ALGOL 68
Also tests for a trival solution only (if n is a perfect square only 1, 0 is solution).
<lang algol68>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</lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
C
For n = 277, the x value is not correct because 64-bits is not enough to represent the value. <lang c>#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;
}</lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 11576121209304062921 and y = 9562401173878027020
C++
As with the C solution, the output for n = 277 is not correct because 64-bits is not enough to represent the value. <lang cpp>#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;
}</lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 11576121209304062921 and y = 9562401173878027020
C#
<lang csharp>using System; using System.Numerics;
static class Program {
static void Fun(ref BigInteger a, ref BigInteger b, int c) { BigInteger t = a; a = b; b = b * c + t; }
static void SolvePell(int n, ref BigInteger a, ref BigInteger b) { int x = (int)Math.Sqrt(n), y = x, z = 1, r = x << 1; BigInteger e1 = 1, e2 = 0, f1 = 0, f2 = 1; while (true) { y = r * z - y; z = (n - y * y) / z; r = (x + y) / z; Fun(ref e1, ref e2, r); Fun(ref f1, ref f2, r); a = f2; b = e2; Fun(ref b, ref a, x); if (a * a - n * b * b == 1) return; } }
static void Main() { BigInteger x, y; foreach (int n in new[] { 61, 109, 181, 277 }) { SolvePell(n, ref x, ref y); Console.WriteLine("x^2 - {0,3} * y^2 = 1 for x = {1,27:n0} and y = {2,25:n0}", n, x, y); } }
}</lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 x^2 - 109 * y^2 = 1 for x = 158,070,671,986,249 and y = 15,140,424,455,100 x^2 - 181 * y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 x^2 - 277 * y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020
D
<lang d>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); }
}</lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Delphi
<lang Delphi> 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.</lang>
Factor
<lang factor>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</lang>
- Output:
x^2 - 61*y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109*y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
FreeBASIC
for n = 277 the result is wrong, I do not know if you can represent such large numbers in FreeBasic! <lang freebasic> Sub Fun(Byref a As LongInt, Byref b As LongInt, c As Integer)
Dim As LongInt t t = a : a = b : b = b * c + t
End Sub
Sub SolvePell(n As Integer, Byref a As LongInt, Byref b As LongInt)
Dim As Integer z, r Dim As LongInt x, y, e1, e2, f1, f2 x = Sqr(n) : y = x : z = 1 : r = 2 * x e1 = 1 : e2 = 0 : f1 = 0 : f2 = 1 While True y = r * z - y : z = (n - y * y) / z : r = (x + y) / z Fun(e1, e2, r) : Fun(f1, f2, r) : a = f2 : b = e2 : Fun(b, a, x) If a * a - n * b * b = 1 Then Exit Sub Wend
End Sub
Dim As Integer i Dim As LongInt x, y Dim As Integer n(0 To 3) = {61, 109, 181, 277} For i = 0 To 3 n In {61, 109, 181, 277}
SolvePell(n(i), x, y) Print Using "x^2 - ### * y^2 = 1 for x = ##################### and y = #####################"; n(i); x; y
Next i </lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = -6870622864405488695 and y = -8884342899831524596
Go
<lang go>package main
import (
"fmt" "math/big"
)
var big1 = new(big.Int).SetUint64(1)
func solvePell(nn uint64) (*big.Int, *big.Int) {
n := new(big.Int).SetUint64(nn) x := new(big.Int).Set(n) x.Sqrt(x) y := new(big.Int).Set(x) z := new(big.Int).SetUint64(1) r := new(big.Int).Lsh(x, 1)
e1 := new(big.Int).SetUint64(1) e2 := new(big.Int) f1 := new(big.Int) f2 := new(big.Int).SetUint64(1)
t := new(big.Int) u := new(big.Int) a := new(big.Int) b := new(big.Int) for { t.Mul(r, z) y.Sub(t, y) t.Mul(y, y) t.Sub(n, t) z.Quo(t, z) t.Add(x, y) r.Quo(t, z) u.Set(e1) e1.Set(e2) t.Mul(r, e2) e2.Add(t, u) u.Set(f1) f1.Set(f2) t.Mul(r, f2) f2.Add(t, u) t.Mul(x, f2) a.Add(e2, t) b.Set(f2) t.Mul(a, a) u.Mul(n, b) u.Mul(u, b) t.Sub(t, u) if t.Cmp(big1) == 0 { return a, b } }
}
func main() {
ns := []uint64{61, 109, 181, 277} for _, n := range ns { x, y := solvePell(n) fmt.Printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", n, x, y) }
}</lang>
- Output:
x^2 - 61*y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109*y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Haskell
<lang haskell>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')</lang>
λ> mapM_ print $ pell <$> [61,109,181,277] (1766319049,226153980) (158070671986249,15140424455100) (2469645423824185801,183567298683461940) (159150073798980475849,9562401173878027020)
J
<lang 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. ) </lang>
Check that task is actually solved <lang J>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 ) </lang>
- Output:
verify '' ┌──┬────────────────────┐ │61│1766319049 226153980│ └──┴────────────────────┘ ┌───┬──────────────────────────────┐ │109│158070671986249 15140424455100│ └───┴──────────────────────────────┘ ┌───┬──────────────────────────────────────┐ │181│2469645423824185801 183567298683461940│ └───┴──────────────────────────────────────┘ ┌───┬─────────────────────────────────────────┐ │277│159150073798980475849 9562401173878027020│ └───┴─────────────────────────────────────────┘
Java
<lang 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; }
}
</lang>
- 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
<lang julia>function pell(n)
x = BigInt(floor(sqrt(n))) y, z, r = x, BigInt(1), x << 1 e1, e2, f1, f2 = BigInt(1), BigInt(0), BigInt(0), BigInt(1) while true y = r * z - y z = div(n - y * y, z) r = div(x + y, z) e1, e2 = e2, e2 * r + e1 f1, f2 = f2, f2 * r + f1 a, b = f2, e2 b, a = a, a * x + b if a * a - n * b * b == 1 return a, b end end
end
for target in BigInt[61, 109, 181, 277]
x, y = pell(target) println("x\u00b2 - $target", "y\u00b2 = 1 for x = $x and y = $y")
end
</lang>
- Output:
x² - 61y² = 1 for x = 1766319049 and y = 226153980 x² - 109y² = 1 for x = 158070671986249 and y = 15140424455100 x² - 181y² = 1 for x = 2469645423824185801 and y = 183567298683461940 x² - 277y² = 1 for x = 159150073798980475849 and y = 9562401173878027020
Kotlin
<lang scala>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)) }
}</lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 x^2 - 109 * y^2 = 1 for x = 158,070,671,986,249 and y = 15,140,424,455,100 x^2 - 181 * y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 x^2 - 277 * y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020
langur
Prior to 0.10, multi-variable declaration/assignment would use parentheses around variable names and values.
<lang langur>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"
} </lang>
- Output:
x² - 61y² = 1 for: x = 1,766,319,049 y = 226,153,980 x² - 109y² = 1 for: x = 158,070,671,986,249 y = 15,140,424,455,100 x² - 181y² = 1 for: x = 2,469,645,423,824,185,801 y = 183,567,298,683,461,940 x² - 277y² = 1 for: x = 159,150,073,798,980,475,849 y = 9,562,401,173,878,027,020 x² - 8941y² = 1 for: x = 2,565,007,112,872,132,129,669,406,439,503,954,211,359,492,684,749,762,901,360,167,370,740,763,715,001,557,789,090,674,216,330,243,703,833,040,774,221,628,256,858,633,287,876,949,448,689,668,281,446,637,464,359,482,677,366,420,261,407,112,316,649,010,675,881,349,744,201 y = 27,126,610,172,119,035,540,864,542,981,075,550,089,190,381,938,849,116,323,732,855,930,990,771,728,447,597,698,969,628,164,719,475,714,805,646,913,222,890,277,024,408,337,458,564,351,161,990,641,948,210,581,361,708,373,955,113,191,451,102,494,265,278,824,127,994,180
Nim
<lang Nim>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})"</lang>
- 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
<lang perl>sub solve_pell {
my ($n) = @_;
use bigint try => 'GMP';
my $x = int(sqrt($n)); my $y = $x; my $z = 1; my $r = 2 * $x;
my ($e1, $e2) = (1, 0); my ($f1, $f2) = (0, 1);
for (; ;) {
$y = $r * $z - $y; $z = int(($n - $y * $y) / $z); $r = int(($x + $y) / $z);
($e1, $e2) = ($e2, $r * $e2 + $e1); ($f1, $f2) = ($f2, $r * $f2 + $f1);
my $A = $e2 + $x * $f2; my $B = $f2;
if ($A**2 - $n * $B**2 == 1) { return ($A, $B); } }
}
foreach my $n (61, 109, 181, 277) {
my ($x, $y) = solve_pell($n); printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", $n, $x, $y);
}</lang>
- Output:
x^2 - 61*y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109*y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Phix
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
Python
<lang python>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))</lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Raku
(formerly Perl 6)
<lang perl6>use Lingua::EN::Numbers;
sub pell (Int $n) {
my $y = my $x = Int(sqrt $n); my $z = 1; my $r = 2 * $x;
my ($e1, $e2) = (1, 0); my ($f1, $f2) = (0, 1);
loop { $y = $r * $z - $y; $z = Int(($n - $y²) / $z); $r = Int(($x + $y) / $z);
($e1, $e2) = ($e2, $r * $e2 + $e1); ($f1, $f2) = ($f2, $r * $f2 + $f1);
my $A = $e2 + $x * $f2; my $B = $f2;
if ($A² - $n * $B² == 1) { return ($A, $B); } }
}
for 61, 109, 181, 277, 8941 -> $n {
next if $n.sqrt.narrow ~~ Int; my ($x, $y) = pell($n); printf "x² - %sy² = 1 for:\n\tx = %s\n\ty = %s\n\n", $n, |($x, $y)».,
}</lang>
- Output:
x² - 61y² = 1 for: x = 1,766,319,049 y = 226,153,980 x² - 109y² = 1 for: x = 158,070,671,986,249 y = 15,140,424,455,100 x² - 181y² = 1 for: x = 2,469,645,423,824,185,801 y = 183,567,298,683,461,940 x² - 277y² = 1 for: x = 159,150,073,798,980,475,849 y = 9,562,401,173,878,027,020 x² - 8941y² = 1 for: x = 2,565,007,112,872,132,129,669,406,439,503,954,211,359,492,684,749,762,901,360,167,370,740,763,715,001,557,789,090,674,216,330,243,703,833,040,774,221,628,256,858,633,287,876,949,448,689,668,281,446,637,464,359,482,677,366,420,261,407,112,316,649,010,675,881,349,744,201 y = 27,126,610,172,119,035,540,864,542,981,075,550,089,190,381,938,849,116,323,732,855,930,990,771,728,447,597,698,969,628,164,719,475,714,805,646,913,222,890,277,024,408,337,458,564,351,161,990,641,948,210,581,361,708,373,955,113,191,451,102,494,265,278,824,127,994,180
REXX
A little extra code was added to align and commatize the gihugeic numbers for readability. <lang 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= 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</lang>
- 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
<lang ruby>def solve_pell(n)
x = Integer.sqrt(n) y = x z = 1 r = 2*x e1, e2 = 1, 0 f1, f2 = 0, 1
loop do y = r*z - y z = (n - y*y) / z r = (x + y) / z e1, e2 = e2, r*e2 + e1 f1, f2 = f2, r*f2 + f1 a, b = e2 + x*f2, f2 break a, b if a*a - n*b*b == 1 end
end
[61, 109, 181, 277].each {|n| puts "x*x - %3s*y*y = 1 for x = %-21s and y = %s" % [n, *solve_pell(n)]} </lang>
- Output:
x*x - 61*y*y = 1 for x = 1766319049 and y = 226153980 x*x - 109*y*y = 1 for x = 158070671986249 and y = 15140424455100 x*x - 181*y*y = 1 for x = 2469645423824185801 and y = 183567298683461940 x*x - 277*y*y = 1 for x = 159150073798980475849 and y = 9562401173878027020
Rust
<lang 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);
} </lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109 * y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181 * y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Sidef
<lang ruby>func solve_pell(n) {
var x = n.isqrt var y = x var z = 1 var r = 2*x
var (e1, e2) = (1, 0) var (f1, f2) = (0, 1)
loop {
y = (r*z - y) z = floor((n - y*y) / z) r = floor((x + y) / z)
(e1, e2) = (e2, r*e2 + e1) (f1, f2) = (f2, r*f2 + f1)
var A = (e2 + x*f2) var B = f2
if (A**2 - n*B**2 == 1) { return (A, B) } }
}
for n in [61, 109, 181, 277] {
var (x, y) = solve_pell(n) printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", n, x, y)
}</lang>
- Output:
x^2 - 61*y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109*y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
Swift
<lang swift>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)")
}</lang>
- Output:
x² - 61y² = 1 for x = 1766319049 and y = 226153980 x² - 109y² = 1 for x = 158070671986249 and y = 15140424455100 x² - 181y² = 1 for x = 2469645423824185801 and y = 183567298683461940 x² - 277y² = 1 for x = 159150073798980475849 and y = 9562401173878027020
Visual Basic .NET
<lang vbnet>Imports System.Numerics
Module Module1
Sub Fun(ByRef a As BigInteger, ByRef b As BigInteger, c As Integer) Dim t As BigInteger = a : a = b : b = b * c + t End Sub
Sub SolvePell(n As Integer, ByRef a As BigInteger, ByRef b As BigInteger) Dim x As Integer = Math.Sqrt(n), y As Integer = x, z As Integer = 1, r As Integer = x << 1, e1 As BigInteger = 1, e2 As BigInteger = 0, f1 As BigInteger = 0, f2 As BigInteger = 1 While True y = r * z - y : z = (n - y * y) / z : r = (x + y) / z Fun(e1, e2, r) : Fun(f1, f2, r) : a = f2 : b = e2 : Fun(b, a, x) If a * a - n * b * b = 1 Then Exit Sub End While End Sub
Sub Main() Dim x As BigInteger, y As BigInteger For Each n As Integer In {61, 109, 181, 277} SolvePell(n, x, y) Console.WriteLine("x^2 - {0,3} * y^2 = 1 for x = {1,27:n0} and y = {2,25:n0}", n, x, y) Next End Sub
End Module</lang>
- Output:
x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980 x^2 - 109 * y^2 = 1 for x = 158,070,671,986,249 and y = 15,140,424,455,100 x^2 - 181 * y^2 = 1 for x = 2,469,645,423,824,185,801 and y = 183,567,298,683,461,940 x^2 - 277 * y^2 = 1 for x = 159,150,073,798,980,475,849 and y = 9,562,401,173,878,027,020
Wren
<lang ecmascript>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])
}</lang>
- Output:
x² - 61y² = 1 for x = 1766319049 and y = 226153980 x² - 109y² = 1 for x = 158070671986249 and y = 15140424455100 x² - 181y² = 1 for x = 2469645423824185801 and y = 183567298683461940 x² - 277y² = 1 for x = 159150073798980475849 and y = 9562401173878027020
zkl
GNU Multiple Precision Arithmetic Library
<lang zkl>var [const] BI=Import("zklBigNum"); // libGMP
fcn solve_pell(n){
x,y,z,r := BI(n).root(2), x.copy(), BI(1), x*2; e1,e2, f1,f2 := BI(1), BI(0), BI(0), BI(1); reg t; // a,b = c,d is a=c; b=d do(30_000){ // throttle this in case of screw up y,z,r = (r*z - y), (n - y*y)/z, (x + y)/z; t,e2,e1 = e2, r*e2 + e1, t; t,f2,f1 = f2, r*f2 + f1, t;
A,B := e2 + x*f2, f2;
if (A*A - B*B*n == 1) return(A,B); }
}</lang> <lang zkl>foreach n in (T(61, 109, 181, 277)){
x,y:=solve_pell(n); println("x^2 - %3d*y^2 = 1 for x = %-21d and y = %d".fmt(n,x,y));
}</lang>
- Output:
x^2 - 61*y^2 = 1 for x = 1766319049 and y = 226153980 x^2 - 109*y^2 = 1 for x = 158070671986249 and y = 15140424455100 x^2 - 181*y^2 = 1 for x = 2469645423824185801 and y = 183567298683461940 x^2 - 277*y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020