code 34601;
procedure QZI(N,A,B,X,ALFR,ALFI,BETA,ITER,EM);
value N;integer N;array A,B,X,ALFR,ALFI,BETA,EM;
integer array ITER;
begin real DWARF,EPS,EPSA,EPSB;
procedure QZIT(N,A,B,X,EPS,EPSA,EPSB,ITER);value N,EPS;
real EPS,EPSA,EPSB;integer N;integer array ITER;array A,B,X;
begin real ANORM,BNORM,ANI,BNI,CONST,A10,A20,A30,B11,
    B22,B33,B44,A11,A12,A21,A22,A33,A34,A43,A44,B12,B34,OLD1,OLD2;
    integer I,Q,M,M1,Q1,J,K,K1,K2,K3,KM1;boolean STATIONARY;
    ANORM:=BNORM:=0;for I:=1 step 1 until N do 
    begin BNI:=0;ITER[I]:=0;ANI:=if I>1then ABS(A[I,I-1])else 0;
        for J:=I step 1 until N do 
        begin ANI:=ANI+ABS(A[I,J]);BNI:=BNI+ABS(B[I,J])
        end;if ANI>ANORM then ANORM:=ANI;if BNI>BNORM"THEN"
        BNORM:=BNI
    end;if ANORM=0 then ANORM:=EPS;if BNORM=0 then BNORM:=EPS;
    EPSA:=EPS*ANORM;EPSB:=EPS*BNORM;
    for M:=N,M while M>=3 do 
  begin 
    for I:=M+1,I-1 while (if I>1 then ABS(A[I,I-1])>EPSA else 
    false ) do Q:=I-1;
    if Q>1 then A[Q,Q-1]:=0;
 L: if Q>=M-1 then M:=Q-1 else 
    begin 
        if ABS(B[Q,Q])<=EPSB then 
        begin B[Q,Q]:=0;Q1:=Q+1;
            HSH2COL(Q,Q,N,Q,A[Q,Q],A[Q1,Q],A,B);A[Q1,Q]:=0;
            Q:=Q1;goto L
        end else M1:=M-1;Q1:=Q+1;CONST:=0.75;ITER[M]:=ITER[M]+1;
        STATIONARY:=if ITER[M]=1 then true else 
        ABS(A[M,M-1])>=CONST*OLD1"AND"ABS(A[M-1,M-2])>=CONST*OLD2;
        if ITER[M]>30and STATIONARY then 
        begin for I:=1 step 1 until M do ITER[I]:=-1;
            goto OUT
        end;
        if ITER[M]=10and STATIONARY then 
        begin A10:=0;A20:=1;A30:=1.1605
        end else 
        begin B11:=B[Q,Q];B22:=if ABS(B[Q1,Q1])<EPSB then EPSB
            else B[Q1,Q1];
            B33:=if ABS(B[M1,M1])<EPSB then EPSB else B[M1,M1];
            B44:=if ABS(B[M,M])<EPSB then EPSB else B[M,M] ;
            A11:=A[Q,Q]/B11;A12:=A[Q,Q1]/B22;A21:=A[Q1,Q]/B11;
            A22:=A[Q1,Q1]/B22;A33:=A[M1,M1]/B33;A34:=A[M1,M]/B44;
            A43:=A[M,M1]/B33;A44:=A[M,M]/B44;B12:=B[Q,Q1]/B22;
            B34:=B[M1,M]/B44;
            A10:=((A33-A11)*(A44-A11)-A34*A43+A43*B34*A11)/A21
                 +A12-A11*B12;
            A20:=(A22-A11-A21*B12)-(A33-A11)-(A44-A11)+A43*B34;
            A30:=A[Q+2,Q1]/B22
        end;OLD1:=ABS(A[M,M-1]);OLD2:=ABS(A[M-1,M-2]);
        for K:=Q step 1 until M1 do 
        begin K1:=K+1;K2:=K+2;K3:=if K+3>M then M else K+3;
            KM1:=if K-1<Q then Q else K-1;
            if K^=M1 then 
            begin if K=Q then 
                HSH3COL(KM1,KM1,N,K,A10,A20,A30,A,B) else 
                begin 
                 HSH3COL(KM1,KM1,N,K,A[K,KM1],A[K1,KM1],A[K2,KM1],A,B);
                  A[K1,KM1]:=A[K2,KM1]:=0
                end;
                HSH3ROW3(1,K3,N,K,B[K2,K2],B[K2,K1],B[K2,K],A,B,X);
                B[K2,K]:=B[K2,K1]:=0 ;
            end else 
            begin HSH2COL(KM1,KM1,N,K,A[K,KM1],A[K1,KM1],A,B);
                A[K1,KM1]:=0
            end;
            HSH2ROW3(1,K3,K3,N,K,B[K1,K1],B[K1,K],A,B,X);B[K1,K]:=0
        end 
    end 
  end; OUT:
end QZIT;
procedure QZVAL(N,A,B,X,EPSA,EPSB,ALFR,ALFI,BETA);value N;
real EPSA,EPSB;integer N;array ALFR,ALFI,BETA,A,B,X;
begin integer M,L,J;real AN,BN,A11,A12,A21,A22,B11,B12,B22,E,C,D,
    ER,EI,A11R,A11I,A12R,A12I,A21R,A21I,A22R,A22I,CZ,SZR,SZI,
    CQ,SQR,SQI,SSR,SSI,TR,TI,BDR,BDI,R;
    for M:=N,M while M>0 do 
    if (if M>1 then A[M,M-1]=0 else true )then 
    begin ALFR[M]:=A[M,M];BETA[M]:=B[M,M];ALFI[M]:=0;M:=M-1
    end else 
    begin L:=M-1;if ABS(B[L,L])<=EPSB then 
        begin B[L,L]:=0;HSH2COL(L,L,N,L,A[L,L],A[M,L],A,B);
            A[M,L]:=B[M,L]:=0;ALFR[L]:=A[L,L];ALFR[M]:=A[M,M];
            BETA[L]:=B[L,L];BETA[M]:=B[M,M];ALFI[M]:=ALFI[L]:=0;
        end else if ABS(B[M,M])<=EPSB then 
        begin B[M,M]:=0;HSH2ROW3(1,M,M,N,L,A[M,M],A[M,L],A,B,X);
            A[M,L]:=B[M,L]:=0;ALFR[L]:=A[L,L];ALFR[M]:=A[M,M];
            BETA[L]:=B[L,L];BETA[M]:=B[M,M];ALFI[M]:=ALFI[L]:=0;
        end else 
      begin 
        AN:=ABS(A[L,L])+ABS(A[L,M])+ABS(A[M,L])+ABS(A[M,M]);
        BN:=ABS(B[L,L])+ABS(B[L,M])+ABS(B[M,M]);
        A11:=A[L,L]/AN;A12:=A[L,M]/AN;A21:=A[M,L]/AN;A22:=A[M,M]/AN;
        B11:=B[L,L]/BN;B12:=B[L,M]/BN;B22:=B[M,M]/BN;
        E:=A11/B11;
        C:=((A22-E*B22)/B22-(A21*B12)/(B11*B22))/2;
        D:=C*C+(A21*(A12-E*B12))/(B11*B22);
        if D>=0 then 
        begin E:=E+(if C<0then C-SQRT(D)else C+SQRT(D));
            A11:=A11-E*B11;A12:=A12-E*B12;A22:=A22-E*B22;
            if ABS(A11)+ABS(A12)>=ABS(A21)+ABS(A22) then 
            HSH2ROW3(1,M,M,N,L,A12,A11,A,B,X)else 
            HSH2ROW3(1,M,M,N,L,A22,A21,A,B,X);
            if AN>=ABS(E)*BN"THEN"
            HSH2COL(L,L,N,L,B[L,L],B[M,L],A,B) else 
            HSH2COL(L,L,N,L,A[L,L],A[M,L],A,B);
            A[M,L]:=B[M,L]:=0;
            ALFR[L]:=A[L,L];ALFR[M]:=A[M,M];BETA[L]:=B[L,L];
            BETA[M]:=B[M,M];ALFI[M]:=ALFI[L]:=0;
        end else 
        begin 
            ER:=E+C;EI:=SQRT(-D);A11R:=A11-ER*B11;A11I:=EI*B11;
            A12R:=A12-ER*B12;A12I:=EI*B12;A21R:=A21;A21I:=0;
            A22R:=A22-ER*B22;A22I:=EI*B22;
            if ABS(A11R)+ABS(A11I)+ABS(A12R)+ABS(A12I)>=
            ABS(A21R)+ABS(A22R)+ABS(A22I)then 
            CHSH2(A12R,A12I,-A11R,-A11I,CZ,SZR,SZI)else 
            CHSH2(A22R,A22I,-A21R,-A21I,CZ,SZR,SZI);
            if AN>=(ABS(ER)+ABS(EI))*BN"THEN"
            CHSH2(CZ*B11+SZR*B12,SZI*B12,SZR*B22,SZI*B22,CQ,SQR,SQI)
            else CHSH2(CZ*A11+SZR*A12,SZI*A12,CZ*A21+SZR*A22,SZI*A22,
            CQ,SQR,SQI);SSR:=SQR*SZR+SQI*SZI;SSI:=SQR*SZI-SQI*SZR;
            TR:=CQ*CZ*A11+CQ*SZR*A12+SQR*CZ*A21+SSR*A22;
            TI:=CQ*SZI*A12-SQI*CZ*A21+SSI*A22;
            BDR:=CQ*CZ*B11+CQ*SZR*B12+SSR*B22;
            BDI:=CQ*SZI*B12+SSI*B22;
            R:=SQRT(BDR*BDR+BDI*BDI);BETA[L]:=BN*R;
            ALFR[L]:=AN*(TR*BDR+TI*BDI)/R;
            ALFI[L]:=AN*(TR*BDI-TI*BDR)/R;
            TR:=SSR*A11-SQR*CZ*A12-CQ*SZR*A21+CQ*CZ*A22;
            TI:=-SSI*A11-SQI*CZ*A12+CQ*SZI*A21;
            BDR:=SSR*B11-SQR*CZ*B12+CQ*CZ*B22;
            BDI:=-SSI*B11-SQI*CZ*B12;
            R:=SQRT(BDR*BDR+BDI*BDI);BETA[M]:=BN*R;
            ALFR[M]:=AN*(TR*BDR+TI*BDI)/R;
            ALFI[M]:=AN*(TR*BDI-TI*BDR)/R;
        end 
      end;M:=M-2
    end 
end QZVAL

procedure QZVEC(N,A,B,X,EPSA,EPSB,ALFR,ALFI,BETA);value N,EPSA,EPSB;
integer N;real EPSA,EPSB;array A,B,ALFR,ALFI,BETA,X;
begin integer M,MR,MI,L,L1,J,K;real BETM,ALFM,SL,SK,D,TKK,TKL,TLK,
    TLL,ALMI,ALMR,TR,TI,SLR,SLI,SKR,SKI,DR,DI,TKKR,TKKI,TKLR,TKLI,TLKR,
    TLKI,TLLR,TLLI,S,R;
    for M:=N step -1 until 1 do 
    if ALFI[M]=0 then 
    begin comment M-TH REAL VECTOR;
        ALFM:=ALFR[M];BETM:=BETA[M];B[M,M]:=1;L1:=M;
        for L:=M-1 step -1 until 1 do 
        begin SL:=0;
            for J:=L1 step 1 until M do 
            SL:=SL+(BETM*A[L,J]-ALFM*B[L,J])*B[J,M];
            if (if L^=1then BETM*A[L,L-1]=0else true )then 
            begin comment 1-1 BLOCK;
                D:=BETM*A[L,L]-ALFM*B[L,L];
                if D=0 then D:=(EPSA+EPSB)/2;B[L,M]:=-SL/D
            end else 
            begin comment 2-2 BLOCK;K:=L-1;SK:=0;
                for J:=L1 step 1 until M do 
                SK:=SK+(BETM*A[K,J]-ALFM*B[K,J])*B[J,M];
                TKK:=BETM*A[K,K]-ALFM*B[K,K];
                TKL:=BETM*A[K,L]-ALFM*B[K,L];
                TLK:=BETM*A[L,K];
                TLL:=BETM*A[L,L]-ALFM*B[L,L];
                D:=TKK*TLL-TKL*TLK;if D=0 then D:=(EPSA+EPSB)/2;
                B[L,M]:=(TLK*SK-TKK*SL)/D;
               B[K,M]:=if ABS(TKK)>=ABS(TLK)then -(SK+TKL*B[L,M])/TKK
                else -(SL+TLL*B[L,M])/TLK;L:=L-1
            end;L1:=L
        end 
    end else 
    begin comment COMPLEX VECTOR;
        ALMR:=ALFR[M-1];ALMI:=ALFI[M-1];BETM:=BETA[M-1];MR:=M-1;MI:=M;
        B[M-1,MR]:=ALMI*B[M,M]/(BETM*A[M,M-1]);
        B[M-1,MI]:=(BETM*A[M,M]-ALMR*B[M,M])/(BETM*A[M,M-1]);
        B[M,MR]:=0;B[M,MI]:=-1;L1:=M-1;
        for L:=M-2 step -1 until 1 do 
        begin SLR:=SLI:=0;
            for J:=L1 step 1 until M do 
            begin TR:=BETM*A[L,J]-ALMR*B[L,J];
                TI:=-ALMI*B[L,J];
                SLR:=SLR+TR*B[J,MR]-TI*B[J,MI];
                SLI:=SLI+TR*B[J,MI]+TI*B[J,MR]
            end;
            if (if L^=1then BETM*A[L,L-1]=0else true )then 
            begin DR:=BETM*A[L,L]-ALMR*B[L,L];
                DI:=-ALMI*B[L,L];
                COMDIV(-SLR,-SLI,DR,DI,B[L,MR],B[L,MI]);
            end else 
            begin K:=L-1;SKR:=SKI:=0;
                for J:=L1 step 1 until M do 
                begin TR:=BETM*A[K,J]-ALMR*B[K,J];
                    TI:=-ALMI*B[K,J];
                    SKR:=SKR+TR*B[J,MR]-TI*B[J,MI];
                    SKI:=SKI+TR*B[J,MI]+TI*B[J,MR]
                end;
                TKKR:=BETM*A[K,K]-ALMR*B[K,K];
                TKKI:=-ALMI*B[K,K];
                TKLR:=BETM*A[K,L]-ALMR*B[K,L];
                TKLI:=-ALMI*B[K,L];
                TLKR:=BETM*A[L,K];TLKI:=0;
                TLLR:=BETM*A[L,L]-ALMR*B[L,L];
                TLLI:=-ALMI*B[L,L];
                DR:=TKKR*TLLR-TKKI*TLLI-TKLR*TLKR;
                DI:=TKKR*TLLI+TKKI*TLLR-TKLI*TLKR;
                if DR=0and DI=0then DR:=(EPSA+EPSB)/2;
                COMDIV(TLKR*SKR-TKKR*SLR+TKKI*SLI,TLKR*SKI-TKKR*SLI-
                TKKI*SLR,DR,DI,B[L,MR],B[L,MI]);
                if ABS(TKKR)+ABS(TKKI)>=ABS(TLKR)then 
                COMDIV(-SKR-TKLR*B[L,MR]+TKLI*B[L,MI],-SKI-TKLR*B[L,MI]
                -TKLI*B[L,MR],TKKR,TKKI,B[K,MR],B[K,MI])else 
                COMDIV(-SLR-TLLR*B[L,MR]+TLLI*B[L,MI],-SLI-TLLR*B[L,MI]
                -TLLI*B[L,MR],TLKR,TLKI,B[K,MR],B[K,MI]);L:=L-1
            end;L1:=L
        end;M:=M-1
    end;
    for M:=N step -1 until 1 do 
    for K:=1 step 1 until N do 
    X[K,M]:=MATMAT(1,M,K,M,X,B);
    for M:=N step -1 until 1 do 
    begin S:=0;if ALFI[M]=0 then 
        begin for K:=1 step 1 until N do 
            begin R:=ABS(X[K,M]);
                if R>=S then begin S:=R;D:=X[K,M] end 
            end;for K:=1 step 1 until N do 
            X[K,M]:=X[K,M]/D
        end else 
        begin for K:=1 step 1 until N do 
            begin R:=ABS(X[K,M-1])+ABS(X[K,M]);
                R:=R*SQRT((X[K,M-1]/R)**2+(X[K,M]/R)**2);
                if R>=S"THEN"
                begin S:=R;DR:=X[K,M-1];DI:=X[K,M] end 
            end;for K:=1 step 1 until N do 
            COMDIV(X[K,M-1],X[K,M],DR,DI,X[K,M-1],X[K,M]);M:=M-1
        end 
    end 
end QZVEC;
DWARF:=EM[0];EPS:=EM[1];
HSHDECMUL(N,A,B,DWARF);
HESTGL3(N,A,B,X);
QZIT(N,A,B,X,EPS,EPSA,EPSB,ITER);
QZVAL(N,A,B,X,EPSA,EPSB,ALFR,ALFI,BETA);
QZVEC(N,A,B,X,EPSA,EPSB,ALFR,ALFI,BETA)
end QZI;
        eop