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)7    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)7  =  466  =   (567)9

Example2:   Convert  (16267)12  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)12  =  31471  =  (7AEF)16

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

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)7    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)7  =  9.510204082  =   (10.45284032)9

Example2:    Convert  (4,123.092,046)128    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)128  =  635.7215576  =  (15,20.29,23,38)41
 

Bell Numbers
 

-Bell numbers are defined by    B(0) = 1  and  B(n) = Cn-10 B(0) + Cn-11 B(1) + ........ + Cn-1n-1 B(n-1)     if  n > 1
  where Cnk = 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 1046  is obtained in 8mn32s.
-Execution time is approximately proportional to n2
 

-Bell numbers may also be computed by the series:  B(n) = e-1 ( 1n/1! + 2n/2! + ...... + kn/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 1099            ( 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) + Cn+11 B(1) +  Cn+12 B(2) + ...... +  Cn+1n B(n) = 0   where   Cnk = 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)n  ( 1/1n + 1/2n + ...... + 1/kn + ...... )   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 1098
 

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   Cn = (2n)! / [ n! (n+1)! ] = [1/(n+1)] Cn2n       where     Ckn = 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"  >>>>   C10  = 16796
  84       R/S          >>>>   C84  = 2.705574509 1047
 172      R/S          >>>>  C172  = 8.904665844 1099
 

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) + SUMk=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) = SUMk=m, m+1,.....,n  ( Cn-kk  Cn-k+1k ) / k      for  n > 0         where  m = ceil [(n+1)/2]  and   Cnk = binomial coefficients
 

Data Registers:     R00 = k        R02 = 1 , 2 , 3 , .......
                                  R01 = n        R03 = ( Cn-kk  Cn-k+1k ) / 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) = Cn2n Cm2m / [ 2 Cnm+n ]    where   Ckn = 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 1099

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

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  Cn
 

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  Cnk =  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"  >>>>   C10084 = 1.345860628 1018
 

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(2k) = 2k-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   (pk)  &  (qk)   by  q -1 = 0 ,  q0 = 1 and  qk = yk.qk-1 + qk-2

   where         yk = INT(xk)             with  x0 = x ,  xk+1 = 1/FRC(xk)     we have then   pk = RND(x.qk)   [ RND must be executed in FIX 0 ]
 

-This program displays the successive   pk/qk  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 <= 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: /
 
 
      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
 

Euler Numbers
 

-These integers  may be computed by the formula:

     E(n) = (-1)n/2 2n+2 n! pi -n-1 ( 1/1n+1 - 1/3n+1 + 1/5n+1 - ..... + (-1)k / (2k+1)n+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 1099
 

Eulerian Numbers  A(n,k)
 

-The integers  A(n;k) are computed by:   A(n;k) =  Cn+10  kn  -  Cn+11  (k-1)n  + Cn+12  (k-2)n - .......... + (-1)k  Cn+1k  (k-k)n    ( 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 1012
 

Farey Sequences
 

-The Farey sequence Fn is the list of  N rational number  { a1/b1 = 0/1 ; a2/b2 = 1/n ; ....... ; aN/bN = 1/1 } sorted in increasing order
  whith  0 <= ak <= bk <= n  ;  GCD(ak;bk) = 1

-The following program uses the formulae:  ak+1 =  ak.mk - ak-1  ;  bk+1 =  bk.mk - bk-1    where  mk = INT( (bk-1+ n)/bk )

  F1 = { 0/1 ; 1/1 }
  F2 = { 0/1 ; 1/2 ; 1/1 }
  F3 = { 0/1 ; 1/3 ; 1/2 ; 2/3 ; 1/1 }
  F4 = { 0/1 ; 1/4 ; 1/3 ; 1/2 ; 2/3 ; 3/4 ; 1/1 }
  .............................. etc .............................
 

Data Registers:  R00 = n ;  R01 = ak-1 ; R02 = bk-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

  u0 = 0  ,   u1 = 1  &   un+1 = un + un-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"  >>>>   u49 = 7778742049   X<>Y   u48 = 4807526976
 480      R/S        >>>>   u480 = 9.216845715 E99   X<>Y   u48 = 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:

     am = the greatest integer n such that  G( n ) = m
     bm = the greatest integer n such that  G(G( n )) = m
     cm = the greatest integer n such that  G(G(G( n ))) = m

-We have   a1 = b1 = c1 = 1 and   am = am-1 + gm ;   bm = bm-1 + m.gm ;  cm = cm-1 + m.gm( am-( gm-1 )/2 )
  and let  ccm,k  defined by   ccm,0 =  cm  ;    ccm,k+1 =  ccm,k + (m+1)(am+k+1)

-First, we find the greatest integer m such that  cm < n
-Then, --------------------------  k  --------   ccm,k < n
-Finally,  G(n) = bm + k(m+1) + INT(( n - ccm,k + am +k )/( am + k +1 ))

-For n = 1010 , m = 330 and fortunately, a very simple formula yields G(m) for m < 423   i-e  G(m) = the nearest integer to 1.2 m0.61832
 

Data Registers:   R00 = n  ;   R01 = am and  am + k + 1  ;  R02 = bm and  bm + k(m+1)(k+1)  ; R03 = gm
                                                 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  N0 = 286
-Sort the digits in increasing order:  268  and in decreasing order:  862
-Calculate the difference  862 - 268 = 594 = N1

-Repeat this operation:  954 - 459 = 495 = N2
-Again:  954 - 459 = 495 = N3

-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  .....................  =  N0  N1  N2  .............................
Flags: /
Subroutine:   "SORT"
 
 
      STACK         INPUT       OUTPUT
           X            N0        bbb.eee

   where  bbb.eee  is the control number of the cycle

Example1:     n = 6  &   N0 = 918682

     6       STO 00
918682  XEQ "KAPR"  >>>>  the HP-41 displays  N0 ,  N1 ,  N2 , ...... 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  &   N0 = 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 ( n! / m! )  Cm-1n-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:  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: /
 
 
 
      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  Nn,k = (1/n) Cnk  Cnk-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"  >>>>   N128,67 = 3.663413794 E72
 

Unrestricted Partitions:   p(n)
 

Formula:      p(n) = 1/(pi) SUMk=1,2,3,..... k1/2 Ak(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     Ak(n) = SUMh=1,2,...,k ; GCD(h,k)=1  cos( pi.S(h,k) - 2(pi)n.h/k )    where  S(h,k) = SUMj=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  1026   ( the last digits should be 57 instead of 78 )
 

Partitions Into Distinct Parts:   q(n)
 

Formula:      q(n) =  SUMk=1,2,3,.....   A2k-1(n).  I1(pi((n+1/24)/3)1/2/(2k-1)) (1/3)1/2 pi/(2k-1) )/(2(n+1/24)1/2)

        where     Ak(n) and S(h,k)  are defined as above ( 1°)b) )  and  I1(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  1021
 

Pell Equation:  x2 - d y2 = m
 

Data Registers: /
Flag:  F29
Subroutines: /
 
 
      STACK        INPUTS      OUTPUTS
           T             /            m
           Z             /             y
           Y             d             d
           X             m             0
           L              /             x

 
Example:                  x2 - 7 y2 = 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:  x2 - d y2 = 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 = 35 74     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 = 24 35 72 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 = | p2 - q2 |  ,  z =  p2 + q2
-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 = ( 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: /
 
 
      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 sk(n) = Sumd | n d k

  Formula:     if  n = p1r1 ....... p1mrm  ,  sk(n) = Producti [ pik(ri+1) - 1 ] / ( pik - 1 ) = Producti ( 1 + pik + pi2k + ...... + pik.ri )
 

Data Registers:   R00: temp   R01 = k   R02 = sk(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"  >>>>  s1(n) = 11577384

  4.007  ENTER^
      2         R/S            >>>>  s2(n) = 1.610126288 1013

Note:

 sk(1) = 1 for all k
 

Stirling Numbers
 

-Stirling numbers of the first kind Snm are defined by the recurrence relation: ( Sn0 = 0 ) ;  Snm = Sn-1m-1 - ( n-1 ) Sn-1m      ( 1 <= m <= n )
-Stirling numbers of the second kind snm are defined by the recurrence relation:  ( sn0 = 0 )  ;   snm =  sn-1m-1  +  m sn-1m      ( 1 <= m <= n )

    (-1)n-m Snm is the number of permutations of n symbols which have exactly m cycles.
               snm  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 = Sn1 , R02 = Sn2 , ...... , Rmm = Snm       ( or snj  )

Flag:   F02
Subroutines:  /
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             /
           X         m <= n         S(n;m) 

 
Example:   Calculate  S127   and    s127

        12  ENTER^
          7  XEQ "SNM1"  >>>>   S127 =  -2637558

     ( and R01 = -39916800 = S121 ; R02 = 120543840 = S122    .........................................      R07 = -2637558 = S127  )

        12  ENTER^
         7  XEQ "SNM2"  >>>>   s127 = 627396

     ( and R01 = 1 =  s121 ; R02 = 2047 =  s122 ; ........ ;  R07 = 627396 = s127 )
 

K-Stirling Numbers
 

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

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

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

                 S2(k;0,0) = 0     ;      S2(k;n,m)  =  [ (k-1).(n-1) + m ] S2(k;n-1,m) + S2(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"   >>>>   S1(3;10,7) = -188370
    3  ENTER^  10  ENTER^  7   XEQ "SKNM2"  >>>>   S2(3;10,7) =  220500

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

Note:    S1(1;n,m) = Snm     ;      S2(1;n,m) = snm
 

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-x2).(1-x3)....... ]24  = Tau(1).x + Tau(2).x2 + Tau(3).x3 + ..... + Tau(n).xn + ......

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

    (n-1).Tau(n) = SUM1<=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) = n4.s(n) - 24 SUM0<k<n k2.(35.k2-52.k.n+18.n2).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 1013  = 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 1010 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  a0 = a  &  a1 = b , we perform the successive Euclidean divisions

    a0 = a1 q1 + a2                     and we                u0 = 1  ,  u1 = 0  ,  uk+1 = uk-1 - qk uk
    a1 = a2 q2 + a3                      define                 v0 = 0  ,  v1 = 1  ,  vk+1 = vk-1 - qk vk
    ......................

  ak-1 = ak qk + ak+1
   .......................

  an-1 = an qn + 0                      it yields                   an = gcd(a,b)        u = c un/ an         v = c vn/ an
 
 

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