Create a free forum in seconds.
zIFBoards - Free Forum Hosting
Welcome to Dozensonline. We hope you enjoy your visit.
You're currently viewing our forum as a guest. This means you are limited to certain areas of the board and there are some features you can't use. If you join our community, you'll be able to access member-only sections, and use many member-only features such as customizing your profile, and sending personal messages. Registration is simple, fast, and completely free. (You will be asked to confirm your email address before we sign you on.)
Join our community!
If you're already a member please log in to your account to access all of our features:

Name:   Password:


 

 Solutions For Pell's Equation, Find X(n) where X^2 - nY^2 = 4.
wendy.krieger
Posted: Jun 3 2017, 12:06 PM


Dozens Demigod


Group: Members
Posts: 2,432
Member No.: 655
Joined: 11-July 12



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
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 + (-c + dx) / f = a; f/(-c+dx)
                   = a; (c+dx)f/e     Where e = d*d*x - c*c
 Since f/e can be reduced, (usually because f divides e)

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 = (c1 + d1 * x ) % 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
/**/
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(x)
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 = (gx + c1) % 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() = 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:
if f2 \= fx then return 1
return 0

sayline:
outcard = format(a2, 6) format(c2, 6) format(e2, 6) format(f2, 6) ':' b2
say outcard
return

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

isocube:
r1 = 1; do while r1 * r1 < b3; r1 = r1 + r1; end
do while r1 * r1 > b3; r1 = (r1 + b3/r1)/2; end
r1 = 1 + trunc(r1)
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 = (r1 + b3/r1/r1) /2; end
 r1 = 1 + trunc(r1)
 if b3 = r1 * (r1 * r1 - 3) then b3 = r1
 end
return

Top
« Next Oldest | General | Next Newest »
zIFBoards - Free Forum Hosting
Free Forums. Reliable service with over 8 years of experience.
Learn More · Sign-up for Free

Topic Options



Hosted for free by zIFBoards* (Terms of Use: Updated 2/10/2010) | Powered by Invision Power Board v1.3 Final © 2003 IPS, Inc.
Page creation time: 0.1084 seconds · Archive