Number Theory Manual

 Overview
 
 
 
 

XROM	Function	Desciption
16,00  16,01  16,02  16,03  16,04  16,05  16,06  16,07  16,08  16,09  16,10  16,11  16,12  16,13  16,14  16,15  16,16  16,17  16,18  16,19  16,20  16,21  16,22  16,23  16,24  16,25  16,26  16,27  16,28  16,29  16,30  16,31  16,32  16,33  16,34  16,35  16,36  16,37  16,38  16,39  16,40  16,41  16,42  16,43  16,44  16,45  16,46  16,47  16,48  16,49  16,50  16,51  16,52  16,53	-NUMBR THRY   "b-X"   "B-X"   "BELL1"   "BELL2"   "BERN"   "BMN"   "CAT"   "CNK"   "CNW"   "DF"   "DF2"   "DNK"   "EULER"   "EULRN"   "FAREY"   "FIB"   "GCAT1"   "GCAT2"   "GCD"   "GOL"   "HOF"   "KAPR"   "LAH"   "MARKOV"   "MARKOV2"   "MLW"   "NARAY"   "PARTPN"   "PARTQN"   "PELL"   "PELL1"   "PHI"   "PMF"   "PMF2"    PRIME?   "PYTH"   "PYTH2"   "PYTH3"   "SCTL"   "SKN"   "SKNM1"   "SKNM2"   "SNM1"   "SNM2"   "SORT"   "TAU"   "TAUL"   "UV"   "UVW"   "X-b"   "X-B"    XROOT    SDGT	Section Header  Base Conversion ( b < 101 , Integer Arguments )  Base Conversion ( Fractional Arguments )  Bell Numbers  Bell Numbers ( program#2 )  Bernoulli Numbers  Beta Function ( Integer Arguments )  Catalan Numbers  Binomial Coefficients  Conway's Batrachions  Decimal to Fraction  Decimal to Fraction ( program#2 )  Decomposition in Sum of Powers  Euler Numbers  Eulerian Numbers  Farey Sequences  Fibonacci Numbers  Generalized Catalan Numbers  Generalized Catalan Numbers ( program#2 )  GCD , LCM & Simplified Fractions  Golomb Sequence  Hofstadter's Batrachions  Kaprekar Series  Lah Numbers  Markov Numbers  Markov Numbers ( program#2 )  Mallow's Batrachions  Narayana's Numbers  Unrestricted partitions p(n)  Partitions into distinct parts q(n)  Pell Equation  x^2 - d.y^2 = m  Pell Equation  x^2 - d.y^2 = 1  Euler Totient Function ( small arguments )  Prime Factorization, Euler & Moebius Functions  Idem ( displayed factorization )  Is a given number a prime ?  Pythagorean Tripples  Pythagorean Tripples ( program#2 )  Pythagorean Quadruples  Super Catalan Numbers of order m  Divisor Functions  K-Stirling Numbers ( 1st Kind )  K-Stirling Numbers ( 2nd Kind )  Stirling Numbers ( 1st Kind )  Stirling Numbers ( 2nd Kind )  Sorting Numbers ( & Alpha Strings )  Ramanujan Tau Numbers  Ramanujan Tau Numbers ( Multiprecision )  Bezout Identity  a.u + b.v = c  Bezout Identity a.u + b.v + c.w = d  Base Conversion ( b < 101 , Integer Arguments )  Base Conversion ( Fractional Arguments )  X-th Root of Y ( non-orthodox definition... )  Sum of Digits

 

Base Conversion - Integer Arguments and b < 101
 

  "b-X" & "X-b" only deal with ( positive ) integers and bases b < 101
  If  10 < b < 101, each "digit" is represented by 2 decimal digits from the decimal point:
  For instance  N = 1A2EF ( base 16 ) must be keyed in   110021415  ( N = 107247 in decimal )
 

Data Registers: /
Flags: /
Subroutines:  /
 

1-Base b >>> Base 10 ( LBL "b-X" )
 
 

STACK	INPUTS	OUTPUTS
Y	Nb	b
X	b	N10

 
2-Base 10 >>> Base b ( LBL "X-b" )
 
 

STACK	INPUTS	OUTPUTS
Y	N10	b
X	b	Nb

 
Example1:    Convert   (1234)₇    to the decimal scale and then to the nonary scale

  1234  ENTER^
        7  XEQ "b-X"   >>>   466
                                              9   XEQ "X-b"  ( or simply R/S )  >>>   567     whence   (1234)₇  =  466  =   (567)₉

Example2:   Convert  (16267)₁₂  to the decimal scale and then to the hexadecimal scale and the septenary scale.

   106020607  ENTER^
                 12  XEQ "b-X"  >>>  31471
                                                         16   R/S  yields  7101415     whence   (16267)₁₂  =  31471  =  (7AEF)₁₆

-Similarly,    31471  ENTER^   7   XEQ "X-b"  gives   160516     whence   (16267)₁₂  =   (160516)₇
 

Base Conversion - Fractional Arguments
 

 "B-X" & "X-B"  deal with ( positive ) fractional arguments
  If   10 < b < 101,   each "digit" is represented by 2 decimal digits from the decimal point
  If 100 < b < 1001, each "digit" is represented by 3 decimal digits from the decimal point ...  and so on  ...

  For instance, in base 234,   R = 5,233,001.122,   is a number with 3 "digits" in its integer part ( 5  233  1 ) and 1 "digit" in its fractional part ( 122 )
  ( in decimal, R = 328303.5214 )

-Roundoff errors are often unavoidable for fractional numbers
 

Data Registers: /
Flags: /
Subroutines:  /
 

1-Base B >>> Base 10 ( LBL "B-X" )
 
 

STACK	INPUTS	OUTPUTS
Y	Nb	b
X	b	N10

 
2-Base 10 >>> Base B ( LBL "X-B" )
 
 

STACK	INPUTS	OUTPUTS
Y	N10	b
X	b	Nb

 
Example1:    Convert   (12.34)₇    to the decimal scale and then to the nonary scale

  12.34  ENTER^
         7  XEQ "B-X"   >>>   9.510204082
                                                     9           XEQ "X-B"  ( or simply R/S )  >>>   10.45284032

             whence   (12.34)₇  =  9.510204082  =   (10.45284032)₉

Example2:    Convert  (4,123.092,046)₁₂₈    to the decimal scale and then to the scale B = 41

   4,123.092,046  ENTER^
           128           XEQ "B-X"  >>>   635.7215576
                                                                 41            R/S  >>>   15,20.29,23,38

           whence     (4,123.092,046)₁₂₈  =  635.7215576  =  (15,20.29,23,38)₄₁
 

Bell Numbers
 

-Bell numbers are defined by    B(0) = 1  and  B(n) = C_n-1⁰ B(0) + C_n-1¹ B(1) + ........ + C_n-1^n-1 B(n-1)     if  n > 1
  where C_n^k = n!/(k!(n-k)!) are the binomial coefficients.
 

Data Registers:             R00 = B(0) ; R01 = B(1) ; ....... ; Rnn = B(n)
Flags: /
Subroutines: /

-Synthetic registers M , N , O are used by this program
 
 

STACK	INPUTS	OUTPUTS
X	n	B(n)
L	/	n

 
Example:     Calculate  B(10)

     10  XEQ "BELL1"  >>>  B(10) = 115975                        We also have the following results in registers R00 thru R10:
 
 

n	0	1	2	3	4	5	6	7	8	9	10
B(n)	1	1	2	5	15	52	203	877	4140	21147	115975

-If  n > 89 there is an "OUT OF RANGE"
-B(49) = 1.072613714 10⁴⁶  is obtained in 8mn32s.
-Execution time is approximately proportional to n²
 

-Bell numbers may also be computed by the series:  B(n) = e^-1 ( 1ⁿ/1! + 2ⁿ/2! + ...... + kⁿ/k! + ..... )
 

Data Registers: /
Flags: /
Subroutines: /

-Synthetic registers M & N are used by "BELL2"
 
 

STACK	INPUTS	OUTPUTS
X	n	B(n)
L	/	n

 
Example:     Evaluate  B(89)

     89  XEQ "BELL2"  >>>  B(89) = 5.225728472 10⁹⁹            ( the last 3 digits should be 505 instead of 472 )

-Roundoff errors are often greater than with "BELL1"
-On the other hand, "BELL2" is much faster than "BELL1" for large n values.
 

Bernoulli Numbers
 

-The Bernoulli numbers could be computed by the relations:

  B(0) = 1 ;  B(0) + C_n+1¹ B(1) +  C_n+1² B(2) + ...... +  C_n+1ⁿ B(n) = 0   where   C_n^k = n!/(k!(n-k)!)  are the binomial coefficients

-However, this recurrence relation is unstable and the results are quite inaccurate for large n  ( for n = 20 , only 4 digits are correct! )
-"BERN" uses a series expansion instead:

  B(n) = (-1)^-1+n/2 2n!/(2pi)ⁿ  ( 1/1ⁿ + 1/2ⁿ + ...... + 1/kⁿ + ...... )   if  n is even  and  B(0) = 1 ; B(1) = -1/2 ; B(2n+1) = 0  if n > 0

-Actually, B(2) = 1/6 ; B(4) = -1/30 ; B(6) = -1/42  are given directly
 

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

STACK	INPUTS	OUTPUTS
X	n	B(n)

 
Example:     116  XEQ "BERN"  gives  B(116) = -1.748892190 10⁹⁸
 

Beta Function - Integer Arguments
 

 "BMN" employs the relation:

 B(m,n) = (m-1)!(n-1)!/(m+n-1)!
 

Data Registers: /
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Y	m	/
X	n	B(m,n)

 
Example:    100  ENTER^   200  XEQ "BMN"  >>>>   B(100,200) = 3.607285453 10 ^-84
 

Catalan Numbers
 
 

-The Catalan numbers are defined by   C_n = (2n)! / [ n! (n+1)! ] = [1/(n+1)] C_n²ⁿ       where     C_kⁿ = n! / [ k! (n-k)! ]
-The first few are   1  1  2  5  14  42  132  429  .....................
 

Data Registers: /
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Y	/	n
X	n	Cn

 
Examples:

  10  XEQ "CAT"  >>>>   C₁₀  = 16796
  84       R/S          >>>>   C₈₄  = 2.705574509 10⁴⁷
 172      R/S          >>>>  C₁₇₂  = 8.904665844 10⁹⁹
 

Generalized Catalan Numbers
 

-Catalan numbers have many generalizations...
-The 2 following programs calculate the sequence  1  1  1  2  4  8  17  37  82  185  423  ...............   cf  Sloane's  A004148

-The coefficients are given by the recurrence relation:   a(0) = 1 ,  a(n+1) = a(n) + SUM_{k=1,2,...,n-1}  a(k) a(n-1-k)
 

Data Registers:   R00 = a(0)   R01 = a(1)   R02 = a(2)  .........................  Rnn = a(n)
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
X	n	a(n)
L	/	n

 
Examples:

  12  XEQ "GCAT1"  >>>>   a(12) = 2283
  29        R/S              >>>>   a(29) = 8622571758
  41        R/S              >>>>   a(41) = 5.442885873 E14

Notes:

-Execute SIZE at least 004 even if n = 0 , 1 , 2
-All the a(k) are stored in Rkk for k <= n
-Execution time becomes prohibitive for large arguments and another approach is used hereunder:
 

>>>  If you only want one term,  a(n) may also be computed by
 

      a(n) = SUM_{k=m, m+1,.....,n}  ( C_n-k^k  C_n-k+1^k ) / k      for  n > 0         where  m = ceil [(n+1)/2]  and   C_n^k = binomial coefficients
 

Data Registers:     R00 = k        R02 = 1 , 2 , 3 , .......
                                  R01 = n        R03 = ( C_n-k^k  C_n-k+1^k ) / k
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
X	n	a(n)

 
Examples:

    12  XEQ "GCAT2"  >>>>  a(12)  = 2283
    29         R/S             >>>>  a(29)  = 8622571758
    41         R/S             >>>>  a(41)  = 5.442885872 E14
   247        R/S             >>>> a(247) = 4.894614237 E99
 

Super-Catalan Numbers of order m
 

-"SCTL" computes "Super-Catalan numbers" defined as    C(n,m) = C_n²ⁿ C_m^2m / [ 2 C_n^m+n ]    where   C_kⁿ = n! / [ k! (n-k)! ]
-All these terms are integers except  C(0,0) = 1/2
 

Data Register:  R00 = n+m , n+m-1 , ............... , 1 , 0
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	m	C(n,m)

 
Examples:

   103  ENTER^
   208  XEQ "SCTL"  >>>>  C(103,208) = 6.694759032 10⁹⁹

   168  ENTER^  R/S  >>>>  C(168,168) = 3.044358792 10⁹⁹

Notes:

-Though C(n,m) = C(m,n) , small differences in the last places are possible - due to roundoff errors.
-For instance 208  ENTER^  103  R/S  gives  6.694759074 E99  whereas an HP-48 returns  6.69475904935 E99
-If you replace R00 by synthetic register M , "SCTL" becomes a SIZE 000 program ( slightly faster ).
-The first few terms are:

    n\m     0    1    2    3    4    5

     0      .5    1    3   10  35  126
     1       1    1    2    5   14  132
     2       3    2    3    6   14   36
     3      10   5    6   10  20   45
     4      35  14  14  20  35   70
     5     126 42  36  45  79  126

-  C(1,n) = C(n,1)  are "the" Catalan numbers  C_n
 

Binomial Coefficients
 

-This routine is only a slight modification of a program by Keith Jarret listed in "Synthetic Programming Made Easy"
-I've simply added 3 lines to avoid a data error if k = 0 or k = n

-"CNK" computes  C_n^k =  n! / [ k! (n-k)! ]
 

Data Registers: /
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
T	T	n
Z	Z	T
Y	n	Z
X	k	Cnk
L	/	k

 
Example:      n = 100  ,  k = 84

    100  ENTER^
     84   XEQ "CNK"  >>>>   C₁₀₀⁸⁴ = 1.345860628 10¹⁸
 

Batrachions
 

   1°) Conway's Batrachion
 

-The sequence is defined by     C(1) = C(2) = 1  ;   C(n) = C(C(n-1)) + C(n-C(n-1))        ( n > 2 )
 

Data Registers:     R00  unused  ;  R01 = C(1)  ;  R02 = C(2)  ; .......... ;  Rnn = C(n)
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Y	/	1+n.nnn*
X	n	C(n)

   * if n > 2

Example:    100  XEQ "CNW"  >>>>   C(100) = 57                             and  C(n) is in register Rnn  for  n < 101

-The first few terms are
 

n	1	2	3	4	5	6	7	8	9	10
C(n)	1	1	2	2	3	4	4	4	5	6

-We have   C(2^k) = 2^k-1
-The graph of the ratio  C(n)/n  looks like a hopping frog and  C(n)/n tends to 1/2 as n tends to infinity.

C(n)/n
   |       *
   |     *  *          *
   |   *     *      *    *              *
   |  *       *   *        *       *       *
   | *         **            *  *             *
--|*-------*----------*------------*------------------------- n
 

   2°) Hofstadter's Batrachion
 

-The sequence is defined by     H(1) = H(2) = 1  ;   H(n) = H(n-H(n-1)) + H(n-H(n-2))
 

Data Registers:     R00 = 0  ;  R01 = H(1)  ;  R02 = H(2)  ; .......... ;  Rnn = H(n)
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Z	/	1+n.nnn
Y	/	H(n-1)
X	n	H(n)

 
Example:     100  XEQ "HOF"  >>>>  H(100) = 56                            and  H(n) is in register Rnn  for  n < 101

-The firs few terms are:
 

n	1	2	3	4	5	6	7	8	9	10
H(n)	1	1	2	3	3	4	5	5	6	6

 

   3°) Mallows' Batrachion
 

-The sequence is defined by     M(1) = M(2) = 1  ;   M(n) = M(M(n-2)) + M(n-M(n-2))
 

Data Registers:     R00 = 0  ;  R01 = M(1)  ;  R02 = M(2)  ; .......... ;  Rnn = M(n)
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Z	/	1+n.nnn
Y	/	M(n-1)
X	n	M(n)

 
Example:     260  XEQ "MLW"  >>>>  M(260) = 180                      and  M(n) is in register Rnn  for  n < 261

-The first few terms are
 

n	1	2	3	4	5	6	7	8	9	10
M(n)	1	1	2	3	3	4	5	6	6	7

 

Decimal to Fraction
 

-"DF" & "DF2" find the best approximations of a decimal number x by a fraction p/q.
-They are based on continued fractions:

-We define   (p_k)  &  (q_k)   by  q _-1 = 0 ,  q₀ = 1 and  q_k = y_k.q_k-1 + q_k-2

   where         y_k = INT(x_k)             with  x₀ = x ,  x_k+1 = 1/FRC(x_k)     we have then   p_k = RND(x.q_k)   [ RND must be executed in FIX 0 ]
 

-This program displays the successive   p_k/q_k  and stops when the fraction is exactly equal to x
-Synthetic register M may be replaced by any data register.
 

Data Registers: /
Flag:   F29
Subroutines: /
 
 

STACK	INPUTS	when AVIEWs	OUTPUTS
T	/	qk-1	qn-1
Z	/	x	x
Y	/	pk	qn
X	x	qk	pn

 
Examples:

a)  PI  XEQ "DF"  the hp-41 successively displays  "22/7"  "333/106"   "355/113"   "104348/33215"

-When the program stops,  X-register = 104348  &  Y-register = 33215

b) 1  E^X  R/S  >>>  "3/1"  "8/3"  "11/4"  "19/7"  "87/32"  "106/39"  "193/71"  "1264/465"  "1457/536"
                                  "2721/1001"  "23225/8544"  "25946/9545"  "49171/18089"  "173459/63812"

  and eventually,  X = 173459  &  Y =  63812

-Some numbers produce many continued fractions:
-A famous example is the golden ratio phi = (1+sqrt(5))/2   "DF" displays 22 fractions before it stops!
-The last one is 75025/46368

-If you set flag F21, the calculator will stop at each AVIEW
-The numerators are then in Y-register and the denominators in X-register
-At the end, X = the numerator, Y = the denominator.
 

-"DF2" stops when the rounded  difference ( x - p/q )  equals 0 in the current display format.
 

Data Registers: /
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
T	/	qn-1
Z	/	pn/qn - x
Y	/	qn
X	x	pn
L	/	x

 
Examples:

  FIX 6   PI    XEQ "DF2"  >>>>   355   RDN  113   RDN   2.66 E-7   whence  PI ~ 355/113  and   355/113 - PI  =  2.66 E-7
  FIX 6   1    E^X    R/S    >>>>  2721  RDN  1001  RDN  -1.10 E-7   whence  e ~  2721/1001 and  2721/1001 - e = -1.1 E-7

  T = the next to last denominator
  The decimal number x is saved in L-register.
 

Decomposition in Sums of Powers
 

-Given 3 positive integers N , n , k ,  "DNK"  searches n integers  1 <= x₁ <= x₂ <= .......... <= x_n  such that

    N = x₁^k + ............. + x_n^k

-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 = x₁ ,  R06 = x₂ , .....  , R3+n = x_n-1  and   x_n is in register Y

Flag:   F07
Subroutines: /
 
 

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 x_i , remember that R05 = x₁ ,  R06 = x₂ , .....  , R3+n = x_n-1 and   x_n is in register Y
 

Euler Numbers
 

-These integers  may be computed by the formula:

     E(n) = (-1)^n/2 2ⁿ⁺² n! pi ^-n-1 ( 1/1ⁿ⁺¹ - 1/3ⁿ⁺¹ + 1/5ⁿ⁺¹ - ..... + (-1)^k / (2k+1)ⁿ⁺¹ + .....  )     if  n is even
     E(n) = 0  if  n is odd

-"EULER" gives E(n) directly if n < 7 or n odd.
-The series expansion is used for n > 7 ; n even.
-The first values are:
 

n	0	2	4	6	8	10
E(n)	1	-1	5	-61	1385	-50521

Data Registers:   R00 = n ; R01 to R02: temp
Flag: /
Subroutine: /
 
 

STACK	INPUTS	OUTPUTS
X	n	E(n)

 
Example:      Find   E(78)

    78   XEQ "EULER"  yields  E(78) =  -7.270601736 10⁹⁹
 

Eulerian Numbers  A(n,k)
 

-The integers  A(n;k) are computed by:   A(n;k) =  C_n+1⁰  kⁿ  -  C_n+1¹  (k-1)ⁿ  + C_n+1²  (k-2)ⁿ - .......... + (-1)^k  C_n+1^k  (k-k)ⁿ    ( 0 < k < n+1 )
-The first values are:

        1
        1   1
        1   4     1
        1  11   11     1
        1  26   66    26    1
        1  57  302  302  57  1
        .....................................

-Note that some authors define these numbers differently ( with k+1 instead of k )
 

Data Registers: /
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Y	n	n
X	k	A(n;k)
L	/	k

 
Example:     Calculate  A(16;7)

    16  ENTER^
     7   XEQ "EULN3"  >>>   A(16;7) = 3.207483180 10¹²
 

Farey Sequences
 

-The Farey sequence F_n is the list of  N rational number  { a₁/b₁ = 0/1 ; a₂/b₂ = 1/n ; ....... ; a_N/b_N = 1/1 } sorted in increasing order
  whith  0 <= a_k <= b_k <= n  ;  GCD(a_k;b_k) = 1

-The following program uses the formulae:  a_k+1 =  a_k.m_k - a_k-1  ;  b_k+1 =  b_k.m_k - b_k-1    where  m_k = INT( (b_k-1+ n)/b_k )

  F₁ = { 0/1 ; 1/1 }
  F₂ = { 0/1 ; 1/2 ; 1/1 }
  F₃ = { 0/1 ; 1/3 ; 1/2 ; 2/3 ; 1/1 }
  F₄ = { 0/1 ; 1/4 ; 1/3 ; 1/2 ; 2/3 ; 3/4 ; 1/1 }
  .............................. etc .............................
 

Data Registers:  R00 = n ;  R01 = a_k-1 ; R02 = b_k-1 ; R03 = k ; When the routine stops, R03 = X-register = N = the number of elements in the list.
Flag:   F29
Subroutines: /
 
 
 

STACK	INPUTS	... AVIEW ...	OUTPUTS
Y	/	bk	/
X	n	ak	N(n)

  where N(n) = the number of fractions

Example:

   5   XEQ "FAREY"  >>>>   0/1  1/5  1/4  1/3  2/5  1/2  3/5  2/3  3/4  4/5  1/1

  and     X-register = 11 = the number of terms in this Farey sequence = R03.
 

Fibonacci Numbers
 

-This famous sequence is defined by

  u₀ = 0  ,   u₁ = 1  &   u_n+1 = u_n + u_n-1
 

Data Registers: /
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Y	/	un-1
X	n	un
L	/	n

  where n is a non-negative integer

Example:

  49  XEQ "FIB"  >>>>   u₄₉ = 7778742049   X<>Y   u₄₈ = 4807526976
 480      R/S        >>>>   u₄₈₀ = 9.216845715 E99   X<>Y   u₄₈ = 5.696323921 E99

Note:

-This routine uses a dummy  u _-1 = 1
 

Greatest Common Divisor & Lowest Common Multiple
 

-"GCD" computes GCD(a,b) & LCM(a,b)  and  calculates  a' = a/GCD(a,b)  & b' = b/GCD(a,b)
-Thus,  a'/b' = a/b  and  a'/b'  is irreductible.
-The Euclidean algorithm is used.
 

Data Registers: /
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
T	/	a'= a/gcd(a,b)
Z	/	b'= b/gcd(a,b)
Y	a	lcm(a,b)
X	b	gcd(a,b)

 
Example:

  950796  ENTER^
  107016  XEQ "GCD"  >>>>  GCD(950796,107016) =   4116
                                      RDN   LCM(950796,107016) = 24720696
                                     RDN       b' = 26
                                      RDN       a' = 231

-Thus,   950796/107016 = 231/26
 

Golomb Sequence
 

-"GOL" uses the following method, let:

     a_m = the greatest integer n such that  G( n ) = m
     b_m = the greatest integer n such that  G(G( n )) = m
     c_m = the greatest integer n such that  G(G(G( n ))) = m

-We have   a₁ = b₁ = c₁ = 1 and   a_m = a_m-1 + g_m ;   b_m = b_m-1 + m.g_m ;  c_m = c_m-1 + m.g_m( a_m-( g_m-1 )/2 )
  and let  cc_m,k  defined by   cc_m,0 =  c_m  ;    cc_m,k+1 =  cc_m,k + (m+1)(a_m+k+1)

-First, we find the greatest integer m such that  c_m < n
-Then, --------------------------  k  --------   cc_m,k < n
-Finally,  G(n) = b_m + k(m+1) + INT(( n - cc_m,k + a_m +k )/( a_m + k +1 ))

-For n = 10¹⁰ , m = 330 and fortunately, a very simple formula yields G(m) for m < 423   i-e  G(m) = the nearest integer to 1.2 m^0.61832
 

Data Registers:   R00 = n  ;   R01 = a_m and  a_m + k + 1  ;  R02 = b_m and  b_m + k(m+1)(k+1)  ; R03 = g_m
                                                 R04 = m+1  ;  R05 = 0.61832  ;  R06 = 1.2  ;  R07 = 0.5
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
X	n	G(n)

 
Examples:

   9,999,999,999  XEQ "GOL"  >>>>  1,820,598
      1,000,000             R/S         >>>>       6,137
 

Kaprekar Series
 

-Take for example a 3-digit number, say  N₀ = 286
-Sort the digits in increasing order:  268  and in decreasing order:  862
-Calculate the difference  862 - 268 = 594 = N₁

-Repeat this operation:  954 - 459 = 495 = N₂
-Again:  954 - 459 = 495 = N₃

-The Kaprekar series in this example is  286 , 594 , 495 , 495 , .....
-So we end up in a 1-number cycle:  { 495 }

-The following programs allows to find such cycles.
 
 

Data Registers:           •  R00 = n = number of digits  ( 1 <= n <= 10 )                   ( Register R00 is to be initialized before executing "KAPR"  )

                                         R01 thru R10 = digits of the numbers    R11 & R12: temp

                                     -  R13  R14  R15  .....................  =  N₀  N₁  N₂  .............................
Flags: /
Subroutine:   "SORT"
 
 

STACK	INPUT	OUTPUT
X	N0	bbb.eee

   where  bbb.eee  is the control number of the cycle

Example1:     n = 6  &   N₀ = 918682

     6       STO 00
918682  XEQ "KAPR"  >>>>  the HP-41 displays  N₀ ,  N₁ ,  N₂ , ...... and finally returns  24.030

-So we eventually end up in a 7-number cycle in registers  R24 , ..... , R30  namely:

       { 851742 , 750843 , 840852 , 860832 , 862632 , 642654 , 420876 }
 

Example2:     n = 5  &   N₀ = 4    which will be read  "00004"  since  n = 5

    5   STO 00
    4     R/S        >>>>   15.018

-The 4-number cycle is in registers  R15 thru R18:

       { 62964 , 71973 , 83952 , 74943 }
 

Lah Numbers
 

-"LAH" computes the integers defined by  L(n,m) = (-1)ⁿ ( n! / m! )  C_m-1^n-1        for  0 < m <= n
 

Data Registers: /
Flags: /
Subroutine:  "CNK"  ( Binomial Coefficients )
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	m	L(n,m)

 
Example:

   41  ENTER^
   24  XEQ "LAH"  >>>>   L(41,24) = -4.784156505 E36
 

Markov Numbers
 

-Markov numbers are integer solutions of the equation:  x² + y² + z² = 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: /
 
 
 

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 variant 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: /
 
 

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

 
Example:

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

Narayana's Numbers
 

-They are defined as  N_n,k = (1/n) C_n^k  C_n^k-1    with  1 <= k <= n
 

Data Registers: /
Flags: /
Subroutine:  "CNK"  ( Binomial Coefficients )
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	k	Nn,k

 
Example:

  128   ENTER^
   67    XEQ "NARAY"  >>>>   N_128,67 = 3.663413794 E72
 

Unrestricted Partitions:   p(n)
 

Formula:      p(n) = 1/(pi) SUM_{k=1,2,3,.....} k^1/2 A_k(n).( ( cosh(pi(2(n-1/24)/3)^1/2/k) (2/3)^1/2 pi/k )/(n-1/24) -  sinh(pi(2(n-1/24)/3)^1/2/k)/(n-1/24)^3/2 )/2

        with     A_k(n) = SUM_{h=1,2,...,k ; GCD(h,k)=1}  cos( pi.S(h,k) - 2(pi)n.h/k )    where  S(h,k) = SUM_{j=0,1,...,k-1} ( frc(h.j/k) -1/2 )( j/k - 1/2 )
 

Data Registers:             R00 = n ; R01 = partial sums and p(n) ;  R02 to R08 = temp
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
X	n	p(n)

 
  Examples:

     100   XEQ "PARTPN"  >>>>   p(100) =  190569292
      721             R/S            >>>>   p(721) = 1.610617578  10²⁶   ( the last digits should be 57 instead of 78 )
 

Partitions Into Distinct Parts:   q(n)
 

Formula:      q(n) =  SUM_{k=1,2,3,.....}   A_2k-1(n).  I₁(pi((n+1/24)/3)^1/2/(2k-1)) (1/3)^1/2 pi/(2k-1) )/(2(n+1/24)^1/2)

        where     A_k(n) and S(h,k)  are defined as above ( 1°)b) )  and  I₁(x) = the modified Bessel function of order 1.
 

Data Registers:             R00 = n ; R01 = partial sums and q(n) ;  R02 to R06 = temp
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
X	n	q(n)

 
Examples:

     200  XEQ "PARTQN"  >>>>  q(200) = 487067755                                ( exact value is 487067746 )
    1000           R/S             >>>>  q(1000) = 8.635565946  10²¹
 

Pell Equation:  x² - d y² = m
 

Data Registers: /
Flag:  F29
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
T	/	m
Z	/	y
Y	d	d
X	m	0
L	/	x

 
Example:                  x² - 7 y² = 2

  7  ENTER^
  2  XEQ "PELL"  >>>>   "3^2-7*1^2=2"
                              R/S     "45^2-7*17^2=2"
                              R/S     "717^2-7*271^2=2"
                              R/S     ..............................

Notes:

-Do not disturb the stack before pressing R/S again.
-Execution time is often prohibitive and the following approach is better for m = 1.
 

Pell Equation:  x² - d y² = 1
 

-Here, d must not be a perfect square.
-This routine uses the continued fraction of sqrt(d).

-"PELL1" takes d in X-register and returns the smallest solution x , y in registers X , Y.
 

Data Registers:    R00 = d     R01 thru R10: temp
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Y	/	y
X	d	x

 
Examples:

   41  XEQ "PELL1"   >>>>   x = 2049
                                  X<>Y   y = 320

   61       R/S       >>>>   x = 1766319049
                          X<>Y   y =  226153980
 

Euler Totient Function
 

-phi(n) is the number of integers not exceeding and relatively prime to n.
-"PHI" can be used if n is small
 

Data Registers:   R00 = phi(n)   R01 = n    R02 = n-1 , n-2 , .... , 0
Flags: /
Subroutines: /
 
 

STACK	INPUT	OUTPUT
X	n	phi(n)

 
Examples:

   1    XEQ "PHI"  >>>>   phi(1)    =  1
 360        R/S        >>>>  phi(360) = 96
 

-The programs hereunder are of course much faster.
 

Prime Factorization, Euler & Moebius Functions
 

-"PMF" stores the factors into contiguous registers from register R04 and it computes phi(n) & the Moebius function µ(n).

   µ(0) = 0   µ(1) = 1

 otherwise,

         µ(n) =   0      if n is a multiple of the square of a prime
         µ(n) = (-1)^k  if n is the product of k distinct primes.

-The factors themselves are stored as follows: if, for instance n = 3⁵ 7⁴     R04 = 3.005   R05 = 7.004
-In other words, the exponents are divided by 1000 and added to the corresponding prime factor p.
-If p > 9999999 , the fractional part will be zero but should be read 0.001
 

Registers:          R00:  temp
                           R01 = phi(N)    R02 = µ(n)    R03 = 4.eee   R04 , R05 , ...  the factors
Flags: /
Subroutines:  the M-code routine  PRIME?
 
 

STACK	INPUTS	OUTPUTS
Z	/	µ(n)
Y	/	phi(n)
X	n	4.eee

  Exception:   if n = 1 ,  X-output = 1

Example:

  n =  3238704  XEQ "PMF"  >>>>       4.007
                                               RDN    phi(n) =  870912
                                               RDN     µ(n)  =   0

  { R04  R05  R06  R07 } = { 2.004  3.005  7.002  17.001 }   whence  3238704 = 2⁴ 3⁵ 7² 17
 

-"PMF2"  calculates phi(n) & µ(n) too but the factorization is gradually displayed.
 

Registers:          R00&R03:  temp
                           R01 = phi(N)    R02 = µ(n)
Flags: /
Subroutines:  the M-code routine  PRIME?
 
 

STACK	INPUTS	OUTPUTS
Y	/	µ(n)
X	n	phi(n)

  where  n  is a positive integer

Example:

  3238704  XEQ "PMF2"  displays  ( successively )

           3238704=2^4*
           3238704=2^4*3^5*
           3238704=2^4*3^5*7*2*
           3238704=2^4*3*5*7^2*17^1        ( if there are more than 24 characters, the left part of the alpha string will be gradually truncated )
 

     phi(3238704) = 870912 is in X-register
       µ(3238704) =      0       is in Y-register
 

Primality Test
 

-PRIME? is an M-code function written by Jason DeLooze.
-It comes from Ángel Martin's SandMath.
-In addition to doing the "skip-if-false" in PRGM mode,
  this routine returns a prime factor if the number isn't prime, or 1 if prime.
 

Examples:

  997  XEQ "PRIME?"  gives  "YES"  and  X-register = 1
  245  XEQ "PRIME?"  gives   "NO"  and  X-register = 5

-The number itself is saved in L-register.
 

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 "PYTH" returns in X , Y , Z  a solution of the Pythagoras equation
 

Data Registers: /
Flags: /
Subroutines: /
 
 

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, "PYHT2" takes 2 positive integers p &q  in X & Y and returns a Pythagorean tripple.
 

Data Registers: /
Flags: /
Subroutines: /
 
 

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 = | p² - q² |  ,  z =  p² + q²
-It can be proved that this method gives all the solutions too !
 

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

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

    x = ( m² + n² - p² ) / p   ,  y = 2 m  ,  z = 2 n ,  t = ( m² + n² + p² ) / p

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

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

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

Divisor Functions
 

-After factorizing n > 1 with "PMF", "SKN" can be used to compute s_k(n) = Sum_{d | n} d ^k

  Formula:     if  n = p₁^r1 ....... p_1m^rm  ,  s_k(n) = Product_i [ p_i^k(ri+1) - 1 ] / ( p_i^k - 1 ) = Product_i ( 1 + p_i^k + p_i^2k + ...... + p_i^k.ri )
 

Data Registers:   R00: temp   R01 = k   R02 = s_k(n)
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Y	bbb.eee	/
X	k	sk(n)

   where  bbb.eee  is the control number of the factorization         bbb > 002

Example:

  n = 3238704 XEQ "PMF" gives  4.007

  4.007  ENTER^
      1     XEQ "SKN"  >>>>  s₁(n) = 11577384

  4.007  ENTER^
      2         R/S            >>>>  s₂(n) = 1.610126288 10¹³

Note:

 s_k(1) = 1 for all k
 

Stirling Numbers
 

-Stirling numbers of the first kind S_n^m are defined by the recurrence relation: ( S_n⁰ = 0 ) ;  S_n^m = S_n-1^m-1 - ( n-1 ) S_n-1^m      ( 1 <= m <= n )
-Stirling numbers of the second kind s_n^m are defined by the recurrence relation:  ( s_n⁰ = 0 )  ;   s_n^m =  s_n-1^m-1  +  m s_n-1^m      ( 1 <= m <= n )

  •  (-1)^n-m S_n^m is the number of permutations of n symbols which have exactly m cycles.
  •             s_n^m  is the number of ways of partitioning a set of n elements into m non-empty subsets.
 

Data Registers:       When the program stops:   R00 = 0 , R01 = S_n¹ , R02 = S_n² , ...... , Rmm = S_n^m       ( or s_n^j  )

Flag:   F02
Subroutines:  /
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	m <= n	S(n;m)

 
Example:   Calculate  S₁₂⁷   and    s₁₂⁷

        12  ENTER^
          7  XEQ "SNM1"  >>>>   S₁₂⁷ =  -2637558

     ( and R01 = -39916800 = S₁₂¹ ; R02 = 120543840 = S₁₂²    .........................................      R07 = -2637558 = S₁₂⁷  )

        12  ENTER^
         7  XEQ "SNM2"  >>>>   s₁₂⁷ = 627396

     ( and R01 = 1 =  s₁₂¹ ; R02 = 2047 =  s₁₂² ; ........ ;  R07 = 627396 = s₁₂⁷ )
 

K-Stirling Numbers
 

  •    k-Stirling numbers of the first kind S₁(k;n,m) are defined by the recurrence relation:

                 S₁(k;0,0) = 0     ;      S₁(k;n,m)  =  [ (1-k).m + 1 - n ] S₁(k;n-1,m) + S₁(k;n-1,m-1)    ( 1 <= m <= n )

  •    k-Stirling numbers of the second kind S₂(k;n,m) are defined by:

                 S₂(k;0,0) = 0     ;      S₂(k;n,m)  =  [ (k-1).(n-1) + m ] S₂(k;n-1,m) + S₂(k;n-1,m-1)    ( 1 <= m <= n )
 

Data Registers:      When the program stops:   R00 = 0 , R01 = S(k;n,1) , R02 = S(k;n,2) , ...... , Rmm = S(k;n,m)

Flags:  F02
Subroutines:  /
 
 

STACK	INPUTS	OUTPUTS
Z	k	/
Y	n	/
X	m <= n	S(k;n,m)

      ( k may be positive, zero or negative )

Example:

    3  ENTER^  10  ENTER^  7   XEQ "SKNM1"   >>>>   S₁(3;10,7) = -188370
    3  ENTER^  10  ENTER^  7   XEQ "SKNM2"  >>>>   S₂(3;10,7) =  220500

-In both cases,      R01 = S_i(3;10,1)   R02 = S_i(3;10,2)  ..................  R07 = S_i(3;10,7)    ( i = 1 or 2 )

Note:    S₁(1;n,m) = S_n^m     ;      S₂(1;n,m) = s_n^m
 

Sorting Numbers & Alpha Strings
 

-"SORT" is called as a subroutine by "KAPR"
 

Data Registers:   Only those of the block of registers.
Flags: /
Subroutines: /
 
 

STACK	INPUT	OUTPUT
X	bbb.eeeii	/

      where  bbb.eeeii  is the control number of the array.

-For example, to sort numbers in registers R01  R03  R05  R07 , key in

  1.00702  XEQ "SORT"

-This program uses a selection sort.
-The numbers are sorted in increasing order.

-It also works to sort alpha strings.
 

Ramanujan Tau Numbers
 

-Ramanujan numbers Tau(n) are defined by

   x.[ (1-x).(1-x²).(1-x³)....... ]²⁴  = Tau(1).x + Tau(2).x² + Tau(3).x³ + ..... + Tau(n).xⁿ + ......

-"TAU" calculates Tau(n) by the recurrence relation:

    (n-1).Tau(n) = SUM_1<=m<=b(n) (-1)^m+1 (2m+1).(n-1-9m(m+1)/2).Tau(n-m(m+1)/2)     where   b(n) = INT ((8n+1)^1/2-1)/2

-This program calculates and stores Tau(1) , Tau(2) , ...... , Tau(n) into registers R01 , R02 , ....... , Rnn
 

Data Registers:  R00 = 0 , R01 = Tau(1) , ........... , Rnn = Tau(n)
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
X	n	tau(n)
L	/	n

 
Example:    Evaluate Tau(10)

   10  XEQ "TAU"  >>>>  Tau(10) = -115920  in X-register and in R10

-We also have   1    , -24 ,  252 , -1472 , 4830 , -6048 , -16744 , 84480 , -113643  ( i-e Tau(1) , .......... , Tau(9) )
    in registers   R01 , R02 , R03 ,   R04  ,  R05 ,    R06 ,     R07  ,    R08  ,     R09       respectively

Notes:

-This program yields Tau(84) = 6,188,510,373 whereas the correct value is Tau(84) = 6,211,086,336
  and 300 XEQ "TAU" produces a completely meaningless result !
-Obviously, roundoff errors overwhelm the wanted function.

-Results are exact up to n = 43 only.
-Therefore, another approach is necessary to compute Tau(n) for larger arguments.

-"TAUL" employs the formula:

    Tau(n) = n⁴.s(n) - 24 SUM_0<k<n k².(35.k²-52.k.n+18.n²).s(k).s(n-k)     ( n > 1 )  where  s(k) is the sum of the divisors of k.

-Multiple precision is used and Tau(n) may be exactly computed up to n = 13868.

-"TAUL" computes Tau(n) by groups of 4 digits in registers R15 R16 R17 R18 R19 R20
 

Data Registers:        R01 = n   ;   R00 and R02 thru R20 are used for temporary data storage and when the program stops:

                                   R03 = Tau(n) = X-register ( approximate value )
                                   R15 to R20 contain the digits of Tau(n) ( exact value )
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Y	/	15.020
X	n	tau(n)

  15.020  = control number of the exact result ( in registers R15 thru R20 )
  X-output is only an approximate value

Example1:     Calculate  Tau(84)

   84  XEQ "TAUL"  >>>>  Tau(84) =  6211086336

In this case, the approximate value is quite exact!

which is confirmed by   R15 = 0 , R16 = 0 , R17 = 0 , R18 = 62 , R19 = 1108 , R20 = 6336

Example2:    Compute Tau(314)

  314  XEQ "TAUL"  >>>>  Tau(314) = -3.156280211 10¹³  = approximate value

and we have:  R15 = 0 , R16 = 0 , R17 = -31 , R18 = -5628 , R19 = -210 ( which is to be read -0210 ) , R20 = -5744

 whence   Tau(314) = -31562802105744

-Remark that registers R15 to R20 are inferior or equal to zero when Tau(n) < 0

Notes:

-With V41 in turbo mode  13868 XEQ "TAUL"  yields:

   R15 = -336 ; R16 = -1948 ; R17 = -0538 ; R18 = -4247 ; R19 = -3535 ; R20 = -1552   ( in about 1h23mn on my PC )

  whence  Tau(13868) = -33619480538424735351552      with a "real" HP-41, execution time is probably of the order of  8 or 9 days!

-For n > 13868 , the calculations involve numbers greater than 10¹⁰ and the results would be only approximate.
 

Bezout Identity  a.u + b.v = c
 

-The equation  a u + b v = c   where  a , b , c  are known integers and u , v are 2 unknowns
  has integer solutions iff c is a multiple of GCD(a,b)
-"UV" computes 2 numbers u & v that satisfy this equation.
-GCD(a,b) is also returned in Z-register & T-register.

-The extended Euclidean algorithm is used:  let  a₀ = a  &  a₁ = b , we perform the successive Euclidean divisions

    a₀ = a₁ q₁ + a₂                     and we                u₀ = 1  ,  u₁ = 0  ,  u_k+1 = u_k-1 - q_k u_k
    a₁ = a₂ q₂ + a₃                      define                 v₀ = 0  ,  v₁ = 1  ,  v_k+1 = v_k-1 - q_k v_k
    ......................

  a_k-1 = a_k q_k + a_k+1
   .......................

  a_n-1 = a_n q_n + 0                      it yields                   a_n = gcd(a,b)        u = c u_n/ a_n         v = c v_n/ a_n
 
 

Data Registers:      R00 thru R04: temp           When the program stops,   R02 = -b & R04 = a  or   R02 = b & R04 = -a
Flags: /
Subroutines: /
 
 

STACK	INPUTS	OUTPUTS
Z	a	gcd(a,b)
Y	b	v
X	c	u

 
Example:

  125  ENTER^
   41   ENTER^
    1    XEQ "UV"  >>>>       u = -20                            -20x125 + 61x41 = 1
                             RDN        v =  61
                             RDN   gcd(125,41) = 1

-If c is not a multiple of gcd(a,b)  u and v are not integers.
-When c = gcd(a,b)  , u and v are minimal solutions. Otherwise, they are not!

-For instance, with  a = 125 , b = 41 , c = 7   "UV" gives  u = -140 & v = 427
 

Bezout Identity a.u + b.v + c.w = d
 

-"UVW" 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: /
 
 

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

-When the program stops, we are inside a local subroutine.
 

X-th Root
 

-This M-code routine calculates the x-th root of y if y is non-negative.
-Otherwise, it computes sign(y) [abs(y)^(1/x)]
-So, it's not an orthodox XROOT !

-It will return correctly cube root ( -64 ) = -4   exactly   ( 13-digit routines allows to get exact integers in these cases )
-But it will also give square root ( -9 ) = -3 !!!

-This function is actually a continuation of  y^(1/x)  by a symmetry with respect to the origine O.
-It's used by "DNK" to test if a positive integer is really an integer power of an integer.
 

Sum of Digits
 

 SDGT  is an M-Code routine written by Ángel Martin.
 It returns the sum of the digits of the mantissa of a real number in X-register.
 
 

STACK	INPUTS	OUTPUTS
Y	/	R
X	R	Sum of Digits

 
Example:

      123456789      XEQ "SDGT"  yields   45
  1.23456789 E15  XEQ "SDGT"  yields   45   too.
 

References:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Philippe Descamps & Jean-Jacques Dhenin , "Programmer HP-41" - PSI - ISBN 2-86595-056-5  ( in French )
[3]  John H. Conway  & Richard K. Guy - "The Book of Numbers"  - Springer Verlag -  ISBN  0-387-97993-X
[4]  Clifford A. Pickover - "Keys to Infinity" - John Wiley & Sons -  ISBN  0-471-11857-5