How to resolve the algorithm Pell's equation step by step in the Pascal programming language
Published on 12 May 2024 09:40 PM
How to resolve the algorithm Pell's equation step by step in the Pascal programming language
Table of Contents
Problem Statement
Pell's equation (also called the Pell–Fermat equation) is a Diophantine equation of the form: with integer solutions for x and y, where n is a given non-square positive integer.
Let's start with the solution:
Step by Step solution about How to resolve the algorithm Pell's equation step by step in the Pascal programming language
Source code in the pascal programming language
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.
You may also check:How to resolve the algorithm Sort an integer array step by step in the Eiffel programming language
You may also check:How to resolve the algorithm String matching step by step in the GDScript programming language
You may also check:How to resolve the algorithm Boolean values step by step in the Elixir programming language
You may also check:How to resolve the algorithm Population count step by step in the Lua programming language
You may also check:How to resolve the algorithm Strip comments from a string step by step in the ALGOL W programming language