(49/50) Pellian Equations

[John Keith]

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
\>>

[Gerald H]

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
;

[Albert Chan]

(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