Solutions For Pell's Equation

Solutions For Pell's Equation

wendy.krieger
wendy.krieger

Jun 3 2017, 12:06 PM #1

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(21) = 4.58257569496

  4.58257569496 = 4 + 1/1.71651513898 = 4 + ( -4 + x)        2
  1.71651513898 = 1 + 1/1.39564392376 = 1 + ( -1 + x) / 5    2
  1.39564392376 = 1 + 1/2.52752523152 = 1 + ( -3 + x) / 4    4 = 1 *  2 +  2
  2.52752523152 = 2 + 1/1.89564392421 = 2 + ( -3 + x) / 3   10 = 2 *  4 +  2
  1.89564392421 = 1 + 1/1.11651513840 = 1 + ( -1 + x) / 4   14 = 1 * 10 +  4
  1.11651513840 = 1 + 1/8.58257573850 = 1 + ( -4 + x) / 5   24 = 1 * 14 + 10
  8.58257573850 = 4 + x               = 8 + ( -4 + x)      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

Quote
Share