Pell's equation: Difference between revisions
Content deleted Content added
Langurmonkey (talk | contribs) |
|||
(13 intermediate revisions by 6 users not shown) | |||
Line 19:
{{trans|Python}}
<
V x = Int(sqrt(n))
V (y, z, r) = (x, 1, x << 1)
Line 40:
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))</
{{out}}
Line 52:
=={{header|Ada}}==
{{Trans|C}}{{works with|Ada 2022}}
<
with Ada.Numerics.Elementary_Functions;
with Ada.Numerics.Big_Numbers.Big_Integers;
Line 116:
Test (181);
Test (277);
end Pells_Equation;</
{{out}}
<pre>
Line 128:
{{Trans|Sidef}} Also tests for a trival solution only (if n is a perfect square only 1, 0 is solution).
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
<
# find solutions to Pell's eqauation: x^2 - ny^2 = 1 for integer x, y, n #
MODE BIGINT = LONG LONG INT;
Line 166:
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</
{{out}}
<pre>
Line 178:
{{trans|Python}}
<
x: to :integer sqrt n
[y, z, r]: @[x, 1, shl x 1]
Line 200:
[x, y]: solvePell n
print ["x² -" n "* y² = 1 for (x,y) =" x "," y]
]</
{{out}}
Line 213:
{{trans|ALGOL 68}}
For n = 277, the x value is not correct because 64-bits is not enough to represent the value.
<
#include <stdbool.h>
#include <stdint.h>
Line 274:
return 0;
}</
{{out}}
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980
Line 284:
{{trans|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 <iostream>
#include <tuple>
Line 333:
return 0;
}</
{{out}}
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980
Line 342:
=={{header|C sharp|C#}}==
{{trans|Sidef}}
<
using System.Numerics;
Line 372:
}
}
}</
{{out}}
<pre>x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980
Line 381:
=={{header|D}}==
{{trans|C#}}
<
import std.math;
import std.stdio;
Line 421:
writefln("x^2 - %3d * y^2 = 1 for x = %27d and y = %25d", n, x, y);
}
}</
{{out}}
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980
Line 431:
{{libheader| Velthuis.BigIntegers}}
{{Trans|Go}}
<syntaxhighlight lang="delphi">
program Pells_equation;
Line 501:
{$IFNDEF UNIX} readln; {$ENDIF}
end.</
=={{header|Factor}}==
{{trans|Sidef}}
<
:: solve-pell ( n -- a b )
Line 536:
dup solve-pell
"x^2 - %3d*y^2 = 1 for x = %-21d and y = %d\n" printf
] each</
{{out}}
<pre>
Line 548:
{{trans|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
Line 573:
Print Using "x^2 - ### * y^2 = 1 for x = ##################### and y = #####################"; n(i); x; y
Next i
</syntaxhighlight>
{{out}}
<pre>
Line 585:
=={{header|Go}}==
{{trans|Sidef}}
<
import (
Line 646:
fmt.Printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", n, x, y)
}
}</
{{out}}
Line 658:
=={{header|Haskell}}==
{{trans|Julia}}
<
pell n = go (x, 1, x * 2, 1, 0, 0, 1)
where
Line 672:
in if a' * a' - n * b' * b' == 1
then (a', b')
else go (y', z', r', e1', e2', f1', f2')</
<pre>λ> mapM_ print $ pell <$> [61,109,181,277]
Line 682:
=={{header|J}}==
<
sqrt_cf =: 3 : 0
rep=. '' [ 'm d'=. 0 1 [ a =. a0=. <. %: y
Line 698:
end.
)
</syntaxhighlight>
Check that task is actually solved
<
assert. 1 = (*: xy) +/ . * 1 _61 [ echo 61 ; xy =. pell 61
assert. 1 = (*: xy) +/ . * 1 _109 [ echo 109 ; xy =. pell 109
Line 707:
assert. 1 = (*: xy) +/ . * 1 _277 [ echo 277 ; xy =. pell 277
)
</syntaxhighlight>
{{out}}
<pre> verify ''
Line 725:
=={{header|Java}}==
<
import java.math.BigInteger;
import java.text.NumberFormat;
Line 811:
}
</syntaxhighlight>
{{out}}
Line 841:
'''Preliminaries'''
<
# otherwise, assuming infinite-precision integer arithmetic,
# if the input and $j are integers, then the result will be an integer.
Line 870:
then irt
else "isqrt requires a non-negative integer for accuracy" | error
end ;</
'''The Task'''
<
. as $n
| ($n|isqrt) as $x
Line 900:
(61, 109, 181, 277)
| solvePell as $res
| "x² - \(.)y² = 1 for x = \($res[0]) and y = \($res[1])"</
{{out}}
<pre>
Line 911:
=={{header|Julia}}==
{{trans|C#}}
<
x = BigInt(floor(sqrt(n)))
y, z, r = x, BigInt(1), x << 1
Line 933:
println("x\u00b2 - $target", "y\u00b2 = 1 for x = $x and y = $y")
end
</
<pre>
x² - 61y² = 1 for x = 1766319049 and y = 226153980
Line 943:
=={{header|Kotlin}}==
{{trans|C#}}
<
import kotlin.math.sqrt
Line 1,012:
println("x^2 - %3d * y^2 = 1 for x = %,27d and y = %,25d".format(it, x.value, y.value))
}
}</
{{out}}
<pre>x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980
Line 1,018:
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</pre>
=={{header|Lambdatalk}}==
Computing big numbers requires the lib_BN library.
<syntaxhighlight lang="Scheme">
{def pell
{lambda {:n}
{let { {:n :n}
{:x {BN.intPart {BN.sqrt :n}}} // x=int(sqrt(n))
} {pell.r :n :x :x 1 {* 2 :x} 1 0 0 1}
}}}
-> pell
{def pell.r
{lambda {:n :x :y :z :r :e1 :e2 :f1 :f2}
{let { {:n :n} {:x :x} {:z :z} {:r :r} // no closure ->
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2} // must reassign :(
{:y {BN.- {BN.* :r :z} :y}} // y=rz-y
} {let { {:n :n} {:x :x} {:y :y} {:r :r}
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}
{:z {BN.intPart
{BN./ {BN.- :n {BN.* :y :y}} :z}}} // z=(n-y*y)//z
} {let { {:n :n} {:x :x} {:y :y} {:z :z}
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}
{:r {BN.intPart {BN./ {BN.+ :x :y} :z}}} // r= (x+y)//z
} {let { {:n :n} {:x :x} {:y :y} {:z :z} {:r :r}
{:e1 :e2} // e1=e2
{:e2 {BN.+ {BN.* :r :e2} :e1}} // e2=r*e2+e1
{:f1 :f2} // f1=f2
{:f2 {BN.+ {BN.* :r :f2} :f1}} // f2=r*f2+f1
} {let { {:n :n} {:x :x} {:y :y} {:z :z} {:r :r}
{:e1 :e1} {:e2 :e2} {:f1 :f1} {:f2 :f2}
{:a {BN.+ :e2 {BN.* :x :f2}}} // a=e2+x*f2
{:b :f2} // b=f2
} {if {= {BN.compare {BN.- {BN.* :a :a}
{BN.* :n {BN.* :b :b}}}
1}
0} // a*a-n*b*b == 1
then {div}x{sup 2} - n*y{sup 2} = 1 for n=:n, x=:a, y=:b
else {pell.r :n :x :y :z :r :e1 :e2 :f1 :f2} // do it again
}}}}}}}}
-> pell.r
{S.map pell 61 109 181 277}
->
x^2 - n*y^2 = 1 for n=61, x=1766319049, y=226153980
x^2 - n*y^2 = 1 for n=109, x=158070671986249, y=15140424455100
x^2 - n*y^2 = 1 for n=181, x=2469645423824185801, y=183567298683461940
x^2 - n*y^2 = 1 for n=277, x=159150073798980475849, y=9562401173878027020
</syntaxhighlight>
=={{header|langur}}==
{{trans|D}}
<syntaxhighlight lang="langur">
val fun = fn a, b, c: [b, b * c + a]
val solvePell = fn(n) {
val x = trunc(n ^/ 2)
var y, z, r = x, 1, x * 2
for {
val
if
}
}
val
# format number string with commas
var
if
}
}
for
val
writeln
}
</syntaxhighlight>
{{out}}
Line 1,080 ⟶ 1,127:
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<
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]</
{{out}}
<pre>{{x -> 1766319049, y -> 226153980}}
Line 1,093 ⟶ 1,140:
{{trans|Python}}
{{libheader|bignum}}
<
import bignum
Line 1,116 ⟶ 1,163:
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})"</
{{out}}
Line 1,123 ⟶ 1,170:
x² - 181 * y² = 1 for (x, y) = ( 2469645423824185801, 183567298683461940)
x² - 277 * y² = 1 for (x, y) = (159150073798980475849, 9562401173878027020)</pre>
=={{header|Pascal}}==
{{libheader|IntXLib4Pascal}}
A console application in Free Pascal, created with the Lazarus IDE.
Pascal has no built-in support for arbitrarily large integers. The program below requires integers larger than 64 bits, and therefore uses an external library. The code could easily be adapted to solve the negative Pell equation x^2 - n*y^2 = -1, or show that no solution exists.
<syntaxhighlight lang="pascal">
program Pell_console;
uses SysUtils,
uIntX; // uIntX is a unit in the library IntXLib4Pascal.
// uIntX.TIntX is an arbitrarily large integer.
// For the given n: if there are non-trivial solutions of x^2 - n*y^2 = 1
// in non-negative integers (x,y), return the smallest.
// Else return the trivial solution (x,y) = (1,0).
procedure SolvePell( n : integer; out x, y : uIntX.TIntX);
var
m, a, c, d : integer;
p, q, p_next, q_next, p_prev, q_prev : uIntX.TIntX;
evenNrSteps : boolean;
begin
if (n >= 0) then m := Trunc( Sqrt( 1.0*n + 0.5)) // or use Rosetta Code Isqrt
else m := 0;
if n <= m*m then begin // if n is not a positive non-square
x := 1; y := 0; exit; // return a trivial solution
end;
c := m; d := 1;
p := 1; q := 0;
p_prev := 0; q_prev := 1;
a := m;
evenNrSteps := true;
repeat
// Get the next convergent p/q in the continued fraction for sqrt(n)
p_next := a*p + p_prev;
q_next := a*q + q_prev;
p_prev := p; p := p_next;
q_prev := q; q := q_next;
// Get the next term a in the continued fraction for sqrt(n)
Assert((n - c*c) mod d = 0); // optional sanity check
d := (n - c*c) div d;
a := (m + c) div d;
c := a*d - c;
evenNrSteps := not evenNrSteps;
until (c = m) and (d = 1);
{
If the first return to (c,d) = (m,1) occurs after an even number of steps,
then p^2 - n*q^2 = 1, and there is no solution to x^2 - n*y^2 = -1.
Else p^2 - n*q^2 = -1, and to get a solution to x^2 - n*y^2 = 1 we can
either continue until we return to (c,d) = (m,1) for the second time,
or use the short cut below.
}
if evenNrSteps then begin
x := p; y := q;
end
else begin
x := 2*p*p + 1; y := 2*p*q
end;
end;
// For the given n: show the Pell solution on the console.
procedure ShowPellSolution( n : integer);
var
x, y : uIntX.TIntX;
lineOut : string;
begin
SolvePell( n, x, y);
lineOut := SysUtils.Format( 'n = %d --> (', [n]);
lineOut := lineOut + x.ToString + ', ' + y.ToString + ')';
WriteLn( lineOut);
end;
// Main routine
begin
ShowPellSolution( 61);
ShowPellSolution( 109);
ShowPellSolution( 181);
ShowPellSolution( 277);
end.
</syntaxhighlight>
{{out}}
<pre>
n = 61 --> (1766319049, 226153980)
n = 109 --> (158070671986249, 15140424455100)
n = 181 --> (2469645423824185801, 183567298683461940)
n = 277 --> (159150073798980475849, 9562401173878027020)
</pre>
=={{header|Perl}}==
<
my ($n) = @_;
Line 1,159 ⟶ 1,292:
my ($x, $y) = solve_pell($n);
printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", $n, $x, $y);
}</
{{out}}
<pre>
Line 1,173 ⟶ 1,306:
{{libheader|Phix/mpfr}}
This now ignores the nonsquare part of the task spec, returning {1,0}.
<!--<
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
Line 1,236 ⟶ 1,369:
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"x^2 - %3d*y^2 = 1 for x = %27s and y = %25s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">xs</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ys</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</
{{out}}
<pre>
Line 1,256 ⟶ 1,389:
=={{header|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.
<syntaxhighlight lang="prolog">
% Find the square root as a continued fraction
Line 1,305 ⟶ 1,438:
X2 is X0*X1 + N*Y0*Y1,
Y2 is X0*Y1 + Y0*X1.
</syntaxhighlight>
{{Out}}
<pre>
Line 1,335 ⟶ 1,468:
=={{header|Python}}==
{{trans|D}}
<
def solvePell(n):
Line 1,356 ⟶ 1,489:
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))</
{{out}}
<pre>x^2 - 61 * y^2 = 1 for x = 1766319049 and y = 226153980
Line 1,368 ⟶ 1,501:
{{trans|Perl}}
<syntaxhighlight lang="raku"
sub pell (Int $n) {
Line 1,400 ⟶ 1,533:
my ($x, $y) = pell($n);
printf "x² - %sy² = 1 for:\n\tx = %s\n\ty = %s\n\n", $n, |($x, $y)».,
}</
{{out}}
<pre>x² - 61y² = 1 for:
Line 1,424 ⟶ 1,557:
=={{header|REXX}}==
A little extra code was added to align and commatize the gihugeic numbers for readability.
<
numeric digits 2200 /*ensure enough decimal digs for answer*/
parse arg $ /*obtain optional arguments from the CL*/
Line 1,453 ⟶ 1,586:
parse value f2 r*f2 + f1 with f1 f2
end /*until*/
return e2 + x * f2 f2</
{{out|output|text= when using the default inputs:}}
<pre>
Line 1,464 ⟶ 1,597:
=={{header|Ruby}}==
{{trans|Sidef}}
<
x = Integer.sqrt(n)
y = x
Line 1,484 ⟶ 1,617:
[61, 109, 181, 277].each {|n| puts "x*x - %3s*y*y = 1 for x = %-21s and y = %s" % [n, *solve_pell(n)]}
</syntaxhighlight>
{{Out}}
<pre>
Line 1,494 ⟶ 1,627:
</pre>
=={{header|Rust}}==
<
use num_bigint::{ToBigInt, BigInt};
use num_traits::{Zero, One};
Line 1,554 ⟶ 1,687:
println!("x^2 - {} * y^2 = 1 for x = {} and y = {}", n, r.v1, r.v2);
}
</syntaxhighlight>
{{out}}
<pre>
Line 1,562 ⟶ 1,695:
x^2 - 277 * y^2 = 1 for x = 159150073798980475849 and y = 9562401173878027020
</pre>
=={{header|Scala}}==
<syntaxhighlight lang="scala">def pellFermat(n: Int): (BigInt,BigInt) = {
import scala.math.{sqrt, floor}
val x = BigInt(floor(sqrt(n)).toInt)
var i = 0
// Use the Continued Fractions method
def converge(y:BigInt, z:BigInt, r:BigInt, e1:BigInt, e2:BigInt, f1:BigInt, f2:BigInt ) : (BigInt,BigInt) = {
val a = f2 * x + e2
val b = f2
if (a * a - n * b * b == 1) {
return (a, b)
}
val yh = r * z - y
val zh = (n - yh * yh) / z
val rh = (x + yh) / zh
converge(yh,zh,rh,e2,e1 + e2 * rh,f2,f1 + f2 * rh)
}
converge(x,BigInt("1"),x << 1,BigInt("1"),BigInt("0"),BigInt("0"),BigInt("1"))
}
val nums = List(61,109,181,277)
val solutions = nums.map{pellFermat(_)}
{
println("For Pell's Equation, x\u00b2 - ny\u00b2 = 1\n")
(nums zip solutions).foreach{ case (n, (x,y)) => println(s"n = $n, x = $x, y = $y")}
}</syntaxhighlight>
{{out}}
<pre>For Pell's Equation, x² - ny² = 1
n = 61, x = 1766319049, y = 226153980
n = 109, x = 158070671986249, y = 15140424455100
n = 181, x = 2469645423824185801, y = 183567298683461940
n = 277, x = 159150073798980475849, y = 9562401173878027020</pre>
=={{header|Sidef}}==
<
var x = n.isqrt
Line 1,594 ⟶ 1,771:
var (x, y) = solve_pell(n)
printf("x^2 - %3d*y^2 = 1 for x = %-21s and y = %s\n", n, x, y)
}</
{{out}}
<pre>
Line 1,608 ⟶ 1,785:
{{libheader|AttaSwift's BigInt}}
<
func swap(_ a: inout T, _ b: inout T, mul by: T) {
(a, b) = (b, b * by + a)
Line 1,647 ⟶ 1,824:
print("x\u{00b2} - \(n)y\u{00b2} = 1 for x = \(x) and y = \(y)")
}</
{{out}}
Line 1,658 ⟶ 1,835:
=={{header|Visual Basic .NET}}==
{{trans|Sidef}}
<
Module Module1
Line 1,682 ⟶ 1,859:
Next
End Sub
End Module</
{{out}}
<pre>x^2 - 61 * y^2 = 1 for x = 1,766,319,049 and y = 226,153,980
Line 1,693 ⟶ 1,870:
{{libheader|Wren-big}}
{{libheader|Wren-fmt}}
<
import "./fmt" for Fmt
var solvePell = Fn.new { |n|
Line 1,725 ⟶ 1,902:
var res = solvePell.call(n)
Fmt.print("x² - $3dy² = 1 for x = $-21i and y = $i", n, res[0], res[1])
}</
{{out}}
Line 1,738 ⟶ 1,915:
{{libheader|GMP}} GNU Multiple Precision Arithmetic Library
{{trans|Raku}}
<
fcn solve_pell(n){
Line 1,754 ⟶ 1,931:
if (A*A - B*B*n == 1) return(A,B);
}
}</
<
x,y:=solve_pell(n);
println("x^2 - %3d*y^2 = 1 for x = %-21d and y = %d".fmt(n,x,y));
}</
{{out}}
|