hp41programs

Diophantine

A few Diophantine Equations for the HP-41


Overview
 

 1°)  Bezout Identity:  a.u + b.v + c.w = d
 2°)  Decomposition in Sums of Powers
 3°)  X^2 + Y^3 + Z^4 = N
 4°)  Markov Numbers:  X^2 + Y^2 + Z^2 = 3 X.Y.Z
 5°)  Pythagoras Equation

  5-a)  X^2+Y^2=Z^2
  5-b)  X^2+Y^2+Z^2=T^2
 
 

1°)  Bezout Identity
 

-A program is listed in "GCD , LCM & Bezout Identity" to find an integer solution of the Diophantine equation  a u + b v = c
-The folllowing routine tries to solve  a u + b v + c w = d  where a , b , c , d  are given integers and  u , v, w  are unknown integers # 0.

-There will be solutions iff  d is a multiple of gcd ( a , b , c )
 

Data Registers:   R00 = d  , R01 = a , R02 = b , R03 = c , R04 = +/-u , R05 = +/-v ,  R06-R07-R08: temp
Flags: /
Subroutines: /
 
 

  01  LBL "UVW"
  02  STO 00
  03  RDN
  04  STO 03
  05  RDN
  06  STO 02       
  07  X<>Y
  08  STO 01
  09  1
  10  STO 07
  11  LBL 00
  12  RCL 07
  13  STO 04
  14  SIGN
  15  STO 05
  16  ST+ 07
  17  LBL 01
  18  RCL 00
  19  RCL 01
  20  RCL 04
  21  *
  22  STO 06
  23  RCL 02
  24  RCL 05       
  25  *
  26  ST- 06
  27  +
  28  ST- Z
  29  +
  30  RCL 03
  31  ST/ Z
  32  /
  33  FRC
  34  X=0?
  35  XEQ 02
  36  RDN
  37  FRC
  38  X=0?
  39  XEQ 03
  40  X<> 06
  41  ST- Z
  42  +
  43  RCL 03
  44  ST/ Z
  45  /
  46  FRC
  47  X=0?
  48  XEQ 04
  49  X<>Y
  50  FRC
  51  X=0?
  52  XEQ 05
  53  ISG 05
  54  CLX
  55  DSE 04
  56  GTO 01
  57  GTO 00
  58  LBL 02
  59  X<> L
  60  X=0?
  61  RTN
  62  X<>Y
  63  STO 08
  64  CLX
  65  RCL 05       
  66  CHS
  67  RCL 04
  68  CHS
  69  STOP
  70  RCL 00
  71  RCL 00
  72  RCL 08
  73  ENTER^
  74  RTN
  75  LBL 03
  76  X<> L
  77  X=0?
  78  RTN
  79  RCL 05       
  80  RCL 04
  81  STOP
  82  RCL 00
  83  ENTER^
  84  ENTER^
  85  RTN
  86  LBL 04
  87  X<> L
  88  X=0?
  89  RTN
  90  X<>Y
  91  STO 08
  92  CLX
  93  RCL 05
  94  RCL 04       
  95  CHS
  96  STOP
  97  RCL 08
  98  X<>Y
  99  RTN
100  LBL 05
101  LASTX
102  X=0?
103  RTN
104  RCL 05
105  CHS
106  RCL 04
107  STOP
108  END

 
    ( 138 bytes / SIZE 009 )
 
 

      STACK        INPUTS     OUTPUTS1   OUTPUTS2...
           T             a             /            /
           Z             b             w           w'
           Y             c             v            v'
           X             d             u            u'

 
Example:

    41   ENTER^
   178  ENTER^
    12   ENTER^
     7    XEQ "UVW"   >>>>    1  RDN  -1  RDN  12    thus  41 x 1 + 178 x (-1) + 12 x 12 = 7

               R/S             >>>>     -3  RDN  1   RDN  -4    so,   41 x (-3) + 178 x 1 + 12 x (-4) = 7

               R/S             >>>>      3   RDN   -2   RDN   20   i-e   41 x 3 + 178 x (-2) + 12 x 20 = 7   .... and so on ....

Notes:

-If there is one solution, the number of solutions is infinite.
-So the program could run ( almost ) forever.
-If a , b , c  are relatively large numbers, the execution time may be large too...

-We cannot replace the STOPs by RTNs because we are inside a local subroutine there.
 

2°)  Decomposition in Sums of Powers
 

-Given 3 positive integers N , n , k ,  "DNK"  searches n integers  1 <= x1 <= x2 <= .......... <= xn  such that

    N = x1k + ............. + xnk

-Each time a solution is found, it is displayed and the program stops ( line 72 )
-Press R/S to seek another result.
-When all the solutions have been displayed - if any - the program stops with X = 5.
 

Data Registers:    R00 to R04: temp

   >>>  When the routine displays a solution,  R05 = x1 ,  R06 = x2 , .....  , R3+n = xn-1  and   xn is in register Y

Flag:   F07
Subroutines: /
 

-The "append" character is denoted  ~
 
 

  01  LBL "DNK"
  02  STO 02
  03  RDN
  04  STO 00
  05  X<>Y
  06  STO 01
  07  CLX
  08  3
  09  +
  10  .1
  11  %
  12  STO 04         
  13  5
  14  +
  15  SIGN
  16  LBL 00
  17  STO IND L
  18  ISG L
  19  GTO 00
  20  CLX
  21  .004
  22  +
  23  STO 03
  24  FIX 0
  25  CF 29
  26  SF 07
  27  LBL 01
  28  RCL 03
  29  RCL 01
  30  LBL 02
  31  RCL IND Y
  32  RCL 02
  33  Y^X
  34  -
  35  DSE Y
  36  GTO 02
  37  X<0?
  38  GTO 05
  39  STO Y
  40  RCL 02         
  41  1/X
  42  Y^X
  43  RND
  44  STO Z
  45  RCL IND 03
  46  X>Y?
  47  GTO 05
  48  SF 07
  49  CLX
  50  RCL 02
  51  Y^X
  52  X#Y?
  53  GTO 04
  54  R^
  55  CLA
  56  ARCL 01
  57  "~="
  58  RCL 04         
  59  FRC
  60  5
  61  +
  62  LBL 03
  63  ARCL IND X
  64  "~^"
  65  ARCL 02
  66  "~+"
  67  ISG X
  68  GTO 03
  69  ARCL Y
  70  "~^"
  71  ARCL 02
  72  PROMPT
  73  LBL 04
  74  SIGN
  75  ST+ IND 03
  76  GTO 01
  77  LBL 05
  78  RCL 03
  79  INT
  80  DSE X
  81  RCL 04         
  82  FRC
  83  +
  84  FS? 07
  85  STO 04
  86  RCL IND 03
  87  RCL IND 04
  88  X=Y?
  89  GTO 07
  90  SIGN
  91  ST+ IND 04
  92  RCL 04
  93  RCL IND 04
  94  LBL 06
  95  STO IND Y
  96  ISG Y
  97  GTO 06
  98  CF 07
  99  GTO 01
100  LBL 07
10 1 SIGN
102  ST- 04
103  RCL 04         
104  5
105  X<=Y?
106  GTO 01
107  FIX 4
108  SF 29
109  END

 
   ( 172 bytes / SIZE ??? )
 
 

      STACK        INPUTS      OUTPUTS
           Z            N             /
           Y             n            xn
           X             k             /

 
Example1:

    17863  ENTER^
        3      ENTER^
        5      XEQ "DNK"  >>>>  "17863=2^5+4^5+7^5"

                     R/S           >>>>    5.0000      i-e no more solution

Example2:

    100   ENTER^
      4     ENTER^
      2     XEQ 'DNK"   >>>>   "100=1^2+1^2+7^2+7^2"

                    R/S         >>>>    "100=1^2+3^2+3^2+9^2"
                    R/S         >>>>    "100=1^2+5^2+5^2+7^2"
                    R/S         >>>>    "100=2^2+4^2+4^2+8^2"
                    R/S         >>>>    "100=5^2+5^2+5^2+5^2"

                    R/S         >>>>     5.0000      i-e no more solution

Notes:

-If the alpha register cannot contain all the xi , remember that R05 = x1 ,  R06 = x2 , .....  , R3+n = xn-1 and   xn is in register Y
-If you have an M-Code function XROOT, replace lines 39 to 54 by

   RCL 02                 SF 07
   XROOT                RDN
   STO Y                  FRC
   RCL IND 03         X#0?
   X>Y?                    GTO 04
   GTO 05                X<>Y
 

3°)  X^2 + Y^3 + Z^4 = N
 

-This routine just illustrates one variant of the program above
 

Data Registers:  R00 = N , R01-R02: temp
Flags:  F29
Subroutines: /

-The "append" character is denoted  ~
 
 

 01  LBL "D234"
 02  STO 00        
 03  CLX
 04  STO 01
 05  LBL 00
 06  CLX
 07  STO 02
 08  SIGN
 09  ST+ 01
 10  RCL 01
 11  X^2
 12  X^2
 13  STO 03
 14  RCL 00        
 15  X<=Y?
 16  GTO 02
 17  LBL 01
 18  SIGN
 19  ST+ 02
 20  RCL 00
 21  RCL 03
 22  -
 23  RCL 02        
 24  3
 25  Y^X
 26  -
 27  X<=0?
 28  GTO 00
 29  SQRT
 30  FRC
 31  X#0?
 32  GTO 01        
 33  FIX 0
 34  CF 29
 35  CLA
 36  ARCL 00
 37  "~="
 38  ARCL L
 39  "~^2+"
 40  ARCL 02
 41  "~^3+"
 42  ARCL 01
 43  "~^4"
 44  FIX 4
 45  SF 29
 46  PROMPT
 47  GTO 01
 48  LBL 02        
 49  END

 
    ( 85 bytes / SIZE 004 )
 
 

      STACK         INPUT       OUTPUT
           X             N             /
 
 
Example1:       N = 2000

  2000  XEQ "D234"  >>>>  "2000=16^2+12^3+2^4"
                   R/S         >>>>   "2000=4^2+12^3+4^4"
                   R/S         >>>>   "2000=19^2+7^3+6^4"

                   R/S         >>>>    2000.0000                          no other solutions

Example2:      N = 16807 = 7^5

 16807        R/S         >>>>    "16807=68^2+23^3+2^4"
                   R/S         >>>>    "16807=74^2+11^3+10^4"

                   R/S         >>>>    16807.0000                          no other solutions
 

4°)  Markov Numbers
 

-Markov numbers are integer solutions of the equation:  x2 + y2 + z2 = 3 xyz   (E)
-All the solutions may be obtained from ( 1 , 1 , 1 ) by the following rule:

-If ( x , y , z ) is a solution then  ( y , z , x' ) and ( x , z , y' ) are also solutions of (E)

  where  x' = 3 y z - x   and   y' = 3 x z - y

-"MARKOV"  takes an alpha string of "A" and "B" characters and returns a tripple ( x , y , z )  in the stack.
 

Data Registers:   R00: unused  ,  R01 = x , R02 = y , R03 = z
Flags: /
Subroutines: /
 
 

 01  LBL "MARKOV"
 02  1
 03  STO 01              
 04  STO 02
 05  STO 03
 06  LBL 01
 07  ATOX
 08  X=0?
 09  GTO 03
 10  65
 11  X=Y?
 12  GTO 02
 13  RCL 01              
 14  RCL 03
 15  *
 16  3
 17  *
 18  RCL 02
 19  -
 20  X<> 03
 21  STO 02              
 22  GTO 01
 23  LBL 02
 24  RCL 02
 25  RCL 03
 26  *
 27  3
 28  *
 29  RCL 01              
 30  -
 31  X<> 03
 32  X<> 02
 33  STO 01
 34  GTO 01
 35  LBL 03
 36  RCL 03              
 37  RCL 02
 38  RCL 01
 39  END

 
   ( 59 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             /             z
           Y             /             y
           X             /             x

 
Examples:

   "ABBAAAB"   XEQ "MARKOV"   >>>>   x = 194 , y = 4400489 , z = 2561077037
     "AAAAA"                 R/S               >>>>   x =  29  , y = 433 , z = 37666
     "BBBBB"                  R/S               >>>>   x =   1   , y =  34  , z =  89
 

Notes:

-You can continue with the results obtained by placing another string of "A" or "B" in the alpha register and  XEQ 01
-Of course, the numbers will be only approximate if they exceed  10^10

-Instead of an alpha string, the following routine takes a positive integer in X and returns 3 Markov numbers..
 

Data Registers:   R00: unused  ,  R01 = x , R02 = y , R03 = z   R04: temp
Flags: /
Subroutines: /
 
 

 01  LBL "MARKV2"
 02  STO 04              
 03  1
 04  STO 01
 05  STO 02
 06  STO 03
 07  LBL 01
 08  RCL 04
 09  X=0?
 10  GTO 03              
 11  2
 12  MOD
 13  ST- 04
 14  LASTX
 15  ST/ 04
 16  X<>Y
 17  X=0?
 18  GTO 02
 19  RCL 02
 20  RCL 01              
 21  GTO 01
 22  LBL 02
 23  RCL 01
 24  RCL 02
 25  STO 01
 26  LBL 01
 27  RCL 03
 28  STO 02              
 29  *
 30  3
 31  *
 32  X<>Y
 33  -
 34  STO 03
 35  GTO 01
 36  LBL 03
 37  RCL 03              
 38  RCL 02
 39  RCL 01
 40  END

 
   ( 57 bytes / SIZE 005 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             /             z
           Y             /             y
           X            N             x

 
Example:

  157   XEQ "MARKV2"  >>>>    x = 89 , y = 2423525 , z = 647072098
 

5°)  Pythagoras Equation
 

    a)  X^2+Y^2=Z^2
 

-All the solutions of this equations may be computed from the solution ( 3 , 4 , 5 ) and 3 matrices A , B , C where

             1   -2    2                       -1    2    2                         1    2    2
   A =    2   -1    2               B =  -2    1    2                 C =  2    1    2
             2    2    3                       -2    2    3                         2    2    3

-You place an alpha string containing the letters  A  B  C  and the HP-41 returns in X , Y , Z  a solution of the Pythagoras equation
 

Data Registers: /
Flags: /
Subroutines: /
 
 

 01  LBL "PYTH"
 02  5
 03  4
 04  3
 05  LBL 10        
 06  ATOX
 07  X=0?
 08  GTO 12
 09  RDN
 10  XEQ IND T
 11  ENTER^
 12  ST+ X
 13  R^
 14  ST+ Y         
 15  R^
 16  ST- T
 17  +
 18  ST+ Y
 19  X<>Y
 20  ST+ Y
 21  ENTER^
 22  R^
 23  -
 24  GTO 10        
 25  LBL 65
 26  X<>Y
 27  CHS
 28  X<>Y
 29  RTN
 30  LBL 66        
 31  CHS
 32  RTN
 33  LBL 67        
 34  RTN
 35  LBL 12
 36  RDN
 37  END
 

 
    ( 59 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             /             z
           Y             /             y
           X             /             x

 
Example:

   ABBCAABCABCCAAABBB   XEQ "PYTH"   >>>>    x = 7866623169  ,  y = 2205672440  ,  z = 8169990881
 

Notes:

-You can continue with the results obtained by placing another string of  A  B  C   in the alpha register and  XEQ 10
-The results will be only approximate if they exceed  10^10

-Instead of an alpha string, the following variant takes 2 positive integers p &q  in X & Y and returns 3 Pythagoras numbers..
 

Data Registers: /
Flags: /
Subroutines: /
 
 

 01  LBL "PYTH2"
 02  ENTER^
 03  X^2
 04  ENTER^
 05  R^
 06  ST* T
 07  X^2
 08  ST+ Z
 09  -
 10  ABS
 11  R^
 12  ST+ X
 13  END

 
( 26 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             /             z
           Y             p             y
           X             q             x

 
Example:

   41247  ENTER^
   87268  XEQ "PYTH2"  >>>>   x = 7199086392 , y = 5914388815 , z = 9317018833

 You can check that  x^2 + y^2 = z^2

Notes:

-We have  x = 2 p.q  ,  y = | p2 - q2 |  ,  z =  p2 + q2
-It can be proved that this method gives all the solutions !
 

    b)  X^2+Y^2+Z^2=T^2
 

-Here we place 2 positive integers  m , n  in registers X , Y  and the HP-41 returns 1 or several solutions of  x^2 + y^2 + z^2 = t^2
-We have:

    x = ( m2 + n2 - p2 ) / p   ,  y = 2 m  ,  z = 2 n ,  t = ( m2 + n2 + p2 ) / p

-The program tests all the integers p from 1 to sqrt( m2 + n2 )  and only keeps those p that are divisors of  ( m2 + n2 )
 

Data Registers:  R00 = p  ,  R01 = m , R02 = n , R03 = m2 + n2
Flags: /
Subroutines: /
 
 

 01  LBL "PYTH3"
 02  STO 01         
 03  X^2
 04  X<>Y
 05  STO 02
 06  X^2
 07  +
 08  STO 03
 09  CLX
 10  STO 00
 11  LBL 01
 12  RCL 03         
 13  ISG 00
 14  CLX
 15  RCL 00
 16  X^2
 17  X>Y?
 18  GTO 00
 19  X<> L
 20  MOD
 21  X#0?
 22  GTO 01
 23  RCL 03         
 24  RCL 03
 25  RCL 00
 26  X^2
 27  ST+ Z
 28  -
 29  RCL 00
 30  ST/ Z
 31  /
 32  RCL 02         
 33  ST+ X
 34  RCL 01
 35  ST+ X
 36  X<> Z
 37  RTN
 38  GTO 01
 39  LBL 00         
 40  "END"
 41  AVIEW
 42  END
 

 
    ( 65 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /             t
           Z             /             z
           Y             n             y
           X             m             x

 
Example:

   5   ENTER^
   7   XEQ "PYTH3"  >>>>   x = 73 , y = 10 , z = 14 , t = 75
              R/S              >>>>   x = 35 , y = 10 , z = 14 , t = 39
              R/S              >>>>   "END"

-Of course, it may take a long time if m , n are large.
-Use a good emulator and activate the turbo mode...

-The following variant returns pythagorean quadruples in a faster way:  place 4 positive integers m , n , p , q  in the stack and  XEQ "PYTH3"

  it yields:     x =  | m2 + n2 - p2 - q2 |  ,   y = 2 ( m.q + n.p )  ,  z = 2 | n.q - m.p |  ,  t = m2 + n2 + p2 + q2
 

Data Registers:  R00 thru R03: temp
Flags: /
Subroutines: /
 
 

 01  LBL "PYTH3"
 02  STO 00         
 03  STO 02
 04  X^2
 05  X<>Y
 06  STO 01
 07  STO 03
 08  X^2
 09  +
 10  X<>Y
 11  ST* 00         
 12  ST* 01
 13  X^2
 14  R^
 15  ST* 02
 16  ST* 03
 17  X^2
 18  +
 19  X<>Y
 20  RCL 01         
 21  RCL 02
 22  +
 23  ST+ X
 24  RDN
 25  ST- Z
 26  +
 27  RCL 03         
 28  ST- 00
 29  X<> 00
 30  ABS
 31  ST+ X         
 32  R^
 33  R^
 34  ABS
 35  END

 
    ( 54 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           T             m             t
           Z             n             z
           Y             p             y
           X             q             x

 
Example:

   13567  ENTER^
   87168  ENTER^
   24116  ENTER^
   31416  XEQ "PYTH3"   >>>>   6213777201
                                          RDN   5056728720
                                          RDN   4822576232
                                          RDN   9350870225

-You can check that   93508702252 = 62137772012 + 50567287202 + 48225762322 = 87438773964791550625
 

References:

[1]  John H. Conway  & Richard K. Guy - "The Book of Numbers"  - Springer Verlag -  ISBN  0-387-97993-X
[2]  Philippe Descamps & Jean-Jacques Dhenin , "Programmer HP-41" - PSI - ISBN 2-86595-056-5  ( in French )