(49/50) Pellian Equations
12-26-2023, 08:38 PM (This post was last modified: 12-26-2023 08:43 PM by John Keith.)
Post: #1
 John Keith Senior Member Posts: 1,034 Joined: Dec 2013
(49/50) Pellian Equations
The first program listed here finds the basic solution to the Pellian equations
x^2 - d*y^2 = 1 or x^2 - d*y^2 = -1. To solve for the negative equation, enter
d as a negative number. For d to be a valid argument for the negative equation,
the repeating part of the continued fraction for the square root of d must be
of odd length. See A031396 and A003814. If d is a square or if d is negative but not
a valid argument, the trivial solution {1 0} is returned.

This program is inspired by a program by Bob Pulluard in Datafile V19N1,
available on the HHC thumb drive and from Jake Schwartz's PPC Archive.
This program uses a completely different algorithm based on convergents to the square root of d.
It is larger but at least 10 times as fast. The program requires the LongFloat Library because
the computations can require very high precision, up to thousands of digits.

The maximum size of d is limited by the amount of available memory, or the
patience of the user. The largest value that I have solved for is 1026061,
which requires over 7 hours on the physical 50g. The x value returned is
a 2555 digit number. See A033316 for values of d whose solutions are increasingly large.

Some notes on the theory of operation:

The basis of the program is the observation that solutions of the Pellian
equation for d are a subset of the convergents to the square root of d.
After determining that d is neither a square (invalid argument) nor 1 plus a
square (x = floor(sqrt(d)), y = 1), LongFloat functions are used to generate
a list of the repeating part of the continued fraction for the square root d.
If the length of the repeating part is odd, there is a solution for the
negative equation. This solution can be found in about half the time of the
solution for the positive equation, so this is done, and the positive
solution is computed from the negative solution if required.

The main loop of the program computes the convergents from the terms
of the continued fraction, and tests them until a solution is found. This is similar
to my recent Convergents program here. If the negative solution was computed and
the positive solution wanted, the former is converted to the latter.

Code:
 @ Given integer d on the stack, returns the basic solution to the Pellian @ equation x^2 - d*y^2 = 1 if d is positive, or x^2 - d*y^2 = -1 if d is @ negative and if the length of the repeating part of the continued fraction @ for sqrt(d) is odd. Requires LongFloat. @ Local variables @ d: absolute value of d @ b: sign of d @ c: repeating part of continued fraction for sqrt(d) @ n: length of c @ a: -1 if d has solution for -1 else 1 @ s: maximum length of loop, also used to exit loop @ k: looping variable \<< DUP ABS DUP ROT SIGN \-> d b           @ Save absolute value and sign   \<<                                      @ Negative means solve for -1     IF DUP ZSqrt                           @ Is d a square?     THEN DROP2 { 1 0 }                     @ Then return trivial solution     ELSE DROP       IF DUP 1 - ZSqrt                     @ Is d 1 plus a square?       THEN NIP 1 b 1 SAME                  @ Return basic solution for -1       ELSE 2 * DUP                         @ 2 * floor(sqrt(d))         IF PICK3 I\->R 8940. \<=           @ If d <= 8940         THEN 20 MAX                        @ At least 20 digit precision         ELSE 3 * 2 /                       @ 1.5 * precision for large d         END 'DIGITS' STO                   @ Store precision for LongFloat         SWAP FSQRT Fl\->SF SF\->CF         @ Continued fraction for sqrt(d)         DUP ROT POS 1. SWAP SUB            @ List up to end of repeating part         DUP HEAD SWAP TAIL                 @ Separate head         DUP SIZE DUP 2. MOD                @ Parity of length of repeating part         IF DUP NOT b -1 SAME AND           @ Even length but solving for -1         THEN DROP2 DROP2 1 0 0.            @ Invalid d for -1         ELSE -1 1 IFTE 4. ALOG \-> c n a s @ Save local variables           \<< c 1. GET                     @ Get second CF term             1 SWAP PICK3 OVER * 1 + SWAP   @ First two convergents             1. s                           @ Loop to generate and test conv's             FOR k c k n MOD 1. + GET       @ Get next term modulo length of c               PICK3 OVER * 6. ROLL +       @ Next numerator               PICK3 ROT * 5. ROLL +        @ Next denominator               OVER SQ OVER SQ d * - a SAME @ Is this convergent a solution?               s 1. IFTE                    @ Exit loop if true             STEP ROT DROP ROT DROP         @ Drop previous convergent             b 1 SAME a -1 SAME AND         @ Flag for converting -1 solution           \>>         END       END { OVER SQ 2 * 1 + UNROT * 2 *    @ Convert -1 solution to +1           } IFT 2. \->LIST                 @ Save as list     END   \>> \>>

The second program returns a list of solutions to the Pellian equation.
It assumes the existence of a variable called PELIAN containing the program above.
Given integers d on level 2 and n on level 1, returns list of the first n
solutions for x^2 - d*y^2 = 1 if d is positive, or x^2 - d*y^2 = -1 if
d is negative. Once the basic solution is found by PELIAN, the subsequent
solutions can be computed by a simple recurrence. Therefor the program is
very fast, although large values of d and/or n will quickly overflow memory.
If d is a square or if d is negative but not a valid argument, a list of n
copies of the trivial solution {1 0} is returned.

Code:
 @ Given d on level 2 and n on level 1, returns list of the first n solutions @ for x^2 - d*y^2 = 1 if d is positive, or x^2 - d*y^2 = -1 if @ d is negative. d must be in A031396 to return a valid result for -1. \<< OVER PELIAN DUP { 1 0 }   IF SAME                          @ If not a valid argument   THEN ROT DROP SWAP NDUPN \->LIST @ Return list of {1 0}   ELSE ROT SIGN UNROT              @ Save real n and     SWAP I\->R 0 \-> n m           @ dummy value for m     \<< DUP 1. GET                 @ x value to compute coefficient       IF ROT 1 SAME                @ Solving for +1?       THEN 2 * 'm' STO             @ Then 2*x         { 1 0 } SWAP               @ Starting list for recurrence       ELSE SQ 4 * 2 + 'm' STO      @ 4*x^2 + 2         DUP { -1 1 } * SWAP        @ Negate x for starting list       END 2. n                     @ Start main loop       START DUP m * PICK3 -       NEXT n \->LIST NIP           @ Save list and remove starting list.     \>>   END \>>
12-27-2023, 06:30 AM
Post: #2
 Gerald H Senior Member Posts: 1,627 Joined: May 2014
RE: (49/50) Pellian Equations
Here a SYS version for 49G/50g, only deals with the +1 case.

49G (1.19-6) takes

124 s

to process 1026061 correctly.

Code:
pell  Size: 257. CkSum: # A5BAh ::   CK1&Dispatch   # FF   ::     DUP     ID x003     DUP     FPTR2 ^ZABS     ZINT 4     EQUAL     IT     ::       ZINT 4       FPTR2 ^ZQUOText       3PICK       FPTR2 ^ZSQ_       2DUPSWAP       ZINT 3       FPTR2 ^QMul       FPTR2 ^QSub       ZINT 2       FPTR2 ^ZQUOText       5ROLL       FPTR2 ^QMul       SWAP3PICK       FPTR2 ^QSub       ZINT 2       FPTR2 ^ZQUOText       4ROLL       FPTR2 ^QMul       ROT     ;     FPTR2 ^ZIsNeg?     IT     ::       OVER       FPTR2 ^QMul       DUP       FPTR2 ^QAdd       SWAP       FPTR2 ^ZSQ_       DUP       FPTR2 ^QAdd       ZINT 1       FPTR2 ^QAdd       SWAP     ;     2DUP     BEGIN     DUP     7PICK     FPTR2 ^ZMod     ZINT 0     EQUALNOT     WHILE     ID x004     REPEAT     6ROLL     FPTR2 ^ZQUOText     TWO{}N     4UNROLL3DROP   ; ; x003  Size: 456. CkSum: # 4883h ::   FPTR2 ^MZSQFF   #2/   ZINT 1   DUPROT   ZERO_DO   ROT   COERCE   BINT2   #/   5PICK   SWAP   FPTR2 ^PPow#   4ROLL   FPTR2 ^QMul   3UNROLL   4ROLLSWAP   FPTR2 ^PPow#   FPTR2 ^QMul   LOOP   DUP   ZINT 5   EQUALcase   ::     ZINT 1     ZINT 1     ZINT -4   ;   DUP   TRUE   SWAPDUP   ::     DUP     ZINT 2     Z<     ?SEMI     ZINT 2     OVER     FPTR2 ^Z>ZH     FPTR2 ^ZBits     SWAPDROP     #2/     #1+     FPTR2 ^PPow#     BEGIN     DUP     3PICKOVER     FPTR2 ^ZQUOText     FPTR2 ^QAdd     ZINT 2     FPTR2 ^ZQUOText     2DUP     Z>     WHILE     SWAPDROP     REPEAT     DROPSWAPDROP   ;   DUPUNROT   FPTR2 ^ZSQ_   FPTR2 ^QSub   ZINT 1   3PICK   ZINT 1   OVER   ZINT 0   ZINT 1   '   NULLLAM   BINT8   NDUPN   DOBIND   BEGIN   7GETLAM   FPTR2 ^DupZIsOne?   NOT   SWAP   ZINT 4   EQUALNOT   AND   WHILE   ::     5GETLAM     8GETLAM     FPTR2 ^QAdd     7GETLAM     FPTR2 ^ZQUOText     DUP     7GETLAM     FPTR2 ^QMul     5GETLAM     FPTR2 ^QSub     5GETLAM     SWAPDUP     5PUTLAM     FPTR2 ^QSub     OVER     FPTR2 ^QMul     6GETLAM     FPTR2 ^QAdd     7GETLAM     6PUTLAM     7PUTLAM     3GETLAM     2DUPSWAP     FPTR2 ^QMul     4GETLAM     FPTR2 ^QAdd     3PUTLAM     4PUTLAM     1GETLAM     DUPUNROT     FPTR2 ^QMul     2GETLAM     FPTR2 ^QAdd     1PUTLAM     2PUTLAM     NOT   ;   REPEAT   3GETLAM   1GETLAM   7GETLAM   ABND   4ROLL   NOT?SEMI   FPTR2 ^RNEGext ; x004  Size: 69.5 CkSum: # CC0Dh ::   OVER   5PICK   FPTR2 ^QMul   OVER   5PICK   FPTR2 ^QMul   7PICK   FPTR2 ^QMul   FPTR2 ^QAdd   SWAP   5PICK   FPTR2 ^QMul   ROT   4PICK   FPTR2 ^QMul   FPTR2 ^QAdd ;
12-27-2023, 12:56 PM (This post was last modified: 12-27-2023 09:24 PM by Albert Chan.)
Post: #3
 Albert Chan Senior Member Posts: 2,686 Joined: Jul 2018
RE: (49/50) Pellian Equations
(12-26-2023 08:38 PM)John Keith Wrote:  If the negative solution was computed and
the positive solution wanted, the former is converted to the latter.

(x² − d*y²) = -1
(x² − d*y²)² = 1

(x² + d*y²)² − d*(2xy)² = 1

Some number theory book prefer solution (x,y) in the form z = x + y*√d

z * conj(z) = -1
z² * conj(z²) = 1

z² = (x + y*√d)² = (x² + d*y²) + (2xy)*√d
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)