Solutions For Pell's Equation

wendy.krieger
wendy.krieger
This is antique code. the file dates are 1998.07.27. As normal with my code, the commentry is written in a txt file, or in a weave-cmd source, the output is a rexx script, sometimes wrapping around some hamster written in UBASIC.

This is native ANSI rexx code.

Because the code is optimised for the equation in 1, this gives the corresponding isocube (since 2 has a three-place period in half of the systems), and we have to check if a cube has not been found.

This is one of the practical uses of continued fractions. We don't use decimals, just integers. Some of the columns require a BIGNUM addin, which is built into the REXX numeric digits command.

This command successfully worked to 1800 and all square-free numbers whose prime divisors are less than 30.

Home grown: no references used.

Code: Select all

``````A program to determine the isomorphic root.

It is desired to create a program, that the input is an natural number n,
and the output is the smallest z, that z*z - n*y*y = 4, for integer z, y.
The form of the command is 'ISOBASE 5' prints '3'.

This can be done by continued fractions, which we expediate by killing it
off as fast as we can.  The worked example is for base 21.

Let base = 21; and x = sqrt&#40;21&#41; = 4.58257569496

4.58257569496 = 4 + 1/1.71651513898 = 4 + &#40; -4 + x&#41;        2
1.71651513898 = 1 + 1/1.39564392376 = 1 + &#40; -1 + x&#41; / 5    2
1.39564392376 = 1 + 1/2.52752523152 = 1 + &#40; -3 + x&#41; / 4    4 = 1 *  2 +  2
2.52752523152 = 2 + 1/1.89564392421 = 2 + &#40; -3 + x&#41; / 3   10 = 2 *  4 +  2
1.89564392421 = 1 + 1/1.11651513840 = 1 + &#40; -1 + x&#41; / 4   14 = 1 * 10 +  4
1.11651513840 = 1 + 1/8.58257573850 = 1 + &#40; -4 + x&#41; / 5   24 = 1 * 14 + 10
8.58257573850 = 4 + x               = 8 + &#40; -4 + x&#41;      110 = 4 * 24 + 14

Where 110 = the cube of 5, also, 110 * 110 - 24 * 24 * 1 = 4.

On paper, one sets this out thus

a    b    c    d    e    f    g
0         2
1         0    0    1         1
2    4    2    4    1    5    5
3    1    2    1    1   20    4
4    1    4    3    1   12    3
5    2   10    3    1   12    4
6    1   14    1    1   20    5
7    1   24    4    1    5    1   *
8    4  110

For a prime like 113

a   b    c    d    e    f    g
0       2
1       0    0    1         1
2  10   2   10    1   13   13
3   1   2    3    1  104    8
4   1   4    5    1   88   11
5   1   6    6    1   77    7
6   2  16    8    1   49    7
7   2  38    6    1   77   11
8   1  54    5    1   88    8
9   1  92    3    1  104   13
10   1 146   10    1   13    1   *
11  10 1552

1552 * 1552 < 146 * 146 * 113
4817412

For the composite number 94

x =  9.695

a   c   d   e    f         b
2
0   1        1         0
9   9   1   13  13         2
1   4   1   78   6         2
2   8   1   30   5         6
3   7   1   45   9        20
1   2   1   90  10        26
1   8   1   30   3        46
5   7   1   45  15       256
1   8   1   30   2       302
8   8   1   30  15      2672
1   7   1   45   3      2974
5   8   1   30  10     17542
1   2   1   90   9     20516
1   7   1   45   5     38058
3   8   1   30   6    134690
2   4   1   78  13    307438
1   9   1   13   1    442128 *
9                    4286590

This method works on evaluating the continued fraction of x, which forms in
the values a. The values in c, d, and f represent the fractional part of the
next term, ie

a + &#40;-c + dx&#41; / f = a; f/&#40;-c+dx&#41;
= a; &#40;c+dx&#41;f/e     Where e = d*d*x - c*c
Since f/e can be reduced, &#40;usually because f divides e&#41;

Some form of code for this.

We operate on a window of current values.

Let us work on row 2, for example, keeping just those variables required
to move from row to row.  Here, a1 means the value of a in row 1.

a0  b0  c0  d0  e0  f0
a1  b1  c1  d1  e1  f1

a2  b2  c2  d2  e2  f2

b0 = 2; b1 = 0; c1 = 0; d1 = 1; f1 = 1

a2 = &#40;c1 + d1 * x &#41; % f1
b2 = a2 * b1 + b0
c2 = a2 * f2 - c1
d2 = d2
e2 = d2 * d2 * bb - c2 * c2
f2 = e2 / f1
if f2 \= 1 then iterate
return b3 = b2 * c2 + b1
``````

Code: Select all

``````/**/
numeric digits 40
numeric fuzz 2
parse arg bb
/* we wish to find the square root of bb */

if bb < 1 then do
say 'The base must be greater than unity'
exit
end
y = 2; yy = y * y;
do while yy < bb;
if bb // yy = 0
then bb = bb / yy
else do; y = y + 1; yy = y * y; end
end

x = 1; do while x * x < bb; x = x + x; end
do while x * x > bb; x = x + bb / x; x = x / 2; end
gx = trunc&#40;x&#41;
if bb = gx * gx then do
say 2
exit
end

/* We now initialise the value to seek */

b0 = 2; b1 = 0; fx = 1
c1 = 0; f1 = 1      /* the fraction x1 + y1 * i / z1 */

/* Main loop */
do forever
a2 = &#40;gx + c1&#41; % f1
b2 = a2 * b1 + b0
c2 = a2 * f1 - c1
e2 = bb - c2 * c2
f2 = e2 % f1; call sayline
if e2 // f1 \= 0 then call errors 4
if isit&#40;&#41; = 0 then leave
b0 = b1; b1 = b2
c1 = c2; f1 = f2
end
b3 = b2 * c2 + b1

if b3 * b3 > b2 * b2 * bb
then b3 = b3
else b3 = b3 * b3 + 2
call isocube
say b3
exit

isit&#58;
if f2 \= fx then return 1
return 0

sayline&#58;
outcard = format&#40;a2, 6&#41; format&#40;c2, 6&#41; format&#40;e2, 6&#41; format&#40;f2, 6&#41; '&#58;' b2
say outcard
return

errors&#58;
parse arg num
select
when num = 4 then msg = 'Value of d1 excedes 1'
otherwise msg = 'Unidentified error'
end
say msg
return

isocube&#58;
r1 = 1; do while r1 * r1 < b3; r1 = r1 + r1; end
do while r1 * r1 > b3; r1 = &#40;r1 + b3/r1&#41;/2; end
r1 = 1 + trunc&#40;r1&#41;
if b3 = r1 * r1 - 2; then b3 = r1
if bb // 4 = 1 then do
r1 = 1; do while r1 * r1 * r1 < b3; r1 = r1 + r1; end
do while r1 * r1 * r1 > b3; r1 = &#40;r1 + b3/r1/r1&#41; /2; end
r1 = 1 + trunc&#40;r1&#41;
if b3 = r1 * &#40;r1 * r1 - 3&#41; then b3 = r1
end
return

``````