Polynomials Z[X]

# Polynomials in Z[X] for the HP-41

Overview

1°)  Irreducible Polynomials

a)  Eisenstein Criterion
b)  Eisenstein Criterion & Translation
c)  Via Polynomials in Z/pZ[X] ( where p is a prime )
d)  Cohn's Criterion
e)  A better Criterion
f)  2 other Criteria

2°)  Factorization

a)  The Brute Force
b)  Via Factorization over finite Fields
c)  Via Bairstow's Method

c-1) Program#1
c-2) Program#2

3°)  Cyclotomic Polynomials

a)   Program#1
b)   Program#2

4°)  Polynomial Factorization over Gaussian Integers

The following routines deal with polynomials whose coefficients are integers ( complex numbers in §4 )
Most of these polynomials are irreducible over Z, so irreducibity criteria are not useless before trying to factor them!

1°)  Irreducible Polynomials

a) Eisenstein Criterion

-Let a polynomial  f(x) = an xn + an--1 xn-1 + ... + a1 x + a0
-If there is a prime p such that

p is a divisor of  a0 , a1 , a2 , ............ , an-1
p is not a divisor of  an ,  p2 is not a divisor of  a0

then  f(x) is irreducible over Z

-The following program returns 1 in X-register if Eisenstein criterion is satisfied and in this case, f is irreducible in Z[X]
-Otherwise, "EISN" returns 0 and we cannot conclude - except if a0 = 0.

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

•  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0          bbb > 004
Flags: /
Subroutines: /

-LBL 01   this loop calculates  GCD( a0 , a1 , ............ , an-1 )
-Lines 41 to 53 compute the prime divisors of  GCD( a0 , a1 , ............ , an-1 )
they don't use the fastest method, but it is usually acceptable if gcd is not too large.

 01  LBL "EISN" 02  STO 01 03  FRC 04   E3 05  * 06  STO 02 07  RCL 01 08  INT 09  - 10  1 11  XY 30  X#0? 31  GTO 02 32  X<>Y 33  GTO 01 34  LBL 00 35  ABS 36  STO 03 37  2 38  STO 00 39  STO 04 40  GTO 06 41  LBL 03   42  RCL 03  43  RCL 00 44  MOD 45  X=0? 46  RTN 47  RCL 04        48  LASTX 49  X=Y? 50  DSE X 51  + 52  STO 00 53  GTO 03 54  LBL 04 55  STO 03 56  RCL 00 57  / 58  FRC 59  X#0? 60  RTN 61  LASTX 62  GTO 04        63  LBL 05 64  RCL 03 65  XEQ 04 66  LBL 06 67  RCL 03 68  DSE X 69  X<0? 70  RTN 71  XEQ 03 72  RCL IND 01 73  RCL 00 74  MOD 75  X=0? 76  GTO 05 77  RCL IND 02 78  LASTX 79  X^2 80  MOD 81  X=0? 82  GTO 05 83  SIGN 84  END

( 116 bytes )

 STACK INPUT OUTPUT X bbb.eee 0 or 1

where   bbb.eee = control number of f(x)

X-output = 1    if Eisenstein criterion is satisfied,         bbb > 004
X-output = 0    otherwise.

Example:           f(x) = 10 x5 + 42 x3 + 98 x2 + 28 x + 56

-Store   [ 10  0  42  98  28  56 ]   into  R05  R06  R07  R08  R09  R10

5.010  XEQ "EISN"  >>>>   1                        --- Execution time = 4s ---

-So,  f is irreducible in Z[X]

-The prime p is in register R00  ( here, p = 7 )

b)  Eisenstein Criterion & Translation

-If Eisenstein criterion is not satisfied, we cannot conclude but if f is irreducible in Z[X]
one may sometimes find an integer "a" such that f(x+a) satisfies Eisenstein property.
( "f(x) irreducible" & "f(x+a) irreducible" are equivalent )

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

•  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0                bbb > 010            Ree+1 ..... temp
Flags: /
Subroutines:  "PCOMP"  composition of 2 polynomials ( cf "Polynomials for the HP-41" )   &   "EISN"     cf § above

 01  LBL "EISN2"   02  STO 10 03  X<>Y 04  STO 09 05  FRC 06   E3 07  * 08  1 09  ST+ Y 10  STO IND Y     11  + 12  STO Y 13  RCL 10 14  STO IND Y     15  CLX 16  .1 17  % 18  + 19  1 20  ST+ Z 21  - 22  RCL 09 23  X<> Z 24  XEQ "PCOMP" 25  XEQ "EISN"    26  RCL 10 27  SIGN 28  RCL 09 29  RCL Z 30  END

( 60 bytes )

 STACK INPUTS OUTPUTS Y bbb.eee bbb.eee X a 0 or 1 L / a

bbb.eee = control number of f(x)                                                 bbb > 010
a       = an integer

X-output = 1    if Eisenstein criterion is satisfied,
X-output = 0    otherwise

Example:      f(x) = x5 - 7 x4 + 19 x3 - 20 x2 + 5 x + 4

-Store for instance   [ 1  -7  19  -20  5  4 ]    into   R11   R12   R13   R14   R15   R16

11.016   XEQ "EISN"  >>>>   0    so,  f(x) doesn't satisfy Eisenstein criterion. Let's try "EISN2" with a = 1

11.016  ENTER^
1       XEQ "EISN2"  >>>>   0    It still doesn't work. But if we try  a = 2

RDN  2  R/S   >>>>   1                 --- Execution time = 32 seconds ---

-So f(x+2) is irreducible and f(x) too!

c)  Via Polynomials in Z/pZ[X]

-If a polynomial f(x) is irreducible over Z/pZ[X] for some prime p, then it is also irreducible over Z[X]  ( the converse is not true - unfortunately )
-For example,  f(x) = x7 + x + 1  is irreducible over Z/2Z  so f is irreducible over Z too.

-Use "IRRM?"  listed in "Polynomial Factorization over finite Fields for the HP-41" paragraph 10.

d) Cohn's Criterion

-This test concerns a poynomial  f(x) = an xn + an-1 xn-1 + ... + a1 x + a0  in Z[X]   with  ak >= 0  for all k's  ( a0 # 0 )

•  If there exists an integer m > max { ak }  such that  f(m) is a prime, then f(x) is irreducible in Z[X]

Data Registers:              R00 to R02:  are used by "PR?"              ( Registers Rbb thru Ree are to be initialized before executing "COHN" )
R03 = bbb.eee    R04 = m

•  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0 # 0             bbb > 004
Flags: /
Subroutines:   "PL"   cf "Polynomials for the HP-41"  § 1-d)
"PR?"  cf "Primes, divisor functions and Euler function for the HP-41"

-Line 24 produces "OUT OF RANGE" if  f(m) > 10^10

 01  LBL "COHN" 02  STO 03 03  0 04  LBL 01 05  RCL IND Y 06  X<0? 07  LN 08  X>Y? 09  X<>Y 10  RDN 11  ISG Y 12  GTO 01 13  X=0? 14  LN 15  STO 04 16  LBL 02          17  ISG 04 18  CLX 19  RCL 03 20  RCL 04 21  XEQ "PL"      22   E10 23  XY 26  XEQ "PR?"   27  X=0? 28  GTO 02 29  END

( 52 bytes )

 STACK INPUT OUTPUT X bbb.eee 1

where   bbb.eee = control number of f(x)            X-output = 1    if Cohn's criterion is satisfied,                  bbb > 004

Example1:       f(x) = 2 x4 + 3 x2 + 5 x + 1

With  [ 2  0  3  5  1 ]  =  [ R05  R06  R07  R08  R09 ]

5.009  XEQ "COHN"  >>>>   1         f(x) is irreducible in Z[X]               --- Execution time = 6s ---

R04 = 6  and  f(6) = 2731 is a prime

Example2:       f(x) = 2 x6 +  3 x5 + 4 x3 + 5 x2 + 2 x + 1

With  [ 2  3  0  4  5  2  1 ]  =  [ R05  R06  R07  R08  R09  R10  R11 ]

5.011  XEQ "COHN"  >>>>   1         f(x) is irreducible in Z[X]               --- Execution time = 43s ---

R04 = 8  and  f(8) = 624977 is a prime

Note:

-If no prime is found, the routine continues until f(m) exceeds 1010 and the HP-41 displays "OUT OF RANGE" and in this case, we cannot conclude.

e) A better Criterion

-Unlike Cohn's criterion, the following one also works if the polynomial has negative coefficients.
-Let  f(x) = an xn + an-1 xn-1 + ... + a1 x + a0  in Z[X]  deg f > 1    ( a0 # 0 )

m1 = max { | ai/an | ,  i = 0 , 1 , .... , n-1 }
m2 = max { | ai/an |1/(n-i) ,  i = 0 , 1 , .... , n-1 }

H = min { 1+m11/(r+1) , 2 m2 }   where  r is such that  an-1 = an-2 = ..... = an-r = 0 ,     ( r = 0 if an-1 # 0 )

•  If there exists an integer m > H  such that  f(m) is a prime or +/-1, then f(x) is irreducible in Z[X]

Data Registers:              R00 to R02:  temp                                 ( Registers Rbb thru Ree are to be initialized before executing "IRR?" )
R03 = bbb.eee    R04 = m

•  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0 # 0                 bbb > 004
Flags: /
Subroutines:   "PL"   cf "Polynomials for the HP-41"  § 1-d)
"PR?"  cf "Primes, divisor functions and Euler function for the HP-41"

-Lines 49-50-51  There is no XROOT function on the HP-41, and  1/X  Y^X  may produce a small error
for example,  343  ENTER^  3  1/X  Y^X  gives  6.999999999  instead of 7
-So, we must add a very small number to avoid a too small M which could lead to a wrong conclusion.
-Lines 61 to 66 may be replaced by the M-code routine "PR?"

 01  LBL "IRR?" 02  STO 01 03  STO 02 04  STO 03 05  SIGN 06  STO 00 07  STO 04 08  ST+ 01 09  ST+ 02 10  LBL 01 11  RCL IND 01 12  X#0? 13  GTO 00 14  SIGN 15  ST+ 04 16  ISG 01 17  GTO 01 18  LBL 00 19  CLST 20  LBL 02 21  RCL IND 02 22  RCL IND 03 23  / 24  ABS 25  STO 01 26  X>Y? 27  X<>Y 28  X<> Z 29  RCL 01 30  RCL 00 31  1/X 32  Y^X 33  XY 35  R^ 36  ISG 00 37  CLX 38  ISG 02 39  GTO 02        40  RCL 04 41  1/X 42  Y^X 43  1 44  + 45  X<>Y 46  ST+ X 47  X>Y? 48  X<>Y 49   E-5    50  +         51  INT     52  STO 04         53  LBL 03 54  RCL 03 55  ISG 04 56  CLX 57  RCL 04 58  XEQ "PL" 59  X=0? 60  RTN 61  ABS   62   E10 63  XY 66  XEQ "PR?"    67  X=0? 68  GTO 03 69  END

( 105 bytes )

 STACK INPUT OUTPUT X bbb.eee 1 or 0 or E10

where   bbb.eee = control number of f(x)                                     bbb > 004

if  X-output = 1   f is irreducible
if  X-output = 0   f is reducible: we've found an integer root
if  X-output = E10 the HP-41 displays "OUT OF RANGE" and we cannot conclude

Example:       f(x) =  2 x7 -  3 x4 + 4 x3 - 9 x2 + 6 x - 1

With  [ 2  0  0  -3  4  -9  6  -1 ]  =  [ R05  R06  R07  R08  R09  R10  R11  R12 ]

5.012  XEQ "IRR?"   >>>>   1    thus,  f is irreducible       --- Execution time = 50 s ---

R04 = m = 6  and  f(6) = 556559  is a prime

Note:

-In some cases - when all the coefficients are non-negative and an is small - Cohn's criterion may be better ( for instance f(x) = x2 + 1 )
-Both criteria may be combined in a single program.

f) 2 other Criteria

-The 2 criteria above compute f(M) , f(M+1) , .....  until we find a prime or a number larger than 1010 ( and we cannot test its primality in that case )
-This "OUT OF RANGE" may occur almost immediately if an is small and the other coefficients of the polynomial are very large.
-So another criterion is used by "IRR2?"

Let  P(f) = card { m Î Z ,  f(m) = +/-p where p is a prime or 1 }

•  If  P(f) > 2 + deg f  then  f is irreducible

-"IRR2?"  tests  f(0) , f(-1) , f(1) , f(-2) , f(2) , ............ until a prime is found or until f(m) > E10

Data Registers:              R00 to R02:  temp                                 ( Registers Rbb thru Ree are to be initialized before executing "IRR2?" )
R03 = bbb.eee    R04 = m
R05 = -n-3 , -n-2 , ..... , 0  if the criterion is satisfied.

•  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0 # 0                    bbb > 005
Flags: /
Subroutines:   "PL"   cf "Polynomials for the HP-41"  § 1-d)
"PR?"  cf "Primes, divisor functions and Euler function for the HP-41"

 01  LBL "IRR2?" 02  STO 03 03  INT 04  LASTX 05  FRC 06   E3 07  * 08  - 09  3 10  - 11  STO 05       12  CLX 13  STO 04 14  LBL 01 15  RCL 03 16  RCL 04 17  XEQ "PL"    18  ABS 19   E10 20  XY 23  X=0? 24  RTN 25  XEQ "PR?"   26  ST+ 05 27  RCL 04 28  CHS 29  X<=Y? 30  DSE X 31  LBL 01 32  STO 04 33  RCL 05       34  X#0? 35  GTO 01 36  SIGN 37  END

( 60 bytes )

 STACK INPUT OUTPUT X bbb.eee 1 or 0 or E10

where   bbb.eee = control number of f(x)                       bbb > 005

if  X-output = 1   f is irreducible
if  X-output = 0   f is reducible
if  X-output = E10 the HP-41 displays "OUT OF RANGE" and we cannot conclude

Example:       f(x) =  x5 + 101 x4 + 3 x + 1

If  [ 1  101  0  0  3  1 ]  is stored into  [ R06  R07  R08  R09  R10  R11 ]

6.011  XEQ "IRR2?"  >>>>  1           --- Execution time = 8m17s ---

So, f is irreducible

-In this example, P(f) = card { 0  -2  -4  4  -10  -14  14  -20  ...... }  > 7

and the corresponding f(m) are { 1  1579  24821  26893  909971  3342151  4417883  12959941 ... }

-"IRR2?" could be improved to store all these results.

-In fact, f(108) = 28434219589 is a prime as you can check if you have an HP-48, and f is irreducible by Cohn's criterion,
but it's difficult to use this method with an HP-41 in that case...

A similar criterion:

Let  P+(f) = card { m Î Z , f(m) > 0 and f(m) = p where p is a prime or 1 }

•  If  P+(f) > deg f , then  f is irreducible

Data Registers:              R00 to R02:  temp                                 ( Registers Rbb thru Ree are to be initialized before executing "IRR3?" )
R03 = bbb.eee    R04 = m
R05 = -n-1 , -n , ..... , 0  if the criterion is satisfied.

•  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0 # 0                      bbb > 005
Flags: /
Subroutines:   "PL"   cf "Polynomials for the HP-41"  § 1-d)
"PR?"  cf "Primes, divisor functions and Euler function for the HP-41"

 01  LBL "IRR3?" 02  STO 03 03  INT 04  LASTX 05  FRC 06   E3 07  * 08  - 09  STO 05 10  CLX 11  STO 04 12  DSE 05 13  LBL 01 14  RCL 03 15  RCL 04 16  XEQ "PL" 17  X=0? 18  RTN 19  X<0? 20  GTO 01       21   E10 22  XY 25  XEQ "PR?"  26  ST+ 05 27  LBL 01 28  RCL 04 29  CHS 30  X<=0? 31  DSE X 32  LBL 01 33  STO 04 34  RCL 05       35  X#0? 36  GTO 01 37  SIGN 38  END

( 63 bytes )

 STACK INPUT OUTPUT X bbb.eee 1 or 0 or E10

where   bbb.eee = control number of f(x)                                                 bbb > 005

if  X-output = 1   f is irreducible
if  X-output = 0   f is reducible
if  X-output = E10 the HP-41 displays "OUT OF RANGE" and we cannot conclude

Example:    With the same polynomial:     f(x) =  x5 + 101 x4 + 3 x + 1

If  [ 1  101  0  0  3  1 ]  is stored into  [ R06  R07  R08  R09  R10  R11 ]

6.011  XEQ "IRR3?"  >>>>  1           --- Execution time = 3mn55s ---   ( instead of 8m17s with "IRR2?" )

So, f is irreducible

Remarks:

-It may happen that f(m) is never a prime. For instance, f(x) = x3 + 105 x + 12  ( Eisenstein  criterion works obviously )
-In that case,  try to apply the criterion to g(x) = f(a0x)/a0 ( a0 = constant term ) since "f irreducible" & "g irreducible" are equivalent.
-The following routine may be used: it replaces the coefficients of f(x) by those of g(x).

-Lines 28 thru 37 are only useful to check that all the coefficients are smaller than E10 in magnitude

 01  LBL "PC/C" 02  ENTER^ 03  STO 00 04  FRC 05   E3 06  * 07  RCL IND X 08  FS? 00 09  ABS 10  ST/ IND Y 11  X<>Y 12  RCL 00 13  INT 14  - 15  SIGN 16  LBL 01 17  DSE L 18  FS? 30 19  GTO 02 20  CLX 21  RCL Y 22  LASTX 23  Y^X 24  ST* IND Z 25  ISG Z 26  GTO 01 27  LBL 02 28  RCL 00  29  LBL 03 30  RCL IND X 31  ABS 32   E10 33  X<=Y? 34  LN 35  X<> Z 36  ISG X 37  GTO 03       38  RCL 00 39  END

( 65 bytes )

 STACK INPUT OUTPUT X bbb.eee bbb.eee

where   bbb.eee = control number of f(x)

-With f(x) = x3 + 105 x + 12 , it yields g(x) = 144 x3 + 105 x + 1
and "IRR?" proves ( in 4mn05s ) that g ( and f ) are irreducible since g(34) = 5663347 is a prime.

-If flag F00 is set, "PC/C" will replace f(x) by h(x) = f(| a0 |x)/| a0 | , it may sometimes work. For example,

with  f(x) = x7 - 2 x6 + 3 x5 - 4 x4 + 7 x3 - 5 x2 + 6 x - 8                                                 "IRR?" doesn't work.

then  g(x) = 262144 x7 + 65536 x6 + 12288 x5 + 2048 x4 + 448 x3 + 40 x2 + 6 x + 1     "IRR?" still doesn't work.
but  h(x) =  262144 x7 - 65536 x6 + 12288 x5 - 2048 x4 + 448 x3 - 40 x2 + 6 x - 1   and "IRR?" finds that h(2) = 29724011 is a prime

-So, h is irreducible and f is irreducible too.

-Of course, the programs will run much faster if you use the M-code routine "PR?"

-In rare cases, "PL" could produce a wrong f(m) < E10 because an intermediate result is > E10. So, "PL" could be modify as follows:

 01  LBL "PL"  02  0  03  LBL 01  04  RCL Y  05  *  06  ABS  07   E10  08  X L  12  RCL IND Z  13  +  14  ABS  15   E10   16  X L  20  ISG Z  21  GTO 01  22  END

2°)  Factorization

a)  The Brute Force

-Given a polynomial  f(x) = an xn + an-1 xn-1 + ... + a1 x + a0  ( a0 # 0 )  ,  "FCTR" seeks divisors g(x) = bm xm + bm-1 xm-1 + ... + b1 x + b0 such that:

bm is a divisor of  an , bm > 0
b0 is a divisor of  a0
| b1 | , | b2 | , ........... , | bm-1 |  <=  N    where N is a positive integer that you choose.

-The HP-41 emits a BEEP for the last factor, and a TONE 9 for the other ones ( if any ).

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

•  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0 # 0           bb > 14         Ree+1 thru Ree+2+n+Int(n/2)  are also used
Flags: /
Subroutines: /

-If you don't have an HP-41CX, replace line 53 by  SIGN  CLX  LBL 16  STO IND L  ISG L  GTO 16
-If you don't have an X-Functions Module, replace lines 66-68 by    RCL 01   RCL 03   INT   XEQ 11
-Lines 78-157-158-164-175-190-205-210-231-260  are three-byte  GTOs
-Line 232  LBL 11:  This subroutine is simply "LCO" ( cf "Miscellaneous short Routines for the HP-41" )

 01  LBL "FCTR"   02  STO 00   03  X<>Y   04  STO 01   05  INT   06  STO 14   07  LASTX   08  FRC   09   E3   10  *   11  STO 03   12  X<>Y   13  -   14  STO 02   15  STO 13   16  RCL IND 01   17  ABS   18  STO 04   19  RCL IND 03   20  ABS   21  STO 05   22  SIGN   23  ST+ 03   24  RCL 02   25  +   26  STO 09   27  RCL 03   28  +   29  STO 06   30  1   31  STO 08   32  +   33  STO 07   34   E3   35  /   36  ST+ 06   37  RCL 09   38  LASTX   39  /   40  RCL 03   41  +   42   E3   43  /   44  ST+ 14 45  RCL 06   46  INT   47  DSE X   48   E3   49  /   50  ST+ 03   51  LBL 01   52  RCL 06       53  CLRGX     54  SIGN   55  STO IND 06   56  STO IND 07   57  RCL 06   58  ISG X   59   E-3   60  -   61  STO 09   62  STO 10   63  LBL 02   64  RCL 13   65  STO 11   66  RCL 14    67  REGMOVE    68  RCL 03   69  STO 12   70  LBL 03    71  RCL 12     72  RCL 06   73  RCL IND Y   74  RCL IND Y   75  /   76  FRC   77  X#0?   78  GTO 07     79  X<> L   80  STO IND Z   81  ISG Z   82  SIGN   83  ISG Y   84  X=0?   85  GTO 05   86  LBL 04   87  CLX   88  RCL IND Y 89  LASTX   90  *   91  ST- IND Z   92  ISG Z   93  CLX   94  ISG Y   95  GTO 04   96  LBL 05   97  ISG 12   98  CLX   99  DSE 11 100  GTO 03 101  RCL 12 102  LBL 06 103  RCL IND X 104  X#0? 105  GTO 00 106  RDN 107  ISG X 108  GTO 06 109  ENTER^ 110  DSE Y 111  LBL 00 112  RCL IND Y 113  X#0? 114  GTO 07 115  RCL 03   116  INT 117  RCL 12 118  INT 119  DSE X 120   E3 121  / 122  + 123  STO 12 124  RCL 06 125  TONE 9 126  RTN 127  RCL 12 128  RCL 01 129  INT 130  XEQ 11 131  STO 01 132  RCL IND X 133  ABS 134  STO 04 135  X<>Y 136  FRC 137   E3 138  * 139  RCL IND X 140  ABS 141  STO 05 142  CLX 143  RCL 01 144  INT 145  - 146  ENTER^ 147  X<> 02 148  - 149  ST+ 13 150   E6 151  / 152  ST+ 14 153  RCL 02 154  RCL 08 155  ST+ X 156  X>Y? 157  GTO 14   158  GTO 02  159  LBL 07 160  RCL IND 07 161  CHS 162  STO IND 07 163  X<0? 164  GTO 02   165  LBL 08 166  ISG IND 07 167  CLX 168  RCL 05 169  RCL IND 07 170  X>Y? 171  GTO 00 172  / 173  FRC 174  X=0? 175  GTO 02  176  GTO 08 177  LBL 00 178  SIGN 179  STO IND 07 180  LBL 09 181  ISG IND 06 182  CLX 183  RCL 04 184  RCL IND 06 185  X>Y? 186  GTO 00 187  / 188  FRC 189  X=0? 190  GTO 02   191  GTO 09 192  LBL 00 193  SIGN 194  STO IND 06 195  RCL 06 196  + 197  ISG X 198  X=0? 199  GTO 13 200  RCL 00 201  RCL IND 09 202  CHS 203  STO IND 09 204  X<0? 205  GTO 02    206  1 207  + 208  STO IND 09 209  X<=Y? 210  GTO 02    211  LBL 10 212  CLX 213  STO IND 09 214  ISG 09 215  X#0? 216  GTO 13 217  RCL IND 09 218  CHS 219  STO IND 09 220  X<0? 221  GTO 00 222  1 223  + 224  STO IND 09 225  RCL 00 226  X L 245   E3 246  / 247  + 248  RTN 249  LBL 13 250   E-3 251  ST+ 06 252  SIGN 253  ST+ 07 254  ST+ 08 255  ST- 13 256  RCL 02 257  RCL 08 258  ST+ X 259  X<=Y? 260  GTO 01  261  LBL 14 262  RCL 01 263  BEEP 264  END

( 386 bytes )

 STACK INPUTS OUTPUTS Y bbb.eee / X N bbb.eeek

bbb.eee = control number of f(x)                      bbb.eeek = control number of the successive factors.           bbb > 014
N     =  positive integer

Example:            f(x) = 2 x5 + 5 x4 + 15 x3 + 28 x2 + 27 x + 35

Store for instance  [ 2  5  15  28  27  35 ]   into  R15  R16  R17  R18  R19  R20   and with  N = 4

15.020  ENTER^
4      XEQ "FCTR"   >>>>  ( TONE 9 )   27.029      --- Execution time = 4mn06s ---

-So, a factor g1(x) has been found,    R27 = 2   R28 = 1   R29 = 5    whence      g1(x) = 2 x2 + x + 5

Press  R/S   >>>>   almost immediately  ( BEEP )   15.018

-The last factor g2(x) has been found   R15 = 1  R16 = 2   R17 = 4   R18 = 7     and       g2(x) = x3 + 2 x2 + 4 x + 7

-Thus,  f(x) = ( 2 x2 + x + 5 )( x3 + 2 x2 + 4 x + 7 )

-This factorization is complete in Z[X] but of course neither in R[X] nor - a fortiori - in C[X]

Notes:

-The N-value does not influence the search for factors of degree 1
-The problem is that we don't know in advance the best choice for N:
if N is too small, we can miss a factor and large N-values lead to prohibitive execution time.
-So, use "FCTR" with a good emulator like V41 in this case.

-For example,  f(x) = x9 + 6 x8 + 17 x7 + 40 x6 + 65 x5 + 74 x4 + 74 x3 + 36 x2 + 23 x + 6
= ( x4 + 4 x3 + 5 x2 + 7 x + 2 )( x5 + 2 x4 + 4 x3 + 7 x2 + x + 3 )

-This factorization is found ( with N = 7 ) in about 9 minutes with  V41 in turbo mode.
-With a true HP-41, it would probably last almost 4 days!

-The following result gives a bound on the bj's when b(x) is a divisor of a(x)

•  | bj |  <=  Cjm-1 || a || + Cj-1m-1 | an |      where   n = deg a , m = deg b , || a || = ( an2 +  an-12 + ..... + a02 )1/2 and  Cjm = m!/[j!(m-j)!]

-Unfortunately, it's seldom a small number.

b)  Via Factorization over finite Fields

-We can use a factorization over Z/pZ - where p is a prime ( cf Polynomial factorization over finite fields for the HP-41" ) - and deduce the factorization in Z[X]

-For instance:  f(x) = x9 + 6 x8 + 17 x7 + 40 x6 + 65 x5 + 74 x4 + 74 x3 + 36 x2 + 23 x + 6   over  Z/17Z

"SQFREE"  proves that f(x) is squarefree in  Z/17Z[X]

"PRFDK"   then gives no factor of degrees 1 , 2 , 3
exactly ( x4 + 4 x3 + 5 x2 + 7 x + 2 )  for the product of polynomials of degree 4
and  ( x5 + 2 x4 + 4 x3 + 7 x2+ x + 3 ) for the product of polynomials of degree 5

-Here, we must of course check that the product of these 2 irreducible polynomials = f(x) and it works!
( if the term 17 x7 had been omitted, we would have gotten the same result - true in Z/17Z[X] - but false in Z[X] )

-This example gives a misleading impression of easiness and in many cases,
we'll have to choose several p-values ( or one large p-value if we use p > twice the bound mentioned above )
-Several "parasitic" factors will appear, and many tests will have to be made before finding the right factorization in Z[X]

-The following routine may be used to check if a polynomial g(x) is really a divisor of a polynomial f(x) in Z[X]

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

•  Rbb = an    •  Rbb+1 = an-1    ........................  •  Ree = a0             bbb > 003
•  Rbb' = a'm  •  Rbb'+1 = a'm-1  ........................  •  Ree' = a'0           bbb' > 003
Flags: /
Subroutines: /

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

( 107 bytes )

 STACK INPUTS OUTPUTS Y bbb.eee(f) bbb.eee(r) X bbb.eee(g) bbb.eee(q) or 0

(bbb.eee)f = the control number of the dividend                        (bbb.eee)r  = the control number of the remainder       ( all  bbb's > 003 )
(bbb.eee)g = the control number of the divisor                           (bbb.eee)q = the control number of the quotient

>>>    we assume that           deg f > deg g

X-output = 0  as soon as a coefficient of the quotient is not an integer ( in this case, Y-output is meaningless )

-So g is a divisor of f in Z[X] if X-output # 0 , eee(r) = bbb(r)  and  Rbb(r) = 0

-This program is called as a subroutine by the program hereunder.

c)  Via Bairstow's Method

c-1)  Program#1

-Bairstow's method expresses a polynomial f(x) as a product of factors of degree 2 , multiplied by a linear term if deg f is odd.
-"FCTR1" first replaces the quadratic factors whose discrimant in non-negative by 2 polynomials of degree 1.
In other words, we get f(x) as a product of irreducible factors in R[X]
-Then it computes all the - useful - combinations of these products, rounds the coefficients according to the current display settings
and tests ( when we get integer coefficients ) if this candidate is really a divisor of f(x).

-To simplify the program, we assume that the leading coefficient  an = +/-1  though it may sometimes work in other cases too.

Data Registers:       R00 to R13+4n: temp                               ( Registers Rbb thru Ree are to be initialized before executing "FCTR1" )

R04 = counter   R05 = eee+1   R06 = address of the first free register
R07 = control number of the list of the irreducible factors in R[X]
R08 = control number of a potential divisor
R09 = control number of f(x)
R10 = deg f                                  R11 = counter
R12 = control number of the product of k-1 factors while we are testing products of k factors
R13 = counter = 2m-1 - 1 , 2m-1 - 2 ,  .......  , 0

•  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0      bb > 13+4n       and several other registers to store the products of polynomials.

Flags:  F00-F01-F26
Subroutines:   "BRST" ( lines 01 to 89 )   cf "Polynomials for the HP-41"
"P2"  "¨PRO"                     cf "Polynomials for the HP-41"
"LCO"   "LDL"                  cf "Miscellaneous short routines for the HP-41"
"DIVN"                             cf  paragraph above

-XEQ "LCO" may of course be replaced by the M-code routine LCO
-Lines 119 to 132 are not absolutely necessary: they simply rearrange the order of the factors
so that the factors of degree 1 are placed before the factors of degree 2 in the list whose control number = R07
-So, you can delete these lines but if you delete lines 119 to 132, delete lines 284 to 293 too!
-The termination criterion wouldn't work well in all cases.
-You'll have sometimes to wait a little more to get the second factor, but it will obviously save bytes.
-You can also employ the other termination criterion used by "FCTR2" below.

-Lines 212-294-307-313-320-330  are three-byte  GTOs

 01  LBL "FCTR1"   02  STO 09   03  ENTER^   04  ENTER^   05  FRC   06   E3   07  *   08  X<>Y   09  INT   10  -   11  STO 10   12  -   13  INT   14  DSE X   15  XEQ "LCO"    16  PI   17  2   18  CF 26   19  XEQ "BRST"   20  SF 26   21  RCL 10   22  14   23  +   24  STO 05   25  STO 06   26  RCL 10   27  ST+ X   28  LASTX   29  +   30  STO 01   31  RCL 10   32  2   33  MOD   34  SF 01   35  X=0?   36  GTO 08   37  CF 01   38  RCL 08   39  INT   40  .1   41  %   42  +   43   E-3   44  +   45  RCL 01   46  XEQ "LCO"   47  STO IND 06   48  2   49  ST+ 01   50  ST+ 08   51  ISG 06   52  LBL 08   53  1   54  FS? 01   55  RCL IND 08   56  STO IND 01 57  ISG 01   58  CLX   59  FS?C 01   60  ISG 08   61  RCL IND 08   62  ISG 08   63  RCL IND 08   64  XEQ "P2"   65  FS?C 00   66  GTO 08   67  CHS      68  STO IND 01   69  RCL 01   70  .1   71  %   72  +   73  1   74  ST+ 01   75  STO IND 01   76  -   77  STO IND 06   78  LASTX   79  ST+ 01   80  ST+ 06   81  CLX   82  R^   83  CHS   84  GTO 00   85  LBL 08   86  RCL 08   87  1   88  -   89  RCL IND X   90  STO IND 01   91  LASTX   92  ST+ 01   93  CHS   94  RCL IND 08   95  LBL 00   96  STO IND 01   97  CLX   98  RCL 01   99  .1 100  % 101  + 102  + 103  1 104  - 105  STO IND 06 106  LASTX 107  ST+ 01 108  ST+ 06 109  ISG 08 110  GTO 08 111  DSE 06 112  RCL 06 113   E3 114  / 115  RCL 05 116  ST- 06 117  + 118  STO 07 119  STO 01    120  STO 02  121  LBL 09 122  RCL IND 01 123  ISG X     124  ISG X 125  GTO 09    126  RCL IND 01 127  X<> IND 02 128  STO IND 01 129  ISG 02 130  LBL 09 131  ISG 01 132  GTO 09 133  RCL 09 134  INT 135  RCL 10 136  + 137  1 138  + 139  STO 05 140  14.014 141  STO 04 142  2 143  RCL 06    144  Y^X 145  DSE X  146  STO 13 147  RCL 07  148  STO 14 149  GTO 12 150  LBL 10  151   E-3       152  ST+ 04  153  LBL 11  154  RCL 07  155   E-3        156  RCL 04  157  STO 11  158  FRC    159   E3 160  * 161  RCL 04 162  INT 163  - 164  * 165  - 166  LBL 01 167  STO IND 11 168  1.001 169  + 170  ISG 11 171  GTO 01 172  LBL 12 173  RCL 04 174  STO 11 175  RCL 05 176  STO 06 177  SF 01 178  LBL 02 179  STO 12 180  RCL IND 11 181  X<>Y 182  RCL IND Y 183  ISG 11 184  X=0? 185  GTO 03 186  FS?C 01 187  GTO 02 188  RCL 06 189  XEQ "PRO" 190  ENTER^ 191  FRC 192   E3 193  * 194  1 195  + 196  STO 06 197  X<>Y 198  GTO 02 199  LBL 03 200  RCL 06 201  FS? 01 202  XEQ "LCO" 203  FC?C 01 204  XEQ "PRO" 205  ENTER^ 206  LBL 04 207  RCL IND Y 208  RND      209  STO IND Z 210  FRC 211  X#0?   212  GTO 06   213  RDN 214  ISG Y 215  GTO 04 216  STO 08 217  RCL 09 218  RCL 08 219  FRC 220   E3 221  * 222  1 223  + 224  XEQ "LCO" 225  RCL 08 226  XEQ "DIVN" 227  X=0? 228  GTO 06 229  RCL IND Y 230  X#0? 231  GTO 06 232  X<>Y 233  STO 01 234  RCL 08 235  TONE 9 236  RTN      237  RCL 01 238  RCL 09 239  INT 240  XEQ "LCO" 241  STO 09 242  RCL 04 243  FRC 244  LASTX 245  INT 246  DSE X 247   E3 248  ST* Z 249  / 250  + 251  STO 01 252  RCL 07 253  LBL 05    254  RCL IND 01 255  RDN 256  RCL IND T 257  XEQ "LDL" 258  DSE 01 259  GTO 05 260  STO 07 261  STO 14 262  FRC 263   E3 264  * 265  RCL 07 266  INT 267  - 268  2 269  X<>Y 270  Y^X 271  1 272  - 273  STO 13 274  X<=0? 275  GTO 00 276  RCL 09 277  FRC 278   E3 279  * 280  RCL 09 281  INT 282  - 283  STO 10 284  RCL 08 285  FRC 286   E3 287  * 288  RCL 08 289  INT 290  -          291  ST+ X 292  X>Y?  293  GTO 00 294  GTO 11  295  LBL 06 296  DSE 13 297  FS? 30 298  GTO 00 299  DSE 11 300  CLX 301  RCL 04 302  ISG X 303  ABS  304  ISG X 305  GTO 06 306  ISG IND 11 307  GTO 12   308  LBL 07 309  RCL 11 310  INT 311  14 312  X=Y? 313  GTO 10  314  SIGN 315  ST- 11 316  ISG IND 11 317  X=0? 318  GTO 07 319  RCL IND 11 320  GTO 01   321  LBL 06 322  ISG IND 11 323  X=0? 324  GTO 07 325  RCL IND 11 326  RCL IND X 327  RCL 12  328  ISG 11 329  CLX 330  GTO 03    331  LBL 00 332  RCL 09 333  BEEP 334  END

( 531 bytes )

 STACK INPUT OUTPUTk X bbb.eee bbb.eeek

where   bbb.eee = control number of f(x)                      bbb.eeek = control number of the successive factors.           bbb > 013+4n     with  n = deg f

-The HP-41 emits a BEEP for the last factor, and a TONE 9 for the other ones, if any.

Example1:         f(x) = x9 + 6 x8 + 17 x7 + 40 x6 + 65 x5 + 74 x4 + 74 x3 + 36 x2 + 23 x + 6

deg f = 9 so we must choose  bbb > 049, for instance store

[ 1  6  17  40  65  74  74  36  23  6 ]  into  [ R50  R51  R52  R53  R54  R55  R56  R57  R58  R59 ]

FIX 4  ( if you choose this format )

50.059  XEQ "FCTR1"  >>>>  ( TONE 9 )    63.067                         --- Eecution time = 7mn41s ---

[ R63  R64  R65  R66  R67 ]  =  [ 1  4  5  7  2 ]

R/S    >>>>   ( BEEP )    50.055

[ R50  R51  R52  R53  R54  R55 ]  =  [ 1  2  4  7  1  3 ]

-Thus    f(x) =  ( x4 + 4 x3 + 5 x2 + 7 x + 2 )( x5 + 2 x4 + 4 x3 + 7 x2 + x + 3 )

Example2:          f(x) = x8 - 128 x6 - x5 + 5120 x4 + 64 x3 - 65538 x2 - 512 x + 262144

We can choose  bbb = 046

With  [ R46  R47  R48  R49  R50  R51  R52  R53  R54 ]  =  [ 1  0  -128  -1  5120  64  -65538  -512  262144 ]

FIX 4
46.054  XEQ "FCTR1"  >>>>   ( TONE 9 )    62.066                            --- Execution time = 17mn32s ---

[ R62  R63  R64  R65  R66 ]  =  [ 1  0  -64  -2  512 ]

R/S    >>>>  ( BEEP )   46.050

[ R46  R47  R48  R49  R50 ]  =  [ 1  0  -64   1   512 ]

-So,    f(x) = ( x4 - 64 x2 - 2 x + 512 )( x4 - 64 x2 + x + 512 )

-The execution time is much longer because f(x) has actually 8 real roots and many more polynomial multiplications have been performed.

-In this example, we wouldn't have found this factorization in FIX 5:
-Before rounding the coefficients - line 208 - the product was  [ 1  0.000000999  -63.99999993  -2.000006500  511.9999998 ]
and  -2.000006500  would have been rounded to  -2.00001  which is not an integer!

-If deg f is large and/or if there are many factors in R[X], we'll have to execute"FCTR1" in FIX 1 or FIX 0
and even then, we could miss some factors because of roundoff errors...
-So the result is less guaranteed than via factorization over finite fields.

Remarks:

-Multiple roots often lead to slow convergence and inaccurate results with "BRST"
-So, try to use "MSR" ( cf polynomials §1-h ) to get a squarefree polynomial.
-If the process seems to diverge while "BRST" is executing, stop the subroutine, store other guesses in R06 & R07 and press GTO 01  R/S
-But do not press  XEQ 01 :  that would clear the return stack !
-You can also use the display setting for the termination criterion of "BRST" ( if you have replaced lines 67-68 by  RND  X#0? )
-In this case, add  FIX IND 11 after line 19 and  STO 11   X<>Y   after line 01  in the listing above,
place bbb.eee in Y-register and 4 in X-register ( for FIX 4 )  before executing "FCTR1"

-If the leading coefficient  c = an # +/-1 , "FCTR1" can fail.
-Apply the program to h(x) = cn-1 f(x/c)
-After factorizing h(x), the factorization of f(x) will be easy to deduce.

c-2)  Program#2

-Instead of calculating all these products of polynomials, we can also calculate the product of their constant coefficients first,
then rounding it in the current format, and perform the polynomial multiplications only if this coefficient is an integer.
-We can thus expect a much faster program... provided the display setting is not FIX 0 !

Data Registers:       R00 to R12+4n: temp                                             ( Registers Rbb thru Ree are to be initialized before executing "FCTR2" )

R04 = counter   R05 = eee+1   R06 = address of the first free register
R07 = control number of the list of the irreducible factors in R[X]
R08 = control number of a potential divisor
R09 = control number of f(x)
R10 = counter = 2m-1 - 1 , 2m-1 - 2 ,  .......  , 0

•  Rbb = an  •  Rbb+1 = an-1  ........................  •  Ree = a0        bb > 11+4n       and several other registers to store the products of polynomials.

Flags:  F00-F01-F26
Subroutines:   "BRST" ( lines 01 to 89 )   cf "Polynomials for the HP-41"
"P2"  "¨PRO"                     cf "Polynomials for the HP-41"
"LCO"   "LDL"                  cf "Miscellaneous short routines for the HP-41"
"DIVN"                             cf  paragraph 2-b)
"CNK"                              cf "Statistics & Probability for the HP-41"

-Lines 182-219-304-312-318-325  are three-byte  GTOs

 01  LBL "FCTR2"   02  STO 12   03  X<>Y   04  STO 09   05  ENTER^   06  ENTER^   07  FRC   08   E3   09  *   10  X<>Y   11  INT   12  -   13  STO 10   14  -   15  INT   16  DSE X   17  XEQ "LCO"   18  PI   19  2   20  CF 26   21  XEQ "BRST"   22  SF 26   23  FIX IND 12   24  RCL 10   25  12   26  +   27  STO 05   28  STO 06   29  RCL 10   30  ST+ X   31  LASTX   32  +   33  STO 01   34  RCL 10   35  2   36  MOD   37  SF 01   38  X=0?   39  GTO 08   40  CF 01   41  RCL 08   42  INT   43  .1   44  %   45  +   46   E-3   47  +   48  RCL 01   49  XEQ "LCO"   50  STO IND 06   51  2   52  ST+ 01   53  ST+ 08   54  ISG 06   55  LBL 08 56  1   57  FS? 01   58  RCL IND 08   59  STO IND 01   60  ISG 01   61  CLX   62  FS?C 01   63  ISG 08   64  RCL IND 08   65  ISG 08   66  RCL IND 08   67  XEQ "P2"   68  FS?C 00   69  GTO 08   70  CHS   71  STO IND 01   72  RCL 01   73  .1   74  %   75  +   76  1   77  ST+ 01   78  STO IND 01   79  -   80  STO IND 06   81  LASTX   82  ST+ 01   83  ST+ 06   84  CLX   85  R^   86  CHS   87  GTO 00   88  LBL 08   89  RCL 08   90  1   91  -   92  RCL IND X   93  STO IND 01   94  LASTX   95  ST+ 01   96  CHS   97  RCL IND 08   98  LBL 00   99  STO IND 01 100  CLX 101  RCL 01 102  .1 103  % 104  + 105  + 106  1 107  - 108  STO IND 06 109  LASTX 110  ST+ 01 111  ST+ 06 112  ISG 08 113  GTO 08 114  DSE 06 115  RCL 06 116   E3 117  / 118  RCL 05 119  ST- 06 120  + 121  STO 07 122  RCL 09 123  INT 124  RCL 10 125  + 126  1 127  + 128  STO 05 129  12.012 130  STO 04 131  2 132  RCL 06 133  Y^X 134  DSE X 135  STO 10 136  RCL 07 137  STO 12 138  GTO 12 139  LBL 10 140   E-3 141  ST+ 04 142  LBL 11 143  RCL 07 144   E-3 145  RCL 04 146  STO 11 147  FRC 148   E3 149  * 150  RCL 04 151  INT 152  - 153  * 154  - 155  LBL 01 156  STO IND 11 157  1.001 158  + 159  ISG 11 160  GTO 01 161  LBL 12 162  RCL 04 163  STO 11 164  RCL 05 165  STO 06 166  SIGN 167  LBL 02 168  RCL IND 11 169  X<>Y 170  RCL IND Y 171  FRC 172   E3 173  * 174  X<>Y 175  RCL IND Y 176  * 177  ISG 11 178  GTO 02 179  RND 180  FRC 181  X#0? 182  GTO 06  183  RCL 04 184  STO 11 185  SF 01 186  LBL 03 187  RCL IND 11 188  X<>Y 189  RCL IND Y 190  ISG 11 191  X=0? 192  GTO 00 193  FS?C 01 194  GTO 03 195  RCL 06 196  XEQ "PRO" 197  ENTER^ 198  FRC 199   E3 200  * 201  1 202  + 203  STO 06 204  X<>Y 205  GTO 03 206  LBL 00 207  RCL 06 208  FS? 01 209  XEQ "LCO" 210  FC?C 01 211  XEQ "PRO" 212  ENTER^ 213  LBL 04 214  RCL IND Y 215  RND 216  STO IND Z 217  FRC 218  X#0? 219  GTO 06  220  RDN 221  ISG Y 222  GTO 04 223  STO 08 224  RCL 09 225  RCL 08 226  FRC 227   E3 228  * 229  1 230  + 231  XEQ "LCO" 232  RCL 08 233  XEQ "DIVN" 234  X=0? 235  GTO 06 236  RCL IND Y 237  X#0? 238  GTO 06 239  X<>Y 240  STO 01 241  RCL 08 242  TONE 9 243  RTN 244  RCL 01 245  RCL 09 246  INT 247  XEQ "LCO" 248  STO 09 249  RCL 04 250  FRC 251  LASTX 252  INT 253  DSE X 254   E3 255  ST* Z 256  / 257  + 258  STO 01 259  RCL 07 260  LBL 05 261  RCL IND 01 262  RDN 263  RCL IND T 264  XEQ "LDL" 265  DSE 01 266  GTO 05 267  STO 07 268  STO 12 269  FRC 270   E3 271  * 272  RCL 07 273  INT 274  - 275  STO 01 276  2 277  X<>Y 278  Y^X 279  1 280  - 281  STO 10 282  RCL 04 283  FRC 284   E3 285  * 286  RCL 04 287  INT 288  - 289  X=0? 290  GTO 00 291  STO 02 292  LBL 09 293  RCL 01 294  RCL 02 295  XEQ "CNK" 296  RND 297  ST- 10 298  DSE 02 299  GTO 09 300  LBL 00 301  RCL 10 302  X<=0? 303  GTO 00 304  GTO 11 305  LBL 06 306  DSE 10 307  FS? 30 308  GTO 00 309  DSE 11 310  CLX 311  ISG IND 11 312  GTO 12 313  LBL 07 314  RCL 11 315  INT 316  12 317  X=Y? 318  GTO 10  319  SIGN 320  ST- 11 321  ISG IND 11 322  X=0? 323  GTO 07 324  RCL IND 11 325  GTO 01  326  LBL 00 327  RCL 09 328  BEEP 329  END

( 520 bytes )

 STACK INPUT OUTPUTk Y bbb.eee / X m bbb.eeek

where   bbb.eee = control number of f(x)                                      bbb.eeek = control number of the successive factors.
and     m is used to round numbers in FIX m ( line 23 )               bbb > 011+4n     with  n = deg f

Example:          f(x) = x8 - 128 x6 - x5 + 5120 x4 + 64 x3 - 65538 x2 - 512 x + 262144

You can choose  bbb = 044

With  [ R44  R45  R46  R47  R48  R49  R50  R51  R52 ]  =  [ 1  0  -128  -1  5120  64  -65538  -512  262144 ]

>>>  FIX 8  if you want to execute "BRST" in this format, and if you want to execute the rest of "FCTR2" in FIX 4

44.052  ENTER^
4      XEQ "FCTR2"  >>>>   ( TONE 9 )    60.064            --- Execution time = 9mn20s ---          ( instead of 17mn32s with "FCRT1" )

[ R60  R61  R62  R63  R64 ]  =  [ 1  0  -64  -2  512 ]

R/S    >>>>    ( BEEP )   44.048     ( in a few seconds )

[ R44  R45  R46  R47  R48 ]  =  [ 1  0  -64   1   512 ]

-And    f(x) = ( x4 - 64 x2 - 2 x + 512 )( x4 - 64 x2 + x + 512 )

-In this example, we have saved more than 8 minutes!

3°)  Cyclotomic Polynomials

a)  Program#1

-Let n a positive integer, the cyclotomic polynomial Fn(x) = Product gcd(k,n)=1 ( x - rk )       where  rk = exp(2i.k.pi/n)     are the primitive ( complex ) roots of 1
-Though it is defined as a complex polynomial, all its coefficients are real, and more precisely integers.

-"CYCLO" uses this definition to compute its coefficients

Data Registers:  R00 to R08: temp  R09 ..... = the coefficients ... and several temporary data.
Flags: /
Subroutine:  "PRO" multiplication of 2 polynomials ( cf "Polynomials for the HP-41" )
"LCO" copying a list ( cf "Miscellaneous Short Routines for the HP-41" )

 01  LBL "CYCLO" 02  STO 04 03  2 04  XY 37  X#0? 38  GTO 02 39  SIGN 40  X#Y?    41  GTO 03 42  STO IND 07  43  RCL 07 44  + 45  RCL 05 46  360 47  * 48  RCL 04 49  / 50  COS 51  ST+ X 52  CHS 53  STO IND Y 54  1 55  ST+ Z 56  STO IND Z 57  RCL 06 58  RCL 07 59  RCL 08 60  XEQ "PRO"   61  RCL 06 62  INT 63  XEQ "LCO" 64  STO 06 65  2.002 66  ST+ 07 67  INT 68  ST+ 08 69  LBL 03 70  DSE 05 71  GTO 01 72  FIX 0 73  LBL 04  74  CLX 75  RCL IND Y   76  RND 77  STO IND Y 78  ISG Y 79  GTO 04 80  RCL 06 81  FIX 4 82  END

( 138 bytes / SIZE 012+ 2.phi(n) )     where  phi = Euler function

 STACK INPUT OUTPUT X n bbb.eee

where bbb.eee = control number of Fn(x)

Examples:

1   XEQ "CYCLO"  >>>>   1.002    [ R01  R02 ] = [ 1  -1 ]                 so  F1(x) = x - 1
2   XEQ "CYCLO"  >>>>   1.002    [ R01  R02 ] =  [ 1  1 ]                  so  F2(x) = x + 1
3   XEQ "CYCLO"  >>>>   9.011    [ R09  R10  R11 ] = [ 1  1  1 ]      so  F3(x) = x2 + x + 1                 --- Execution time =  7s ---

10           R/S            >>>>   9.013    [ R09  R10  R11  R12  R13 ]  =  [ 1  -1  1  -1  1 ]                          --- Execution time = 16s ---

thus,  F10(x) = x4 - x3 + x2 - x + 1

Note:

-If n is a prime, Fn(x) = xn-1 + xn-2 + ............. + x + 1     all the coefficients = 1
-Deg Fn = phi(n)  where  phi = Euler Totient function
-Roundoff errors become huge as n increases
-The coefficients are still exact for n = 43 but they are wrong for n = 47
-So we must use another method for large n-values.

b)  Program#2

-"CYCLO2"  uses the formula:    Fn(x) = Product d | n ( xd - 1 )µ(n/d)    where  µ = Moebius function

-First, the product is performed when µ(n/d) = +1 ( lines 22 to 65 )
-Then, the result is divided by the terms corresponding to µ(n/d) = -1  ( lines 66 to 99 )
-So, there is no roundoff error at all!

Data Registers:  R00 to R06: temp ( R03 = n )
R07  R08 .....  = coefficients of Fn(x) ... and other temporary data
Flags: /
Subroutine:  "MOEB" or "MOEB2"  ( cf "Primes & Number Theoretic Functions for the HP-41" )

-If you don't have an HP-41CX, replace lines 17 & 48 by  XEQ "LCL"  ( cf "Miscellaneous Short Routines for the HP-41" )
replace line 45 by RCL Y and line 49 by CLX

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

( 151 bytes )

 STACK INPUT OUTPUT X n bbb.eee

where bbb.eee = control number of Fn(x)

Example1:

20    XEQ "CYCLO2"  >>>>  7.015                       --- Execution time = 33s ---

[ R07  R08  ........  R15 ]  =  [ 1  0  -1  0  1  0  -1  0  1 ]       whence     F20(x) = x8 - x6 + x4 - x2 + 1

Example2:

SIZE 128 or greater

105   XEQ "CYCLO2"   >>>>  7.055                       --- Execution time = 3m59s ---           ( 3m47s if you use "MOEB2" instead of "MOEB" )

[ R07 .......  R55 ] = [ 1  1  1  0  0  -1  -1  -2  -1  -1  0  0  1  1  1  1  1  1  0  0  -1  0  -1  0  -1
0  -1  0  -1  0  0  1  1  1  1  1  1  0  0  -1  -1  -2  -1  -1  0  0  1  1  1 ]

whence  F105(x) = x48 + x47 + x46 - x43 -x42 - 2 x41 - x40 - x39 + x36 + x35 + x34 + x33 + x32 + x31 - x28 - x26 - x24
- x22 - x20 + x17 + x16 + x15 + x14 + x13 + x12 - x9 - x8 - 2 x7 - x6 - x5 + x2 + x + 1

Notes:

F105(x) is the first cyclotomic polynomial with a coefficient ( in fact 2 ) different from  -1 , 0 , +1
All the cyclotomic polynomials are irreducible in Z[X]

4°)  Polynomial Factorization over Gaussian Integers

-We assume that f(z) is unitary - leading coefficient = 1 - and that the other coefficients
are Gaussian integers i-e complex numbers a + b.i  where  a and b are integers.

f(z) = cn zn + cn-1 zn-1 + .................. + c1 z + c0          with         ck = ak + i.bk

Data Registers:       R00 to R15+4n: temp                               ( Registers R01 thru Ree are to be initialized before executing "FCTRZ" )

R08 = counter   R09 = eee'+1   R10 = address of the first free register
R11 = control number of the list of the irreducible factors in C[X]
R12 = control number of a potential divisor
R13 = control number of f(z)
R14 = counter = 2m-1 - 1 , 2m-1 - 2 ,  .......  , 0    and several other registers to store the products of polynomials.

•  R01 = an = 1   •  R03 = an-1    ........................  •  Ree-1 = a0
•  R02 = bn = 0   •  R04 = bn-1    ........................  •  Ree = b0

Flags:  /
Subroutines:

"PLRZ"       cf "Polynomials for the HP-41"
"PROZ"      cf "Polynomials for the HP-41"
"DIVZ"       cf "Polynomials for the HP-41"
"LDL"         cf "Miscellaneous short routines for the HP-41"

-If possible, employ the versions that use the M-code routine Z*Z & Z/Z  ( cf "A few M-code Routines for the HP-41" )
-If you use the 1st version of "DIVZ" ( without Z*Z and Z/Z ), delete the tiny leading-coefficients in the remainder  ( "DELZ" )
before testing if the remainder = 0 ( lines 193-194 )
-And  add a loop to round the coefficients of the quotient  (  R-P  P-R  doesn't always return integers )
-For example, add   LBL 14   RCL IND X   RND   STO IND Y   X<>Y   ISG X   GTO 14   after line 207

-Other M-code routines like LCO , DEG F , CNK are used but they may be replaced by standard programs.

LCO  may be replaced by  XEQ "LCO"  cf "Miscellaneous short routines"

-Lines 05-14  may be replaced by  FRC   E3   *   RCL Y   INT   -
-Lines 90-237 may be replaced by  FRC   E3   *   RCL 08   INT   -
-Line 228 may be replaced by   FRC   E3   *   RCL 11   INT   -

-Lines 155 to 160 may be replaced by  RCL IND 15   RCL IND X   RCL 10   XEQ "LCO"
but add  STO 01  after line 147  and add   RCL 01  after line 163

-Lines 169-170 may be replaced  by  FRC   E3   *   1   +

-Line 244 may be replaced by  XEQ "CNK"  cf "Statistics & Probability for the HP-41"

-Lines  122-127-253-261-267-274 are three-byte  GTOs

 01  LBL "FCTRZ"   02  STO 00    03  X<>Y     04  ENTER^   05  DEG F      06  DSE X   07  ST+ X   08  16   09  +   10  LCO   11  STO 13   12  ENTER^   13  ENTER^   14  DEG F     15  STO 14   16  -   17  INT   18  DSE X   19  LCO       20  RCL 00   21  STO 12   22  CLX   23  PI   24  2   25  XEQ "PLRZ"   26  FIX IND 12   27  STO 01   28  STO 02   29  RCL 14   30  DSE X   31  2   32  /   33  STO 03   34  16   35  +   36  STO 04   37  LBL 08   38  RCL 01   39  INT   40  .1   41  %   42  +   43   E-3   44  +   45  STO IND Y   46  SIGN   47  + 48  ISG 01   49  ISG 01   50  GTO 08   51  LASTX   52  -   53   E3   54  /   55  RCL 04   56  +   57  STO 11   58  STO 16   59  SIGN   60  CHS   61  LBL 13   62  ST* IND 02   63  ISG 02   64  GTO 13   65  RCL 13   66  ISG X   67  RCL 14   68  +   69  INT   70  STO 09   71  16.016   72  STO 08   73  2   74  RCL 03   75  1   76  -   77  Y^X   78  1   79  -   80  STO 14   81  GTO 12   82  LBL 10   83   E-3   84  ST+ 08   85  LBL 11   86  RCL 11   87   E-3   88  RCL 08   89  STO 15   90  DEG F   91  *   92  -   93  LBL 01   94  STO IND 15 95  1.001   96  +   97  ISG 15   98  GTO 01   99  LBL 12 100  RCL 08 101  STO 15 102  RCL 09 103  STO 10 104  CLST 105  LBL 02        106  RCL IND 15 107  RDN 108  RCL IND T  109  X<>Y 110  RCL IND Y  111  + 112  ISG Y   113  RCL IND Y 114  R^ 115  + 116  X<>Y 117  ISG 15 118  GTO 02 119  RND 120  FRC 121  X#0? 122  GTO 06    123  X<>Y 124  RND 125  FRC 126  X#0? 127  GTO 06   128  RCL 08 129  STO 15 130  SIGN 131  STO IND 09 132  ST+ 10 133  CLX 134  STO IND 10 135  SIGN 136  ST+ 10 137  RCL IND 15 138  RCL IND X 139  RCL 10 140  LCO    141  2 142  ST+ 10 143  - 144  LBL 03 145  ISG 15 146  X=0? 147  GTO 00 148  1 149  STO IND 10 150  ST+ 10 151  CLX 152  STO IND 10 153  SIGN 154  ST+ 10 155  CLX 156  RCL IND 15 157  X<>Y 158  RCL IND Y 159  RCL 10 160  LCO   161  2         162  ST+ 10 163  - 164  RCL 10 165  XEQ "PROZ" 166  RCL 09 167  LCO   168  ENTER^ 169  DEG F  170  X<> L 171  STO 10 172  X<>Y 173  GTO 03 174  LBL 00 175  ENTER^ 176  LBL 04 177  RCL IND Y 178  RND 179  STO IND Z 180  FRC 181  X#0? 182  GTO 06 183  RDN 184  ISG Y 185  GTO 04 186  STO 12 187  RCL 13 188  RCL 10 189  LCO    190  RCL 12 191  XEQ "DIVZ" 192  STO 01     193  RCL IND Y 194  X#0? 195  GTO 06 196  ISG Z 197  RCL IND Z 198  X#0? 199  GTO 06 200  RCL 12 201  TONE 9 202  RTN 203  RCL 01 204  RCL 13 205  INT 206  LCO  207  STO 13   208  RCL 08   209  FRC 210  LASTX 211  INT 212  DSE X 213   E3 214  ST* Z 215  / 216  + 217  STO 01 218  RCL 11 219  LBL 05 220  RCL IND 01 221  RDN 222  RCL IND T 223  XEQ "LDL" 224  DSE 01 225  GTO 05 226  STO 11 227  STO 16 228  DEG F   229  STO 01 230  2 231  X<>Y 232  Y^X 233  1 234  - 235  STO 14 236  RCL 08 237  DEG F  238  X=0? 239  GTO 00 240  STO 02 241  LBL 09 242  RCL 01 243  RCL 02 244  CNK    245  RND 246  ST- 14 247  DSE 02 248  GTO 09 249  LBL 00 250  RCL 14 251  X<=0? 252  GTO 00 253  GTO 11  254  LBL 06 255  DSE 14 256  FS? 30 257  GTO 00 258  DSE 15 259  CLX 260  ISG IND 15 261  GTO 12   262  LBL 07 263  RCL 15 264  INT 265  16 266  X=Y? 267  GTO 10  268  SIGN 269  ST- 15 270  ISG IND 15 271  X=0? 272  GTO 07 273  RCL IND 15 274  GTO 01  275  LBL 00 276  RCL 13 277  BEEP 278  END

( 431 bytes )

 STACK INPUT OUTPUTk Y bbb.eee / X m bbb.eeek

where   bbb.eee = 001.eee = control number of f(z)                      bbb.eeek = control number of the successive factors.
and     m is used to round numbers in FIX m ( line 26 )

Example:           f(z) = z5 + ( 4+6.i) z4 + (-4+20.i) z3 + (-27+35.i) z2 + (-31+8.i) z + (-31+24.i)

-Store  1  0  4  6  -4  20  -27  35  -31  8  -31  24  into  R01 R02 ............... R12  respectively:  control number = 1.012

FIX 8   if you want to execute "PLRZ" in this format

and if you want to execute the rest of "FCTRZ" in FIX 4

1.012  ENTER^
4     XEQ "FCTRZ"  >>>>     ( TONE 9 )     48.053         --- Execution time = 5mn10s ---

[ R48  R49  R50  R51  R52  R53 ]  =  [ 1  0  1  2  2  7 ]       which corresponds to g(z) = z2 + (1+2.i) z + (2+7.i)

R/S    >>>>    ( BEEP )    36.043

[ R36  R37  R38  R39  R40  R41  R42  R43 ]  =  [ 1  0  3  4  -1  3  2  5 ]       which corresponds to h(z) = z3 + (3+4.i) z2 + (-1+3.i) z + (2+5.i)

-Thus,  f(z) = [ z2 + (1+2.i) z + (2+7.i) ] [ z3 + (3+4.i) z2 + (-1+3.i) z + (2+5.i) ]

Note:

-If the process seems to diverge while "PLRZ" is executing, stop the subroutine, store other guesses in R01 & R02 and press GTO 01  R/S
( XEQ 01  would clear the return stack ! )

Reference:

  R. Thangadurai - Mathematics Newsletter - Vol 17 #2 - "Irreducibility of Polynomials whose Coefficients are Integers"