|
|
!' アーバス法(Aberth法)
OPTION ARITHMETIC COMPLEX
OPTION BASE 0
PUBLIC NUMERIC N,EPS
LET N=10 !'次数
LET EPS=1E-8
DIM Z(N),A(N)
LET A(N)=1 !' A(N)*X^N+A(N-1)*X^(N-1)+...A(1)*X+A(0)=0
LET A(0)=-1
CALL HORNER(A,-A(N-1)/(N*A(N))/A(N),Y,DY)
FOR I=1 TO N
LET Z(I)=-A(N-1)/(N*A(N))+ABS(Y)^(1/N)*EXP(SQR(-1)*2*PI*(I-1)/N+3/(2*N))
NEXT I
DO
FOR I=1 TO N
LET S=0
FOR J=1 TO N
IF I<>J THEN LET S=S+1/(Z(I)-Z(J))
NEXT J
CALL HORNER(A,Z(I),Y,DY)
LET Z(I)=Z(I)-Y/(DY-Y*S)
NEXT I
LET FL=0
FOR I=1 TO N
CALL HORNER(A,Z(I),Y,DY)
IF ABS(RE(Y))>EPS OR ABS(IM(Y))>EPS THEN
LET FL=1
EXIT FOR
END IF
NEXT I
LOOP UNTIL FL=0
FOR I=1 TO N
PRINT "ANSWER";I;":";
CALL DISPLAY(Z(I))
NEXT I
END
EXTERNAL SUB DISPLAY(Z)
OPTION ARITHMETIC COMPLEX
IF ABS(RE(Z))>EPS THEN PRINT RE((Z));
IF ABS(IM(Z))>EPS THEN
IF IM(Z)<0 THEN
PRINT "-";
ELSE
IF ABS(RE(Z))>EPS THEN PRINT "+";
END IF
PRINT ABS(IM(Z));"i";
END IF
PRINT
END SUB
EXTERNAL SUB HORNER(A(),X,Y,DY)
OPTION ARITHMETIC COMPLEX
LET Y=A(N)
LET DY=Y
FOR I=N-1 TO 0 STEP -1
LET Y=Y*X+A(I)
IF I>0 THEN LET DY=DY*X+Y
NEXT I
END SUB
|
|