Schur-Cohn Algorithm

[Eddie W. Shore]

Source:

Henrici, Peter. Computational Analysis With the HP-25 Calculator A Wiley-Interscience Publication. John Wiley & Sons: New York 1977 . ISBN 0-471-02938-6

Introduction:

The Schur-Cohn Algorithm tests whether the roots of a polynomial p(x) lies with in the unit circle. That is for the polynomial p(x):

p(x) = a_0 * x^n + a_1 * x^(n - 1) + a_2 * x^(n - 2) + ... + a_n

For all the roots of p(x), r_0, r_1, r_2, ... , r_n , |r_k| < 1. The test covers both real and complex roots. Keep in mind that this test doesn't tell us what the roots are, only a clue to whether the roots lie in the unit circle or not.

The Schur-Cohn Algorithm returns the test answers b_0, b_1, etc. If |b_k| < 1 holds true for each b_k for k = 0 to n, then we can conclude that |r_k| < 1.

The following program is adopted from Peter Henrici's book Computational Analysis With the HP-25 Pocket Calculator [Henrici, 111] for the HP 41C and the Swiss Micros DM 41L. The program deals with polynomials up to the fourth order, but can be expanded on an RPN calculator with more program space and memory registers.

Instructions:

For the fourth order polynomial:

p(x) = a_0 * x^4 + a_1 * x^3 + a_2 * x^2 + a_1 * x + a_4

To the coefficients in the following registers:

R00 = a_0
R01 = a_1
R02 = a_2
R03 = a_3
R04 = a_4

The test results will pop momentarily during execution. The program is done when after four test values, 4 is in the display. The last test value is stored in R05.

DM 41L and HP 41C Program: SCHUR (77 bytes, 11 registers of memory)

Code:


01  LBL^T SCHUR
02  CLX
03  STO 07
04  LBL 03
05  STO 06
06  RCL 04
07  RCL 00 
08  /
09  STO 05    
10 PSE   // show test values
11 RCL 04
12  *
13 ST- 00
14 RCL 01
15 RCL 02
16 RCL 03
17 LBL 15
18 RCL 06
19 X=0?
20 GTO 24
21 RDN   // R↓
22 RCL 04
23 1
24 ST- 06
25 RDN
26 GTO 15
27 LBL 24
28 RDN
29 RCL 05
30 *
31 ST- 01
32 RDN
33 RCL 05
34 *
35 ST- 02
36 RDN
37 RCL 05
38  *
39 ST- 03
40 RCL 03
41 STO 04
42 RCL 02
43 STO 03
44 RCL 01
45 STO 02
46 RCL 00
47 STO 01
48 1
49 ST+ 07
50 4
51 RCL 07
52 X<Y?
53 GTO 03

Link to blog post entry: https://edspi31415.blogspot.com/2018/12/...rithm.html