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