Just to keep you all updated, the SandMatrix module is close to release. While preparing the manual I’ve been tweaking and improving some of the routines, taking advantage of its extended function set. One of these has been CHRPOL, used to obtain the characteristic polynomial of a matrix.

The new version code is about 1/3rd. of the original one, and does the job in about half the time – if not faster. Such is the improvement derived from the “direct approach”. An image is worth 1,000 words, so here’s the program listing for you.- Stay tuned, the real think will be packed…

** Magic in 60 program steps.**

Uses the Faddevv-Leverrier method - you can almost read the algorithm from the ALPHA register, directly using the abstraction layer of the Matrix functions.

1 LBL "CHRPOL" MNAME in Alpha

2 DIM? n,00n

3 I#J? not square?

4 -ADV MATRX show error

5 ASTO 01 MNAME

6 -CCD MATRIX shows 'RUNNING…"

7 "|-,P" requires auxiliary array

8 MAT= P = A

9 ASWAP swap alpha around comma

10 DIM? n,00n

11 INT n

12 E

13 + n+1

14 MDET independent term

15 STO IND Y stored in Rn+1

16 ASWAP swap again

17 MAT= avoids LU issues

18 DIM?

19 "#" another auxiliary array

20 MATDIM

21 FRC 0,00n

22 2

23 + 2,00n

24 STO 00

25 CF 21 not halting VIEW

26 LBL 00

27 VIEW 00 shows current index

28 "#"

29 MIDN [#] = [I]

30 "P"

31 MTRACE tr(P)

32 RCL 00

33 INT k+1

34 E

35 - k

36 /

37 CHS

38 STO IND 00 pk = -tr(P) / k

39 "X,#,#"

40 MAT* [#] = pk [I]

41 "P,#,#"

42 MAT+ [#] = [P] + p[ I]

43 CLA

44 ARCL 01

45 "|-,#,P"

46 M*M B= A (B - p I)

47 ISG 00 increment counter

48 GTO 00 next coefficient

49 DIM? n,00n

50 FRC 0,00n

51 E

52 STO 01 it's monic (!)

53 E3/E+ 1.001

54 + 1.00(n+1) - cnt'l word

55 "#"

56 PURFL remove auxiliary arrays

57 "P"

58 PURFL

59 MNAME? bring MNAME back

60 END done

The polynomial control word is left in X upon termination - that is, registers containing the n+1 coefficients in the form bbb.eee -from Rbbb to Reee.

Cheers,

ÁM

PS. Here’s an article that may be useful to jog your memory:

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=944

*Edited: 19 Aug 2013, 9:10 a.m. *