--- a
+++ b/SpezFunkt1.mod
@@ -0,0 +1,304 @@
+IMPLEMENTATION MODULE SpezFunkt1;
+
+ (*------------------------------------------------------------------------*)
+ (* Stellt einige "Sonderfunktionen" fuer Polynome und Nullstellen- *)
+ (* berechnungen bereit. *)
+ (* Provides some special procedure for the evaluation of polynominals and *)
+ (* for the calculation of roots. *)
+ (*------------------------------------------------------------------------*)
+ (* Letzte Bearbeitung: *)
+ (* *)
+ (* 17.07.93, MRi: Durchsicht *)
+ (* 07.10.15, MRi: Umbenennen von Funkt nach SpezFunkt1 *)
+ (* 02.01.16, MRi: Prozedur EvalCheb1Coeff eingefuegt. *)
+ (* 04.01.16, MRi: Prozedur ChebPoly1Sum eingefuegt *)
+ (* 29.06.16, MRi: Prozedur Nullstellen erweitert *)
+ (* 29.06.16, MRi: Prozedur Regula um Test auf Ueberlauf erweitert *)
+ (* 02.10.16, MRi: Prozedur QuadGl und KubGl aus ISO COMPLEX Typ umge- *)
+ (* stellt *)
+ (* 19.09.17, MRi: Korrektur in QuadGl fuer X1re = 0 *)
+ (*------------------------------------------------------------------------*)
+ (* Implementation : Michael Riedl *)
+ (* Licence : GNU Lesser General Public License (LGPL) *)
+ (*------------------------------------------------------------------------*)
+
+ (* $Id: SpezFunkt1.mod,v 1.6 2016/10/12 09:22:26 mriedl Exp mriedl $ *)
+
+FROM Errors IMPORT Fehler,Fehlerflag,ErrOut;
+FROM LongMath IMPORT pi,sqrt,cos,arccos,exp,ln;
+FROM LongComplexMath IMPORT zero,scalarMult,conj;
+ IMPORT LongComplexMath;
+FROM LMathLib IMPORT MachEps,Small,sqrt3,
+ sign2,sqr,CardPot,MinCard,Funktion1;
+
+PROCEDURE Horner(VAR x : LONGREAL;
+ VAR Koeff : ARRAY OF LONGREAL;
+ n : CARDINAL) : LONGREAL;
+
+ VAR y : LONGREAL;
+ i : CARDINAL;
+
+BEGIN
+ y:=Koeff[n]; FOR i:=n-1 TO 0 BY -1 DO y:=y*x + Koeff[i]; END;
+ RETURN y;
+END Horner;
+
+PROCEDURE Clenshaw(x : LONGREAL; (* Funktionswert *)
+ n : CARDINAL; (* Grad d. T-Entwicklung *)
+ VAR C : ARRAY OF LONGREAL;
+ U,O : LONGREAL) : LONGREAL;
+
+ VAR i : CARDINAL;
+ y,y2,D2,D3,Zw : LONGREAL;
+
+BEGIN
+ y := (2.0*x - O - U) / (O - U); y2 := 2.0*y;
+ D2:=0.0; D3:=0.0;
+ FOR i:=n TO 1 BY -1 DO Zw:=D2; D2:=y2*D2 - D3 + C[i]; D3:=Zw; END;
+ RETURN y*D2 - D3 + 0.5*C[0];
+END Clenshaw;
+
+PROCEDURE EvalCheb1Coeff( x : LONGREAL;
+ VAR T : ARRAY OF LONGREAL;
+ n : CARDINAL);
+
+ VAR i : CARDINAL;
+BEGIN
+ T[0]:=1.0;
+ IF (n >= 1) THEN T[1]:=x; END;
+ IF (n >= 2) THEN
+ FOR i:=2 TO n DO
+ T[i]:=2.0*x*T[i-1] - T[i-2];
+ END;
+ END;
+END EvalCheb1Coeff;
+
+PROCEDURE ChebPoly1Sum( n : CARDINAL;
+ x : LONGREAL;
+ VAR A : ARRAY OF LONGREAL) : LONGREAL;
+
+ VAR k : CARDINAL;
+ h,r,s,tx : LONGREAL;
+BEGIN
+ IF (n = 0) THEN
+ RETURN A[0];
+ ELSIF (n = 1) THEN
+ RETURN A[0] + A[1]*x;
+ ELSE
+ tx := x + x;
+ r := A[n];
+ h := A[n-1] + r*tx;
+ FOR k:=n-2 TO 1 BY -1 DO
+ s := r;
+ r := h;
+ h := A[k] + r*tx - s;
+ END;
+ RETURN A[0] - r + h*x;
+ END;
+END ChebPoly1Sum;
+
+PROCEDURE RatPoly( x : LONGREAL;
+ VAR A,B : ARRAY OF LONGREAL;
+ n,m : CARDINAL) : LONGREAL;
+
+ VAR i : CARDINAL;
+ Zaehler,Nenner : LONGREAL;
+
+BEGIN
+ Zaehler:=A[n];
+ FOR i:=n-1 TO 0 BY -1 DO Zaehler := x*Zaehler + A[i]; END;
+ Nenner :=B[m];
+ FOR i:=m-1 TO 0 BY -1 DO Nenner := x*Nenner + B[i]; END;
+ IF (ABS(Nenner) > 1.0E-34) THEN
+ RETURN Zaehler / Nenner;
+ ELSE
+ RETURN MAX(LONGREAL);
+ END;
+END RatPoly;
+
+PROCEDURE QuadGl( A,B,C : LONGREAL;
+ VAR X1,X2 : LONGCOMPLEX);
+
+ VAR Nenner,Dis : LONGREAL;
+ X1re,X1im,X2re,X2im : LONGREAL;
+BEGIN
+ Dis := B*B - 4.0*A*C;
+ Nenner := 2.0*A;
+ X1re := - B / Nenner;
+ X2re := X1re;
+ IF (Dis < 0.0) THEN (* Komplexe Wurzeln *)
+ X1im := - sqrt(-Dis) / Nenner;
+ X2im := - X1im;
+ ELSE (* Reelle Wurzeln *)
+ X1im := 0.0;
+ X2im := 0.0;
+ IF (B < 0.0) THEN
+ X1re := X1re + sqrt(Dis) / Nenner;
+ ELSE
+ X1re := X1re - sqrt(Dis) / Nenner;
+ END;
+ IF (X1re = 0.0) THEN
+ X2re:=0.0; (* ??? *)
+ ELSE
+ X2re := C / (A*X1re);
+ END;
+ END;
+ X1 := CMPLX(X1re,X1im);
+ X2 := CMPLX(X2re,X2im);
+END QuadGl;
+
+PROCEDURE KubGl( A,B,C,D : LONGREAL;
+ VAR X1,X2,X3 : LONGCOMPLEX);
+
+ PROCEDURE Root3(x : LONGREAL) : LONGREAL;
+
+ VAR sign : LONGREAL;
+ BEGIN
+ sign:=1.0; IF (x < 0.0) THEN sign:=-1.0; END;
+ RETURN sign*exp(ln(ABS(x))/3.0);
+ END Root3;
+
+ CONST sqrt3 = 1.7320508075688772;
+ W1 = CMPLX(-0.5, 0.5*sqrt3);
+ W2 = CMPLX(-0.5,-0.5*sqrt3);
+
+ VAR P,Q,Z,T,Phi : LONGREAL;
+ Z1,Z2,Cp,Px,A0,A3 : LONGREAL;
+ V,U : LONGCOMPLEX;
+ Y1,Y2,Y3 : LONGCOMPLEX;
+BEGIN
+ A0 := A; A := B/A0; B := C/A0; C := D/A0;
+
+ P := - sqr(A)/9.0 + B/3.0;
+ Q := A*(sqr(A)/27.0 - B/6.0) + C/2.0;
+
+ Z := sqr(Q) + CardPot(P,3);
+ A3 := A / 3.0;
+
+ IF (ABS(Z) < 1.0E-15) THEN
+ T := Root3(-Q);
+ Y1 := CMPLX(2.0*T - 0.0, 0.0);
+ Y2 := CMPLX( - T - 0.0, 0.0);
+ Y3 := CMPLX( - T - 0.0, 0.0);
+ ELSIF (Z > 0.0) THEN
+ Z1 := - Q + sqrt(Z);
+ Z2 := - Q - sqrt(Z);
+ U := CMPLX(Root3(Z1),0.0);
+ V := CMPLX(Root3(Z2),0.0);
+ Y1 := U + V;
+ Y2 := W1*U + W2*V;
+ Y3 := W2*U + W1*V;
+ ELSIF (Z < 0.0) THEN
+ Cp := - Q / sqrt(-CardPot(P,3));
+ Phi := arccos(Cp);
+ Px :=-2.0*sqrt(-P);
+ Y1 := CMPLX(-Px*cos( Phi /3.0), 0.0);
+ Y2 := CMPLX( Px*cos((Phi + pi)/3.0), 0.0);
+ Y3 := CMPLX( Px*cos((Phi - pi)/3.0), 0.0);
+ ELSE (* Unfug wg. Compilerwarnung *)
+ Y1 := zero; Y2 := zero; Y3 := zero;
+ END;
+ X1 := Y1 - CMPLX(A3,0.0);
+ X2 := Y2 - CMPLX(A3,0.0);
+ X3 := Y3 - CMPLX(A3,0.0);
+END KubGl;
+
+PROCEDURE CubicEq( A,B,C,D : LONGCOMPLEX; (* coeffs of the polynomial *)
+ VAR X1,X2,X3 : LONGCOMPLEX);
+
+ CONST Zeta1 = CMPLX(-0.5, 0.5*sqrt3);
+ Zeta2 = CMPLX(-0.5, -0.5*sqrt3);
+ third = 1.0 / 3.0;
+
+ VAR P,Q,P2,A2,A0,s1,s2,delta : LONGCOMPLEX;
+BEGIN
+ A0 := A; A := - B/A0; B := C/A0; C := - D/A0;
+
+ A2 := A*A;
+ P := A*(scalarMult(2.0,A2) - scalarMult(9.0,B) + scalarMult(27.0,C));
+ Q := A2 - scalarMult(3.0,B);
+ P2 := P*P;
+ delta := LongComplexMath.sqrt(P2 - scalarMult(4.0,(Q*Q*Q)));
+ IF (RE(conj(P)*delta) >= 0.0) THEN
+ s1 := LongComplexMath.power(scalarMult(0.5,(P + delta)),third);
+ ELSE
+ s1 := LongComplexMath.power(scalarMult(0.5,(P - delta)),third);
+ END;
+ IF (s1 = zero) THEN
+ s2 := zero;
+ ELSE
+ s2 := Q / s1;
+ END;
+
+ X1 := scalarMult(third, A + s1 + s2);
+ X2 := scalarMult(third, A + s1*Zeta2 + s2*Zeta1);
+ X3 := scalarMult(third, A + s1*Zeta1 + s2*Zeta2);
+END CubicEq;
+
+PROCEDURE Regula(VAR Funkt : Funktion1; (* y:=Funkt(x), mit x,y : LONGREAL *)
+ Xu,Xo : LONGREAL; (* Startwerte fuer die Nullstelle *)
+ fXu,fXo : LONGREAL; (* Funktionswerte von Xu,Xo. *)
+ VAR X0 : LONGREAL; (* Nullstelle von Funkt *)
+ genau : LONGREAL); (* Geforderte Genauigkeit *)
+
+ CONST MaxIter = 32;
+
+ VAR Iter : CARDINAL;
+ Nenner : LONGREAL;
+BEGIN
+ Iter:=0; Fehler:=FALSE;
+ LOOP
+ INC(Iter);
+ IF (Iter > MaxIter) THEN
+ Fehler:=TRUE; ErrOut('Iter > Maxiter (Regula).');
+ RETURN;
+ END;
+ Nenner := (fXo - fXu);
+ IF (ABS(Nenner) < Small) THEN
+ Nenner := sign2(Small,Nenner); (* resultiert in grosser Zahl *)
+ END;
+ X0 := Xo - fXo*(Xo - Xu) / Nenner;
+ IF ((ABS(X0 - Xo)) < genau*ABS(X0)) THEN RETURN END;
+ fXu:=fXo; fXo:=Funkt(X0);
+ Xu:=Xo; Xo:=X0;
+ END;
+END Regula;
+
+PROCEDURE NullStellen(u,o : LONGREAL; (* Intervallgrenzen *)
+ dH : LONGREAL; (* Schrittweite *)
+ f : Funktion1;
+ VAR X0 : ARRAY OF LONGREAL; (* Nullstellen von f *)
+ VAR nx0 : CARDINAL; (* Anzahl der Nullstellen *)
+ i : CARDINAL; (* Anzahl zu suchender Nullst. *)
+ genau : LONGREAL);
+
+ VAR xu,xo,fxu,fxo : LONGREAL;
+ max : CARDINAL;
+
+BEGIN
+ IF (i = 0) THEN max:=HIGH(X0)+1 ELSE max:=MinCard(HIGH(X0)+1,i) END;
+ fxo:=f(u); xo:=u; nx0:=0;
+ IF (ABS(fxo) < genau) THEN
+ X0[0]:=u; INC(nx0);
+ IF (max = 1) THEN RETURN; END;
+ END;
+ REPEAT
+ xu:=xo; xo:=xo + dH;
+ fxu:=fxo; fxo:=f(xo);
+ IF (fxo*fxu < 0.0) THEN (* Vorzeichenwechsel *)
+ Regula(f,xu,xo,fxu,fxo,X0[nx0],genau);
+ INC(nx0);
+ ELSIF ((ABS(fxo)-MachEps) < genau) THEN (* Funktionswert ist Null *)
+ X0[nx0]:=xo;
+ INC(nx0);
+ xo:=xo + dH;
+ END;
+ IF (nx0 = max) THEN RETURN; END;
+ UNTIL (xo >= o);
+ Fehler := (nx0 = 0);
+ IF Fehler THEN
+ Fehlerflag:='Keine Nullstellen im Intervall.';
+ END;
+END NullStellen;
+
+END SpezFunkt1.