hp41programs

Polynomials finite Fields

Polynomial Factorization over Finite Fields for the HP-41


Overview
 

 1°)  Euclidean Division
 2°)  Polynomial modulo p
 3°)  Polynomial modulo f(x) and modulo p
 4°)  Power of a Polynomial modulo f(x) and modulo p
 5°)  Greatest Common Divisor in Z/pZ [X]
 6°)  Coefficients of the First Derivative of a Polynomial
 7°)  Polynomial >>> Squarefree Polynomial
 8°)  Squarefree Polynomial >>> Products of Irreductible Factors of Degree k = 1 , 2 , ............
 9°)  Factorization in Z/pZ [X] ( p > 2 ) , Cantor-Zassenhaus Algorithm
10°) Irreducible Polynomials
 

-In the following programs, p is a positive prime,  and in paragraph 9, p is an odd prime ( p > 2 )
-In most cases, p must be smaller than 10^5

-The coefficients of a polynomial  a(x) = an xn + an-1 xn-1 + ... + a1 x + a0  in Z/pZ[X]    are to be stored
  [ an an-1 ............ a0 ]  into consecutive registers  [ Rbb ............ Ree ] and a(x) is identified by the control number  bbb.eee
 

1°) Euclidean Division
 

-Given  a(x) = an xn + an-1 xn-1 + ... + a1 x + a0  &  b(x) = bm xm + bm-1 xm-1 + ... + b1 x + b0   ( deg a = n , deg b = m )
 "DIVM" finds the 2 polynomials  q(x) and r(x) such that

   a = b.q + r   deg r < deg b
 

Data Registers:   R00 thru R09:  temp

                        and the registers that contain the coefficients of the 2 polynomials            ( which are to be initialized before executing "DIVM" )

Flags: /
Subroutine:  "UV"  ( cf "GCD , LCM and Bezout's Identity for the HP-41" )
 
 

01  LBL "DIVM"
02  STO 09
03  RDN
04  STO 06
05  FRC
06  X<>Y
07  STO 05
08  FRC
09  X<>Y
10  -
11   E3
12  *
13  RCL 05
14  INT
15  STO 07
16  -
17  RCL 06
18  INT
19  +
20  1
21  +
22  STO 08
23  X>0?
24  GTO 00
25  RCL 05
26  0
27  RTN
28  LBL 00
29  RCL IND 06
30  RCL 09
31  1
32  XEQ "UV"
33  STO 04
34  LBL 01
35  RCL 04
36  RCL IND 05
37  *
38  RCL 09
39  MOD
40  STO 01
41  STO IND 05
42  RCL 05
43  STO 02
44  RCL 06
45  STO 03
46  SIGN
47  ST+ 02
48  ISG 03
49  X=0?
50  GTO 00
51  LBL 02
52  RCL IND 02
53  RCL IND 03
54  RCL 01
55  *
56  -
57  RCL 09
58  MOD 
59  STO IND 02
60  ISG 02
61  CLX
62  ISG 03
63  GTO 02
64  LBL 00
65  ISG 05
66  CLX
67  DSE 08
68  GTO 01
69  RCL 05
70  RCL 07
71  RCL 05       
72  INT
73  1
74  -
75   E3
76  /
77  +
78  END

 
  ( 108 bytes )
 
 

      STACK        INPUTS      OUTPUTS
           Y      bbb.eee(a)      bbb.eee(r)
           X      bbb.eee(b)   bbb.eee(q) or 0

  X-OUTPUT = 0  when  deg a < deg b   All  bbb's > 009

Example:   in  Z/41Z [X]   a(x) = 3 x5 + 24 x4 - 49 x3 + 16 x2 - 7 x + 28  &  b(x) =  2 x3 + 18 x2 - 9 x - 37

  Store for instance  [ 3  24  -49  16  -7  28 ]  into  [ R10  R11  R12  R13  R14  R15 ]    control number = 10.015
                     and         [ 2  18  -9  -37 ]         into       [ R16  R17  R18  R19 ]                control number = 16.019

  10.015  ENTER^
  16.019  ENTER^
      41     XEQ "DIVM"  >>>>  10.012                          --- Execution time = 10 seconds ---
                                       RDN   13.015

    we find  [ R10  R11  R12 ] = [ 22  19  6 ]     So,   q(x) = 22 x2 + 19 x + 6
          and  [ R13  R14  R15 ] = [ 32  12  4 ]             r(x) =  32 x2 + 12 x + 4

-Thus,  a(x) = b(x).q(x) + r(x)  ( for all x )  modulo 41  -  the result would be wrong in R[X]

Notes:

-The registers that contain c(x) are undisturbed
-If  b(x) is a constant,  r = 0 and  bbb.eee(r)  verifies  bbb = eee + 1 ( for instance 20.019 )
 

2°) Polynomial Modulo p
 

-After performing a product of 2 polynomials ( or another operation ), it may be useful to calculate all the coefficients modulo p
-If leading coefficients = 0 mod p , they are deleted by "DEL" which is a part of the routine listed in the next paragraph.
 

Data Registers:  Only those of the polynomial
Flags: /
Subroutine:  "DEL"  cf next paragraph
 
 

 01  LBL "PMOD"
 02  SIGN
 03  X<>Y
 04  ENTER^
 05  LBL 01
 06  RCL IND Y
 07  LASTX
 08  MOD
 09  STO IND Z
 10  RDN
 11  ISG Y
 12  GTO 01
 13  GTO "DEL"
 14  END

 
   ( 31 bytes )
 
 

      STACK        INPUTS     OUTPUTS
           Y       bbb.eee(f)            /
           X             p      bbb.eee(f)

  The output will be different from the input if the leading coefficient(s) = 0 mod m  except if f = 0

Example:     f(x) = 42 x3 + 16 x2 - 7 x + 22   in   Z/7Z [X]

  If you've stored  [ 42  16  -7  22 ]   in registers  [ R01  R02  R03  R04 ]

  1.004   ENTER^
      7      XEQ "PMOD"   >>>>  2.004

   [ R02  R03  R04 ]  =  [ 2  0  1 ]    So,  f(x) = 2 x2 + 1
 

Note:

-Here, p is not necessarily a prime.
 

3°) Polynomial Modulo f(x) & Modulo p
 

-In paragraphs 8-9, we have to perform relatively complex calculations that are greatly simplified if we do them in (Z/pZ [X])/ f
-The following routine actually computes the remainder r(x) in the division of a polynomial a(x) by another polynomial f(x) - all modulo p.
 

Data Registers:   R00 thru R09:  temp

                        and the registers that contain the coefficients of the 2 polynomials ( which are to be initialized before executing "PMODF" )

Flags: /
Subroutine:  "DIVM"  ( cf paragraph 1 )

-If you don't have an HP-41CX, replace line 09 by  0  STO IND Y  RDN
 
 

 01  LBL "PMODF"
 02  XEQ "DIVM"
 03  SIGN
 04  -
 05  ISG X
 06  GTO 01
 07  LASTX
 08  -
 09  CLRGX  
 10  LBL "DEL"
 11  LBL 01
 12  RCL IND X
 13  X#0?
 14  GTO 02
 15  X<>Y
 16  ISG X
 17  GTO 01
 18  ENTER^
 19  DSE Y
 20  LBL 02
 21  X<>Y
 22  END

 
  ( 51 bytes )
 
 

      STACK        INPUTS      OUTPUTS
           Z      bbb.eee(a)             /
           Y      bbb.eee(f)             /
           X            p      bbb.eee(r)

       All  bbb's > 009

Example:         in  Z/41Z [X]   a(x) = 3 x5 + 24 x4 - 49 x3 + 16 x2 - 7 x + 28   ,   f(x) = 2 x2 + x + 5

     If you have stored   [ 3  24  -49  16  -7  28 ]  into  [ R10  R11  R12  R13  R14  R15 ]
                 and                      [ 2  1  5 ]                into              [ R21  R22  R23 ]

  10.015  ENTER^
  21.023  ENTER^
      41     XEQ "PMODF"   >>>>   14.015   --- Execution time = 10s ---

  [ R14  R15 ]  =  [ 40  26 ]   So,  r(x) = 40 x + 26 = -x + 26  ( mod 41 )
 

4°) Power of a Polynomial Modulo f(x) & Modulo p
 

-Given a polynomial  a(x) , a positive integer k , a polynomial f(x) and a prime p,  "POWMF" returns  ak mod ( f , p )
-We assume that deg a < deg f  ( otherwise, execute "PMODF" first )
 

Data Registers:   R00 thru R14:  temp                          ( Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "POWMF" )

                           •  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0                           All bbb's > 014
                           •  Rbb' = fm  •  Rbb'+1 = fm-1  ........................  •  Ree' = f0

                        Registers  Rbb ...............  Rbb-2+4 deg f     are also used  and all these blocks of registers cannot overlap ( choose for instance  ee' < bb )
Flags: /
Subroutines:  "PRO"         ( cf "Polynomials for the HP-41" )
                         "PMODF"  ( cf paragraph 3 )
                         "LCO"        ( cf "Miscellaneous short routines for the HP-41" )
 
 

01  LBL "POWMF"
02  STO 09
03  RDN
04  STO 10
05  RDN
06  STO 11
07  X<>Y
08  STO 12
09  INT
10  RCL 10
11  FRC
12   E3
13  *
14  RCL 10
15  INT
16  -
17  STO Z
18  +
19  .1
20  %
21  +
22  STO 13
23  +
24  INT
25  STO 14
26  SIGN
27  STO IND 13   
28  LBL 01
29  RCL 11
30  2
31  /
32  STO 11
33  FRC
34  ST- 11
35  X=0?
36  GTO 02
37  RCL 12
38  RCL 13
39  RCL 14
40  XEQ "PRO"    
41  RCL 10
42  RCL 09
43  XEQ "PMODF"
44  RCL 13
45  INT
46  XEQ "LCO"
47  STO 13
48  LBL 02
49  RCL 11
50  X=0?
51  GTO 03
52  RCL 12
53  RCL 12
54  RCL 14
55  XEQ "PRO"
56  RCL 10
57  RCL 09
58  XEQ "PMODF"
59  RCL 12
60  INT
61  XEQ "LCO"
62  STO 12
63  GTO 01
64  LBL 03
65  RCL 13
66  END
 

 
   ( 112 bytes )
 
 

      STACK        INPUTS      OUTPUTS
           T      bbb.eee(a)             /
           Z            k             /
           Y      bbb.eee(f)             /
           X            p      bbb.eee(r)

       All  bbb's > 014

Example:       in  Z/41Z [X]     a(x) =   6 x4 + 5 x3 + 4 x2 + x + 8     k = 757       f(x) = 3 x5 + 2 x4 - 4 x3 + 16 x2 - 7 x + 28

  If   [ 3  2  -4  16  -7  28 ]  has been stored in  [ R15  R16  R17  R18  R19  R20 ]   control number = 15.020
and       [ 6  5  4  1  8 ]        has been stored in     [ R21  R22  R23  R24  R25 ]         control number = 21.025

   21.025   ENTER^
     757      ENTER^
   15.020   ENTER^
      41       XEQ "POWMF"   >>>>   26.030  = R13                     --- Execution time = 7mn08s ---

  [ R26  R27  R28  R29  R30 ]  =  [ 23  30  38  9  16 ]

-So,  a757 mod ( f , 41 ) = 23 x4 + 30 x3 + 38 x2 + 9 x + 16

Notes:

-SIZE 040 was enough for this example.
-The registers that contain the coefficients of f(x) are unchanged.
-Actually, we are working in ( Z/pZ [X] ) / < f >
-If deg a >= deg f , execute "PMODF" first
 

5°)  Greatest Common Divisor in Z/pZ [X]
 

-"PGCD" uses the Euclidean algorithm to find the greatest common divisor of 2 polynomials

  a(x) = an xn + an-1 xn-1 + ... + a1 x + a0  &  b(x) = bm xm + bm-1 xm-1 + ... + b1 x + b0         ( deg a = n , deg b = m )
 

Data Registers:   R00 thru R11:  temp                          ( Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "PGCD" )

                           •  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0                           All bbb's > 011
                           •  Rbb' = bm  •  Rbb'+1 = bm-1  .....................  •  Ree' = b0

Flags: /
Subroutines:     "DIVM"     ( cf paragraph 1 )
                            "PMOD"    ( cf paragraph 2 )
                            "DEL"        ( cf paragraph 3 )
                            "UV"          ( cf "GCD , LCM and Bezout's Identity for the HP-41" )

-Lines 37 to 49 are useful to produce a unitary polynomial ( leading coefficient = 1 )
 otherwise, they may be replaced by  RCL 11 and the routine will be faster
 
 

01  LBL "PGCD"
02  STO 09
03  RDN
04  STO 11
05  X<> T
06  XEQ "PMOD"
07  X<> 11
08  RCL 09
09  XEQ "PMOD"
10  RCL IND X
11  X=0?
12  GTO 02
13  X<>Y
14  X<> 11
15  RCL IND X  
16  X=0?
17  GTO 02
18  CLX
19  RCL 11
20  LBL 01
21  STO 11
22  RCL 09
23  XEQ "DIVM"
24  SIGN
25  -
26  ISG X
27  X=0?
28  GTO 02
29  XEQ "DEL"
30  RCL IND X  
31  X=0?
32  GTO 02
33  RCL 11
34  RCL Z
35  GTO 01
36  LBL 02
37  RCL IND 11  
38  RCL 09   
39  1
40  XEQ "UV"
41  RCL 11
42  X<>Y
43  LBL 03
44  ST* IND Y
45  ISG Y
46  GTO 03
47  RCL 11
48  RCL 09
49  XEQ "PMOD"
50  END

 
  ( 103 bytes )
 
 

         STACK          INPUTS       OUTPUTS
             Z        bbb.eee(a)              /
             Y        bbb.eee(b)              /
             X              p    bbb.eee(gcd(a,b))

       All  bbb's > 011

Example:      a(x) =  x7 + 2 x6 + x5 + 3 x4 + x2 + 3 x + 2          b(x) = 4 x8 + 4 x6 + x4 + x3 + 2 x2 + x + 6      in  Z/7Z [X]

  Store for instance  [ 1  2  1  3  0  1  3  2 ]  into  [ R12  R13  R14  R15  R16  R17  R18  R19 ]
              and         [ 4  0  4  0  1  1  2  1  6 ] into  [ R20  R21  R22  R23  R24  R25  R26  R27  R28 ]

    12.019   ENTER^
    20.028   ENTER^
        7        XEQ "PGCD"  >>>>   22.028     --- Execution time = 36s ---

  [ R22  R23  R24  R25  R26  R27  R28 ]  =  [ 1  5  2  2  6  5  4 ]

-Thus,   GCD(a,b) =  x6 + 5 x5 + 2 x4 + 2 x3 + 6 x2 + 5 x + 4 = g(x)

-You can check that  a(x)/g(x) = x + 4  and  b(x)/g(x) = 4 x2 + x + 5     modulo 7
 

Note:

-If you have replaced lines 37 to 49 by  RCL 11, you'll get  GCD(a,b) = 2 x6 + 3 x5 + 4 x4 + 4 x3 + 5 x2 + 3 x + 1  in 31 seconds.
-All the coefficients have been mutiplied by 2 ( modulo 7 )
 

6°) Coefficients of the first Derivative of a Polynomial
 

-"dPLC" replaces the coefficients ak of a polynomial  a(x) = an xn + an-1 xn-1 + ... + a1 x + a0  by the coefficients of a'(x)  i-e  k ak  ( except a0 )
-It works for polynomials over Z[X] or R[X] too, but not over C[X].
 

Data Registers:  Only those of the polynomial
Flags: /
Subroutines: /
 
 

 01  LBL "dPLC"
 02  ENTER^
 03  INT
 04  LASTX
 05  FRC
 06   E3
 07  *
 08  X<>Y
 09  -
 10  X=0?
 11  GTO 00
 12  X<>Y
 13   E-3
 14  -
 15  X<>Y
 16  LBL 00
 17  X<>Y
 18  STO Z
 19  X<>Y
 20  LBL 01
 21  ST* IND Z
 22  ISG Z
 23  ""      
 24  DSE X
 25  GTO 01
 26  X<>Y
 27  END

 
  ( 45 bytes )
 
 

      STACK        INPUT      OUTPUT
           X      bbb.eee(a)      bbb.eee(a') 

Example:         a(x) =  x5 + 3 x4 + x2 + 3 x + 2

  If  [ 1  3  0  1  3  2 ]  is stored in  [ R01  R02  R03  R04  R05  R06 ]

  1.006   XEQ "dPLC"  >>>>  1.005

  and  [ R01  R02  R03  R04  R05 ]  =  [ 5  12  0  2  3 ]    So,  a'(x) = 5 x4 + 12 x3 + 2 x + 3
 

7°) Polynomial >>>> Squarefree Polynomial
 

-Given a polynomial a(x) over Z/pZ[X] , the following routine returns the polynomial  b(x) = a(x) / GCD(a(x),a'(x))
-Usually, these polynomials have the same irreducible factors fi(x), but  (fi(x))2 is never a divisor of b(x)

-However, some problems should be kept in mind:  a'(x) may be equal to zero even if a(x) is not constant
-For example, in Z/7Z[X]   a(x) = 2 x14 + 3 x7 + 4  has a derivative  a'(x) = 0   since  7 = 0  ( mod 7 )
-In such cases, we always have:  a(x) = [ g(x) ]p   Here,  a(x) = ( 2 x2 + 3 x + 4 )7   ( Fermat's Theorem )
-"SQFREE" will display a  DATA ERROR message line 25 when a'(x) = 0

-In similar cases, we can also lose factors. For instance,  if  a(x) = ( x7 + 1 )( x + 2 ) = x8 + 2 x7 + x + 2    ( mod 7 )

  a'(x) = 8 x7 + 1  and  b(x) = x + 2    thus, we've lost the factor ( x7 + 1 )

-We will, however, retrieve it by dividing a(x) by b(x)

-These pitfalls cannot occur if  p > deg a or if deg a - deg b < p
 

Data Registers:   R00 thru R11:  temp                                    ( Rbb thru Ree are to be initialized before executing "SQFREE" )

                           •  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0                           bbb > 011

                              Registers  Ree+1 thru Ree+2+2.deg a  are also used
Flags: /
Subroutines:     "DIVM"     ( cf paragraph 1 )
                            "PMOD"    ( cf paragraph 2 )
                            "dPLC"      ( cf paragraph 6 )
                            "LCO"       ( cf "Miscellaneous short routines for the HP-41" )
 
 

01  LBL "SQFREE"
02  STO 09
03  X<>Y
04  STO 10
05  ENTER^
06  FRC
07   E3
08  *
09  1
10  +
11  XEQ "LCO"   
12  STO 11
13  ENTER^
14  FRC
15   E3
16  *
17  1
18  +
19  XEQ "LCO"
20  XEQ "dPLC" 
21  RCL 09
22  XEQ "PMOD"
23  RCL IND X
24  X=0?
25  LN
26  CLX
27  RCL 11
28  X<>Y
29  RCL 09
30  XEQ "PGCD"
31  RCL 10
32  X<>Y
33  RCL 09
34  XEQ "DIVM"
35  END

 
  ( 77 bytes / SIZE eee+3+2.deg a )
 
 

         STACK          INPUTS       OUTPUTS
             Y        bbb.eee(a)              /
             X               p       bbb.eee(b)

         bbb > 011

Example:      a(x) =  x7 + 8 x6 + 27 x5 + 52 x4 + 65 x3 + 54 x2 + 28 x + 8   in  Z/97Z[X]

  If  [ 1  8  27  52  65  54  28  8 ]  is stored in  [ R12  R13  R14  R15  R16  R17  R18  R19 ]

  12.019  ENTER^
      97     XEQ "SQFREE  >>>>   12.015                                           --- Execution time = 69s ---

  [ R12  R13  R14  R15 ]  =  [ 1  3  3  2 ]    So  b(x) =  x3 + 3 x2 + 3 x + 2

Notes:

-Actually,  a(x) = ( x+ 2 )3( x2 + x + 1 )2    and    b(x) = ( x+ 2 )( x2 + x + 1 )
-If you have replaced lines 37 to 40 by RCL 11 in the routine "PGCD" , you get  b(x) = 11 x3 + 33 x2 + 33 x + 22    and execution time = 65 seconds
 

8°)  Squarefree Polynomial >>> Products of Irreductible Factors of Degree k = 1 , 2 , ............
 

-"PRFDK" takes a squarefree polynomial  f(x) = an xn + an-1 xn-1 + ... + a1 x + a0  in Z/pZ[X]
  and returns gradually  g1(x) , g2(x) , ..... , gk(x) , ......    where gk(x)  is a product of irreducible factors of degree k.

-Let  f0 = f ,  the HP-41 computes         fk = fk-1/gk    where  gk =  gcd(fk-1,xp^k-x)      for  k = 1 , 2 , ...     until  fk is a constant

-If  deg fk = k , then  fk  is irreducible over Z/pZ[X]
-Otherwise, we'll have to go further ( cf next paragraph )
 

Data Registers:   R00 thru R18:  temp                                    ( Rbb thru Ree are to be initialized before executing "PRFDK" )

                           •  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0                           bbb > 018

                              Registers  Ree+1 thru Ree+6.deg f  are also used for temporary data storage
Flags: /
Subroutines:     "DIVM"     ( cf paragraph 1 )
                            "POWMF" ( cf paragraph 4 )
                            "PGCD"     ( cf paragraph 5 )
                            "LCO"       ( cf "Miscellaneous short routines for the HP-41" )

-If you don't have an HP-41CX, replace line 37  by  XEQ "LCL"   ( cf "Miscellaneous short routines for the HP-41" )
 
 

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

 
   ( 177 bytes / SIZE eee+1+6.deg f )
 
 

         STACK          INPUTS       OUTPUTS1      OUTPUTS2      ......................      OUTPUTSk
             Z               /       bbb.eee(f1)       bbb.eee(f2)      ......................       bbb.eee(fk)
             Y        bbb.eee(f)       bbb.eee(g1)       bbb.eee(g2)      ......................       bbb.eee(gk)
             X               p             1               2              ...             k

       bbb > 018      The HP-41 emits a BEEP for the last output and a TONE 9 for the other ones
                                       The last bbb.eee(fk)  is sometimes replaced by a single "1"

Example:     f(x) = x6 + 2 x5 + x4 + 3 x3 + 6 x2 + 7 x + 4   in  Z/13Z[X]       ( use "SQFREE" to check that f(x) is squarefree over Z/13Z[X] )

   If  [ 1  2  1  3  6  7  4 ]   is stored in   [ R19  R20  R21  R22  R23  R24  R25 ]

   SIZE 062  at least

   19.025  ENTER^
       13     XEQ "PRFDK"  >>>>       1                 --- Execution time = 2mn32s ---
                                          RDN   41.044
                                          RDN   26.029

        [ R41  R42  R43  R44 ]  =  [ 1  2  2  1 ]         g1(x) = x3 + 2 x2 + 2 x + 1       is a product of factors of degree 1
        [ R26  R27  R28  R29 ]  =  [ 1  0  12  4 ]       f1(x) = x3 + 12 x + 4

                     R/S               >>>>        2                 --- Execution time = 4mn09s ---
                                         RDN     44.044
                                         RDN     26.029

                                 [ R44 ]  =  [ 1 ]                      g2(x) = 1                                  There is no factor of degree 2
       [ R26  R27  R28  R29 ]  =  [ 1  0  12  4 ]         f2(x) = x3 + 12 x + 4

                     R/S              >>>>        3                  --- Execution time = 1s ---        the next step is skipped since deg f2 = 3
                                         RDN     26.029
                                         RDN        1

       [ R26  R27  R28  R29 ]  =  [ 1  0  12  4 ]         g3(x) = x3 + 12 x + 4                is "the" factor of degree 3
                                                                               f3(x) = 1

-So,  f(x) = ( x3 + 2 x2 + 2 x + 1 ).( x3 + 12 x + 4 )     and we still have to factor  ( x3 + 2 x2 + 2 x + 1 )

Notes:

-Incidentally, we have proved that  x3 + 12 x + 4  is an irreducible polynomial over Z/13Z[X]
-In fact, we could have stopped after the first output since f1(x) cannot be a multiple of a product of polynomial(s) of degree 2

-If you have replaced lines 37 to 40 by RCL 11 in the routine "PGCD" , you get

     g1(x) = 7 x3 + x2 + x + 7     f1(x) = 2 x3 + 11 x + 8          whence  f(x) = 10.( 7 x3 + x2 + x + 7 ).( 8 x3 + 5 x + 6 )
     g2(x) = 10                           f2(x) = 8 x3 + 5 x + 6
     g3(x) = 8 x3 + 5 x + 6         f3(x) = 1                                though the 2 factorizations may seem different at first sight, they are equivalent modulo 13.

-The coefficients of f(x) are still in registers R19 .... R25
 

9°)  Factorization in Z/pZ [X] ( p > 2 ) , Cantor-Zassenhaus Algorithm
 

-"CANZAS" takes a squarefree polynomial  f(x) = an xn + an-1 xn-1 + ... + a1 x + a0  in Z/pZ[X]  whose irreducible factors have the same degree k
  ( "SQFREE" & "PRFDK" are to be used first )

-Cantor-Zassenhaus method is a probabilistic algorithm:

-We take a random polynomial a(x) of degree n-1

       •  If gcd( f , a ) is not a constant, then we have found a non-trivial factor g
       •  Otherwise, we compute b = 1 + a(p^k-1)/2   -If  deg [ gcd ( f , b ) ] # 1 and # n , then we have found a non-trivial factor g.

-When a non-trivial factor g is found, the routine continues after replacing f by ( f / g )
-Otherwise, we try another random polynomial a(x)

-If one or several g have a degree > k , we must execute "CANZAS" again for these factors.

-Another choice for b(x) must be made if p = 2, but the following program only works for p > 2.
 
 

Data Registers:   R00 thru R19:  temp                                    ( Rbb thru Ree are to be initialized before executing "CANZAS" )

                           •  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0                           bbb > 020

                              Registers  Ree+1 thru Ree+5.deg f  are also used for temporary data storage
Flags: /
Subroutines:     "DIVM"     ( cf paragraph 1 )
                            "POWMF" ( cf paragraph 4 )
                            "PGCD"     ( cf paragraph 5 )
                            "LCO"       ( cf "Miscellaneous short routines for the HP-41" )

-Line 07,  the control number of f(x) is used as a random seed.
-Line 37 is a three-byte  GTO 04
-Lines 64-65 ,  the random number generator    R-D   FRC   is used to build the random polynomial a(x)
-  p^k  must not exceed  10^10  So you could replace line 93 by   E10   X<Y?   ACOS  SIGN  -
  to produce a DATA ERROR in that case.
-Lines 108-118 are three-byte  GTO 01
 
 

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

 
   ( 205 bytes / SIZE 1+eee+5.deg f )
 
 

         STACK         INPUTS        OUTPUT1       OUTPUT2      ......................       OUTPUTk
             Z        bbb.eee(f)               /               /      ......................                /
             Y               k               /               /      ......................                /
             X            p > 2         bbb.eee1         bbb.eee2              ...         bbb.eeek

       bbb > 018       The HP-41 emits a BEEP for the last output and a TONE 9 for the other ones
 

Example1:      In Z/37Z[X] ,  f(x) = 4 x6 - 4 x5 + 37 x4 + 2 x3 + 118 x2 + 26 x + 105    is squarefree and it is a product of irreductible polynomials of degree 2
                                                                                                                                                   ( use "SQFREE" & "PRFDK" to check )

   If you have stored   [ 4  -4  37  2  118  26  105 ]   into  [ R21  R22  R23  R24  R25  R26  R27 ]

   SIZE 058

   21.027  ENTER^
        2      ENTER^
       37     XEQ "CANZAS"   >>>>   ( TONE 9 )     44.046          --- Execution time = 10mn08s ---

   [ R44  R45  R46 ]  =  [ 1  35  7 ]     So, we have a non-trivial factor    x2 + 35 x + 7 = x2 - 2 x + 7  (mod 37 )

                          R/S              >>>>   ( TONE 9 )     35.037          --- Execution time = 5mn58s ---

   [ R35  R36  R37 ]  =  [ 1  19  20 ]    A second factor is                       x2 + 19 x + 20

                         R/S               >>>>     ( BEEP )        21.023          --- Execution time = 9s ---

   [ R21  R22  R23 ]  =  [ 4  2  10 ]      The last factor is                        4 x2 + 2 x + 10
 

-Thus,   f(x) =  ( x2 - 2 x + 7 ) (  x2 + 19 x + 20 ) (  4 x2 + 2 x + 10 )   which may be re-written:
            f(x) =  ( x2 - 2 x + 7 ) ( 2 x2 + x + 3 ) ( 2 x2 + x + 5 )               modulo 37

-It is also the factorization of f(x) in Z[X]
 

Example2:      In  Z/13Z [X]    f(x) = 2 x4 + 4 x3 + x2 + x + 5     is squarefree and it is a product of irreducible polynomials of degree 1
                                                                                                                       ( use "SQFREE" & "PRFDK" to check )

-If   [ 2  4  1  1  5 ]  =  [ R21  R22  R23  R24  R25 ]

       21.025   ENTER^
            1       ENTER^
           13      XEQ "CANZAS"   >>>>  ( TONE 9 )    28.030          --- Execution time = 2mn05s ---

   [ R28  R29  R30 ]  =  [ 1  12  1 ]     A factor is   g(x) =  x2 + 12 x + 1

                               R/S             >>>>    ( TONE 9 )    29.030          --- Execution time = 26s ---

   [ R29  R30 ]  =  [ 1  4 ]      x + 3    is another factor

                               R/S             >>>>     ( BEEP )       21.022          --- Execution time = 6s ---

  [ R21  R22 ]  =  [ 2  11 ]     2 x + 11 = 2 ( x + 12 )   mod 13    is also factor

-But  deg g = 2 > k = 1   so we execute "CANZAS" again for g(x) = x2 + 12 x + 1     After storing  [ 1  12  1 ]  into  [ R21  R22  R23 ]

      21.023   ENTER^
           1       ENTER^
          13      XEQ "CANZAS"   >>>>  ( TONE 9 )    29.030          --- Execution time = 2mn14s ---

  [ R29  R30 ]  =  [ 1  9 ]      x + 9     is factor

                            R/S                >>>>    ( BEEP )       21.022          --- Execution time = 6s ---

  [ R21  R22 ]  =  [ 1  3 ]      x + 3     is the last factor

-Thus, we have found    f(x) = ( x2 + 12 x + 1 ) ( x + 4 ) 2 ( x + 12 )
         and finally                   = 2 ( x + 9 ) ( x + 3 ) ( x + 4 ) ( x + 12 )
 

Notes:

-Since Cantor-Zassenhaus algorithm is probabilistic, it may be - relatively - fast or extremely slow if we are unlucky !
-One could perhaps write a program that handles all the steps ( "SQFREE" "PRFDK" "CANZAS" )
  but it would be awkward to use all the control numbers without any register conflict.
-Moreover, a lot of data registers would be required, but if you are in the mood ...

Exercise:    Factorize   f(x) =  x15 - 5 x14 + 4 x13 - 2 x12 - x11 - 5 x10 - x9 - x8 + 2 x7 + 6 x6 + x5 - 6 x3 + 4 x2 + 6 x - 4   over  Z/13Z

Answer:      f(x) = ( x + 1 )( x + 6 )3( x2 + x + 5 )2( x3 + x - 6 )( x4 - 3 x2 + 5 x + 1 )

  SIZE 099  for this example.

  SQFREE needs 3mn41s to find a squarefree g(x) with the same irreducible factors as f(x)     ( deg g = 11 )

  PRFDK   needs 3mn58s to find the product of irreductible polynomials of degree 1
                          11mn20s to find the product of irreductible polynomials of degree 2 ( actually one polynomial )
                          10mn51s to find the product of irreductible polynomials of degree 3 ( actually one polynomial )
             and just 1 second to find the product of irreductible polynomials of degree 4 ( actually one polynomial )

  CANZAS needs 72 seconds to find ( x+ 6 ) and 6 seconds to return ( x + 1 )

-After g(x) is factorized, calculate f(x) / g(x) = h(x)  ( deg h = 4 )
   test if ( x + 1 ) is a divisor of h:  it doesn't work!
   test if ( x + 6 ) is a divisor of h:  it works twice, and the last quotient is precisely the last unknown factor i-e ( x2 + x + 5 )
 

10°) Irreducible Polynomials
 

-This routine takes a polynomial f of degree n  and a positive prime p   and returns:

   1 if  f is irreducible in Z/pZ[X]
   0 if  f has a non-trivial factor.
 

Data Registers:   R00 thru R19:  temp                                     ( Rbb thru Ree are to be initialized before executing "IRRM?" )

                           •  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0                           bbb > 019

                              Registers  Ree+1 thru Ree+6.deg f  are also used for temporary data storage

Flag:  F26 is cleared
Subroutines:     "SQFREE"  ( §7 )  &  "PRFDK"   ( §8 )   after adding  LBL "NEXT" after line 84 RTN  of this subroutine

-Line 06 produces a  DATA ERROR message if  an = 0  mod p
-Line 14:  if f is not squarefree, it's not irreducible
-Line 46:  the HP-41 computes gcd (f,xp^k-x) while k <= int(n/2)
 
 

01  LBL "IRRM?"
02  RCL IND Y
03  RCL Y
04  MOD
05  X=0?
06  LN    
07  X<> Z
08  STO 12
09  X<>Y
10  XEQ "SQFREE"
11  RCL 12
12  X=Y?
13  GTO 00           
14  CLX   
15  RTN
16  LBL 00
17  FRC
18   E3
19  *
20  RCL 12
21  INT
22  -
23  2
24  X>Y?
25  GTO 02            
26  /
27  INT
28  STO 19
29  RCL 12
30  RCL 09
31  CF 26
32  XEQ "PRFDK" 
33  GTO 00
34  LBL 01
35  XEQ "NEXT"  
36  LBL 00
37  ISG Y
38  X=0?
39  GTO 00
40  SIGN  
41  -
42  0
43  RTN
44  LBL 00
45  RCL 19             
46  X>Y?    
47  GTO 01
48  LBL 02
49  SIGN
50  END

 
    ( 91 bytes / SIZE eee+1+6.deg f )
 
 

      STACK        INPUTS      OUTPUTS
           Y        bbb.eee    /   or  bbb.eee'
           X             p       1  or   0

  If f is irreducible over Z/pZ[X] , X-output = 1                                                                                                  bbb > 019
  If f is not irreducible X-output = 0 and Y-output is the control number of a divisor of f
  This divisor is non-trivial except if f = product of irreducible polynomials of degree k where k = R16

Example1:        f(x) = x7 + x + 1  over  Z/2Z

-With, for instance  [ 1  0  0  0  0  0  1  1 ]  =  [ R20  R21  R22  R23  R24  R25  R26  R27 ]

   20.027  ENTER^
       2       XEQ "IRRM?"   >>>>   1                          --- Execution time = 4mn05s ---

-Thus,  f is irreducible in Z/2Z[X] and in Z[X] too !

Example2:       f(x) = x5 + 2 x4 + 4 x3 + 7 x2 + x + 3  over  Z/11Z

-If   [ 1  2  4  7  1  3 ]  =  [ R20  R21  R22  R23  R24  R25 ]

   20.025  ENTER^
      11      XEQ "IRRM?"   >>>>      0                       --- Execution time = 6mn35s ---
                                          RDN   39.041

-So, f is not irreducible in Z/11Z[X]  and  [ R39  R40  R41 ]  =  [ 1  6  1 ]   whence  x2 + 6 x + 1  is a non-trivial divisor of f.
 ( moreover, since R16 = 2 ,  x2 + 6 x + 1 is irreducible over Z/11Z )
 

Example3:       f(x) = x5 + 2 x4 + 4 x3 + 7 x2 + x + 3  over  Z/17Z

-If   [ 1  2  4  7  1  3 ]  =  [ R20  R21  R22  R23  R24  R25 ]  again

   20.025  ENTER^
      17      XEQ "IRRM?"   >>>>    1                       --- Execution time = 7mn01s ---

-So,   f is irreducible in Z/17Z[X] and in Z[X] too.