sw/src/pascal/BAIRST.PS

   1 PROGRAM 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
  51 BEGIN 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])
  87 END.