Add README.md
[pdp8.git] / sw / src / pascal / BAIRST.PS
... / ...
CommitLineData
1PROGRAM BAIRSTOW(INPUT,OUTPUT);
2
3 CONST NMAX=20; EPSILON=1E-8;
4
5 TYPE KOMPLEX = RECORD RE,IM: REAL END;
6
7 VAR N,NR,I: INTEGER;
8 A: ARRAY[0..NMAX] OF REAL;
9 B,C: ARRAY[-1..NMAX] OF REAL;
10 P,Q,PDELTA,QDELTA,DISKR: REAL;
11
12 PROCEDURE DRUCKELOESUNG(X: KOMPLEX);
13 BEGIN
14 NR := NR + 1;
15 WRITE("X", NR:1, "=":2, X.RE :11:5);
16 IF X.IM<>0 THEN
17 BEGIN IF X.IM>0 THEN WRITE("+":2)
18 ELSE WRITE("-":2);
19 WRITE(ABS(X.IM):11:5, " * J")
20 END;
21 WRITELN
22 END (* DRUCKELOESUNG *);
23
24
25 PROCEDURE QUADRATFAKTOR(P,Q: REAL);
26 VAR D: REAL;
27 X: KOMPLEX;
28 BEGIN
29 D := P*P/4 - Q;
30 IF D>=0 THEN
31 BEGIN X.RE := -P/2 + SQRT(D); X.IM := 0;
32 DRUCKELOESUNG(X);
33 X.RE := -P/2 - SQRT(D);
34 DRUCKELOESUNG(X)
35 END ELSE
36 BEGIN X.RE := -P/2; X.IM := SQRT(-D);
37 DRUCKELOESUNG(X);
38 X.IM := -X.IM;
39 DRUCKELOESUNG(X)
40 END
41 END (* QUADRATFAKTOR *);
42
43
44 PROCEDURE LINEARFAKTOR(Q: REAL);
45 VAR X: KOMPLEX;
46 BEGIN X.RE := -Q; X.IM := 0;
47 DRUCKELOESUNG(X)
48 END (* LINEARFAKTOR *);
49
50
51BEGIN READ(N);
52 WRITELN;
53 WRITELN("P O L Y N O M", N:3, ". GRADES:");
54 FOR I := 0 TO N DO
55 BEGIN READ(A[I]);
56 WRITELN(A[I]:8:2, " X^", N-I :1)
57 END;
58 WRITELN;
59 IF A[0]<>1 THEN BEGIN Q := A[0];
60 A[0] := 1;
61 FOR I := 1 TO N DO A[I] := A[I]/Q;
62 END;
63
64 NR := 0;
65
66 WHILE N>2 DO
67 BEGIN P := -1; Q := 1;
68 REPEAT
69 B[-1] := 0; B[0] := 1;
70 C[-1] := 0; C[0] := 1;
71 FOR I := 1 TO N DO
72 BEGIN
73 B[I] := A[I] - P*B[I-1] - Q*B[I-2];
74 C[I] := B[I] - P*C[I-1] - Q*C[I-2]
75 END;
76 DISKR := SQR(C[N-2]) - (C[N-1]-B[N-1])*C[N-3];
77 PDELTA:= (B[N-1]*C[N-2]-B[N]*C[N-3])/DISKR;
78 QDELTA:= (C[N-2]*B[N]-(C[N-1]-B[N-1])*B[N-1])/DISKR;
79 P := P + PDELTA; Q := Q + QDELTA
80 UNTIL ABS(PDELTA)+ABS(QDELTA) < EPSILON;
81 QUADRATFAKTOR(P,Q);
82 N := N-2;
83 FOR I := 1 TO N DO A[I] := B[I];
84 END;
85 IF N=2 THEN QUADRATFAKTOR(A[1],A[2])
86 ELSE LINEARFAKTOR(A[1])
87END.