Polynomials2

# Discriminant for the HP-41

Overview

1°)  Polynomials of degree 3
2°)  Polynomials of degree 4
3°)  Polynomials of degree n

1°)  Polynomials of Degree 3

-The discriminant of a polynomial of degree 2:   a x2 + b x + c   is easy to compute:    D = b2 - 4 a c
-With a polynomial of degree 3:    a x3 + b x2 + c x + d  ,   the formula is already more complicated:

D = b2 c2 - 4 a c3 - 4 b3 d - 27 a2 d2 + 18 a b c d

Data Registers:  R00 unused

R01 = a   R02 = b  R03 = c  R04 = d
Flags: /
Subroutines: /

 01 LBL "DIS3"  02 STO 04  03 X<>Y  04 STO 03            05 *  06 X<>Y  07 STO 02  08 * 09 X<>Y  10 STO 01  11 *  12 ST+ X  13 RCL 01            14 RCL 04  15 *  16 X^2 17 3  18 *  19 -  20 9  21 *  22 RCL 02            23 3  24 Y^X 25 RCL 04  26 *  27 RCL 03            28 3  29 Y^X  30 RCL 01  31 *  32 + 33 4  34 *  35 -  36 RCL 02            37 RCL 03  38 *  39 X^2  40 +  41 END

( 51 bytes / SIZE 005 )

 STACK INPUTS OUTPUTS T a Z b Y c / X d D

Examples:

•  2 x3 - 5 x2 - 7 x + 1

2  ENTER^
-5  ENTER^
-7  ENTER^
1  XEQ "DIS3"   >>>>    D = 5621               ---Execution time = 2s---

D > 0   so there are 3 distinct roots

•  2 x3 - 5 x2 + 7 x + 1

2  ENTER^
-5  ENTER^
7  ENTER^
1     R/S     >>>>    D = -2387

D < 0   so there is one real root

•  x3 + x2 - 5 x + 3

1  ENTER^  ENTER^
-5  ENTER^
3     R/S   >>>>    D = 0

D = 0   so there are ( at least ) 2 equal roots

2°)  Polynomials of Degree 4

-Of course, with a polynomial of degree 4:    a x4 + b x3 + c x2 + d x + e   ,  the formula is much more complicated:

D = 256 a3 e3 -192 a2 b d e2 - 128 a2 c2 e2 + 144 a2 c d2 e - 27 a2 d4 + 144 a b2 c e2 - 6 a b2 d2 e
- 80 a b c2 d e + 18 a b c d3 + 16 a c4 e - 4 a c3 d2 - 27 b4 e2 + 18 b3 c d e - 4 b3 d3 - 4 b2 c3 e + b2 c2 d2

Data Registers:  R00: temp                                      ( Registers R01 thru R05 are to be initialized before executing "DIS4" )

• R01 = a   • R02 = b   • R03 = c   • R04 = d   • R05 = e
Flags: /
Subroutines: /

 01 LBL "DIS4"  02 RCL 01  03 RCL 05  04 *  05 STO 00  06 16  07 *  08 RCL 02  09 RCL 04  10 *  11 12  12 *  13 -  14 RCL 03            15 X^2  16 8  17 *  18 -  19 RCL 00  20 X^2  21 *  22 RCL 00  23 RCL 03 24 X^2  25 X^2  26 *  27 +  28 RCL 00  29 RCL 02  30 *  31 RCL 03  32 X^2  33 *  34 RCL 04            35 *  36 5  37 *  38 -  39 16  40 *  41 RCL 01  42 RCL 04  43 *  44 X^2  45 RCL 03  46 * 47 RCL 05  48 *  49 RCL 02  50 LASTX  51 *  52 X^2  53 RCL 01  54 *  55 RCL 03            56 *  57 +  58 144  59 *  60 +  61 RCL 01  62 RCL 04  63 X^2  64 *  65 X^2  66 RCL 02  67 X^2  68 RCL 05  69 * 70 X^2  71 +  72 27  73 *  74 -  75 RCL 02  76 RCL 04            77 *  78 X^2  79 RCL 00  80 *  81 6  82 *  83 -  84 RCL 01  85 RCL 02  86 *  87 RCL 03  88 *  89 RCL 04  90 3  91 Y^X  92 * 93 RCL 02  94 3  95 Y^X  96 RCL 03  97 *  98 RCL 04            99 * 100 RCL 05 101 * 102 + 103 18 104 * 105 + 106 RCL 02 107 RCL 04 108 * 109 3 110 Y^X 111 RCL 03 112 LASTX 113 Y^X 114 STO 00 115 RCL 02 116 X^2 117 * 118 RCL 05 119 * 120 + 121 RCL 00         122 RCL 01 123 * 124 RCL 04 125 X^2 126 * 127 + 128 4 129 * 130 - 131 RCL 02 132 RCL 03 133 * 134 RCL 04 135 * 136 X^2 137 + 138 END

( 154 bytes / SIZE 005 )

 STACK INPUT OUTPUT X / D

Examples:

•  x4 + 3 x3 + x2 - 3 x - 2

1  STO 01  STO 03
3  STO 02  CHS  STO 04
-2  STO 05

XEQ "DIS4"  >>>>   D = 0          ( at least 2 equal roots )            ---Execution time = 6s---

•  x4 - 2 x3 - 13 x2 + 38 x - 24

1   STO 01
-2   STO 02
-13  STO 03
38  STO 04
-24  STO 05   R/S   >>>>   D = 176400   ( All roots are real or all roots are non-real )

•  x4 - 13 x3 + 60 x2 - 112 x + 64

1    STO 01
-13   STO 02
60   STO 03
-112  STO 04
64   STO 05   R/S  >>>>    D = 10

-Unfortunately, this result is wrong because it involves numbers > 1010
-So, roundoff-errors yield  10  wheras the exact result is  D = 0
-However, this program gives the correct discriminant if we use free42 !

-Thus, with an HP41, it's better to employ the following routine which calculates D with a determinant of a square matrix.

3°)  Polynomials of Degree n

-If p(x) = an xn + an-1 xn-1 + ..... + a1 x + a0  is a polynomial of degree n > 1,  the discriminant may be computed by the following formula:

D = (-1)n(n-1)/2 (1/an)  Det (M)    where  M  is the square matrix of order  2n - 1

an    0    0   .....      bn-1    0       0   ........
an-1  an   0   .....      bn-2   bn-1    0   ........                     where  the 1st derivative of p(x) is p'(x) = bn-1 xn-1 + ..... + b1 x + b0
an-2  an-1 an   .....    bn-3   bn-2  bn-1
...............................................................

0  ......... 0 ......a0     0       0       0    ...  b0

with   (n-1) columns      and   n columns
with p(x)               with p'(x)
coefficients           coefficients
and zeros              and zeros

-For instance,  if  p(x) = x4 - 13 x3 + 60 x2 - 112 x + 64  we have to calculate the determinant of the 7x7 matrix:

1      0     0       4      0       0      0
-13     1     0    -39     4       0      0
60   -13    1     120  -39     4      0
-112   60  -13  -112  120  -39     4
64  -112  60     0    -112  120  -39
0     64  -112   0       0   -112  120
0      0     64     0       0      0  -112

Data Registers:  R00: temp                                      ( Registers R01 thru Rn+1 are to be initialized before executing "DISC" )

• R01 = an   • R02 = an-1  ......   • Rn+1 = a0

Flags:     CF 00 = a pivot p is regarded as zero if  | p | < 10-7                   CF 01 = partial pivoting
SF 00 =  a pivot p is regarded as zero if    p = 0                         SF 01 = no pivoting

Subroutines: "DET"  ( cf "Determinant for the HP41"  §4a) )  after adding  RCL M  just after  LBL 08

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

( 145 bytes / SIZE (2n-1)2+2 )

 STACK INPUT OUTPUT X n D

Examples:

•  2 x3 - 5 x2 - 7 x + 1

CF 00   CF 01

2   STO 01
-5  STO 02
-7  STO 03
1  STO 04

3  XEQ "DISC"   >>>>    D = 5621               ---Execution time = 36s---

•  x4 - 13 x3 + 60 x2 - 112 x + 64

CF 00   CF 01

1    STO 01
-13   STO 02
60   STO 03
-112  STO 04
64   STO 05

4    XEQ "DISC"   >>>>    D = 0              ---Execution time = 75s---

Notes:

-We can also use a program that solves linear systems.
-For instance, with "LS3" listed in "Linear & Non-Linear Systems for the HP41" §1c) , replace lines 81 to 92 with

 81 ENTER^  82 XEQ "LS3"  83 R^  84 INT  85 X^2  86 1  87 +  88 X<>Y  89 RCL IND Y  90 /  91 END

-With the examples above, execution times become 48 & 92 seconds.
-We can also insert "DET" in the same program: it saves one data register:

Data Registers:  R00: temp                                      ( Registers R01 thru Rn+1 are to be initialized before executing "DISC" )

• R01 = an   • R02 = an-1  ......   • Rn+1 = a0

Flags:     CF 00 = a pivot p is regarded as zero if  | p | < 10-7
SF 00 =  a pivot p is regarded as zero if    p = 0

Subroutines: /

-Line 176 is a three-byte GTO 07

 01 LBL "DISC"  02 STO 00  03 STO N  04 ST+ X  05 DSE X  06 STO M  07 X^2  08 STO O  09 LBL 01  10 RCL 00            11 1  12 SIGN  13 LBL 02  14 CLX  15 RCL IND Y  16 LASTX  17 *  18 STO IND O  19 DSE O  20 ISG L  21 CLX  22 DSE Y  23 GTO 02  24 DSE N  25 FS? 30  26 GTO 00 27 RCL 00  28 SIGN  29 CLX  30 LBL 03  31 STO IND O  32 DSE O  33 DSE L  34 GTO 03  35 GTO 01  36 LBL 00  37 RCL 00            38 DSE X  39 STO N  40 LBL 04  41 RCL 00  42 1  43 +  44 SIGN  45 DSE N  46 X=0?  47 GTO 00  48 LBL 05  49 RCL IND L  50 STO IND O  51 DSE O  52 DSE L 53 GTO 05  54 RCL 00  55 DSE X  56 SIGN  57 CLX  58 LBL 06  59 STO IND O  60 DSE O  61 DSE L  62 GTO 06  63 GTO 04  64 LBL 00  65 1  66 CHS  67 RCL 00            68 X^2  69 LASTX  70 -  71 2  72 /  73 Y^X  74 RCL 01  75 /  76 STO 00  77 X<> M  78 STO Y 79 .1  80 %  81 +  82 STO N  83 FRC  84 STO O  85 *  86 1  87 ST+ O  88 LASTX  89 %  90 +  91 +  92 LBL 07  93 STO M            94 INT  95 RCL O  96 FRC  97 +  98 ENTER  99 ENTER 100 CLX 101 LBL 08 102 RCL IND Z 103 ABS 104 X>Y? 105 STO Z 106 X>Y? 107 + 108 RDN 109 ISG Z 110 GTO 08 111 RCL M 112 ENTER 113 FRC 114 R^ 115 INT 116 + 117 X=Y? 118 GTO 00 119 LBL 09 120 RCL IND X 121 X<> IND Z 122 STO IND Y 123 ISG Y 124 RDN 125 ISG Y 126 GTO 09 127 RCL 00         128 CHS 129 STO 00 130 LBL 00 131 CLX 132 FC? 00 133  E-7 134 RCL IND M 135 ST* 00 136 ABS 137 X<=Y? 138 GTO 12         139 ISG O 140 X<0? 141 GTO 13 142 RCL O 143 STO P 144 LBL 10 145 RCL M 146 ENTER 147 FRC 148 RCL P 149 INT 150 + 151 RCL IND X 152 RCL IND Z 153 / 154 SIGN 155 ISG Y 156 CLX 157 ISG Z 158 LBL 11 159 CLX 160 RCL IND Z 161 LASTX 162 * 163 ST- IND Y 164 ISG Y 165 CLX 166 ISG Z 167 GTO 11 168 ISG P 169 GTO 10 170 RCL N 171 ST+ O 172 SIGN 173 RCL M           174 + 175 ISG X 176 GTO 07 177 LBL 12 178 CLX 179 STO 00 180 LBL 13 181 RCL 00 182 CLA 183 END

( 275 bytes / SIZE (2n-1)2+1 )

 STACK INPUT OUTPUT X n D

Examples:

•  2 x3 - 5 x2 - 7 x + 1

CF 00

2   STO 01
-5  STO 02
-7  STO 03
1  STO 04

3  XEQ "DISC"   >>>>    D = 5621               ---Execution time = 35s---

•  x4 - 13 x3 + 60 x2 - 112 x + 64

1    STO 01
-13   STO 02
60   STO 03
-112  STO 04
64   STO 05

4    R/S    >>>>    D = 0                                   ---Execution time = 74s---

•  2 x7+ 3 x6 + 4 x5 + 5 x4 + 6 x3 + 7 x2 + 8 x + 9

2  STO 01  3  STO 02  4  STO 03  5  STO 04  6  STO 05  7  STO 06  8  STO 07  9  STO 08

7  R/S   >>>>   D =  -6.302249843  1012           ---Execution time = 5mn59s---

Notes:

-Since there are at most 319 data registers, the maximum degree is 9
-But we can modify this program if you want to use it with free42:

Data Registers:  R00.... R09: temp            R10 = D                          ( Registers R11 thru Rn+11 are to be initialized before executing "DISC" )

• R11 = an   • R12 = an-1  ......   • Rn+11 = a0

Flags:  /
Subroutines: /

-Line 207 is a three-byte GTO 07

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

( 302 bytes / SIZE (2n-1)2+11 )

 STACK INPUT OUTPUT X n D = R10

Examples:    With free42,  line 142 should be replaced by a smaller number ( perhaps E-20 )

•  2 x7+ 3 x6 + 4 x5 + 5 x4 + 6 x3 + 7 x2 + 8 x + 9

2  STO 11  3  STO 12  4  STO 13  5  STO 14  6  STO 15  7  STO 16  8  STO 17  9  STO 18

7  XEQ "DISC"   >>>>   D =  -6302249844736

•  2 x10+ 3 x9 + 4 x8 + 5 x7 + 6 x6 + 7 x5 + 8 x4 + 9 x3 + 10 x2 + 11 x + 12

2  STO 11  3  STO 12  4  STO 13  5  STO 14  6  STO 15  7  STO 16  8  STO 17  9  STO 18  10 STO 19  11 STO 20  12 STO 21

10   R/S   >>>>   D = -3727403657953361647360

Notes:

-I've rounded the results above: in the 2nd example, X-output = -3727403657953361647359,999999999999
-With free42, we can can compute the discriminant of polynomials of degree much larger than 9...

Reference:

[1]  https://en.wikipedia.org/w/index.php?title=Discriminant&oldid=1057323016