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) = a_n xⁿ + a_n--1 x^n-1 + ... + a₁ x + a₀
-If there is a prime p such that

p is a divisor of a₀ , a₁ , a₂ , ............ , a_n-1
p is not a divisor of a_n , p² is not a divisor of a₀

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 a₀ = 0.

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

• Rbb = a_n • Rbb+1 = a_n-1 ........................ • Ree = a₀ bbb > 004
Flags: /
Subroutines: /

-LBL 01 this loop calculates GCD( a₀ , a₁ , ............ , a_n-1 )
-Lines 41 to 53 compute the prime divisors of GCD( a₀ , a₁ , ............ , a_n-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 X<Y?
12 X=0?
13 RTN
14 RCL IND 02
15 X=0?
16 RTN
17 RCL 01

18 ISG X
19 STO 03
20 RCL IND X
21 LBL 01
22 ISG 03
23 FS? 30
24 GTO 00
25 RCL IND 03
26 LBL 02
27 MOD
28 LASTX
29 X<>Y
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 x⁵ + 42 x³ + 98 x² + 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 = a_n • Rbb+1 = a_n-1 ........................ • Ree = a₀ 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) = x⁵ - 7 x⁴ + 19 x³ - 20 x² + 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) = x⁷ + 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) = a_n xⁿ + a_n-1 x^n-1 + ... + a₁ x + a₀ in Z[X] with a_k >= 0 for all k's ( a₀ # 0 )

• If there exists an integer m > max { a_k } 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 = a_n • Rbb+1 = a_n-1 ........................ • Ree = a₀ # 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 X<Y?
24 E^X

25 X<>Y
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 x⁴ + 3 x² + 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 x⁶ + 3 x⁵ + 4 x³ + 5 x² + 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 10¹⁰ 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) = a_n xⁿ + a_n-1 x^n-1 + ... + a₁ x + a₀ in Z[X] deg f > 1 ( a₀ # 0 )

m₁ = max { | a_i/a_n | , i = 0 , 1 , .... , n-1 }
m₂ = max { | a_i/a_n |^1/(n-i) , i = 0 , 1 , .... , n-1 }

H = min { 1+m₁^1/(r+1) , 2 m₂ } where r is such that a_n-1 = a_n-2 = ..... = a_n-r = 0 , ( r = 0 if a_n-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 = a_n • Rbb+1 = a_n-1 ........................ • Ree = a₀ # 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 X<Y?
34 X<>Y
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 X<Y?
64 E^X
65 X<>Y
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 x⁷ - 3 x⁴ + 4 x³ - 9 x² + 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 a_n is small - Cohn's criterion may be better ( for instance f(x) = x² + 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 10¹⁰ ( and we cannot test its primality in that case )
-This "OUT OF RANGE" may occur almost immediately if a_n 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 = a_n • Rbb+1 = a_n-1 ........................ • Ree = a₀ # 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 X<Y?
21 E^X
22 X<>Y
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) = x⁵ + 101 x⁴ + 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 = a_n • Rbb+1 = a_n-1 ........................ • Ree = a₀ # 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 X<Y?
23 E^X
24 X<>Y

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) = x⁵ + 101 x⁴ + 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) = x³ + 105 x + 12 ( Eisenstein criterion works obviously )
-In that case, try to apply the criterion to g(x) = f(a₀x)/a₀ ( a₀ = 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) = x³ + 105 x + 12 , it yields g(x) = 144 x³ + 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(| a₀ |x)/| a₀ | , it may sometimes work. For example,

with f(x) = x⁷ - 2 x⁶ + 3 x⁵ - 4 x⁴ + 7 x³ - 5 x² + 6 x - 8 "IRR?" doesn't work.

then g(x) = 262144 x⁷ + 65536 x⁶ + 12288 x⁵ + 2048 x⁴ + 448 x³ + 40 x² + 6 x + 1 "IRR?" still doesn't work.
but h(x) = 262144 x⁷ - 65536 x⁶ + 12288 x⁵ - 2048 x⁴ + 448 x³ - 40 x² + 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<Y?
09 E^X
10 RDN
11 X<> L
12 RCL IND Z
13 +
14 ABS
15 E10
16 X<Y?
17 E^X
18 RDN
19 X<> L
20 ISG Z
21 GTO 01
22 END

2°) Factorization

a) The Brute Force

-Given a polynomial f(x) = a_n xⁿ + a_n-1 x^n-1 + ... + a₁ x + a₀ ( a₀ # 0 ) , "FCTR" seeks divisors g(x) = b_m x^m + b_m-1 x^m-1 + ... + b₁ x + b₀ such that:

        b_m is a divisor of a_n , b_m > 0
        b₀ is a divisor of a₀
        | b₁ | , | b₂ | , ........... , | b_m-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 = a_n • Rbb+1 = a_n-1 ........................ • Ree = a₀ # 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<Y?
227 GTO 10
228 LBL 00
229 RCL 10
230 STO 09
231 GTO 02
232 LBL 11
233 ENTER^
234 SIGN
235 ST- L
236 LBL 12
237 ISG L
238 CLX
239 CLX
240 RCL IND Z
241 STO IND L
242 ISG Z
243 GTO 12
244 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.eee_k

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

Example: f(x) = 2 x⁵ + 5 x⁴ + 15 x³ + 28 x² + 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 g₁(x) has been found, R27 = 2 R28 = 1 R29 = 5 whence g₁(x) = 2 x² + x + 5

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

-The last factor g₂(x) has been found R15 = 1 R16 = 2 R17 = 4 R18 = 7 and g₂(x) = x³ + 2 x² + 4 x + 7

-Thus, f(x) = ( 2 x² + x + 5 )( x³ + 2 x² + 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) = x⁹ + 6 x⁸ + 17 x⁷ + 40 x⁶ + 65 x⁵ + 74 x⁴ + 74 x³ + 36 x² + 23 x + 6
= ( x⁴ + 4 x³ + 5 x² + 7 x + 2 )( x⁵ + 2 x⁴ + 4 x³ + 7 x² + 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 b_j's when b(x) is a divisor of a(x)

• | b_j | <= C_j^m-1 || a || + C_j-1^m-1 | a_n | where n = deg a , m = deg b , || a || = ( a_n² + a_n-1² + ..... + a₀² )^1/2 and C_j^m = 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) = x⁹ + 6 x⁸ + 17 x⁷ + 40 x⁶ + 65 x⁵ + 74 x⁴ + 74 x³ + 36 x² + 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 ( x⁴ + 4 x³ + 5 x² + 7 x + 2 ) for the product of polynomials of degree 4
                         and ( x⁵ + 2 x⁴ + 4 x³ + 7 x²+ 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 x⁷ 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 = a_n • Rbb+1 = a_n-1 ........................ • Ree = a₀ bbb > 003
• Rbb' = a'_m • Rbb'+1 = a'_m-1 ........................ • Ree' = a'₀ 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 a_n = +/-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 = 2^m-1 - 1 , 2^m-1 - 2 , ....... , 0

• Rbb = a_n • Rbb+1 = a_n-1 ........................ • Ree = a₀ 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	OUTPUT_k
X	bbb.eee	bbb.eee_k

where bbb.eee = control number of f(x) bbb.eee_k = 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) = x⁹ + 6 x⁸ + 17 x⁷ + 40 x⁶ + 65 x⁵ + 74 x⁴ + 74 x³ + 36 x² + 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) = ( x⁴ + 4 x³ + 5 x² + 7 x + 2 )( x⁵ + 2 x⁴ + 4 x³ + 7 x² + x + 3 )

Example2: f(x) = x⁸ - 128 x⁶ - x⁵ + 5120 x⁴ + 64 x³ - 65538 x² - 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) = ( x⁴ - 64 x² - 2 x + 512 )( x⁴ - 64 x² + 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 = a_n # +/-1 , "FCTR1" can fail.
-Apply the program to h(x) = c^n-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" )

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

-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	OUTPUT_k
Y	bbb.eee	/
X	m	bbb.eee_k

where bbb.eee = control number of f(x) bbb.eee_k = 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) = x⁸ - 128 x⁶ - x⁵ + 5120 x⁴ + 64 x³ - 65538 x² - 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) = ( x⁴ - 64 x² - 2 x + 512 )( x⁴ - 64 x² + 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 F_n(x) = Product_gcd(k,n)=1 ( x - r_k ) where r_k = 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 X<Y?
05 GTO 00
06 SIGN
07 STO 01
08 X=Y?
09 CHS
10 STO 02
11 1.002
12 RTN
13 LBL 00
14 /
15 ENTER^
16 FRC
17 X=0?

18 SIGN
19 -
20 STO 05
21 9.009
22 STO 06
23 10.012
24 STO 07
25 13
26 STO 08
27 SIGN
28 STO 09
29 DEG
30 LBL 01
31 RCL 04
32 RCL 05
33 LBL 02
34 MOD

35 LASTX
36 X<>Y
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 F_n(x)

Examples:

   1   XEQ "CYCLO" >>>>   1.002    [ R01 R02 ] = [ 1 -1 ]                 so F₁(x) = x - 1
   2   XEQ "CYCLO" >>>>   1.002    [ R01 R02 ] = [ 1 1 ]                 so F₂(x) = x + 1
   3   XEQ "CYCLO" >>>>   9.011    [ R09 R10 R11 ] = [ 1 1 1 ]      so F₃(x) = x² + x + 1                 --- Execution time = 7s ---

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

thus, F₁₀(x) = x⁴ - x³ + x² - x + 1

Note:

-If n is a prime, F_n(x) = x^n-1 + x^n-2 + ............. + x + 1 all the coefficients = 1
-Deg F_n = 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: F_n(x) = Product_{d | n} ( x^d - 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 F_n(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 F_n(x)

Example1:

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

[ R07 R08 ........ R15 ] = [ 1 0 -1 0 1 0 -1 0 1 ] whence F₂₀(x) = x⁸ - x⁶ + x⁴ - x² + 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 F₁₀₅(x) = x⁴⁸ + x⁴⁷ + x⁴⁶ - x⁴³ -x⁴² - 2 x⁴¹ - x⁴⁰ - x³⁹ + x³⁶ + x³⁵ + x³⁴ + x³³ + x³² + x³¹ - x²⁸ - x²⁶ - x²⁴
- x²² - x²⁰ + x¹⁷ + x¹⁶ + x¹⁵ + x¹⁴ + x¹³ + x¹² - x⁹ - x⁸ - 2 x⁷ - x⁶ - x⁵ + x² + x + 1

Notes:

F₁₀₅(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) = c_n zⁿ + c_n-1 z^n-1 + .................. + c₁ z + c₀ with c_k = a_k + i.b_k

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 = 2^m-1 - 1 , 2^m-1 - 2 , ....... , 0    and several other registers to store the products of polynomials.

• R01 = a_n = 1 • R03 = a_n-1 ........................ • Ree-1 = a₀
• R02 = b_n = 0 • R04 = b_n-1 ........................ • Ree = b₀

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	OUTPUT_k
Y	bbb.eee	/
X	m	bbb.eee_k

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

Example: f(z) = z⁵ + ( 4+6.i) z⁴ + (-4+20.i) z³ + (-27+35.i) z² + (-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) = z² + (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) = z³ + (3+4.i) z² + (-1+3.i) z + (2+5.i)

-Thus, f(z) = [ z² + (1+2.i) z + (2+7.i) ] [ z³ + (3+4.i) z² + (-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:

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

hp41programs

Polynomials in Z[X] for the HP-41