#include <perms.h>
// Edinburgh IMP77 Compiler - Version 8.4
//
// 1 !This is an old ERCC IMP library (advisory have documentation). It contains
// 2 ! a variety of useful procedures
// 3
// 4 ! IMP DFLEPO IMP DGAHER IMP DGALAG IMP DGALEG IMP DIGVEC IMP DKUSCR
// 5 ! IMP DSCRCH IMP FLEPO IMP GAHER IMP GALAG IMP GALEG IMP LOGFCT
// 6 ! IMP DAPRO IMP DCHLIN IMP DELI3 IMP PLEX IMP DPLEX IMP VPLEX
// 7 ! IMP DVPLEX IMP DEXPAN IMP DERFIN IMP CHINT IMP DCHINT IMP DKUMER
// 8 ! IMP DAVDEN IMP DEXPFT IMP DIVINT IMP DCHCOL IMP DRUNGE IMP DRUMOD
// 9 ! IMP DBUTCH IMP DNORSK IMP DGROTA IMP DLINPR IMP INFNT IMP DINFNT
// 10 ! IMP GRAPH9 IMP DCURFT IMP DROOTS IMP DSQRTM IMP DIFSYS IMP CDEVLR
// 11 ! IMP DRATNT IMP ELI3 IMP SELECT IMP FILON IMP DFILON IMP DLEQSP
// 12 ! IMP QADMN IMP DETSP IMP INSORT IMP RLSORT IMP CDHESS
// 13
// 14 !*TOP*IMPDFLEPO
// 15 %external %routine imp dflepo(%routine funct(%long %real %name f, %long %real %array %name x, g), %integer n, limit, iprin,
void IMPDFLEPO(void (*FUNCT)(double *F, double *X, double *G), int N, int LIMIT, int IPRIN, int *CONV, double EPS,
double *F, double *X, double *H) {
// 16+ %integer %name conv, %long %real eps, %long %real %name f,
// 17+ %long %real %array %name x, h)
// 18 %long %real oldf, sg, ghg, z, step, ita, fa, fb, ga, gb, w, lambda
double OLDF;
double SG;
double GHG;
double Z;
double STEP;
double ITA;
double FA;
double FB;
double GA;
double GB;
double W;
double LAMBDA;
// 19 %integer i, j, k, count, fnevalns
int I;
int J;
int K;
int COUNT;
int FNEVALNS;
// 20 %long %real %array g, s, gamma, sigma, y(1:n)
double G[N - 1 + 1];
double S[N - 1 + 1];
double GAMMA[N - 1 + 1];
double SIGMA[N - 1 + 1];
double Y[N - 1 + 1];
// 21
// 22 %long %real %function dot(%long %real %array %name a, b)
double DOT(double *A, double *B) {
// 23 %integer i
int I;
// 24 %long %real s
double S;
// 25 s = 0
S = 0;
// 26 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0001:
if (I == _final) goto L_0002;
I += 1;
// 27 s = s+a(i)*b(i)
S = S + [/*BUG*/ 0] * [/*BUG*/ 0];
// 28 %repeat
goto L_0001;
}
L_0002:;
// 29 %result = s
return S;
// 30 %end
} /* end of DOT */
// 31
// 32 %long %real %function updot(%integer longi, %long %real %array %name a, b)
double UPDOT(int LONGI, double *A, double *B) {
// 33 %comment MULTIPLY BY THE ITH ROW OF THE SYMMETRIC MATRIX A, %c
// 34 !WHOSE UPPER TRIANGLE IS STORED BY ROWS.
// 35 %integer j, k; %long %real s
int J;
int K;
double S;
// 36 k = i; s = 0
K = I;
S = 0;
// 37 ->nrl1 %if i=1
if (I != 1) goto L_0004;
goto U_1000;
L_0004:;
// 38 %cycle j = 1, 1, i-1
// %for J = 1, 1, I - 1 %cycle
{
int _final = I - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, I - 1");
J = 0;
L_0005:
if (J == _final) goto L_0006;
J += 1;
// 39 s = s+a(k)*b(j)
S = S + [/*BUG*/ 0] * [/*BUG*/ 0];
// 40 k = k+n-j
K = K + N - J;
// 41 %repeat; !STEPS TO DIAGONAL:NOW GO ALONG ROW.
goto L_0005;
}
L_0006:;
// 42 nrl1: %cycle j = i, 1, n
U_1000:;
// %for J = I, 1, N %cycle
{
int _initial = I, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = I, 1, N");
J = _initial - (1);
L_0007:
if (J == _final) goto L_0008;
J += 1;
// 43 s = s+a(k+j-i)*b(j)
S = S + [/*BUG*/ 0] * [/*BUG*/ 0];
// 44 %repeat
goto L_0007;
}
L_0008:;
// 45 %result = s
return S;
// 46 %end;
} /* end of UPDOT */
// ?LONGI unused
// 47 ; !OF FN UPDOT
// 47
// 48 nrl1:
U_1001:;
// 49 ->nrl2 %unless conv=1
if (*CONV == 1) goto L_0009;
goto U_1002;
L_0009:;
// 50 k = 1
K = 1;
// 51 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_000A:
if (I == _final) goto L_000B;
I += 1;
// 52 h(k) = 1;
[/*BUG*/ 0] = 1;
// 53 ->nrl30 %if i=n
if (I != N) goto L_000C;
goto U_1003;
L_000C:;
// 54 %cycle j = 1, 1, n-i
// %for J = 1, 1, N - I %cycle
{
int _final = N - I;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, N - I");
J = 0;
L_000D:
if (J == _final) goto L_000E;
J += 1;
// 55 h(k+j) = 0
[/*BUG*/ 0] = 0;
// 56 %repeat
goto L_000D;
}
L_000E:;
// 57 nrl30:k = k+n-i+1
U_1003:;
K = K + N - I + 1;
// 58 %repeat; !FORMATION OF UNIT MATRIX IN H.
goto L_000A;
}
L_000B:;
// 59 nrl2:conv = 1; step = 1; fnevalns = 0
U_1002:;
*CONV = 1;
STEP = 1;
FNEVALNS = 0;
// 60 funct(f, x, g); fnevalns = fnevalns+1
FUNCT(F, X, G);
FNEVALNS = FNEVALNS + 1;
// 61 oldf = f+1
OLDF = *F + 1;
// 62 %cycle count = 1, 1, limit
// %for COUNT = 1, 1, LIMIT %cycle
{
int _final = LIMIT;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle COUNT = 1, 1, LIMIT");
COUNT = 0;
L_0010:
if (COUNT == _final) goto L_0011;
COUNT += 1;
// 63 %return %unless oldf>f
if (OLDF > *F) goto L_0012;
return;
L_0012:;
// 64 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0013:
if (I == _final) goto L_0014;
I += 1;
// 65 sigma(i) = x(i); gamma(i) = g(i); s(i) = -updot(i, h, g)
SIGMA[I - 1] = [/*BUG*/ 0];
GAMMA[I - 1] = G[I - 1];
S[I - 1] = -UPDOT(I, H, G);
// 66 %repeat
goto L_0013;
}
L_0014:;
// 67 fb = f; gb = dot(g, s)
FB = *F;
GB = DOT(G, S);
// 68 ->nrl21 %unless gb>=0
if (GB >= 0) goto L_0015;
goto U_1004;
L_0015:;
// 69 conv = -1; %return
*CONV = -1;
return;
// 70 nrl21:oldf = f; ita = step
U_1004:;
OLDF = *F;
ITA = STEP;
// 71 ! ITA IS DETERMINED FROM ITS VALUE AT THE PREVIOUS ITERATION (STEP)
// 72 ! AND IS INCREASED OR DECREASED BY 4 WHERE NECESSARY. IT SHOULD
// 73 ! TEND TO 1 AT THE MINIMUM.
// 74 nrl4: fa = fb; ga = gb
U_1005:;
FA = FB;
GA = GB;
// 75 %cycle i = 1, 1, n; x(i) = x(i)+ita*s(i); %repeat
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0016:
if (I == _final) goto L_0017;
I += 1;
[/*BUG*/ 0] = [/*BUG*/ 0] + ITA * S[I - 1];
goto L_0016;
}
L_0017:;
// 76 funct(f, x, g); fb = f; gb = dot(g, s); fnevalns = fnevalns+1
FUNCT(F, X, G);
FB = *F;
GB = DOT(G, S);
FNEVALNS = FNEVALNS + 1;
// 77 ->nrl3 %unless (gb<0 %and fb<fa)
if (GB >= 0) goto L_000F;
if (FB < FA) goto L_0018;
L_000F:;
goto U_1006;
L_0018:;
// 78 ita = 4*ita; step = 4*step;
ITA = 4 * ITA;
STEP = 4 * STEP;
// 79 ->nrl4
goto U_1005;
// 80 nrl3: %if ita<=0 %then %start; conv = -1; %return; %finish
U_1006:;
if (ITA > 0) goto L_0019;
*CONV = -1;
return;
L_0019:;
// 81 %cycle i = 1, 1, n; y(i) = g(i); %repeat
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_001A:
if (I == _final) goto L_001B;
I += 1;
Y[I - 1] = G[I - 1];
goto L_001A;
}
L_001B:;
// 82 z = 3*(fa-fb)/ita+ga+gb
Z = (double)((3 * (FA - FB)) / ITA) + GA + GB;
// 83 w = sqrt(z*z-ga*gb)
W = _imp_SQRT(Z * Z - GA * GB);
// 84 lambda = ita*(gb+w-z)/(gb-ga+2*w)
LAMBDA = (double)((ITA * (GB + W - Z)) / (GB - GA + 2 * W));
// 85 %cycle i = 1, 1, n; x(i) = x(i)-lambda*s(i); %repeat
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_001D:
if (I == _final) goto L_001E;
I += 1;
[/*BUG*/ 0] = [/*BUG*/ 0] - LAMBDA * S[I - 1];
goto L_001D;
}
L_001E:;
// 86 funct(f, x, g); fnevalns = fnevalns+1
FUNCT(F, X, G);
FNEVALNS = FNEVALNS + 1;
// 87 ->nrl6 %unless (f>fa %or f>fb)
if (*F > FA) goto L_001F;
if (*F > FB) goto L_001F;
goto U_1000;
L_001F:;
// 88 step = step/4
STEP = (double)(STEP / 4);
// 89 ->nrl5 %unless fb<fa
if (FB < FA) goto L_0020;
goto U_1007;
L_0020:;
// 90 %cycle i = 1, 1, n; x(i) = x(i)+lambda*s(i); g(i) = y(i)
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0021:
if (I == _final) goto L_0022;
I += 1;
[/*BUG*/ 0] = [/*BUG*/ 0] + LAMBDA * S[I - 1];
G[I - 1] = Y[I - 1];
// 91 %repeat
goto L_0021;
}
L_0022:;
// 92 f = fb;
*F = FB;
// 93 ->nrl6
goto U_1000;
// 94 nrl5: gb = dot(g, s)
U_1007:;
GB = DOT(G, S);
// 95 ->nrl7 %unless (gb<0 %and count>n %and step<9.999999999999971@-7)
if (GB >= 0) goto L_001C;
if (COUNT <= N) goto L_001C;
if (STEP < 9.999999999999971E-7) goto L_0023;
L_001C:;
goto U_1008;
L_0023:;
// 96 conv = -1; %return
*CONV = -1;
return;
// 97 nrl7: fb = f; ita = ita-lambda;
U_1008:;
FB = *F;
ITA = ITA - LAMBDA;
// 98 ->nrl3
goto U_1006;
// 99 nrl6: %cycle i = 1, 1, n
U_1000:;
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0024:
if (I == _final) goto L_0025;
I += 1;
// 100 sigma(i) = x(i)-sigma(i)
SIGMA[I - 1] = [/*BUG*/ 0] - SIGMA[I - 1];
// 101 gamma(i) = g(i)-gamma(i)
GAMMA[I - 1] = G[I - 1] - GAMMA[I - 1];
// 102 %repeat
goto L_0024;
}
L_0025:;
// 103 sg = dot(sigma, gamma)
SG = DOT(SIGMA, GAMMA);
// 104 ->nrl13 %unless count>=n
if (COUNT >= N) goto L_0027;
goto U_1009;
L_0027:;
// 105 %return %if sqrt(dot(s, s))<eps %and sqrt(dot(sigma, sigma))<eps
if (_imp_SQRT(DOT(S, S)) >= EPS) goto L_0028;
if (_imp_SQRT(DOT(SIGMA, SIGMA)) >= EPS) goto L_0028;
return;
L_0028:;
// 106 %comment %end OFTESTFORVANISHINGDERIVATIVE.
// 107 nrl13:%cycle i = 1, 1, n
U_1009:;
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0029:
if (I == _final) goto L_002A;
I += 1;
// 108 s(i) = updot(i, h, gamma)
S[I - 1] = UPDOT(I, H, GAMMA);
// 109 %repeat
goto L_0029;
}
L_002A:;
// 110 ghg = dot(s, gamma)
GHG = DOT(S, GAMMA);
// 111 k = 1
K = 1;
// 112 ->nrl8 %if mod(sg)<9.999999999999987@-31 %or %c
if (_imp_MOD(SG) < 9.999999999999987E-31) goto L_0026;
if (_imp_MOD(GHG) >= 9.999999999999987E-31) goto L_002B;
L_0026:;
goto U_100A;
L_002B:;
// 113+ mod(ghg)<9.999999999999987@-31
// 114 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_002C:
if (I == _final) goto L_002D;
I += 1;
// 115 %cycle j = i, 1, n
// %for J = I, 1, N %cycle
{
int _initial = I, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = I, 1, N");
J = _initial - (1);
L_002F:
if (J == _final) goto L_0030;
J += 1;
// 116 h(k) = h(k)+sigma(i)*sigma(j)/sg-s(i)*s(j)/ghg
[/*BUG*/ 0] =
[/*BUG*/ 0] + (double)((SIGMA[I - 1] * SIGMA[J - 1]) / SG) - (double)((S[I - 1] * S[J - 1]) / GHG);
// 117 k = k+1; %repeat
K = K + 1;
goto L_002F;
}
L_0030:;
// 118 %repeat
goto L_002C;
}
L_002D:;
// 119 ->nrl20 %if iprin=0
if (IPRIN != 0) goto L_0032;
goto U_100B;
L_0032:;
// 120 ->nrl20 %unless count//iprin*iprin=count
if (((long)COUNT / (long)IPRIN) * IPRIN == COUNT) goto L_0033;
goto U_100B;
L_0033:;
// 121 nrl8: newlines(2);
U_100A:;
_imp_NEWLINES(2);
// 122 printstring("
_imp_PRINTSTRING(_imp_str_literal("\nITERATION "));
// 123+ ITERATION ")
// 124 print(count, 2, 0)
_imp_PRINT(COUNT, 2, 0);
// 125 printstring(":
_imp_PRINTSTRING(_imp_str_literal(":\n NO. OF FN EVALNS = "));
// 126+ NO. OF FN EVALNS = ")
// 127 write(fnevalns, 2)
_imp_WRITE(FNEVALNS, 2);
// 128 newlines(2);
_imp_NEWLINES(2);
// 129 printstring(" X G
_imp_PRINTSTRING(_imp_str_literal(" X G\n"));
// 130+ ")
// 131 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0034:
if (I == _final) goto L_0035;
I += 1;
// 132 printfl(x(i), 6); spaces(4); printfl(g(i), 6); newline
_imp_PRINTFL([/*BUG*/ 0], 6);
_imp_SPACES(4);
_imp_PRINTFL(G[I - 1], 6);
_imp_NEWLINE();
// 133 %repeat;
goto L_0034;
}
L_0035:;
// 134 printstring("
_imp_PRINTSTRING(_imp_str_literal("\n PRESENT MINIMUM VALUE OF FN ="));
// 135+ PRESENT MINIMUM VALUE OF FN =")
// 136 printfl(f, 6)
_imp_PRINTFL(*F, 6);
// 137 newlines(2)
_imp_NEWLINES(2);
// 138 nrl20:%repeat; ! COUNT LOOP; CONV=-1 ;
U_100B:;
goto L_0010;
}
L_0011:;
// 139 printstring("
_imp_PRINTSTRING(_imp_str_literal("\nPROBLEM DID NOT CONVERGE IN SPECIFIED NO OF ITERATIONS\nROUTINE TERMINATED."));
// 140+ PROBLEM DID NOT CONVERGE IN SPECIFIED NO OF ITERATIONS
// 141+ ROUTINE TERMINATED.")
// 142 newlines(2)
_imp_NEWLINES(2);
// 143 %end;
return;
} /* end of IMPDFLEPO */
// ?NRL1 unused
// 144 ; !OF ROUTINE IMP DFLEPO
// 144 !*END*IMPDFLEPO
// 145
// 146 !*TOP*IMPDGAHER
// 147 %external %routine imp dgaher(%long %real %function f(%long %real x), %long %real eps,
void IMPDGAHER(double (*F)(double X), double EPS, double *ANTGRL, double *TEST) {
// 148+ %long %real %name antgrl, test)
// 149 %long %real i4, i6, i9, i12, i15, c
double I4;
double I6;
double I9;
double I12;
double I15;
double C;
// 150 i4 = .813128354472@-1*(f(1.65068012388)+f(-1.65068012388))+.804914090005* %c
I4 = .813128354472E-1 * (F(1.65068012388) + F(-1.65068012388)) +
.804914090005 * (F(.524647623275) + F(-.524647623275));
// 151+ (f(.524647623275)+f(-.524647623275))
// 152 i6 = .453000990551@-2*(f(2.35060497367)+f(-2.35060497367))+.157067320323* %c
I6 = .453000990551E-2 * (F(2.35060497367) + F(-2.35060497367)) +
.157067320323 * (F(1.33584907401) + F(-1.33584907401)) + .724629595224 * (F(.436077411928) + F(-.436077411928));
// 153+ (f(1.33584907401)+f(-1.33584907401))+.724629595224*(f(.436077411928)+f %c
// 154+ (-.436077411928))
// 155 c = mod(mod(i6)-mod(i4))
C = _imp_MOD(_imp_MOD(I6) - _imp_MOD(I4));
// 156 %if (mod(i6))>1 %then c = c/(mod(i6))
if (_imp_MOD(I6) <= 1) goto L_0036;
C = (double)(C / _imp_MOD(I6));
L_0036:;
// 157 ->nrl1 %unless c<eps
if (C < EPS) goto L_0037;
goto U_100C;
L_0037:;
// 158 antgrl = i6
*ANTGRL = I6;
// 159 test = c
*TEST = C;
// 160 %return
return;
// 161 nrl1:i9 = .396069772632@-4*(f(3.19099320178)+f(-3.19099320178))+ %c
U_100C:;
I9 = .396069772632E-4 * (F(3.19099320178) + F(-3.19099320178)) + .494362427553 -
2 * (F(2.26658058453) + F(-2.26658058453));
// 162+ .494362427553-2*(f(2.26658058453)+f(-2.26658058453))
// 163 i9 = i9+.884745273943-1*(f(1.46855328921)+f(-1.46855328921))+ %c
I9 = I9 + .884745273943 - 1 * (F(1.46855328921) + F(-1.46855328921)) +
.432651559002 * (F(.723551018752) + F(-.723551018752)) + .720235215606 * F(0);
// 164+ .432651559002*(f(.723551018752)+f(-.723551018752))+.720235215606*f(0)
// 165 c = mod(mod(i9)-mod(i6))
C = _imp_MOD(_imp_MOD(I9) - _imp_MOD(I6));
// 166 %if (mod(i9))>1 %then c = c/(mod(i9))
if (_imp_MOD(I9) <= 1) goto L_0038;
C = (double)(C / _imp_MOD(I9));
L_0038:;
// 167 ->nrl2 %unless c<eps
if (C < EPS) goto L_0039;
goto U_100D;
L_0039:;
// 168 antgrl = i9
*ANTGRL = I9;
// 169 test = c
*TEST = C;
// 170 %return
return;
// 171 nrl2:i12 = .265855168436@-6*(f(3.88972489787)+f(-3.88972489787))+ %c
U_100D:;
I12 = .265855168436E-6 * (F(3.88972489787) + F(-3.88972489787)) +
.857368704358E-4 * (F(3.02063702512) + F(-3.02063702512)) +
.390539058463E-2 * (F(2.27950708050) + F(-2.27950708050));
// 172+ .857368704358@-4*(f(3.02063702512)+f(-3.02063702512))+.390539058463@-2* %c
// 173+ (f(2.27950708050)+f(-2.27950708050))
// 174 i12 = i12+.516079856158@-1*(f(1.59768263515)+f(-1.59768263515))+ %c
I12 = I12 + .516079856158E-1 * (F(1.59768263515) + F(-1.59768263515)) +
.260492310264 * (F(.947788391240) + F(-.947788391240)) + .570135236262 * (F(.314240376254) + F(-.314240376254));
// 175+ .260492310264*(f(.947788391240)+f(-.947788391240))+.570135236262*(f( %c
// 176+ .314240376254)+f(-.314240376254))
// 177 c = mod(mod(i12)-mod(i9))
C = _imp_MOD(_imp_MOD(I12) - _imp_MOD(I9));
// 178 %if (mod(i12))>1 %then c = c/(mod(i12))
if (_imp_MOD(I12) <= 1) goto L_003A;
C = (double)(C / _imp_MOD(I12));
L_003A:;
// 179 ->nrl3 %unless c<eps
if (C < EPS) goto L_003B;
goto U_100E;
L_003B:;
// 180 antgrl = i12
*ANTGRL = I12;
// 181 test = c
*TEST = C;
// 182 %return
return;
// 183 nrl3:i15 = .152247580425@-8*(f(4.49999070731)+f(-4.49999070731))+ %c
U_100E:;
I15 = .152247580425E-8 * (F(4.49999070731) + F(-4.49999070731)) +
.105911554771E-5 * (F(3.66995037340) + F(-3.66995037340));
// 184+ .105911554771@-5*(f(3.66995037340)+f(-3.66995037340))
// 185 i15 = i15+.100004441232@-3*(f(2.96716692791)+f(-2.96716692791))+ %c
I15 = I15 + .100004441232E-3 * (F(2.96716692791) + F(-2.96716692791)) +
.277806884291E-2 * (F(2.32573248617) + F(-2.32573248617));
// 186+ .277806884291@-2*(f(2.32573248617)+f(-2.32573248617))
// 187 i15 = i15+.307800338725@-1*(f(1.71999257519)+f(-1.71999257519))+ %c
I15 = I15 + .307800338725E-1 * (F(1.71999257519) + F(-1.71999257519)) +
.158488915796 * (F(1.13611558521) + F(-1.13611558521));
// 188+ .158488915796*(f(1.13611558521)+f(-1.13611558521))
// 189 i15 = i15+.412028687499*(f(.565069583256)+f(-.565069583256))+ %c
I15 = I15 + .412028687499 * (F(.565069583256) + F(-.565069583256)) + .564100308726 * F(0);
// 190+ .564100308726*f(0)
// 191 c = mod(mod(i15)-mod(i12))
C = _imp_MOD(_imp_MOD(I15) - _imp_MOD(I12));
// 192 %if (mod(i15))>1 %then c = c/(mod(i15))
if (_imp_MOD(I15) <= 1) goto L_003C;
C = (double)(C / _imp_MOD(I15));
L_003C:;
// 193 antgrl = i15
*ANTGRL = I15;
// 194 test = c
*TEST = C;
// 195 %return %if c<eps
if (C >= EPS) goto L_003D;
return;
L_003D:;
// 196 printstring("
_imp_PRINTSTRING(_imp_str_literal("\n\nMETHOD FAILED TO SATISFY GIVEN ACCURACY. LAST"));
// 197+
// 198+ METHOD FAILED TO SATISFY GIVEN ACCURACY. LAST")
// 199 printstring(" ESTIMATE OF INTEGRAL IS ")
_imp_PRINTSTRING(_imp_str_literal(" ESTIMATE OF INTEGRAL IS "));
// 200 printfl(antgrl, 10)
_imp_PRINTFL(*ANTGRL, 10);
// 201 printstring("
_imp_PRINTSTRING(_imp_str_literal("\n\nTHE ERROR OF THIS ANSWER IS "));
// 202+
// 203+ THE ERROR OF THIS ANSWER IS ")
// 204 printfl(test, 10)
_imp_PRINTFL(*TEST, 10);
// 205 %return
return;
// 206 %end; !OF IMP DGAHER
} /* end of IMPDGAHER */
// 207 !*END*IMPDGAHER
// 208
// 209 !*TOP*IMPDGALAG
// 210 %external %routine imp dgalag(%long %real %function f(%long %real x), %long %real eps,
void IMPDGALAG(double (*F)(double X), double EPS, double *ANTGRL, double *TEST) {
// 211+ %long %real %name antgrl, test)
// 212 %long %real s4, s6, s9, s12, s15, c
double S4;
double S6;
double S9;
double S12;
double S15;
double C;
// 213 s4 = 6.03154104342@-1*f(0.322547689619)+3.57418692438@-1*f(1.74576110116) %c
S4 = 6.03154104342E-1 * F(0.322547689619) + 3.57418692438E-1 * F(1.74576110116) +
3.88879085150E-2 * F(4.53662029692) + 5.39294705561E-4 * F(9.395070912301);
// 214+ +3.88879085150@-2*f(4.53662029692)+5.39294705561@-4*f(9.395070912301)
// 215 s6 = 4.58964673950@-1*f(.222846604179)+4.17000830772@-1*f(1.18893210167)
S6 = 4.58964673950E-1 * F(.222846604179) + 4.17000830772E-1 * F(1.18893210167);
// 216 s6 = s6+1.13373382074@-1*f(2.99273632606)+1.03991974531@-2*f(5 %c
S6 = S6 + 1.13373382074E-1 * F(2.99273632606) + 1.03991974531E-2 * F(5.77514356910) +
2.61017202815E-4 * F(9.83746741838) + 8.98547906430E-7 * F(15.9828739806);
// 217+ .77514356910)+2.61017202815@-4*f(9.83746741838)+8.98547906430@-7*f(15 %c
// 218+ .9828739806)
// 219 c = mod(mod(s6)-mod(s4))
C = _imp_MOD(_imp_MOD(S6) - _imp_MOD(S4));
// 220 %if (mod(s6))>1 %then c = c/(mod(s6))
if (_imp_MOD(S6) <= 1) goto L_003E;
C = (double)(C / _imp_MOD(S6));
L_003E:;
// 221 ->nrl1 %unless c<eps
if (C < EPS) goto L_003F;
goto U_100D;
L_003F:;
// 222 antgrl = s6
*ANTGRL = S6;
// 223 test = c
*TEST = C;
// 224 %return
return;
// 225 nrl1:s9 = 3.36126421798@-1*f(.152322227732)+4.11213980424@-1*f(.807220022742)
U_100D:;
S9 = 3.36126421798E-1 * F(.152322227732) + 4.11213980424E-1 * F(.807220022742);
// 226 s9 = s9+1.99287525371@-1*f(2.00513515562)+4.74605627657@-2*f(3.78347397333)
S9 = S9 + 1.99287525371E-1 * F(2.00513515562) + 4.74605627657E-2 * F(3.78347397333);
// 227 s9 = s9+5.59962661079@-3*f(6.20495677788)+3.05249767093@-4*f(9.37298525169)
S9 = S9 + 5.59962661079E-3 * F(6.20495677788) + 3.05249767093E-4 * F(9.37298525169);
// 228 s9 = s9+6.59212302608@-6*f(13.4662369111)+4.11076933035@-8*f(18.8335977890)
S9 = S9 + 6.59212302608E-6 * F(13.4662369111) + 4.11076933035E-8 * F(18.8335977890);
// 229 s9 = s9+3.29087403035@-11*f(26.3740718909)
S9 = S9 + 3.29087403035E-11 * F(26.3740718909);
// 230 c = mod(mod(s9)-mod(s6))
C = _imp_MOD(_imp_MOD(S9) - _imp_MOD(S6));
// 231 %if (mod(s9))>1 %then c = c/(mod(s9))
if (_imp_MOD(S9) <= 1) goto L_0040;
C = (double)(C / _imp_MOD(S9));
L_0040:;
// 232 ->nrl2 %unless c<eps
if (C < EPS) goto L_0041;
goto U_100E;
L_0041:;
// 233 antgrl = s9
*ANTGRL = S9;
// 234 test = c
*TEST = C;
// 235 %return
return;
// 236 nrl2:s12 = 2.64731371055@-1*f(.115722117358)+3.77759275873@-1*f(.61175748452) %c
U_100E:;
S12 = 2.64731371055E-1 * F(.115722117358) + 3.77759275873E-1 * F(.61175748452) + 9.04492222117E-2 * F(2.83375133774);
// 237+ +9.04492222117@-2*f(2.83375133774)
// 238 s12 = s12+2.44082011320@-1*f(1.51261026977)
S12 = S12 + 2.44082011320E-1 * F(1.51261026977);
// 239 s12 = s12+2.01023811546@-2*f(4.59922763942)+2.66397354187@-3*f(6 %c
S12 = S12 + 2.01023811546E-2 * F(4.59922763942) + 2.66397354187E-3 * F(6.84452545312);
// 240+ .84452545312)
// 241 s12 = s12+2.03231592663@-4*f(9.62131684246)+8.36505585682@-6*f(13 %c
S12 = S12 + 2.03231592663E-4 * F(9.62131684246) + 8.36505585682E-6 * F(13.0060549933);
// 242+ .0060549933)
// 243 s12 = s12+1.66849387654@-7*f(17.1168551874)+1.34239103051@-9*f(22 %c
S12 = S12 + 1.66849387654E-7 * F(17.1168551874) + 1.34239103051E-9 * F(22.1510903794);
// 244+ .1510903794)
// 245 s12 = s12+3.06160163503@-12*f(28.4879672510)+8.14807746743@-16*f(37 %c
S12 = S12 + 3.06160163503E-12 * F(28.4879672510) + 8.14807746743E-16 * F(37.0991210445);
// 246+ .0991210445)
// 247 c = mod(mod(s12)-mod(s9))
C = _imp_MOD(_imp_MOD(S12) - _imp_MOD(S9));
// 248 %if (mod(s12))>1 %then c = c/(mod(s12))
if (_imp_MOD(S12) <= 1) goto L_0042;
C = (double)(C / _imp_MOD(S12));
L_0042:;
// 249 ->nrl3 %unless c<eps
if (C < EPS) goto L_0043;
goto U_100F;
L_0043:;
// 250 antgrl = s12
*ANTGRL = S12;
// 251 test = c
*TEST = C;
// 252 %return
return;
// 253 nrl3:s15 = 2.18234885940@-1*f(.093307812017)+3.42210177923@-1*f(.492691740302)
U_100F:;
S15 = 2.18234885940E-1 * F(.093307812017) + 3.42210177923E-1 * F(.492691740302);
// 254 s15 = s15+2.63027577942@-1*f(1.21559541207)+1.26425818106@-1*f(2 %c
S15 = S15 + 2.63027577942E-1 * F(1.21559541207) + 1.26425818106E-1 * F(2.26994952620);
// 255+ .26994952620)
// 256 s15 = s15+4.02068649210@-2*f(3.66762272175)+8.56387780361@-3*f(5 %c
S15 = S15 + 4.02068649210E-2 * F(3.66762272175) + 8.56387780361E-3 * F(5.42533662741);
// 257+ .42533662741)
// 258 s15 = s15+1.2124361472@-3*f(7.56591622661)+1.11674392344@-4*f(10.1202285680)
S15 = S15 + 1.2124361472E-3 * F(7.56591622661) + 1.11674392344E-4 * F(10.1202285680);
// 259 s15 = s15+6.45992676202@-6*f(13.1302824822)+2.22631690710@-7*f(16 %c
S15 = S15 + 6.45992676202E-6 * F(13.1302824822) + 2.22631690710E-7 * F(16.6544077083);
// 260+ .6544077083)
// 261 s15 = s15+4.22743038498@-9*f(20.7764788994)+3.92189726704@-11*f(25 %c
S15 = S15 + 4.22743038498E-9 * F(20.7764788994) + 3.92189726704E-11 * F(25.6238942267);
// 262+ .6238942267)
// 263 s15 = s15+1.45651526407@-13*f(31.4075191697)+1.48302705111@-16*f(38 %c
S15 = S15 + 1.45651526407E-13 * F(31.4075191697) + 1.48302705111E-16 * F(38.5306833065);
// 264+ .5306833065)
// 265 s15 = s15+1.60059490621@-20*f(48.0260855727)
S15 = S15 + 1.60059490621E-20 * F(48.0260855727);
// 266 c = mod(mod(s15)-mod(s12))
C = _imp_MOD(_imp_MOD(S15) - _imp_MOD(S12));
// 267 %if (mod(s15))>1 %then c = c/(mod(s15))
if (_imp_MOD(S15) <= 1) goto L_0044;
C = (double)(C / _imp_MOD(S15));
L_0044:;
// 268 antgrl = s15
*ANTGRL = S15;
// 269 test = c
*TEST = C;
// 270 ->nrl4 %unless c<eps
if (C < EPS) goto L_0045;
goto U_1010;
L_0045:;
// 271 %return
return;
// 272 nrl4:
U_1010:;
// 273 printstring("
_imp_PRINTSTRING(_imp_str_literal("\n\nMETHOD FAILED. LAST APPROXIMATION TO VALUE OF INTEGRAL IS "));
// 274+
// 275+ METHOD FAILED. LAST APPROXIMATION TO VALUE OF INTEGRAL IS ")
// 276 printfl(antgrl, 10)
_imp_PRINTFL(*ANTGRL, 10);
// 277 printstring("
_imp_PRINTSTRING(_imp_str_literal("\n\nTHE ERROR OF THIS ANSWER IS "));
// 278+
// 279+ THE ERROR OF THIS ANSWER IS ")
// 280 printfl(test, 10)
_imp_PRINTFL(*TEST, 10);
// 281 %return
return;
// 282 %end; !OF IMP DGALAG
} /* end of IMPDGALAG */
// 283 !*END*IMPDGALAG
// 284
// 285 !*TOP*IMPDGALEG
// 286 %external %routine imp dgaleg(%long %real %function f(%long %real x), %long %real xup, xlo,
void IMPDGALEG(double (*F)(double X), double XUP, double XLO, double *TIN, double *EPS, int LIM) {
// 287+ %long %real %name tin, eps, %integer lim)
// 288 %long %real i5, i6, int5, int6, h, a, xu, xl
double I5;
double I6;
double INT5;
double INT6;
double H;
double A;
double XU;
double XL;
// 289 %integer cntr, i
int CNTR;
int I;
// 290
// 291 %routine comp(%long %real a, %long %real %name i5, i6,
void COMP(double A, double *I5, double *I6, double (*F)(double X)) {
// 292+ %long %real %function f(%long %real x))
// 293 %long %real c
double C;
// 294 c = 0.453089922969332*h
C = 0.453089922969332 * H;
// 295 i5 = 0.118463442528095*(f(a+c)+f(a-c))
*I5 = 0.118463442528095 * (F(A + C) + F(A - C));
// 296 c = 0.269234655052842*h
C = 0.269234655052842 * H;
// 297 i5 = i5+0.239314335249683*(f(a+c)+f(a-c))
*I5 = *I5 + 0.239314335249683 * (F(A + C) + F(A - C));
// 298 i5 = h*(i5+0.284444444444445*f(a))
*I5 = H * (*I5 + 0.284444444444445 * F(A));
// 299 c = 0.466234757101576*h
C = 0.466234757101576 * H;
// 300 i6 = 0.85662246189585@-1*(f(a+c)+f(a-c))
*I6 = 0.85662246189585E-1 * (F(A + C) + F(A - C));
// 301 c = 0.330604693233133*h
C = 0.330604693233133 * H;
// 302 i6 = i6+0.180380786524070*(f(a+c)+f(a-c))
*I6 = *I6 + 0.180380786524070 * (F(A + C) + F(A - C));
// 303 c = 0.119309593041599*h
C = 0.119309593041599 * H;
// 304 i6 = (i6+0.233956967286346*(f(a+c)+f(a-c)))*h
*I6 = (*I6 + 0.233956967286346 * (F(A + C) + F(A - C))) * H;
// 305 %end; ! OF COMP
return;
} /* end of COMP */
// 306
// 307 cntr = 0
CNTR = 0;
// 308 a = .5*(xup+xlo)
A = .5 * (XUP + XLO);
// 309 h = xup-xlo
H = XUP - XLO;
// 310 comp(a, i5, i6, f)
COMP(A, /*?*/ &I5, /*?*/ &I6, F);
// 311 int5 = i5
INT5 = I5;
// 312 int6 = i6
INT6 = I6;
// 313 nrl1:
U_1011:;
// 314 ->nrl3 %unless mod(int6-int5)<eps
if (_imp_MOD(INT6 - INT5) < *EPS) goto L_0046;
goto U_1012;
L_0046:;
// 315 tin = int6
*TIN = INT6;
// 316 eps = mod(int6-int5)
*EPS = _imp_MOD(INT6 - INT5);
// 317 %return
return;
// 318 nrl3:cntr = cntr+1
U_1012:;
CNTR = CNTR + 1;
// 319 ->nrl2 %if cntr<=lim
if (CNTR > LIM) goto L_0047;
goto U_1013;
L_0047:;
// 320 newlines(3)
_imp_NEWLINES(3);
// 321 printstring("
_imp_PRINTSTRING(_imp_str_literal("\nMETHOD FAILED: REQUIRED ACCURACY NOT ATTAINED AFTER"));
// 322+ METHOD FAILED: REQUIRED ACCURACY NOT ATTAINED AFTER")
// 323 write(lim, 5)
_imp_WRITE(LIM, 5);
// 324 printstring(" BISECTIONS OF THE INTERVAL")
_imp_PRINTSTRING(_imp_str_literal(" BISECTIONS OF THE INTERVAL"));
// 325 newlines(3)
_imp_NEWLINES(3);
// 326 %return
return;
// 327 nrl2:xu = xup
U_1013:;
XU = XUP;
// 328 int5 = 0
INT5 = 0;
// 329 int6 = 0
INT6 = 0;
// 330 h = h/2
H = (double)(H / 2);
// 331 %cycle i = 1, 1, 2^^cntr
// %for I = 1, 1, _imp_IEXP(2, CNTR) %cycle
{
int _final = _imp_IEXP(2, CNTR);
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, _imp_IEXP(2, CNTR)");
I = 0;
L_0048:
if (I == _final) goto L_0049;
I += 1;
// 332 xl = xu-h
XL = XU - H;
// 333 a = .5*(xu+xl)
A = .5 * (XU + XL);
// 334 comp(a, i5, i6, f)
COMP(A, /*?*/ &I5, /*?*/ &I6, F);
// 335 int5 = int5+i5
INT5 = INT5 + I5;
// 336 int6 = int6+i6
INT6 = INT6 + I6;
// 337 xu = xl
XU = XL;
// 338 %repeat
goto L_0048;
}
L_0049:;
// 339 ->nrl1
goto U_1011;
// 340 %end; ! OF ROUTINE IMP DGALEG
} /* end of IMPDGALEG */
// 341 !*END*IMPDGALEG
// 342
// 343 !*TOP*IMPDIGVEC
// 344 %external %routine imp digvec(%integer n, r, e, %long %real %array (2) %name a,
void IMPDIGVEC(int N, int R, int E, double *A, double *L, double *M, double *X, double *Y, int *IT, int *VEC,
double EPS) {
auto void DIRHESSE(int N, double *A, int *INT);
auto void QR(int N, double T, double SUM, double PROD, int *ITS, double *A, double *V, double *COM);
auto void REALLU(double *A, int *PERM, int N, double EPS, int *D);
auto void REALSOLVE(double *A, double *B, double *X, int *PERM, int N);
auto void COMLU(double *AX, double *AY, int *PERM, int N, double EPS, int *D);
auto void COMSOLVE(double *AX, double *AY, double *B, double *X, int *PERM, int N);
// 345+ %long %real %array %name l, m, %long %real %array (2) %name x, y, %integer %array %name it, vec, %long %real eps)
// 346 %routine %spec dirhesse(%integer n, %long %real %array (2) %name a,
// 347+ %integer %array %name int)
// 348 %routine %spec qr(%integer n, %long %real t, sum, prod,
// 349+ %integer %array %name its, %long %real %array (2) %name a, %long %real %array %name v, com)
// 350 %routine %spec reallu(%long %real %array (2) %name a,
// 351+ %integer %array %name perm, %integer n, %long %real eps,
// 352+ %integer %name d)
// 353 %routine %spec realsolve(%long %real %array (2) %name a, %long %real %array %name b, x,
// 354+ %integer %array %name perm, %integer n)
// 355 %routine %spec comlu(%long %real %array (2) %name ax, ay,
// 356+ %integer %array %name perm, %integer n, %long %real eps,
// 357+ %integer %name d)
// 358 %routine %spec comsolve(%long %real %array (2) %name ax, ay, %long %real %array %name b, x,
// 359+ %integer %array %name perm, %integer n)
// 360 %long %real s, t, p, q, sm, max, abs, norm, tol
double S;
double T;
double P;
double Q;
double SM;
double MAX;
double ABS;
double NORM;
double TOL;
// 361 %integer i, k, d, rho, mult, sh, j, n2, ctr
int I;
int K;
int D;
int RHO;
int MULT;
int SH;
int J;
int N2;
int CTR;
// 362 %long %real %array b, z(1:2*n), hx, hy(1:n, 1:n)
double B[2 * N - 1 + 1];
double Z[2 * N - 1 + 1];
double HX[N - 1 + 1][N - 1 + 1];
double HY[N - 1 + 1];
// 363 %integer %array inter, sort, perm(1:n)
int INTER[N - 1 + 1];
int SORT[N - 1 + 1];
int PERM[N - 1 + 1];
// 364 %if n<1 %then ->er
if (N >= 1) goto L_004A;
goto U_1014;
L_004A:;
// 365 %if n=1 %then %start
if (N != 1) goto L_004B;
// 366 l(1) = a(1, 1)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 367 m(1) = 0
[/*BUG*/ 0] = 0;
// 368 x(1, 1) = 1
[/*BUG*/ 0] = 1;
// 369 y(1, 1) = 0
[/*BUG*/ 0] = 0;
// 370 it(1) = 0
[/*BUG*/ 0] = 0;
// 371 vec(1) = 0
[/*BUG*/ 0] = 0;
// 372 %return
return;
// 373 %finish
L_004B:;
// 374 e = e//2
E = ((long)E / (long)2);
// 375 abs = 0
ABS = 0;
// 376 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_004C:
if (I == _final) goto L_004D;
I += 1;
// 377 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_004E:
if (K == _final) goto L_004F;
K += 1;
// 378 %if abs<mod(a(i, k)) %then abs = mod(a(i, k))
if (ABS >= _imp_MOD([/*BUG*/ 0])) goto L_0050;
ABS = _imp_MOD([/*BUG*/ 0]);
L_0050:;
// 379 %repeat
goto L_004E;
}
L_004F:;
// 380 vec(i) = -1
[/*BUG*/ 0] = -1;
// 381 it(i) = -n
[/*BUG*/ 0] = -N;
// 382 sort(i) = i
SORT[I - 1] = I;
// 383 %repeat
goto L_004C;
}
L_004D:;
// 384 %if n=2 %then %start
if (N != 2) goto L_0051;
// 385 s = (a(1, 1)+a(2, 2))/2
S = (double)(([/*BUG*/ 0] + [/*BUG*/ 0]) / 2);
// 386 t = (a(1, 1)-a(2, 2))^2+4*a(1, 2)*a(2, 1)
T = _imp_REXP([/*BUG*/ 0] - [/*BUG*/ 0], 2) + 4 * [/*BUG*/ 0] * [/*BUG*/ 0];
// 387 p = sqrt(mod(t))/2
P = (double)(_imp_SQRT(_imp_MOD(T)) / 2);
// 388 %if t<0 %then %start
if (T >= 0) goto L_0052;
// 389 l(1) = s
[/*BUG*/ 0] = S;
// 390 l(2) = s
[/*BUG*/ 0] = S;
// 391 m(1) = p
[/*BUG*/ 0] = P;
// 392 m(2) = -p
[/*BUG*/ 0] = -P;
// 393 %finish %else %start
goto L_0053;
L_0052:;
// 394 l(1) = s+p
[/*BUG*/ 0] = S + P;
// 395 l(2) = s-p
[/*BUG*/ 0] = S - P;
// 396 m(1) = 0
[/*BUG*/ 0] = 0;
// 397 m(2) = 0
[/*BUG*/ 0] = 0;
// 398 %if l(1)<l(2) %then %start
if ([/*BUG*/ 0] >= [/*BUG*/ 0]) goto L_0054;
// 399 t = l(1)
T = [/*BUG*/ 0];
// 400 l(1) = l(2)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 401 l(2) = t
[/*BUG*/ 0] = T;
// 402 %finish
L_0054:;
// 403 %finish
L_0053:;
// 404 it(1) = 0
[/*BUG*/ 0] = 0;
// 405 it(2) = 0
[/*BUG*/ 0] = 0;
// 406 %finish %else %start
goto L_0055;
L_0051:;
// 407 dirhesse(n, a, inter)
DIRHESSE(N, A, INTER);
// 408 s = a(n, n)+a(n-1, n-1)
S = [/*BUG*/ 0] + [/*BUG*/ 0];
// 409 t = a(n, n)*a(n-1, n-1)-a(n, n-1)*a(n-1, n)
T = [/*BUG*/ 0] * [/*BUG*/ 0] - [/*BUG*/ 0] * [/*BUG*/ 0];
// 410 qr(n, eps, s, t, it, a, b, l)
QR(N, EPS, S, T, IT, A, B, L);
// 411 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_0056:
if (K == _final) goto L_0057;
K += 1;
// 412 m(k) = b(k)^2+l(k)^2
[/*BUG*/ 0] = _imp_REXP(B[K - 1], 2) + _imp_REXP([/*BUG*/ 0], 2);
// 413 %repeat
goto L_0056;
}
L_0057:;
// 414 %cycle i = 1, 1, n-1
// %for I = 1, 1, N - 1 %cycle
{
int _final = N - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N - 1");
I = 0;
L_0058:
if (I == _final) goto L_0059;
I += 1;
// 415 s = m(i)
S = [/*BUG*/ 0];
// 416 %cycle k = i+1, 1, n
// %for K = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = I + 1, 1, N");
K = _initial - (1);
L_005A:
if (K == _final) goto L_005B;
K += 1;
// 417 %if s<m(k) %then %start
if (S >= [/*BUG*/ 0]) goto L_005C;
// 418 m(i) = m(k)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 419 m(k) = s
[/*BUG*/ 0] = S;
// 420 s = m(i)
S = [/*BUG*/ 0];
// 421 d = sort(i)
D = SORT[I - 1];
// 422 sort(i) = sort(k)
SORT[I - 1] = SORT[K - 1];
// 423 sort(k) = d
SORT[K - 1] = D;
// 424 %finish
L_005C:;
// 425 %repeat
goto L_005A;
}
L_005B:;
// 426 %repeat
goto L_0058;
}
L_0059:;
// 427 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_005D:
if (K == _final) goto L_005E;
K += 1;
// 428 m(k) = l(sort(k))
[/*BUG*/ 0] = [/*BUG*/ 0];
// 429 %repeat
goto L_005D;
}
L_005E:;
// 430 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_005F:
if (K == _final) goto L_0060;
K += 1;
// 431 l(k) = b(sort(k))
[/*BUG*/ 0] = B[SORT[K - 1] - 1];
// 432 %repeat
goto L_005F;
}
L_0060:;
// 433 %finish
L_0055:;
// 434 %if r=0 %then %return
if (R != 0) goto L_0061;
return;
L_0061:;
// 435 %if r>0 %then d = 1 %else d = -1
if (R <= 0) goto L_0062;
D = 1;
goto L_0063;
L_0062:;
D = -1;
L_0063:;
// 436 r = imod(r)
R = _imp_IMOD(R);
// 437 %if r>n %then r = n
if (R <= N) goto L_0064;
R = N;
L_0064:;
// 438 ctr = 0
CTR = 0;
// 439 n2 = 2*n
N2 = 2 * N;
// 440 tol = (n*eps*100)^(-2)
TOL = _imp_REXP(N * EPS * 100, -2);
// 441 %cycle rho = 1, 1, r
// %for RHO = 1, 1, R %cycle
{
int _final = R;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle RHO = 1, 1, R");
RHO = 0;
L_0065:
if (RHO == _final) goto L_0066;
RHO += 1;
// 442 sh = 1
SH = 1;
// 443 %if d<0 %then j = n+1-rho %else j = rho
if (D >= 0) goto L_0067;
J = N + 1 - RHO;
goto L_0068;
L_0067:;
J = RHO;
L_0068:;
// 444 %if it(j)<0 %then ->v3
if ([/*BUG*/ 0] >= 0) goto L_0069;
goto U_1015;
L_0069:;
// 445 p = l(j)
P = [/*BUG*/ 0];
// 446 q = m(j)
Q = [/*BUG*/ 0];
// 447 %if ctr=1 %then %start
if (CTR != 1) goto L_006A;
// 448 vec(j) = 0
[/*BUG*/ 0] = 0;
// 449 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_006B:
if (K == _final) goto L_006C;
K += 1;
// 450 x(k, j) = x(k, j-d)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 451 y(k, j) = -y(k, j-d)
[/*BUG*/ 0] = -[/*BUG*/ 0];
// 452 %repeat
goto L_006B;
}
L_006C:;
// 453 ctr = 0
CTR = 0;
// 454 ->v3
goto U_1015;
// 455 %finish
L_006A:;
// 456 %if mod(q)<=eps*abs %then ctr = 0 %else ctr = 1
if (_imp_MOD(Q) > EPS * ABS) goto L_006D;
CTR = 0;
goto L_006E;
L_006D:;
CTR = 1;
L_006E:;
// 457 v1: %cycle i = 1, 1, n
U_1016:;
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_006F:
if (I == _final) goto L_0070;
I += 1;
// 458 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_0071:
if (K == _final) goto L_0072;
K += 1;
// 459 hx(i, k) = a(i, k)
HX[K - 1] = [/*BUG*/ 0];
// 460 %repeat
goto L_0071;
}
L_0072:;
// 461 hx(i, i) = hx(i, i)-p
HX[I - 1] = HX[I - 1] - P;
// 462 z(i+n) = 0
Z[I + N - 1] = 0;
// 463 %repeat
goto L_006F;
}
L_0070:;
// 464 sm = 1/2^intpt(e/sh)
SM = (double)(1 / (_imp_REXP(2, _imp_INTPT((double)(E / SH)))));
// 465 %if ctr=1 %then %start
if (CTR != 1) goto L_0073;
// 466 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0074:
if (I == _final) goto L_0075;
I += 1;
// 467 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_0076:
if (K == _final) goto L_0077;
K += 1;
// 468 hy(i, k) = 0
HY[K - 1] = 0;
// 469 %repeat
goto L_0076;
}
L_0077:;
// 470 hy(i, i) = -q
HY[I - 1] = -Q;
// 471 %repeat
goto L_0074;
}
L_0075:;
// 472 comlu(hx, hy, perm, n, sm, sh)
COMLU(HX, HY, PERM, N, SM, /*?*/ &SH);
// 473 %finish %else reallu(hx, perm, n, sm, sh)
goto L_0079;
L_0073:;
REALLU(HX, PERM, N, SM, /*?*/ &SH);
L_0079:;
// 474 %if sh<0 %then %start
if (SH >= 0) goto L_007A;
// 475 sh = imod(sh)+1
SH = _imp_IMOD(SH) + 1;
// 476 ->v1
goto U_1016;
// 477 %finish
L_007A:;
// 478 %cycle k = 1, 1, n2
// %for K = 1, 1, N2 %cycle
{
int _final = N2;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N2");
K = 0;
L_007B:
if (K == _final) goto L_007C;
K += 1;
// 479 b(k) = 0
B[K - 1] = 0;
// 480 %repeat
goto L_007B;
}
L_007C:;
// 481 b(sort(j)) = 1
B[SORT[J - 1] - 1] = 1;
// 482 %cycle i = 1, 1, 4
// %for I = 1, 1, 4 %cycle
{
I = 0;
L_007D:
if (I == 4) goto L_007E;
I += 1;
// 483 %if ctr=1 %then %start
if (CTR != 1) goto L_007F;
// 484 comsolve(hx, hy, b, z, perm, n)
COMSOLVE(HX, HY, B, Z, PERM, N);
// 485 %finish %else %start
goto L_0080;
L_007F:;
// 486 realsolve(hx, b, z, perm, n)
REALSOLVE(HX, B, Z, PERM, N);
// 487 %finish
L_0080:;
// 488 max = 0
MAX = 0;
// 489 norm = 0
NORM = 0;
// 490 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_0081:
if (K == _final) goto L_0082;
K += 1;
// 491 p = z(k)^2+z(k+n)^2
P = _imp_REXP(Z[K - 1], 2) + _imp_REXP(Z[K + N - 1], 2);
// 492 norm = norm+p
NORM = NORM + P;
// 493 %if max<p %then %start
if (MAX >= P) goto L_0083;
// 494 max = p
MAX = P;
// 495 mult = k
MULT = K;
// 496 %finish
L_0083:;
// 497 %repeat
goto L_0081;
}
L_0082:;
// 498 %if norm>=tol %then ->v2
if (NORM < TOL) goto L_0084;
goto U_1017;
L_0084:;
// 499 p = z(mult)
P = Z[MULT - 1];
// 500 q = z(mult+n)
Q = Z[MULT + N - 1];
// 501 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_0085:
if (K == _final) goto L_0086;
K += 1;
// 502 b(k) = (z(k)*p+z(k+n)*q)/max
B[K - 1] = (double)((Z[K - 1] * P + Z[K + N - 1] * Q) / MAX);
// 503 b(k+n) = (z(k+n)*p-z(k)*q)/max
B[K + N - 1] = (double)((Z[K + N - 1] * P - Z[K - 1] * Q) / MAX);
// 504 %repeat
goto L_0085;
}
L_0086:;
// 505 %if norm<1@4 %then %start
if (NORM >= 1E4) goto L_0087;
// 506 mult = -1
MULT = -1;
// 507 %cycle sh = 1, 1, n
// %for SH = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle SH = 1, 1, N");
SH = 0;
L_0088:
if (SH == _final) goto L_0089;
SH += 1;
// 508 k = sh
K = SH;
// 509 v5:
U_1018:;
// 510 %if mod(b(k))>9.999999999999989@-9 %then %start
if (_imp_MOD(B[K - 1]) <= 9.999999999999989E-9) goto L_008A;
// 511 b(k) = mult*b(k)
B[K - 1] = MULT * B[K - 1];
// 512 %if k=sh %then mult = -mult
if (K != SH) goto L_008B;
MULT = -MULT;
L_008B:;
// 513 %finish %else b(k) = b(k)+1000*eps
goto L_008C;
L_008A:;
B[K - 1] = B[K - 1] + 1000 * EPS;
L_008C:;
// 514 %if k=sh %and ctr=1 %then %start
if (K != SH) goto L_008D;
if (CTR != 1) goto L_008D;
// 515 k = sh+n
K = SH + N;
// 516 ->v5
goto U_1018;
// 517 %finish
L_008D:;
// 518 %repeat
goto L_0088;
}
L_0089:;
// 519 %finish
L_0087:;
// 520 %repeat
goto L_007D;
}
L_007E:;
// 521 i = -4
I = -4;
// 522 v2: vec(j) = i
U_1017:;
[/*BUG*/ 0] = I;
// 523 p = sqrt(norm)
P = _imp_SQRT(NORM);
// 524 %cycle k = 1, 1, n2
// %for K = 1, 1, N2 %cycle
{
int _final = N2;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N2");
K = 0;
L_008E:
if (K == _final) goto L_008F;
K += 1;
// 525 b(k) = z(k)/p
B[K - 1] = (double)(Z[K - 1] / P);
// 526 %repeat
goto L_008E;
}
L_008F:;
// 527 z(1) = b(1)
Z[0] = B[0];
// 528 z(2) = b(2)
Z[1] = B[1];
// 529 z(n+1) = b(n+1)
Z[N + 1 - 1] = B[N + 1 - 1];
// 530 z(n+2) = b(n+2)
Z[N + 2 - 1] = B[N + 2 - 1];
// 531 %if n=2 %then ->v4
if (N != 2) goto L_0090;
goto U_1019;
L_0090:;
// 532 %cycle sh = 3, 1, n
// %for SH = 3, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(3, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle SH = 3, 1, N");
SH = 2;
L_0091:
if (SH == _final) goto L_0092;
SH += 1;
// 533 i = 0
I = 0;
// 534 v6: s = b(sh+i)
U_101A:;
S = B[SH + I - 1];
// 535 %cycle k = 2, 1, sh-1
// %for K = 2, 1, SH - 1 %cycle
{
int _final = SH - 1;
if (invalid_for_loop(2, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 2, 1, SH - 1");
K = 1;
L_0093:
if (K == _final) goto L_0094;
K += 1;
// 536 s = s+a(sh, k-1)*b(k+i)
S = S + [/*BUG*/ 0] * B[K + I - 1];
// 537 %repeat
goto L_0093;
}
L_0094:;
// 538 z(sh+i) = s
Z[SH + I - 1] = S;
// 539 %if i=0 %and ctr=1 %then %start
if (I != 0) goto L_0095;
if (CTR != 1) goto L_0095;
// 540 i = n
I = N;
// 541 ->v6
goto U_101A;
// 542 %finish
L_0095:;
// 543 %repeat
goto L_0091;
}
L_0092:;
// 544 %cycle k = n-1, -1, 1
// %for K = N - 1, -1, 1 %cycle
{
int _initial = N - 1;
if (invalid_for_loop(_initial, -1, 1)) _imp_signal(5, 1, 0, "Illegal cycle K = N - 1, -1, 1");
K = _initial - (-1);
L_0096:
if (K == 1) goto L_0097;
K -= 1;
// 545 mult = inter(k)
MULT = INTER[K - 1];
// 546 %if mult#k+1 %then %start
if (MULT == K + 1) goto L_0098;
// 547 i = k+1
I = K + 1;
// 548 v7: s = z(mult)
U_101B:;
S = Z[MULT - 1];
// 549 z(mult) = z(i)
Z[MULT - 1] = Z[I - 1];
// 550 z(i) = s
Z[I - 1] = S;
// 551 %if mult<=n %then %start
if (MULT > N) goto L_0099;
// 552 mult = mult+n
MULT = MULT + N;
// 553 i = i+n
I = I + N;
// 554 ->v7
goto U_101B;
// 555 %finish
L_0099:;
// 556 %finish
L_0098:;
// 557 %repeat
goto L_0096;
}
L_0097:;
// 558 v4: norm = 0
U_1019:;
NORM = 0;
// 559 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_009A:
if (K == _final) goto L_009B;
K += 1;
// 560 norm = norm+z(k)^2+z(k+n)^2
NORM = NORM + _imp_REXP(Z[K - 1], 2) + _imp_REXP(Z[K + N - 1], 2);
// 561 %repeat
goto L_009A;
}
L_009B:;
// 562 p = sqrt(norm)
P = _imp_SQRT(NORM);
// 563 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_009C:
if (K == _final) goto L_009D;
K += 1;
// 564 x(k, j) = z(k)/p
[/*BUG*/ 0] = (double)(Z[K - 1] / P);
// 565 y(k, j) = z(k+n)/p
[/*BUG*/ 0] = (double)(Z[K + N - 1] / P);
// 566 %repeat
goto L_009C;
}
L_009D:;
// 567 v3:%repeat
U_1015:;
goto L_0065;
}
L_0066:;
// 568 %return
return;
// 569 er:newline
U_1014:;
_imp_NEWLINE();
// 570 printstring("N TOO SMALL, NAMELY N = ")
_imp_PRINTSTRING(_imp_str_literal("N TOO SMALL, NAMELY N = "));
// 571 write(n, 3)
_imp_WRITE(N, 3);
// 572 printstring(" THIS ERROR OCCURRED IN ROUTINE == IMP DIGVEC ==")
_imp_PRINTSTRING(_imp_str_literal(" THIS ERROR OCCURRED IN ROUTINE == IMP DIGVEC =="));
// 573 newline
_imp_NEWLINE();
// 574 %routine reallu(%long %real %array (2) %name a, %integer %array %name perm,
void REALLU(double *A, int *PERM, int N, double EPS, int *D) {
// 575+ %integer n, %long %real eps, %integer %name d)
// 576 %long %real c
double C;
// 577 %integer i, j, k, co
int I;
int J;
int K;
int CO;
// 578 co = 0
CO = 0;
// 579 %cycle i = 1, 1, n-1
// %for I = 1, 1, N - 1 %cycle
{
int _final = N - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N - 1");
I = 0;
L_009E:
if (I == _final) goto L_009F;
I += 1;
// 580 j = i+1
J = I + 1;
// 581 %if mod(a(j, i))>mod(a(i, i)) %then %start
if (_imp_MOD([/*BUG*/ 0]) <= _imp_MOD([/*BUG*/ 0])) goto L_00A0;
// 582 perm(i) = 1
[/*BUG*/ 0] = 1;
// 583 %cycle k = i, 1, n
// %for K = I, 1, N %cycle
{
int _initial = I, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = I, 1, N");
K = _initial - (1);
L_00A1:
if (K == _final) goto L_00A2;
K += 1;
// 584 c = a(i, k)
C = [/*BUG*/ 0];
// 585 a(i, k) = a(j, k)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 586 a(j, k) = c
[/*BUG*/ 0] = C;
// 587 %repeat
goto L_00A1;
}
L_00A2:;
// 588 %finish %else perm(i) = 0
goto L_00A3;
L_00A0:;
[/*BUG*/ 0] = 0;
L_00A3:;
// 589 %if mod(a(i, i))<=eps %then %start
if (_imp_MOD([/*BUG*/ 0]) > EPS) goto L_00A4;
// 590 co = co+1
CO = CO + 1;
// 591 %if co>d %then %start
if (CO <= *D) goto L_00A5;
// 592 d = -d
*D = -*D;
// 593 %return
return;
// 594 %finish
L_00A5:;
// 595 a(i, i) = eps
[/*BUG*/ 0] = EPS;
// 596 %finish
L_00A4:;
// 597 c = a(j, i)/a(i, i)
C = (double)([/*BUG*/ 0] / [/*BUG*/ 0]);
// 598 %cycle k = j, 1, n
// %for K = J, 1, N %cycle
{
int _initial = J, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = J, 1, N");
K = _initial - (1);
L_00A6:
if (K == _final) goto L_00A7;
K += 1;
// 599 a(j, k) = a(j, k)-c*a(i, k)
[/*BUG*/ 0] = [/*BUG*/ 0] - C * [/*BUG*/ 0];
// 600 %repeat
goto L_00A6;
}
L_00A7:;
// 601 a(j, i) = c
[/*BUG*/ 0] = C;
// 602 %repeat
goto L_009E;
}
L_009F:;
// 603 %if mod(a(n, n))<=eps %then %start
if (_imp_MOD([/*BUG*/ 0]) > EPS) goto L_00A8;
// 604 co = co+1
CO = CO + 1;
// 605 %if co>d %then %start
if (CO <= *D) goto L_00A9;
// 606 d = -d
*D = -*D;
// 607 %return
return;
// 608 %finish
L_00A9:;
// 609 a(n, n) = eps
[/*BUG*/ 0] = EPS;
// 610 %finish
L_00A8:;
// 611 %end
return;
} /* end of REALLU */
// 612 %routine realsolve(%long %real %array (2) %name a, %long %real %array %name b, x,
void REALSOLVE(double *A, double *B, double *X, int *PERM, int N) {
// 613+ %integer %array %name perm, %integer n)
// 614 %long %real s
double S;
// 615 %integer i, j, k
int I;
int J;
int K;
// 616 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_00AA:
if (I == _final) goto L_00AB;
I += 1;
// 617 x(i) = b(i)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 618 %repeat
goto L_00AA;
}
L_00AB:;
// 619 %cycle i = 1, 1, n-1
// %for I = 1, 1, N - 1 %cycle
{
int _final = N - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N - 1");
I = 0;
L_00AC:
if (I == _final) goto L_00AD;
I += 1;
// 620 j = i+1
J = I + 1;
// 621 %if perm(i)#0 %then %start
if ([/*BUG*/ 0] == 0) goto L_00AE;
// 622 s = x(j)
S = [/*BUG*/ 0];
// 623 x(j) = x(i)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 624 x(i) = s
[/*BUG*/ 0] = S;
// 625 %finish
L_00AE:;
// 626 x(j) = x(j)-a(j, i)*x(i)
[/*BUG*/ 0] = [/*BUG*/ 0] - [/*BUG*/ 0] * [/*BUG*/ 0];
// 627 %repeat
goto L_00AC;
}
L_00AD:;
// 628 x(n) = x(n)/a(n, n)
[/*BUG*/ 0] = (double)([/*BUG*/ 0] / [/*BUG*/ 0]);
// 629 %cycle i = n-1, -1, 1
// %for I = N - 1, -1, 1 %cycle
{
int _initial = N - 1;
if (invalid_for_loop(_initial, -1, 1)) _imp_signal(5, 1, 0, "Illegal cycle I = N - 1, -1, 1");
I = _initial - (-1);
L_00AF:
if (I == 1) goto L_00B0;
I -= 1;
// 630 s = x(i)
S = [/*BUG*/ 0];
// 631 %cycle k = i+1, 1, n
// %for K = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = I + 1, 1, N");
K = _initial - (1);
L_00B1:
if (K == _final) goto L_00B2;
K += 1;
// 632 s = s-a(i, k)*x(k)
S = S - [/*BUG*/ 0] * [/*BUG*/ 0];
// 633 %repeat
goto L_00B1;
}
L_00B2:;
// 634 x(i) = s/a(i, i)
[/*BUG*/ 0] = (double)(S / [/*BUG*/ 0]);
// 635 %repeat
goto L_00AF;
}
L_00B0:;
// 636 %end
return;
} /* end of REALSOLVE */
// 637 %routine comlu(%long %real %array (2) %name ax, ay,
void COMLU(double *AX, double *AY, int *PERM, int N, double EPS, int *D) {
// 638+ %integer %array %name perm, %integer n, %long %real eps,
// 639+ %integer %name d)
// 640 %long %real cx, cy, mod
double CX;
double CY;
double MOD;
// 641 %integer i, j, k, co
int I;
int J;
int K;
int CO;
// 642 co = 0
CO = 0;
// 643 %cycle i = 1, 1, n-1
// %for I = 1, 1, N - 1 %cycle
{
int _final = N - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N - 1");
I = 0;
L_00B3:
if (I == _final) goto L_00B4;
I += 1;
// 644 j = i+1
J = I + 1;
// 645 cx = ax(j, i)^2+ay(j, i)^2
CX = _imp_REXP([/*BUG*/ 0], 2) + _imp_REXP([/*BUG*/ 0], 2);
// 646 cy = ax(i, i)^2+ay(i, i)^2
CY = _imp_REXP([/*BUG*/ 0], 2) + _imp_REXP([/*BUG*/ 0], 2);
// 647 %if cx>cy %then %start
if (CX <= CY) goto L_00B5;
// 648 perm(i) = 1
[/*BUG*/ 0] = 1;
// 649 %cycle k = i, 1, n
// %for K = I, 1, N %cycle
{
int _initial = I, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = I, 1, N");
K = _initial - (1);
L_00B6:
if (K == _final) goto L_00B7;
K += 1;
// 650 cx = ax(i, k)
CX = [/*BUG*/ 0];
// 651 ax(i, k) = ax(j, k)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 652 ax(j, k) = cx
[/*BUG*/ 0] = CX;
// 653 cx = ay(i, k)
CX = [/*BUG*/ 0];
// 654 ay(i, k) = ay(j, k)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 655 ay(j, k) = cx
[/*BUG*/ 0] = CX;
// 656 %repeat
goto L_00B6;
}
L_00B7:;
// 657 %finish %else perm(i) = 0
goto L_00B8;
L_00B5:;
[/*BUG*/ 0] = 0;
L_00B8:;
// 658 mod = ax(i, i)^2+ay(i, i)^2
MOD = _imp_REXP([/*BUG*/ 0], 2) + _imp_REXP([/*BUG*/ 0], 2);
// 659 %if mod<=eps^2 %then %start
if (MOD > _imp_REXP(EPS, 2)) goto L_00B9;
// 660 co = co+1
CO = CO + 1;
// 661 %if co>d %then %start
if (CO <= *D) goto L_00BA;
// 662 d = -d
*D = -*D;
// 663 %return
return;
// 664 %finish
L_00BA:;
// 665 ax(i, i) = eps
[/*BUG*/ 0] = EPS;
// 666 ay(i, i) = 0
[/*BUG*/ 0] = 0;
// 667 mod = eps^2
MOD = _imp_REXP(EPS, 2);
// 668 %finish
L_00B9:;
// 669 cx = (ax(j, i)*ax(i, i)+ay(j, i)*ay(i, i))/mod
CX = (double)(([/*BUG*/ 0] * [/*BUG*/ 0] + [/*BUG*/ 0] * [/*BUG*/ 0]) / MOD);
// 670 cy = (ay(j, i)*ax(i, i)-ax(j, i)*ay(i, i))/mod
CY = (double)(([/*BUG*/ 0] * [/*BUG*/ 0] - [/*BUG*/ 0] * [/*BUG*/ 0]) / MOD);
// 671 %cycle k = j, 1, n
// %for K = J, 1, N %cycle
{
int _initial = J, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = J, 1, N");
K = _initial - (1);
L_00BB:
if (K == _final) goto L_00BC;
K += 1;
// 672 mod = ay(j, k)-cx*ay(i, k)-cy*ax(i, k)
MOD = [/*BUG*/ 0] - CX * [/*BUG*/ 0] - CY * [/*BUG*/ 0];
// 673 ax(j, k) = ax(j, k)-cx*ax(i, k)+cy*ay(i, k)
[/*BUG*/ 0] = [/*BUG*/ 0] - CX * [/*BUG*/ 0] + CY * [/*BUG*/ 0];
// 674 ay(j, k) = mod
[/*BUG*/ 0] = MOD;
// 675 %repeat
goto L_00BB;
}
L_00BC:;
// 676 ax(j, i) = cx
[/*BUG*/ 0] = CX;
// 677 ay(j, i) = cy
[/*BUG*/ 0] = CY;
// 678 %repeat
goto L_00B3;
}
L_00B4:;
// 679 mod = ax(n, n)^2+ay(n, n)^2
MOD = _imp_REXP([/*BUG*/ 0], 2) + _imp_REXP([/*BUG*/ 0], 2);
// 680 %if mod<=eps^2 %then %start
if (MOD > _imp_REXP(EPS, 2)) goto L_00BD;
// 681 co = co+1
CO = CO + 1;
// 682 %if co>d %then %start
if (CO <= *D) goto L_00BE;
// 683 d = -d
*D = -*D;
// 684 %return
return;
// 685 %finish
L_00BE:;
// 686 ax(n, n) = eps
[/*BUG*/ 0] = EPS;
// 687 ay(n, n) = 0
[/*BUG*/ 0] = 0;
// 688 %finish
L_00BD:;
// 689 %end
return;
} /* end of COMLU */
// 690 %routine comsolve(%long %real %array (2) %name ax, ay, %long %real %array %name b, x,
void COMSOLVE(double *AX, double *AY, double *B, double *X, int *PERM, int N) {
// 691+ %integer %array %name perm, %integer n)
// 692 %long %real sx, sy, mod
double SX;
double SY;
double MOD;
// 693 %integer i, j, k
int I;
int J;
int K;
// 694 %cycle i = 1, 1, 2*n
// %for I = 1, 1, 2 * N %cycle
{
int _final = 2 * N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, 2 * N");
I = 0;
L_00BF:
if (I == _final) goto L_00C0;
I += 1;
// 695 x(i) = b(i)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 696 %repeat
goto L_00BF;
}
L_00C0:;
// 697 %cycle i = 1, 1, n-1
// %for I = 1, 1, N - 1 %cycle
{
int _final = N - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N - 1");
I = 0;
L_00C1:
if (I == _final) goto L_00C2;
I += 1;
// 698 j = i+1
J = I + 1;
// 699 %if perm(i)#0 %then %start
if ([/*BUG*/ 0] == 0) goto L_00C3;
// 700 mod = x(j)
MOD = [/*BUG*/ 0];
// 701 x(j) = x(i)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 702 x(i) = mod
[/*BUG*/ 0] = MOD;
// 703 mod = x(j+n)
MOD = [/*BUG*/ 0];
// 704 x(j+n) = x(i+n)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 705 x(i+n) = mod
[/*BUG*/ 0] = MOD;
// 706 %finish
L_00C3:;
// 707 x(j) = x(j)-ax(j, i)*x(i)+ay(j, i)*x(i+n)
[/*BUG*/ 0] = [/*BUG*/ 0] - [/*BUG*/ 0] * [/*BUG*/ 0] + [/*BUG*/ 0] * [/*BUG*/ 0];
// 708 x(j+n) = x(j+n)-ax(j, i)*x(i+n)-ay(j, i)*x(i)
[/*BUG*/ 0] = [/*BUG*/ 0] - [/*BUG*/ 0] * [/*BUG*/ 0] - [/*BUG*/ 0] * [/*BUG*/ 0];
// 709 %repeat
goto L_00C1;
}
L_00C2:;
// 710 mod = ax(n, n)^2+ay(n, n)^2
MOD = _imp_REXP([/*BUG*/ 0], 2) + _imp_REXP([/*BUG*/ 0], 2);
// 711 sx = (x(n)*ax(n, n)+x(2*n)*ay(n, n))/mod
SX = (double)(([/*BUG*/ 0] * [/*BUG*/ 0] + [/*BUG*/ 0] * [/*BUG*/ 0]) / MOD);
// 712 x(2*n) = (x(2*n)*ax(n, n)-x(n)*ay(n, n))/mod
[/*BUG*/ 0] = (double)(([/*BUG*/ 0] * [/*BUG*/ 0] - [/*BUG*/ 0] * [/*BUG*/ 0]) / MOD);
// 713 x(n) = sx
[/*BUG*/ 0] = SX;
// 714 %cycle i = n-1, -1, 1
// %for I = N - 1, -1, 1 %cycle
{
int _initial = N - 1;
if (invalid_for_loop(_initial, -1, 1)) _imp_signal(5, 1, 0, "Illegal cycle I = N - 1, -1, 1");
I = _initial - (-1);
L_00C4:
if (I == 1) goto L_00C5;
I -= 1;
// 715 sx = x(i)
SX = [/*BUG*/ 0];
// 716 sy = x(i+n)
SY = [/*BUG*/ 0];
// 717 %cycle k = i+1, 1, n
// %for K = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = I + 1, 1, N");
K = _initial - (1);
L_00C6:
if (K == _final) goto L_00C7;
K += 1;
// 718 sx = sx-ax(i, k)*x(k)+ay(i, k)*x(k+n)
SX = SX - [/*BUG*/ 0] * [/*BUG*/ 0] + [/*BUG*/ 0] * [/*BUG*/ 0];
// 719 sy = sy-ax(i, k)*x(k+n)-ay(i, k)*x(k)
SY = SY - [/*BUG*/ 0] * [/*BUG*/ 0] - [/*BUG*/ 0] * [/*BUG*/ 0];
// 720 %repeat
goto L_00C6;
}
L_00C7:;
// 721 mod = ax(i, i)^2+ay(i, i)^2
MOD = _imp_REXP([/*BUG*/ 0], 2) + _imp_REXP([/*BUG*/ 0], 2);
// 722 x(i) = (sx*ax(i, i)+sy*ay(i, i))/mod
[/*BUG*/ 0] = (double)((SX * [/*BUG*/ 0] + SY * [/*BUG*/ 0]) / MOD);
// 723 x(i+n) = (sy*ax(i, i)-sx*ay(i, i))/mod
[/*BUG*/ 0] = (double)((SY * [/*BUG*/ 0] - SX * [/*BUG*/ 0]) / MOD);
// 724 %repeat
goto L_00C4;
}
L_00C5:;
// 725 %end
return;
} /* end of COMSOLVE */
// 726 %routine dirhesse(%integer n, %long %real %array (2) %name a,
void DIRHESSE(int N, double *A, int *INT) {
// 727+ %integer %array %name int)
// 728 %long %real m
double M;
// 729 %integer i, j, k, l
int I;
int J;
int K;
int L;
// 730 %if n<2 %then ->er
if (N >= 2) goto L_00C8;
goto U_101C;
L_00C8:;
// 731 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_00C9:
if (I == _final) goto L_00CA;
I += 1;
// 732 int(i) = i+1
[/*BUG*/ 0] = I + 1;
// 733 %repeat
goto L_00C9;
}
L_00CA:;
// 734 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_00CB:
if (K == _final) goto L_00CC;
K += 1;
// 735 %if k=1 %then ->l1
if (K != 1) goto L_00CD;
goto U_101D;
L_00CD:;
// 736 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_00CE:
if (I == _final) goto L_00CF;
I += 1;
// 737 m = a(i, k)
M = [/*BUG*/ 0];
// 738 l = k-1
L = K - 1;
// 739 %if k=n %then ->l2
if (K != N) goto L_00D0;
goto U_101E;
L_00D0:;
// 740 %cycle j = k+1, 1, n
// %for J = K + 1, 1, N %cycle
{
int _initial = K + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = K + 1, 1, N");
J = _initial - (1);
L_00D1:
if (J == _final) goto L_00D2;
J += 1;
// 741 m = m+a(i, j)*a(j, l)
M = M + [/*BUG*/ 0] * [/*BUG*/ 0];
// 742 %repeat
goto L_00D1;
}
L_00D2:;
// 743 l2: %if i<=2 %then ->l4
U_101E:;
if (I > 2) goto L_00D3;
goto U_101F;
L_00D3:;
// 744 %if i>k %then ->l3
if (I <= K) goto L_00D4;
goto U_1020;
L_00D4:;
// 745 %cycle j = i-1, -1, 2
// %for J = I - 1, -1, 2 %cycle
{
int _initial = I - 1;
if (invalid_for_loop(_initial, -1, 2)) _imp_signal(5, 1, 0, "Illegal cycle J = I - 1, -1, 2");
J = _initial - (-1);
L_00D5:
if (J == 2) goto L_00D6;
J -= 1;
// 746 m = m-a(i, j-1)*a(j, k)
M = M - [/*BUG*/ 0] * [/*BUG*/ 0];
// 747 %repeat
goto L_00D5;
}
L_00D6:;
// 748 ->l4
goto U_101F;
// 749 l3: %cycle j = k, -1, 2
U_1020:;
// %for J = K, -1, 2 %cycle
{
int _initial = K;
if (invalid_for_loop(_initial, -1, 2)) _imp_signal(5, 1, 0, "Illegal cycle J = K, -1, 2");
J = _initial - (-1);
L_00D7:
if (J == 2) goto L_00D8;
J -= 1;
// 750 m = m-a(i, j-1)*a(j, k)
M = M - [/*BUG*/ 0] * [/*BUG*/ 0];
// 751 %repeat
goto L_00D7;
}
L_00D8:;
// 752 l4: a(i, k) = m
U_101F:;
[/*BUG*/ 0] = M;
// 753 %repeat
goto L_00CE;
}
L_00CF:;
// 754 l1: l = k+1
U_101D:;
L = K + 1;
// 755 i = l
I = L;
// 756 %if i>=n %then ->l5
if (I < N) goto L_00D9;
goto U_1021;
L_00D9:;
// 757 m = mod(a(i, k))
M = _imp_MOD([/*BUG*/ 0]);
// 758 %cycle j = i+1, 1, n
// %for J = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = I + 1, 1, N");
J = _initial - (1);
L_00DA:
if (J == _final) goto L_00DB;
J += 1;
// 759 %if mod(a(j, k))<=m %then ->l8
if (_imp_MOD([/*BUG*/ 0]) > M) goto L_00DC;
goto U_1022;
L_00DC:;
// 760 m = mod(a(j, k))
M = _imp_MOD([/*BUG*/ 0]);
// 761 i = j
I = J;
// 762 l8: %repeat
U_1022:;
goto L_00DA;
}
L_00DB:;
// 763 int(k) = i
[/*BUG*/ 0] = I;
// 764 %if i=l %then ->l7
if (I != L) goto L_00DD;
goto U_1023;
L_00DD:;
// 765 %cycle j = 1, 1, n
// %for J = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, N");
J = 0;
L_00DE:
if (J == _final) goto L_00DF;
J += 1;
// 766 m = a(l, j)
M = [/*BUG*/ 0];
// 767 a(l, j) = a(i, j)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 768 a(i, j) = m
[/*BUG*/ 0] = M;
// 769 %repeat
goto L_00DE;
}
L_00DF:;
// 770 %cycle j = 1, 1, n
// %for J = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, N");
J = 0;
L_00E0:
if (J == _final) goto L_00E1;
J += 1;
// 771 m = a(j, l)
M = [/*BUG*/ 0];
// 772 a(j, l) = a(j, i)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 773 a(j, i) = m
[/*BUG*/ 0] = M;
// 774 %repeat
goto L_00E0;
}
L_00E1:;
// 775 l7: m = a(l, k)
U_1023:;
M = [/*BUG*/ 0];
// 776 %if m=0 %or k>n-2 %then ->l5
if (M == 0) goto L_002E;
if (K <= N - 2) goto L_00E2;
L_002E:;
goto U_1021;
L_00E2:;
// 777 %cycle i = k+2, 1, n
// %for I = K + 2, 1, N %cycle
{
int _initial = K + 2, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = K + 2, 1, N");
I = _initial - (1);
L_00E3:
if (I == _final) goto L_00E4;
I += 1;
// 778 a(i, k) = a(i, k)/m
[/*BUG*/ 0] = (double)([/*BUG*/ 0] / M);
// 779 %repeat
goto L_00E3;
}
L_00E4:;
// 780 l5: %repeat
U_1021:;
goto L_00CB;
}
L_00CC:;
// 781 %return
return;
// 782 er: newline
U_101C:;
_imp_NEWLINE();
// 783 printstring("N TOO SMALL, NAMELY N = ")
_imp_PRINTSTRING(_imp_str_literal("N TOO SMALL, NAMELY N = "));
// 784 print(n, 3, 0)
_imp_PRINT(N, 3, 0);
// 785 printstring(" THIS ERROR OCCURRED IN ROUTINE == DIRHESSE ==")
_imp_PRINTSTRING(_imp_str_literal(" THIS ERROR OCCURRED IN ROUTINE == DIRHESSE =="));
// 786 newline
_imp_NEWLINE();
// 787 %end
return;
} /* end of DIRHESSE */
// 788 %routine qr(%integer n, %long %real t, sum, prod,
void QR(int N, double T, double SUM, double PROD, int *ITS, double *A, double *V, double *COM) {
// 789+ %integer %array %name its, %long %real %array (2) %name a, %long %real %array %name v, com)
// 790 %long %real z, norm, s, x, p, q
double Z;
double NORM;
double S;
double X;
double P;
double Q;
// 791 %integer i, j, k, l, m, u, preset
int I;
int J;
int K;
int L;
int M;
int U;
int PRESET;
// 792 %long %real %array aa(1:n, 1:n), c(1:6)
double AA[N - 1 + 1][N - 1 + 1];
double C[6];
// 793 %switch branch(1:3)
int BRANCH_index;
int BRANCH_line;
char *BRANCH_file;
static const void *BRANCH[3] = {&&case_BRANCH_1, &&case_BRANCH_2, &&case_BRANCH_3};
// 794 %if n<3 %then ->er
if (N >= 3) goto L_00E6;
goto U_1024;
L_00E6:;
// 795 preset = 1
PRESET = 1;
// 796 %if sum*prod#0 %then ->q1
if (SUM * PROD == 0) goto L_00E7;
goto U_1025;
L_00E7:;
// 797 preset = -1
PRESET = -1;
// 798 q1: norm = 0
U_1025:;
NORM = 0;
// 799 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_00E8:
if (I == _final) goto L_00E9;
I += 1;
// 800 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_00EA:
if (K == _final) goto L_00EB;
K += 1;
// 801 aa(i, k) = a(i, k)
AA[I - 1][K - 1] = [/*BUG*/ 0];
// 802 %if i<k+2 %then ->q2
if (I >= K + 2) goto L_00EC;
goto U_1026;
L_00EC:;
// 803 aa(i, k) = 0
AA[I - 1][K - 1] = 0;
// 804 q2:
U_1026:;
// 805 norm = norm+mod(aa(i, k))
NORM = NORM + _imp_MOD(AA[I - 1][K - 1]);
// 806 %repeat
goto L_00EA;
}
L_00EB:;
// 807 its(i) = 0
[/*BUG*/ 0] = 0;
// 808 %repeat
goto L_00E8;
}
L_00E9:;
// 809 z = norm*t
Z = NORM * T;
// 810 l = 1
L = 1;
// 811 u = n
U = N;
// 812 l2: %cycle i = u, -1, l+1
U_1027:;
// %for I = U, -1, L + 1 %cycle
{
int _initial = U, _final = L + 1;
if (invalid_for_loop(_initial, -1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = U, -1, L + 1");
I = _initial - (-1);
L_00ED:
if (I == _final) goto L_00EE;
I -= 1;
// 813 %if mod(aa(i, i-1))<z %then ->l1
if (_imp_MOD(AA[I - 1][I - 1 - 1]) >= Z) goto L_00EF;
goto U_1028;
L_00EF:;
// 814 %repeat
goto L_00ED;
}
L_00EE:;
// 815 %if its(u)>=40 %then %start
if ([/*BUG*/ 0] < 40) goto L_00F0;
// 816 sum = 0
SUM = 0;
// 817 %cycle i = 1, 1, u
// %for I = 1, 1, U %cycle
{
int _final = U;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, U");
I = 0;
L_00F1:
if (I == _final) goto L_00F2;
I += 1;
// 818 sum = sum+mod(aa(l, i))
SUM = SUM + _imp_MOD(AA[L - 1][I - 1]);
// 819 %repeat
goto L_00F1;
}
L_00F2:;
// 820 %cycle i = l+1, 1, u
// %for I = L + 1, 1, U %cycle
{
int _initial = L + 1, _final = U;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = L + 1, 1, U");
I = _initial - (1);
L_00F3:
if (I == _final) goto L_00F4;
I += 1;
// 821 prod = 0
PROD = 0;
// 822 %cycle j = i-1, 1, u
// %for J = I - 1, 1, U %cycle
{
int _initial = I - 1, _final = U;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = I - 1, 1, U");
J = _initial - (1);
L_00F5:
if (J == _final) goto L_00F6;
J += 1;
// 823 prod = prod+mod(aa(i, j))
PROD = PROD + _imp_MOD(AA[I - 1][J - 1]);
// 824 %repeat
goto L_00F5;
}
L_00F6:;
// 825 %if sum>prod %then ->q3
if (SUM <= PROD) goto L_00F7;
goto U_1029;
L_00F7:;
// 826 sum = prod
SUM = PROD;
// 827 q3: %repeat
U_1029:;
goto L_00F3;
}
L_00F4:;
// 828 prod = sum*sum
PROD = SUM * SUM;
// 829 sum = 1.5*sum
SUM = 1.5 * SUM;
// 830 preset = 1
PRESET = 1;
// 831 %if its(u)>100*n %then %start
if ([/*BUG*/ 0] <= 100 * N) goto L_00F8;
// 832 %cycle i = 1, 1, u
// %for I = 1, 1, U %cycle
{
int _final = U;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, U");
I = 0;
L_00F9:
if (I == _final) goto L_00FA;
I += 1;
// 833 its(i) = -u
[/*BUG*/ 0] = -U;
// 834 v(i) = 0
[/*BUG*/ 0] = 0;
// 835 com(i) = 0
[/*BUG*/ 0] = 0;
// 836 %repeat
goto L_00F9;
}
L_00FA:;
// 837 %return
return;
// 838 %finish
L_00F8:;
// 839 %finish
L_00F0:;
// 840 %if preset=1 %then ->l9
if (PRESET != 1) goto L_00FB;
goto U_102A;
L_00FB:;
// 841 sum = aa(u-1, u-1)+aa(u, u)
SUM = AA[U - 1 - 1][U - 1 - 1] + AA[U - 1][U - 1];
// 842 prod = aa(u-1, u-1)*aa(u, u)-aa(u-1, u)*aa(u, u-1)
PROD = AA[U - 1 - 1][U - 1 - 1] * AA[U - 1][U - 1] - AA[U - 1 - 1][U - 1] * AA[U - 1][U - 1 - 1];
// 843 l9: %if u-2<l+1 %then ->m3
U_102A:;
if (U - 2 >= L + 1) goto L_00FC;
goto U_102B;
L_00FC:;
// 844 %cycle i = u-2, -1, l+1
// %for I = U - 2, -1, L + 1 %cycle
{
int _initial = U - 2, _final = L + 1;
if (invalid_for_loop(_initial, -1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = U - 2, -1, L + 1");
I = _initial - (-1);
L_00FD:
if (I == _final) goto L_00FE;
I -= 1;
// 845 c(1) = aa(i, i)*(aa(i, i)-sum)+prod+aa(i, i+1)*aa(i+1, i)
C[0] = AA[I - 1][I - 1] * (AA[I - 1][I - 1] - SUM) + PROD + AA[I - 1][I + 1 - 1] * AA[I + 1 - 1][I - 1];
// 846 %if c(1)#0 %then %start
if (C[0] == 0) goto L_00FF;
// 847 c(2) = aa(i+1, i)*(aa(i, i)+aa(i+1, i+1)-sum)
C[1] = AA[I + 1 - 1][I - 1] * (AA[I - 1][I - 1] + AA[I + 1 - 1][I + 1 - 1] - SUM);
// 848 c(3) = 0
C[2] = 0;
// 849 %if i+2>u %then ->q4
if (I + 2 <= U) goto L_0100;
goto U_102C;
L_0100:;
// 850 c(3) = aa(i+2, i+1)*aa(i+1, i)
C[2] = AA[I + 2 - 1][I + 1 - 1] * AA[I + 1 - 1][I - 1];
// 851 q4:
U_102C:;
// 852 %if mod(aa(i, i-1)*(mod(c(2))+mod(c(3)))/c(1))<z %then %start
if (_imp_MOD((double)((AA[I - 1][I - 1 - 1] * (_imp_MOD(C[1]) + _imp_MOD(C[2]))) / C[0])) >= Z) goto L_0101;
// 853 m = i
M = I;
// 854 ->l8
goto U_102D;
// 855 %finish
L_0101:;
// 856 %finish
L_00FF:;
// 857 %repeat
goto L_00FD;
}
L_00FE:;
// 858 m3: m = l
U_102B:;
M = L;
// 859 c(1) = aa(l, l)*(aa(l, l)-sum)+prod+aa(l, l+1)*aa(l+1, l)
C[0] = AA[L - 1][L - 1] * (AA[L - 1][L - 1] - SUM) + PROD + AA[L - 1][L + 1 - 1] * AA[L + 1 - 1][L - 1];
// 860 c(2) = aa(l+1, l)*(aa(l, l)+aa(l+1, l+1)-sum)
C[1] = AA[L + 1 - 1][L - 1] * (AA[L - 1][L - 1] + AA[L + 1 - 1][L + 1 - 1] - SUM);
// 861 c(3) = aa(l+2, l+1)*aa(l+1, l)
C[2] = AA[L + 2 - 1][L + 1 - 1] * AA[L + 1 - 1][L - 1];
// 862 l8: %if m>u-1 %then ->m4
U_102D:;
if (M <= U - 1) goto L_0102;
goto U_102E;
L_0102:;
// 863 %cycle k = m, 1, u-1
// %for K = M, 1, U - 1 %cycle
{
int _initial = M, _final = U - 1;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = M, 1, U - 1");
K = _initial - (1);
L_0103:
if (K == _final) goto L_0104;
K += 1;
// 864 %if k#m %then %start
if (K == M) goto L_0105;
// 865 c(1) = aa(k, k-1)
C[0] = AA[K - 1][K - 1 - 1];
// 866 c(2) = aa(k+1, k-1)
C[1] = AA[K + 1 - 1][K - 1 - 1];
// 867 c(3) = 0
C[2] = 0;
// 868 %if k+2>u %then ->q5
if (K + 2 <= U) goto L_0106;
goto U_102F;
L_0106:;
// 869 c(3) = aa(k+2, k-1)
C[2] = AA[K + 2 - 1][K - 1 - 1];
// 870 ->q5
goto U_102F;
// 871 %finish
L_0105:;
// 872 its(u) = its(u)+1
[/*BUG*/ 0] = [/*BUG*/ 0] + 1;
// 873 q5: s = c(1)*c(1)+c(2)*c(2)+c(3)*c(3)
U_102F:;
S = C[0] * C[0] + C[1] * C[1] + C[2] * C[2];
// 874 s = sqrt(s)
S = _imp_SQRT(S);
// 875 %if c(1)>0 %then ->q6
if (C[0] <= 0) goto L_0107;
goto U_1030;
L_0107:;
// 876 s = -s
S = -S;
// 877 q6:
U_1030:;
// 878 %if mod(s)<z %then %start
if (_imp_MOD(S) >= Z) goto L_0108;
// 879 %if k=m %then %start
if (K != M) goto L_0109;
// 880 p = aa(u-1, u-1)-aa(u, u)
P = AA[U - 1 - 1][U - 1 - 1] - AA[U - 1][U - 1];
// 881 s = sum
S = SUM;
// 882 q = aa(u-1, u)*aa(u, u-1)
Q = AA[U - 1 - 1][U - 1] * AA[U - 1][U - 1 - 1];
// 883 %if p*p+4*q>0 %and preset=-1 %then %start
if (P * P + 4 * Q <= 0) goto L_010A;
if (PRESET != -1) goto L_010A;
// 884 q = sqrt(p*p+4*q)
Q = _imp_SQRT(P * P + 4 * Q);
// 885 x = (s+q)/2
X = (double)((S + Q) / 2);
// 886 %if mod(s+q)<=mod(s-q) %then ->q7
if (_imp_MOD(S + Q) > _imp_MOD(S - Q)) goto L_010B;
goto U_1031;
L_010B:;
// 887 x = (s-q)/2
X = (double)((S - Q) / 2);
// 888 %finish
L_010A:;
// 889 x = s/2
X = (double)(S / 2);
// 890 q7: s = x
U_1031:;
S = X;
// 891 %cycle i = m, 1, u
// %for I = M, 1, U %cycle
{
int _initial = M, _final = U;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = M, 1, U");
I = _initial - (1);
L_010C:
if (I == _final) goto L_010D;
I += 1;
// 892 aa(i, i) = aa(i, i)-s
AA[I - 1][I - 1] = AA[I - 1][I - 1] - S;
// 893 %repeat
goto L_010C;
}
L_010D:;
// 894 %cycle k = m, 1, u-1
// %for K = M, 1, U - 1 %cycle
{
int _initial = M, _final = U - 1;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = M, 1, U - 1");
K = _initial - (1);
L_010E:
if (K == _final) goto L_010F;
K += 1;
// 895 x = sqrt(aa(k, k)*aa(k, k)+aa(k+1, k)*aa(k+1, k))
X = _imp_SQRT(AA[K - 1][K - 1] * AA[K - 1][K - 1] + AA[K + 1 - 1][K - 1] * AA[K + 1 - 1][K - 1]);
// 896 v(k) = aa(k, k)/x
[/*BUG*/ 0] = (double)(AA[K - 1][K - 1] / X);
// 897 com(k) = aa(k+1, k)/x
[/*BUG*/ 0] = (double)(AA[K + 1 - 1][K - 1] / X);
// 898 aa(k, k) = x
AA[K - 1][K - 1] = X;
// 899 aa(k+1, k) = 0
AA[K + 1 - 1][K - 1] = 0;
// 900 %cycle j = k+1, 1, u
// %for J = K + 1, 1, U %cycle
{
int _initial = K + 1, _final = U;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = K + 1, 1, U");
J = _initial - (1);
L_0110:
if (J == _final) goto L_0111;
J += 1;
// 901 x = v(k)*aa(k, j)+com(k)*aa(k+1, j)
X = [/*BUG*/ 0] * AA[K - 1][J - 1] + [/*BUG*/ 0] * AA[K + 1 - 1][J - 1];
// 902 aa(k+1, j) = v(k)*aa(k+1, j)-com(k)*aa(k, j)
AA[K + 1 - 1][J - 1] = [/*BUG*/ 0] * AA[K + 1 - 1][J - 1] - [/*BUG*/ 0] * AA[K - 1][J - 1];
// 903 aa(k, j) = x
AA[K - 1][J - 1] = X;
// 904 %repeat
goto L_0110;
}
L_0111:;
// 905 %repeat
goto L_010E;
}
L_010F:;
// 906 %cycle k = m, 1, u-1
// %for K = M, 1, U - 1 %cycle
{
int _initial = M, _final = U - 1;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = M, 1, U - 1");
K = _initial - (1);
L_0113:
if (K == _final) goto L_0114;
K += 1;
// 907 %cycle i = k+1, -1, l
// %for I = K + 1, -1, L %cycle
{
int _initial = K + 1, _final = L;
if (invalid_for_loop(_initial, -1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = K + 1, -1, L");
I = _initial - (-1);
L_0115:
if (I == _final) goto L_0116;
I -= 1;
// 908 x = v(k)*aa(i, k)+com(k)*aa(i, k+1)
X = [/*BUG*/ 0] * AA[I - 1][K - 1] + [/*BUG*/ 0] * AA[I - 1][K + 1 - 1];
// 909 aa(i, k+1) = v(k)*aa(i, k+1)-com(k)*aa(i, k)
AA[I - 1][K + 1 - 1] = [/*BUG*/ 0] * AA[I - 1][K + 1 - 1] - [/*BUG*/ 0] * AA[I - 1][K - 1];
// 910 aa(i, k) = x
AA[I - 1][K - 1] = X;
// 911 %repeat
goto L_0115;
}
L_0116:;
// 912 %repeat
goto L_0113;
}
L_0114:;
// 913 %cycle i = m, 1, u
// %for I = M, 1, U %cycle
{
int _initial = M, _final = U;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = M, 1, U");
I = _initial - (1);
L_0117:
if (I == _final) goto L_0118;
I += 1;
// 914 aa(i, i) = aa(i, i)+s
AA[I - 1][I - 1] = AA[I - 1][I - 1] + S;
// 915 %repeat
goto L_0117;
}
L_0118:;
// 916 preset = -1
PRESET = -1;
// 917 ->l2
goto U_1027;
// 918 %finish
L_0109:;
// 919 aa(k+1, k-1) = 0
AA[K + 1 - 1][K - 1 - 1] = 0;
// 920 %if k+2>u %then ->l7
if (K + 2 <= U) goto L_0119;
goto U_1032;
L_0119:;
// 921 aa(k+2, k-1) = 0
AA[K + 2 - 1][K - 1 - 1] = 0;
// 922 ->l7
goto U_1032;
// 923 %finish
L_0108:;
// 924 c(1) = c(1)+s
C[0] = C[0] + S;
// 925 c(4) = c(1)/s
C[3] = (double)(C[0] / S);
// 926 c(5) = c(2)/s
C[4] = (double)(C[1] / S);
// 927 c(6) = c(3)/s
C[5] = (double)(C[2] / S);
// 928 c(2) = c(2)/c(1)
C[1] = (double)(C[1] / C[0]);
// 929 c(3) = c(3)/c(1)
C[2] = (double)(C[2] / C[0]);
// 930 %if k#m %then %start
if (K == M) goto L_011A;
// 931 aa(k, k-1) = -s
AA[K - 1][K - 1 - 1] = -S;
// 932 aa(k+1, k-1) = 0
AA[K + 1 - 1][K - 1 - 1] = 0;
// 933 %if k+2>u %then ->q8
if (K + 2 <= U) goto L_011B;
goto U_1033;
L_011B:;
// 934 aa(k+2, k-1) = 0
AA[K + 2 - 1][K - 1 - 1] = 0;
// 935 ->q8
goto U_1033;
// 936 %finish
L_011A:;
// 937 %if m=l %then ->q8
if (M != L) goto L_011C;
goto U_1033;
L_011C:;
// 938 aa(k, k-1) = -aa(k, k-1)
AA[K - 1][K - 1 - 1] = -AA[K - 1][K - 1 - 1];
// 939 q8: %cycle j = k, 1, u
U_1033:;
// %for J = K, 1, U %cycle
{
int _initial = K, _final = U;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = K, 1, U");
J = _initial - (1);
L_011D:
if (J == _final) goto L_011E;
J += 1;
// 940 x = aa(k, j)+c(2)*aa(k+1, j)
X = AA[K - 1][J - 1] + C[1] * AA[K + 1 - 1][J - 1];
// 941 %if k+2>u %then ->q9
if (K + 2 <= U) goto L_011F;
goto U_1034;
L_011F:;
// 942 x = x+c(3)*aa(k+2, j)
X = X + C[2] * AA[K + 2 - 1][J - 1];
// 943 aa(k+2, j) = aa(k+2, j)-x*c(6)
AA[K + 2 - 1][J - 1] = AA[K + 2 - 1][J - 1] - X * C[5];
// 944 q9: aa(k, j) = aa(k, j)-x*c(4)
U_1034:;
AA[K - 1][J - 1] = AA[K - 1][J - 1] - X * C[3];
// 945 aa(k+1, j) = aa(k+1, j)-x*c(5)
AA[K + 1 - 1][J - 1] = AA[K + 1 - 1][J - 1] - X * C[4];
// 946 %repeat
goto L_011D;
}
L_011E:;
// 947 j = u
J = U;
// 948 %if k+3>u %then ->m1
if (K + 3 <= U) goto L_0120;
goto U_1035;
L_0120:;
// 949 j = k+3
J = K + 3;
// 950 m1: %cycle i = 1, 1, j
U_1035:;
// %for I = 1, 1, J %cycle
{
int _final = J;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, J");
I = 0;
L_0121:
if (I == _final) goto L_0122;
I += 1;
// 951 x = c(4)*aa(i, k)+c(5)*aa(i, k+1)
X = C[3] * AA[I - 1][K - 1] + C[4] * AA[I - 1][K + 1 - 1];
// 952 %if k+2>u %then ->m2
if (K + 2 <= U) goto L_0123;
goto U_1036;
L_0123:;
// 953 x = x+c(6)*aa(i, k+2)
X = X + C[5] * AA[I - 1][K + 2 - 1];
// 954 aa(i, k+2) = aa(i, k+2)-x*c(3)
AA[I - 1][K + 2 - 1] = AA[I - 1][K + 2 - 1] - X * C[2];
// 955 m2: aa(i, k) = aa(i, k)-x
U_1036:;
AA[I - 1][K - 1] = AA[I - 1][K - 1] - X;
// 956 aa(i, k+1) = aa(i, k+1)-x*c(2)
AA[I - 1][K + 1 - 1] = AA[I - 1][K + 1 - 1] - X * C[1];
// 957 %repeat
goto L_0121;
}
L_0122:;
// 958 l7: %repeat
U_1032:;
goto L_0103;
}
L_0104:;
// 959 m4: preset = -1
U_102E:;
PRESET = -1;
// 960 ->l2
goto U_1027;
// 961 l1: %if i=u %then %start
U_1028:;
if (I != U) goto L_0124;
// 962 branch(2):v(u) = aa(u, u)
case_BRANCH_2:;
[/*BUG*/ 0] = AA[U - 1][U - 1];
// 963 com(u) = 0
[/*BUG*/ 0] = 0;
// 964 u = u-1
U = U - 1;
// 965 ->l3
goto U_1037;
// 966 %finish
L_0124:;
// 967 %if i=u-1 %then %start
if (I != U - 1) goto L_0126;
// 968 branch(3):p = aa(i, i)-aa(u, u)
case_BRANCH_3:;
P = AA[I - 1][I - 1] - AA[U - 1][U - 1];
// 969 x = aa(i, i)+aa(u, u)
X = AA[I - 1][I - 1] + AA[U - 1][U - 1];
// 970 q = aa(u, i)*aa(i, u)
Q = AA[U - 1][I - 1] * AA[I - 1][U - 1];
// 971 s = aa(i, i)*aa(u, u)
S = AA[I - 1][I - 1] * AA[U - 1][U - 1];
// 972 %if p*p+4*q<=0 %then %start
if (P * P + 4 * Q > 0) goto L_0127;
// 973 com(u) = (sqrt(-p*p-4*q))/2
[/*BUG*/ 0] = (double)(_imp_SQRT(-(P * P) - 4 * Q) / 2);
// 974 com(u-1) = -com(u)
[/*BUG*/ 0] = -[/*BUG*/ 0];
// 975 v(u) = x/2
[/*BUG*/ 0] = (double)(X / 2);
// 976 v(u-1) = v(u)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 977 %finish %else %start
goto L_0128;
L_0127:;
// 978 j = 1
J = 1;
// 979 %if x>=0 %then ->q0
if (X < 0) goto L_0129;
goto U_1038;
L_0129:;
// 980 j = -1
J = -1;
// 981 q0: v(u) = (x+j*sqrt(p*p+4*q))/2
U_1038:;
[/*BUG*/ 0] = (double)((X + J * _imp_SQRT(P * P + 4 * Q)) / 2);
// 982 v(u-1) = (s-q)/v(u)
[/*BUG*/ 0] = (double)((S - Q) / [/*BUG*/ 0]);
// 983 com(u) = 0
[/*BUG*/ 0] = 0;
// 984 com(u-1) = 0
[/*BUG*/ 0] = 0;
// 985 %finish
L_0128:;
// 986 u = u-2
U = U - 2;
// 987 ->l3
goto U_1037;
// 988 %finish
L_0126:;
// 989 l = i
L = I;
// 990 ->l9
goto U_102A;
// 991 l3: %if u-l>=2 %then ->l2
U_1037:;
if (U - L < 2) goto L_012A;
goto U_1027;
L_012A:;
// 992 i = l
I = L;
// 993 ->branch(u-l+2)
BRANCH_line = __LINE__;
BRANCH_file = __FILE__;
BRANCH_index = U - L + 2;
goto case_BRANCH_dispatch;
// 994 branch(1):%if l#1 %then %start
case_BRANCH_1:;
if (L == 1) goto L_012B;
// 995 l = 1
L = 1;
// 996 %if u>2 %then ->l2
if (U <= 2) goto L_012C;
goto U_1027;
L_012C:;
// 997 i = 2
I = 2;
// 998 %if u=2 %and mod(aa(2, 1))<z %then ->l1
if (U != 2) goto L_012D;
if (_imp_MOD(AA[1][0]) >= Z) goto L_012D;
goto U_1028;
L_012D:;
// 999 i = 1
I = 1;
// 1000 ->l1
goto U_1028;
// 1001 %finish
L_012B:;
// 1002 %return
return;
// 1003 er: newline
U_1024:;
_imp_NEWLINE();
// 1004 printstring("N TOO SMALL, NAMELY N = ")
_imp_PRINTSTRING(_imp_str_literal("N TOO SMALL, NAMELY N = "));
// 1005 print(n, 3, 0)
_imp_PRINT(N, 3, 0);
// 1006 printstring(" THIS ERROR OCCURRED IN ROUTINE == QR ==")
_imp_PRINTSTRING(_imp_str_literal(" THIS ERROR OCCURRED IN ROUTINE == QR =="));
// 1007 newline
_imp_NEWLINE();
// 1008 %end
return;
case_BRANCH_dispatch:
if (BRANCH_index < 1 || BRANCH_index > 3)
_imp_signalf(6, 3, BRANCH_index, "->BRANCH(%d) out of bounds at %s line %d", BRANCH_index, BRANCH_file,
BRANCH_line);
goto *BRANCH[BRANCH_index - 1];
} /* end of QR */
// 1009 %end
return;
} /* end of IMPDIGVEC */
// 1010 !*END*IMPDIGVEC
// 1011
// 1012 !*TOP*IMPDKUSCR
// 1013 %external %routine imp dkuscr(%routine aux(%long %real %array %name f, y, %long %real x), %long %real x0, x1, e,
void IMPDKUSCR(void (*AUX)(double *F, double *Y, double X), double X0, double X1, double E, double *Y, int N, int K) {
auto void INTSTEP(double *Y, double X, double H, int N, double *E, void (*AUX)(double *F, double *Y, double X));
// 1014+ %long %real %array %name y, %integer n, k)
// 1015 %routine %spec int step(%long %real %array %name y, %long %real x, h,
// 1016+ %integer n, %long %real %name e, %routine aux(%long %real %array %name f, y, %long %real x))
// 1017 %integer m, i, j, p, q, r
int M;
int I;
int J;
int P;
int Q;
int R;
// 1018 %long %real h, e1, e2, e3, x
double H;
double E1;
double E2;
double E3;
double X;
// 1019 %long %real %array z(1:n)
double Z[N - 1 + 1];
// 1020 m = 1
M = 1;
// 1021 nrl2:
U_1039:;
// 1022 %if m>=k %then ->nrl1
if (M < K) goto L_012E;
goto U_103A;
L_012E:;
// 1023 m = 2*m
M = 2 * M;
// 1024 ->nrl2
goto U_1039;
// 1025 nrl1:h = (x1-x0)/m
U_103A:;
H = (double)((X1 - X0) / M);
// 1026 p = 1; q = m; e2 = 0; e3 = 0; x = x0
P = 1;
Q = M;
E2 = 0;
E3 = 0;
X = X0;
// 1027 nrl8:r = 1
U_103B:;
R = 1;
// 1028 nrl10:%cycle i = p, 1, q
U_103C:;
// %for I = P, 1, Q %cycle
{
int _initial = P, _final = Q;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = P, 1, Q");
I = _initial - (1);
L_012F:
if (I == _final) goto L_0130;
I += 1;
// 1029 %cycle j = 1, 1, n
// %for J = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, N");
J = 0;
L_0131:
if (J == _final) goto L_0132;
J += 1;
// 1030 z(j) = y(j)
Z[J - 1] = [/*BUG*/ 0];
// 1031 %repeat
goto L_0131;
}
L_0132:;
// 1032 int step(y, x, h, n, e1, aux)
INTSTEP(Y, X, H, N, /*?*/ &E1, AUX);
// 1033 %if e1>e %then ->nrl3
if (E1 <= E) goto L_0133;
goto U_103D;
L_0133:;
// 1034 %if e1<1.000000000000000@-2*e %and r=1 %and i#q %then ->nrl4
if (E1 >= 1.000000000000000E-2 * E) goto L_0134;
if (R != 1) goto L_0134;
if (I == Q) goto L_0134;
goto U_103E;
L_0134:;
// 1035 nrl9: e2 = 0; e3 = 0
U_103F:;
E2 = 0;
E3 = 0;
// 1036 x = x+h
X = X + H;
// 1037 %repeat
goto L_012F;
}
L_0130:;
// 1038 %return
return;
// 1039 nrl3:
U_103D:;
// 1040 %if e2=0 %and e3=0 %then ->nrl5
if (E2 != 0) goto L_0135;
if (E3 != 0) goto L_0135;
goto U_1040;
L_0135:;
// 1041 %if e3=0 %then ->nrl6
if (E3 != 0) goto L_0136;
goto U_1001;
L_0136:;
// 1042 %if e1>e2 %and e1>e3 %then ->nrl7
if (E1 <= E2) goto L_0137;
if (E1 <= E3) goto L_0137;
goto U_1002;
L_0137:;
// 1043 nrl6:e3 = e2
U_1001:;
E3 = E2;
// 1044 nrl5:e2 = e1
U_1040:;
E2 = E1;
// 1045 p = 1;
P = 1;
// 1046 q = 2*q-2*i+2
Q = 2 * Q - 2 * I + 2;
// 1047 h = h/2
H = (double)(H / 2);
// 1048 %cycle j = 1, 1, n
// %for J = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, N");
J = 0;
L_0138:
if (J == _final) goto L_0139;
J += 1;
// 1049 y(j) = z(j)
[/*BUG*/ 0] = Z[J - 1];
// 1050 %repeat
goto L_0138;
}
L_0139:;
// 1051 ->nrl8
goto U_103B;
// 1052 nrl7:
U_1002:;
// 1053 e = 2*e1
E = 2 * E1;
// 1054 ->nrl9
goto U_103F;
// 1055 nrl4:
U_103E:;
// 1056 %if fracpt((q-i)/2)>0 %then ->nrl9
if (_imp_FRACPT((double)((Q - I) / 2)) <= 0) goto L_013A;
goto U_103F;
L_013A:;
// 1057 p = 1; q = int((q-i)/2); x = x+h
P = 1;
Q = _imp_INT((double)((Q - I) / 2));
X = X + H;
// 1058 h = 2*h
H = 2 * H;
// 1059 ->nrl8
goto U_103B;
// 1060 %routine int step(%long %real %array %name y, %long %real x, h,
void INTSTEP(double *Y, double X, double H, int N, double *E, void (*AUX)(double *F, double *Y, double X)) {
// 1061+ %integer n, %long %real %name e, %routine aux(%long %real %array %name f, y, %long %real x))
// 1062 %long %real %array a, b, c, d, q, r(1:n)
double A[N - 1 + 1];
double B[N - 1 + 1];
double C[N - 1 + 1];
double D[N - 1 + 1];
double Q[N - 1 + 1];
double R[N - 1 + 1];
// 1063 %long %real u, m; %integer k, i
double U;
double M;
int K;
int I;
// 1064 k = n
K = N;
// 1065 aux(y, a, x)
AUX(Y, A, X);
// 1066 %cycle i = 1, 1, k
// %for I = 1, 1, K %cycle
{
int _final = K;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, K");
I = 0;
L_013B:
if (I == _final) goto L_013C;
I += 1;
// 1067 a(i) = h*a(i)
A[I - 1] = H * A[I - 1];
// 1068 q(i) = -5.555555559999999@-2*a(i)
Q[I - 1] = -(5.555555559999999E-2 * A[I - 1]);
// 1069 r(i) = 7.916666666999999@-1*a(i)
R[I - 1] = 7.916666666999999E-1 * A[I - 1];
// 1070 y(i) = z(i)+2.222222221999999@-1*a(i)
[/*BUG*/ 0] = Z[I - 1] + 2.222222221999999E-1 * A[I - 1];
// 1071 %repeat
goto L_013B;
}
L_013C:;
// 1072 u = x+2.222222221999999@-1*h
U = X + 2.222222221999999E-1 * H;
// 1073 aux(y, b, u)
AUX(Y, B, U);
// 1074 %cycle i = 1, 1, k
// %for I = 1, 1, K %cycle
{
int _final = K;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, K");
I = 0;
L_013D:
if (I == _final) goto L_013E;
I += 1;
// 1075 b(i) = h*b(i)
B[I - 1] = H * B[I - 1];
// 1076 r(i) = r(i)-3.375000000000000@0*b(i)
R[I - 1] = R[I - 1] - 3.375000000000000 * B[I - 1];
// 1077 y(i) = z(i)+8.333333333999993@-2*a(i)+2.500000000000000@-1*b(i)
[/*BUG*/ 0] = Z[I - 1] + 8.333333333999993E-2 * A[I - 1] + 2.500000000000000E-1 * B[I - 1];
// 1078 %repeat
goto L_013D;
}
L_013E:;
// 1079 u = x+3.333333332999999@-1*h
U = X + 3.333333332999999E-1 * H;
// 1080 aux(y, c, u)
AUX(Y, C, U);
// 1081 %cycle i = 1, 1, k
// %for I = 1, 1, K %cycle
{
int _final = K;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, K");
I = 0;
L_013F:
if (I == _final) goto L_0140;
I += 1;
// 1082 c(i) = h*c(i)
C[I - 1] = H * C[I - 1];
// 1083 q(i) = q(i)+1.588235293999999@-1*c(i)
Q[I - 1] = Q[I - 1] + 1.588235293999999E-1 * C[I - 1];
// 1084 r(i) = r(i)+2.850000000000000@0*c(i)
R[I - 1] = R[I - 1] + 2.850000000000000 * C[I - 1];
// 1085 y(i) = z(i)+5.390625000000000@-1*a(i)-1.898437499999999@0*b(i)+2 %c
[/*BUG*/ 0] =
Z[I - 1] + 5.390625000000000E-1 * A[I - 1] - 1.898437499999999 * B[I - 1] + 2.109374999999999 * C[I - 1];
// 1086+ .109374999999999@0*c(i)
// 1087 %repeat
goto L_013F;
}
L_0140:;
// 1088 u = x+7.500000000000000@-1*h
U = X + 7.500000000000000E-1 * H;
// 1089 aux(y, d, u)
AUX(Y, D, U);
// 1090 %cycle i = 1, 1, k
// %for I = 1, 1, K %cycle
{
int _final = K;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, K");
I = 0;
L_0141:
if (I == _final) goto L_0142;
I += 1;
// 1091 d(i) = h*d(i)
D[I - 1] = H * D[I - 1];
// 1092 q(i) = q(i)-2.666666666999999@-1*d(i)
Q[I - 1] = Q[I - 1] - 2.666666666999999E-1 * D[I - 1];
// 1093 r(i) = r(i)-2.666666666999999@-1*d(i)
R[I - 1] = R[I - 1] - 2.666666666999999E-1 * D[I - 1];
// 1094 y(i) = z(i)-3.104999999999999@-1*a(i)+1.822499999999999@0*b(i)-1 %c
[/*BUG*/ 0] = Z[I - 1] - 3.104999999999999E-1 * A[I - 1] + 1.822499999999999 * B[I - 1] -
1.101600000000000 * C[I - 1] + 4.895999999999999E-1 * D[I - 1];
// 1095+ .101600000000000@0*c(i)+4.895999999999999@-1*d(i)
// 1096 %repeat
goto L_0141;
}
L_0142:;
// 1097 u = x+9.000000000000000@-1*h
U = X + 9.000000000000000E-1 * H;
// 1098 aux(y, b, u)
AUX(Y, B, U);
// 1099 e = 0
*E = 0;
// 1100 %cycle i = 1, 1, k
// %for I = 1, 1, K %cycle
{
int _final = K;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, K");
I = 0;
L_0143:
if (I == _final) goto L_0144;
I += 1;
// 1101 b(i) = h*b(i)
B[I - 1] = H * B[I - 1];
// 1102 q(i) = q(i)+1.633986927999999@-1*b(i)
Q[I - 1] = Q[I - 1] + 1.633986927999999E-1 * B[I - 1];
// 1103 y(i) = z(i)+1.049382715999999@-1*a(i)+4.764705881999999@-1*c(i)+2 %c
[/*BUG*/ 0] = Z[I - 1] + 1.049382715999999E-1 * A[I - 1] + 4.764705881999999E-1 * C[I - 1] +
2.370370370000000E-1 * D[I - 1] + 1.815541030999999E-1 * B[I - 1];
// 1104+ .370370370000000@-1*d(i)+1.815541030999999@-1*b(i)
// 1105 %if (d(i)-a(i))=0 %then ->nrl1
if (D[I - 1] - A[I - 1] != 0) goto L_0145;
goto U_1019;
L_0145:;
// 1106 m = -(q(i)*(i))/(d(i)-a(i))
M = -((double)((Q[I - 1] * I) / (D[I - 1] - A[I - 1])));
// 1107 %if mod(m)<=e %then ->nrl1
if (_imp_MOD(M) > *E) goto L_0146;
goto U_1019;
L_0146:;
// 1108 e = mod(m)
*E = _imp_MOD(M);
// 1109 nrl1: %repeat
U_1019:;
goto L_0143;
}
L_0144:;
// 1110 %end
return;
} /* end of INTSTEP */
// 1111 %end
return;
} /* end of IMPDKUSCR */
// ?NRL10 unused
// 1112 !*END*IMPDKUSCR
// 1113
// 1114 !*TOP*IMPDSCRCH
// 1115 %external %routine imp dscrch(%long %real %function p(%long %real x), %long %real %function f(%long %real x), %long %real x,
void IMPDSCRCH(double (*P)(double X), double (*F)(double X), double X, double *Y, double *A, int N, int *IND) {
auto void COPIV(double *A, double *B, int N, int M, double EPS, double *DET, int *SING);
// 1116+ %long %real %name y, %long %real %array %name a, %integer n,
// 1117+ %integer %name ind)
// 1118 %routine %spec copiv(%long %real %array (2) %name a, b, %integer n, m,
// 1119+ %long %real eps, %long %real %name det, %integer %name sing)
// 1120 %long %real pi, val, th, sum
double PI;
double VAL;
double TH;
double SUM;
// 1121 %integer k, r, s, d
int K;
int R;
int S;
int D;
// 1122 %long %real %array pv(0:2*n+1), fv(0:n+1), aux(0:2*n+1, 1:1), ch(0:2*n+1,
double PV[2 * N + 1 + 1];
double FV[N + 1 + 1];
double AUX[2 * N + 1 + 1][1];
double CH[2 * N + 1 + 1][2 * N + 1 + 1];
auto void EXPAND(int D, double *C, double (*FCT)(double X));
// 1123+ 0:2*n+1)
// 1124 %routine %spec expand(%integer d, %long %real %array %name c,
// 1125+ %long %real %function fct(%long %real x))
// 1126 %if ind#0 %then ->nrl1
if (*IND == 0) goto L_0147;
goto U_103B;
L_0147:;
// 1127 %if n<0 %then ->nrl9
if (N >= 0) goto L_0148;
goto U_103C;
L_0148:;
// 1128 %if n=0 %then ->nrl7
if (N != 0) goto L_0149;
goto U_103D;
L_0149:;
// 1129 %if n=1 %then ->nrl8
if (N != 1) goto L_014A;
goto U_103E;
L_014A:;
// 1130 pi = 4*x-2
PI = 4 * X - 2;
// 1131 sum = a(n)
SUM = [/*BUG*/ 0];
// 1132 val = a(n-1)+pi*sum
VAL = [/*BUG*/ 0] + PI * SUM;
// 1133 %cycle k = n-2, -1, 0
// %for K = N - 2, -1, 0 %cycle
{
int _initial = N - 2;
if (invalid_for_loop(_initial, -1, 0)) _imp_signal(5, 1, 0, "Illegal cycle K = N - 2, -1, 0");
K = _initial - (-1);
L_014B:
if (K == 0) goto L_014C;
K -= 1;
// 1134 th = sum
TH = SUM;
// 1135 sum = val
SUM = VAL;
// 1136 val = a(k)-th+pi*sum
VAL = [/*BUG*/ 0] - TH + PI * SUM;
// 1137 %repeat
goto L_014B;
}
L_014C:;
// 1138 y = (val-th)/2
*Y = (double)((VAL - TH) / 2);
// 1139 %return
return;
// 1140 nrl7:y = a(0)/2
U_103D:;
*Y = (double)([/*BUG*/ 0] / 2);
// 1141 %return
return;
// 1142 nrl8:y = a(0)/2+(2*x-1)*a(1)
U_103E:;
*Y = (double)([/*BUG*/ 0] / 2) + (2 * X - 1) * [/*BUG*/ 0];
// 1143 %return
return;
// 1144 nrl9:ind = 1
U_103C:;
*IND = 1;
// 1145 ->nrl10
goto U_103F;
// 1146 nrl1:
U_103B:;
// 1147 %if n<2 %then ->nrl10
if (N >= 2) goto L_014D;
goto U_103F;
L_014D:;
// 1148 pi = 3.14159265359/2
PI = (double)(3.14159265359 / 2);
// 1149 d = 2*n+1
D = 2 * N + 1;
// 1150 %cycle r = 0, 1, d
// %for R = 0, 1, D %cycle
{
int _final = D;
if (invalid_for_loop(0, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle R = 0, 1, D");
R = -1;
L_014E:
if (R == _final) goto L_014F;
R += 1;
// 1151 %cycle k = 0, 1, d
// %for K = 0, 1, D %cycle
{
int _final = D;
if (invalid_for_loop(0, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 0, 1, D");
K = -1;
L_0150:
if (K == _final) goto L_0151;
K += 1;
// 1152 ch(r, k) = 0
CH[K] = 0;
// 1153 %repeat
goto L_0150;
}
L_0151:;
// 1154 %repeat
goto L_014E;
}
L_014F:;
// 1155 ch(0, 0) = 1
CH[0] = 1;
// 1156 ch(1, 1) = 2
CH[1] = 2;
// 1157 ch(0, 1) = -1
CH[1] = -1;
// 1158 %cycle r = 2, 1, d
// %for R = 2, 1, D %cycle
{
int _final = D;
if (invalid_for_loop(2, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle R = 2, 1, D");
R = 1;
L_0152:
if (R == _final) goto L_0153;
R += 1;
// 1159 %cycle k = 0, 1, r
// %for K = 0, 1, R %cycle
{
int _final = R;
if (invalid_for_loop(0, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 0, 1, R");
K = -1;
L_0154:
if (K == _final) goto L_0155;
K += 1;
// 1160 ch(k, r) = -ch(k, r-2)-2*ch(k, r-1)
CH[R] = -CH[R - 2] - 2 * CH[R - 1];
// 1161 %if k=0 %then ->nrl2
if (K != 0) goto L_0156;
goto U_1040;
L_0156:;
// 1162 ch(k, r) = ch(k, r)+4*ch(k-1, r-1)
CH[R] = CH[R] + 4 * CH[R - 1];
// 1163 nrl2: %repeat
U_1040:;
goto L_0154;
}
L_0155:;
// 1164 %repeat
goto L_0152;
}
L_0153:;
// 1165 expand(d, pv, p)
EXPAND(D, PV, P);
// 1166 expand(n+1, fv, f)
EXPAND(N + 1, FV, F);
// 1167 a(-1) = (mod(pv(d)))
[/*BUG*/ 0] = _imp_MOD(PV[D]);
// 1168 a(-2) = (mod(fv(n+1)))
[/*BUG*/ 0] = _imp_MOD(FV[N + 1]);
// 1169 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_0157:
if (K == _final) goto L_0158;
K += 1;
// 1170 val = ch(k, k)
VAL = CH[K];
// 1171 %cycle r = k-1, -1, 0
// %for R = K - 1, -1, 0 %cycle
{
int _initial = K - 1;
if (invalid_for_loop(_initial, -1, 0)) _imp_signal(5, 1, 0, "Illegal cycle R = K - 1, -1, 0");
R = _initial - (-1);
L_0159:
if (R == 0) goto L_015A;
R -= 1;
// 1172 val = val*x+ch(r, k)
VAL = VAL * X + CH[K];
// 1173 %repeat
goto L_0159;
}
L_015A:;
// 1174 a(k) = val
[/*BUG*/ 0] = VAL;
// 1175 %repeat
goto L_0157;
}
L_0158:;
// 1176 %cycle r = 1, 1, n
// %for R = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle R = 1, 1, N");
R = 0;
L_015B:
if (R == _final) goto L_015C;
R += 1;
// 1177 sum = pv(r-1)-pv(r+1)
SUM = PV[R - 1] - PV[R + 1];
// 1178 aux(r, 1) = 2*(fv(r-1)-fv(r+1)-sum*y)
AUX[1] = 2 * (FV[R - 1] - FV[R + 1] - SUM * *Y);
// 1179 %cycle s = 1, 1, n
// %for S = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle S = 1, 1, N");
S = 0;
L_015D:
if (S == _final) goto L_015E;
S += 1;
// 1180 %if r=s %then k = 1
if (R != S) goto L_015F;
K = 1;
L_015F:;
// 1181 %if r#s %then k = 0
if (R == S) goto L_0160;
K = 0;
L_0160:;
// 1182 pi = 8*k*r-2*sum*a(s)
PI = 8 * K * R - 2 * SUM * [/*BUG*/ 0];
// 1183 th = pv(r+s-1)-pv(r+s+1)
TH = PV[R + S - 1] - PV[R + S + 1];
// 1184 %if r=s %then ->nrl6
if (R != S) goto L_0161;
goto U_1001;
L_0161:;
// 1185 k = (imod(r-s))
K = _imp_IMOD(R - S);
// 1186 val = pv(k-1)-pv(k+1)
VAL = PV[K - 1] - PV[K + 1];
// 1187 %if r<s %then val = -val
if (R >= S) goto L_0162;
VAL = -VAL;
L_0162:;
// 1188 th = th+val
TH = TH + VAL;
// 1189 nrl6: ch(r, s) = pi+th
U_1001:;
CH[S] = PI + TH;
// 1190 %repeat
goto L_015D;
}
L_015E:;
// 1191 %repeat
goto L_015B;
}
L_015C:;
// 1192 copiv(ch, aux, n, 1, 1@-10, th, ind)
COPIV(CH, AUX, N, 1, 1E-10, /*?*/ &TH, IND);
// 1193 %if ind=1 %then %return
if (*IND != 1) goto L_0163;
return;
L_0163:;
// 1194 sum = y
SUM = *Y;
// 1195 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_0164:
if (K == _final) goto L_0165;
K += 1;
// 1196 sum = sum-a(k)*aux(k, 1)
SUM = SUM - [/*BUG*/ 0] * AUX[1];
// 1197 a(k) = aux(k, 1)
[/*BUG*/ 0] = AUX[1];
// 1198 %repeat
goto L_0164;
}
L_0165:;
// 1199 a(0) = 2*sum
[/*BUG*/ 0] = 2 * SUM;
// 1200 %return
return;
// 1201 nrl10:
U_103F:;
// 1202 printstring("
_imp_PRINTSTRING(_imp_str_literal("\nN TOO SMALL, NAMELY N = "));
// 1203+ N TOO SMALL, NAMELY N = ")
// 1204 print(n, 3, 0)
_imp_PRINT(N, 3, 0);
// 1205 printstring(" THIS ERROR OCCURRED IN ROUTINE ")
_imp_PRINTSTRING(_imp_str_literal(" THIS ERROR OCCURRED IN ROUTINE "));
// 1206 printstring("== IMP DSCRSH ==
_imp_PRINTSTRING(_imp_str_literal("== IMP DSCRSH ==\n"));
// 1207+ ")
// 1208 %routine expand(%integer d, %long %real %array %name c,
void EXPAND(int D, double *C, double (*FCT)(double X)) {
// 1209+ %long %real %function fct(%long %real x))
// 1210 %long %real %array z(0:2*n+1)
double Z[2 * N + 1 + 1];
// 1211 %cycle k = 0, 1, d
// %for K = 0, 1, D %cycle
{
int _final = D;
if (invalid_for_loop(0, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 0, 1, D");
K = -1;
L_0166:
if (K == _final) goto L_0167;
K += 1;
// ?Non-local
// 1212 th = (2*k+1)*pi/(d+1)
TH = (double)(((2 * K + 1) * PI) / (D + 1));
// 1213 th = cos(th/2)^2
TH = _imp_REXP(_imp_COS((double)(TH / 2)), 2);
// 1214 z(k) = th
Z[K] = TH;
// 1215 aux(k, 1) = fct(th)
AUX[1] = FCT(TH);
// 1216 %repeat
goto L_0166;
}
L_0167:;
// 1217 %cycle k = 0, 1, d
// %for K = 0, 1, D %cycle
{
int _final = D;
if (invalid_for_loop(0, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 0, 1, D");
K = -1;
L_0168:
if (K == _final) goto L_0169;
K += 1;
// ?Non-local
// 1218 sum = 0
SUM = 0;
// 1219 %cycle r = 0, 1, d
// %for R = 0, 1, D %cycle
{
int _final = D;
if (invalid_for_loop(0, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle R = 0, 1, D");
R = -1;
L_016A:
if (R == _final) goto L_016B;
R += 1;
// ?Non-local
// 1220 val = ch(k, k)
VAL = CH[K];
// 1221 %if k=0 %then ->nrl104
if (K != 0) goto L_016C;
goto U_1006;
L_016C:;
// 1222 %cycle s = k-1, -1, 0
// %for S = K - 1, -1, 0 %cycle
{
int _initial = K - 1;
if (invalid_for_loop(_initial, -1, 0)) _imp_signal(5, 1, 0, "Illegal cycle S = K - 1, -1, 0");
S = _initial - (-1);
L_016D:
if (S == 0) goto L_016E;
S -= 1;
// ?Non-local
// 1223 val = val*z(r)+ch(s, k)
VAL = VAL * Z[R] + CH[K];
// 1224 %repeat
goto L_016D;
}
L_016E:;
// 1225 nrl104: sum = sum+val*aux(r, 1)
U_1006:;
SUM = SUM + VAL * AUX[1];
// 1226 %repeat
goto L_016A;
}
L_016B:;
// 1227 c(k) = 2*sum/(d+1)
[/*BUG*/ 0] = (double)((2 * SUM) / (D + 1));
// 1228 %repeat
goto L_0168;
}
L_0169:;
// 1229 s = -1
S = -1;
// 1230 %cycle k = 0, 1, d
// %for K = 0, 1, D %cycle
{
int _final = D;
if (invalid_for_loop(0, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 0, 1, D");
K = -1;
L_016F:
if (K == _final) goto L_0170;
K += 1;
// ?Non-local
// 1231 %if s>=0 %then ->nrl102
if (S < 0) goto L_0171;
goto U_1000;
L_0171:;
// 1232 %if (mod(c(k)))<=4.999999999999999@-5 %then s = 0
if (_imp_MOD([/*BUG*/ 0]) > 4.999999999999999E-5) goto L_0172;
S = 0;
L_0172:;
// 1233 ->nrl101
goto U_1007;
// 1234 nrl102: s = k
U_1000:;
S = K;
// 1235 %if (mod(c(k)))<=(mod(c(k-1))) %then ->nrl101
if (_imp_MOD([/*BUG*/ 0]) > _imp_MOD([/*BUG*/ 0])) goto L_0173;
goto U_1007;
L_0173:;
// 1236 %if (mod(c(k)))<=4.999999999999999@-5 %then ->nrl103
if (_imp_MOD([/*BUG*/ 0]) > 4.999999999999999E-5) goto L_0174;
goto U_1008;
L_0174:;
// 1237 s = -1
S = -1;
// 1238 nrl101:%repeat
U_1007:;
goto L_016F;
}
L_0170:;
// 1239 s = -1
S = -1;
// 1240 nrl103:%if s<=0 %then %return
U_1008:;
if (S > 0) goto L_0175;
return;
L_0175:;
// 1241 %cycle k = s, 1, d
// %for K = S, 1, D %cycle
{
int _initial = S, _final = D;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = S, 1, D");
K = _initial - (1);
L_0176:
if (K == _final) goto L_0177;
K += 1;
// ?Non-local
// 1242 c(k) = 0
[/*BUG*/ 0] = 0;
// 1243 %repeat
goto L_0176;
}
L_0177:;
// 1244 %end
return;
} /* end of EXPAND */
// 1245 %routine copiv(%long %real %array (2) %name a, b, %integer n, m,
void COPIV(double *A, double *B, int N, int M, double EPS, double *DET, int *SING) {
// 1246+ %long %real eps, %long %real %name det, %integer %name sing)
// 1247 %long %real c, sum
double C;
double SUM;
// 1248 %integer i, j, k, ri, ci, con
int I;
int J;
int K;
int RI;
int CI;
int CON;
// 1249
// 1250 !========================================================================
// 1251 ! FILE CORRUPTION! Guessed missing lines here:
// 1252
// 1253 %integerarray s(1:n)
int S[N - 1 + 1];
// 1254 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0178:
if (I == _final) goto L_0179;
I += 1;
// 1255 s(i) = i
S[I - 1] = I;
// 1256 %repeat
goto L_0178;
}
L_0179:;
// 1257
// 1258 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_017A:
if (I == _final) goto L_017B;
I += 1;
// 1259 ! Complete pivoting: search for maximum element in the submatrix
// 1260 !========================================================================
// 1261
// 1262
// 1263 %cycle j = i+1, 1, n
// %for J = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = I + 1, 1, N");
J = _initial - (1);
L_017C:
if (J == _final) goto L_017D;
J += 1;
// 1264 c = a(j, i)/a(i, i)
C = (double)([/*BUG*/ 0] / [/*BUG*/ 0]);
// 1265 %cycle k = i+1, 1, n
// %for K = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = I + 1, 1, N");
K = _initial - (1);
L_017E:
if (K == _final) goto L_017F;
K += 1;
// 1266 a(j, k) = a(j, k)-c*a(i, k)
[/*BUG*/ 0] = [/*BUG*/ 0] - C * [/*BUG*/ 0];
// 1267 %repeat
goto L_017E;
}
L_017F:;
// 1268 %cycle k = 1, 1, m
// %for K = 1, 1, M %cycle
{
int _final = M;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, M");
K = 0;
L_0180:
if (K == _final) goto L_0181;
K += 1;
// 1269 b(j, k) = b(j, k)-c*b(i, k)
[/*BUG*/ 0] = [/*BUG*/ 0] - C * [/*BUG*/ 0];
// 1270 %repeat
goto L_0180;
}
L_0181:;
// 1271 a(j, i) = c
[/*BUG*/ 0] = C;
// 1272 %repeat
goto L_017C;
}
L_017D:;
// 1273 %repeat
goto L_017A;
}
L_017B:;
// 1274 nrl4:%cycle k = 1, 1, m
U_1018:;
// %for K = 1, 1, M %cycle
{
int _final = M;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, M");
K = 0;
L_0182:
if (K == _final) goto L_0183;
K += 1;
// 1275 b(n, k) = b(n, k)/a(n, n)
[/*BUG*/ 0] = (double)([/*BUG*/ 0] / [/*BUG*/ 0]);
// 1276 %repeat
goto L_0182;
}
L_0183:;
// 1277 %cycle i = n-1, -1, 1
// %for I = N - 1, -1, 1 %cycle
{
int _initial = N - 1;
if (invalid_for_loop(_initial, -1, 1)) _imp_signal(5, 1, 0, "Illegal cycle I = N - 1, -1, 1");
I = _initial - (-1);
L_0184:
if (I == 1) goto L_0185;
I -= 1;
// 1278 %cycle j = 1, 1, m
// %for J = 1, 1, M %cycle
{
int _final = M;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, M");
J = 0;
L_0186:
if (J == _final) goto L_0187;
J += 1;
// 1279 sum = 0
SUM = 0;
// 1280 %cycle k = i+1, 1, n
// %for K = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = I + 1, 1, N");
K = _initial - (1);
L_0188:
if (K == _final) goto L_0189;
K += 1;
// 1281 sum = sum+a(i, k)*b(k, j)
SUM = SUM + [/*BUG*/ 0] * [/*BUG*/ 0];
// 1282 %repeat
goto L_0188;
}
L_0189:;
// 1283 b(i, j) = (b(i, j)-sum)/a(i, i)
[/*BUG*/ 0] = (double)(([/*BUG*/ 0] - SUM) / [/*BUG*/ 0]);
// 1284 %repeat
goto L_0186;
}
L_0187:;
// 1285 %repeat
goto L_0184;
}
L_0185:;
// 1286 c = 1
C = 1;
// 1287 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_018A:
if (I == _final) goto L_018B;
I += 1;
// 1288 c = c*a(i, i)
C = C * [/*BUG*/ 0];
// 1289 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_018C:
if (K == _final) goto L_018D;
K += 1;
// 1290 %if s(k)=i %then ->nrl5
if (S[K - 1] != I) goto L_018E;
goto U_1019;
L_018E:;
// 1291 %repeat
goto L_018C;
}
L_018D:;
// 1292 nrl5: %cycle j = 1, 1, m
U_1019:;
// %for J = 1, 1, M %cycle
{
int _final = M;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, M");
J = 0;
L_018F:
if (J == _final) goto L_0190;
J += 1;
// 1293 sum = b(k, j)
SUM = [/*BUG*/ 0];
// 1294 b(k, j) = b(i, j)
[/*BUG*/ 0] = [/*BUG*/ 0];
// 1295 b(i, j) = sum
[/*BUG*/ 0] = SUM;
// 1296 %repeat
goto L_018F;
}
L_0190:;
// 1297 con = s(k)
CON = S[K - 1];
// 1298 s(k) = s(i)
S[K - 1] = S[I - 1];
// 1299 s(i) = con
S[I - 1] = CON;
// 1300 %repeat
goto L_018A;
}
L_018B:;
// 1301 sing = 0
*SING = 0;
// 1302 det = (mod(c))
*DET = _imp_MOD(C);
// 1303 %return
return;
// 1304 nrl10:
U_101A:;
// 1305 printstring("
_imp_PRINTSTRING(_imp_str_literal("\nN OR M TOO SMALL, NAMELY N = "));
// 1306+ N OR M TOO SMALL, NAMELY N = ")
// 1307 print(n, 3, 0)
_imp_PRINT(N, 3, 0);
// 1308 printstring(" M = ")
_imp_PRINTSTRING(_imp_str_literal(" M = "));
// 1309 print(m, 3, 0)
_imp_PRINT(M, 3, 0);
// 1310 printstring("
_imp_PRINTSTRING(_imp_str_literal("\nTHIS ERROR OCCURRED IN ROUTINE "));
// 1311+ THIS ERROR OCCURRED IN ROUTINE ")
// 1312 printstring("== COPIV ==
_imp_PRINTSTRING(_imp_str_literal("== COPIV ==\n"));
// 1313+ ")
// 1314 %end
return;
} /* end of COPIV */
// ?NRL10 unused
// ?NRL4 unused
// ?CI unused
// ?RI unused
// ?EPS unused
// 1315 %end
return;
} /* end of IMPDSCRCH */
// 1316 !*END*IMPDSCRCH
// 1317
// 1318 !*TOP*IMPFLEPO
// 1319 %external %routine imp flepo(%routine funct(%real %name f, %real %array %name x, g), %integer n, limit, iprin,
void IMPFLEPO(void (*FUNCT)(float *F, float *X, float *G), int N, int LIMIT, int IPRIN, int *CONV, float EPS, float *F,
float *X, float *H) {
// 1320+ %integer %name conv, %real eps, %real %name f, %real %array %name x, h)
// 1321 %real oldf, sg, ghg, z, step, ita, fa, fb, ga, gb, w, lambda
float OLDF;
float SG;
float GHG;
float Z;
float STEP;
float ITA;
float FA;
float FB;
float GA;
float GB;
float W;
float LAMBDA;
// 1322 %integer i, j, k, count, fnevalns
int I;
int J;
int K;
int COUNT;
int FNEVALNS;
// 1323 %real %array g, s, gamma, sigma, y(1:n)
float G[N - 1 + 1];
float S[N - 1 + 1];
float GAMMA[N - 1 + 1];
float SIGMA[N - 1 + 1];
float Y[N - 1 + 1];
// 1324 %real %function dot(%real %array %name a, b)
float DOT(float *A, float *B) {
// 1325 %integer i
int I;
// 1326 %real s
float S;
// 1327 s = 0
S = 0;
// 1328 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0191:
if (I == _final) goto L_0192;
I += 1;
// 1329 s = s+a(i)*b(i)
S = S + [/*BUG*/ 0] * [/*BUG*/ 0];
// 1330 %repeat
goto L_0191;
}
L_0192:;
// 1331 %result = s
return S;
// 1332 %end
} /* end of DOT */
// 1333 %real %function updot(%integer i, %real %array %name a, b)
float UPDOT(int I, float *A, float *B) {
// 1334 %comment MULTIPLYBYTHEITHROWOFTHESYMMETRICMATRIXA, %c
// 1335 !WHOSE UPPER TRIANGLE IS STORED BY ROWS.
// 1336 %integer j, k; %real s
int J;
int K;
float S;
// 1337 k = i; s = 0
K = I;
S = 0;
// 1338 ->nrl1 %if i=1
if (I != 1) goto L_0193;
goto U_1014;
L_0193:;
// 1339 %cycle j = 1, 1, i-1
// %for J = 1, 1, I - 1 %cycle
{
int _final = I - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, I - 1");
J = 0;
L_0194:
if (J == _final) goto L_0195;
J += 1;
// 1340 s = s+a(k)*b(j)
S = S + [/*BUG*/ 0] * [/*BUG*/ 0];
// 1341 k = k+n-j
K = K + N - J;
// 1342 %repeat; !STEPS TO DIAGONAL:NOW GO ALONG ROW.
goto L_0194;
}
L_0195:;
// 1343 nrl1: %cycle j = i, 1, n
U_1014:;
// %for J = I, 1, N %cycle
{
int _initial = I, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = I, 1, N");
J = _initial - (1);
L_0196:
if (J == _final) goto L_0197;
J += 1;
// 1344 s = s+a(k+j-i)*b(j)
S = S + [/*BUG*/ 0] * [/*BUG*/ 0];
// 1345 %repeat
goto L_0196;
}
L_0197:;
// 1346 %result = s
return S;
// 1347 %end; !OF FN UPDOT
} /* end of UPDOT */
// 1348 nrl1:
U_1007:;
// 1349 ->nrl2 %unless conv=1
if (*CONV == 1) goto L_0198;
goto U_1008;
L_0198:;
// 1350 k = 1
K = 1;
// 1351 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_0199:
if (I == _final) goto L_019A;
I += 1;
// 1352 h(k) = 1;
[/*BUG*/ 0] = 1;
// 1353 ->nrl30 %if i=n
if (I != N) goto L_019B;
goto U_1009;
L_019B:;
// 1354 %cycle j = 1, 1, n-i
// %for J = 1, 1, N - I %cycle
{
int _final = N - I;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, N - I");
J = 0;
L_019C:
if (J == _final) goto L_019D;
J += 1;
// 1355 h(k+j) = 0
[/*BUG*/ 0] = 0;
// 1356 %repeat
goto L_019C;
}
L_019D:;
// 1357 nrl30:k = k+n-i+1
U_1009:;
K = K + N - I + 1;
// 1358 %repeat; !FORMATION OF UNIT MATRIX IN H.
goto L_0199;
}
L_019A:;
// 1359 nrl2:conv = 1; step = 1; fnevalns = 0
U_1008:;
*CONV = 1;
STEP = 1;
FNEVALNS = 0;
// 1360 funct(f, x, g); fnevalns = fnevalns+1
FUNCT(F, X, G);
FNEVALNS = FNEVALNS + 1;
// 1361 oldf = f+1
OLDF = *F + 1;
// 1362 %cycle count = 1, 1, limit
// %for COUNT = 1, 1, LIMIT %cycle
{
int _final = LIMIT;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle COUNT = 1, 1, LIMIT");
COUNT = 0;
L_019E:
if (COUNT == _final) goto L_019F;
COUNT += 1;
// 1363 %return %unless oldf>f
if (OLDF > *F) goto L_01A0;
return;
L_01A0:;
// 1364 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_01A1:
if (I == _final) goto L_01A2;
I += 1;
// 1365 sigma(i) = x(i); gamma(i) = g(i); s(i) = -updot(i, h, g)
SIGMA[I - 1] = [/*BUG*/ 0];
GAMMA[I - 1] = G[I - 1];
S[I - 1] = -UPDOT(I, H, G);
// 1366 %repeat
goto L_01A1;
}
L_01A2:;
// 1367 fb = f; gb = dot(g, s)
FB = *F;
GB = DOT(G, S);
// 1368 ->nrl21 %unless gb>=0
if (GB >= 0) goto L_01A3;
goto U_100A;
L_01A3:;
// 1369 conv = -1; %return
*CONV = -1;
return;
// 1370 nrl21:oldf = f; ita = step
U_100A:;
OLDF = *F;
ITA = STEP;
// 1371 ! ITA IS DETERMINED FROM ITS VALUE AT THE PREVIOUS ITERATION (STEP)
// 1372 ! AND IS INCREASED OR DECREASED BY 4 WHERE NECESSARY. IT SHOULD
// 1373 ! TEND TO 1 AT THE MINIMUM.
// 1374 nrl4: fa = fb; ga = gb
U_100B:;
FA = FB;
GA = GB;
// 1375 %cycle i = 1, 1, n; x(i) = x(i)+ita*s(i); %repeat
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_01A4:
if (I == _final) goto L_01A5;
I += 1;
[/*BUG*/ 0] = [/*BUG*/ 0] + ITA * S[I - 1];
goto L_01A4;
}
L_01A5:;
// 1376 funct(f, x, g); fb = f; gb = dot(g, s); fnevalns = fnevalns+1
FUNCT(F, X, G);
FB = *F;
GB = DOT(G, S);
FNEVALNS = FNEVALNS + 1;
// 1377 ->nrl3 %unless (gb<0 %and fb<fa)
if (GB >= 0) goto L_00E5;
if (FB < FA) goto L_01A6;
L_00E5:;
goto U_1041;
L_01A6:;
// 1378 ita = 4*ita; step = 4*step;
ITA = 4 * ITA;
STEP = 4 * STEP;
// 1379 ->nrl4
goto U_100B;
// 1380 nrl3: %if ita<=0 %then %start; conv = -1; %return; %finish
U_1041:;
if (ITA > 0) goto L_01A7;
*CONV = -1;
return;
L_01A7:;
// 1381 %cycle i = 1, 1, n; y(i) = g(i); %repeat
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_01A8:
if (I == _final) goto L_01A9;
I += 1;
Y[I - 1] = G[I - 1];
goto L_01A8;
}
L_01A9:;
// 1382 z = 3*(fa-fb)/ita+ga+gb
Z = (double)((3 * (FA - FB)) / ITA) + GA + GB;
// 1383 w = sqrt(z*z-ga*gb)
W = _imp_SQRT(Z * Z - GA * GB);
// 1384 lambda = ita*(gb+w-z)/(gb-ga+2*w)
LAMBDA = (double)((ITA * (GB + W - Z)) / (GB - GA + 2 * W));
// 1385 %cycle i = 1, 1, n; x(i) = x(i)-lambda*s(i); %repeat
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_01AB:
if (I == _final) goto L_01AC;
I += 1;
[/*BUG*/ 0] = [/*BUG*/ 0] - LAMBDA * S[I - 1];
goto L_01AB;
}
L_01AC:;
// 1386 funct(f, x, g); fnevalns = fnevalns+1
FUNCT(F, X, G);
FNEVALNS = FNEVALNS + 1;
// 1387 ->nrl6 %unless (f>fa %or f>fb)
if (*F > FA) goto L_01AD;
if (*F > FB) goto L_01AD;
goto U_1014;
L_01AD:;
// 1388 step = step/4
STEP = (double)(STEP / 4);
// 1389 ->nrl5 %unless fb<fa
if (FB < FA) goto L_01AE;
goto U_1015;
L_01AE:;
// 1390 %cycle i = 1, 1, n; x(i) = x(i)+lambda*s(i); g(i) = y(i)
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_01AF:
if (I == _final) goto L_01B0;
I += 1;
[/*BUG*/ 0] = [/*BUG*/ 0] + LAMBDA * S[I - 1];
G[I - 1] = Y[I - 1];
// 1391 %repeat
goto L_01AF;
}
L_01B0:;
// 1392 f = fb;
*F = FB;
// 1393 ->nrl6
goto U_1014;
// 1394 nrl5: gb = dot(g, s)
U_1015:;
GB = DOT(G, S);
// 1395 ->nrl7 %unless (gb<0 %and count>n %and step<9.999999999999971@-7)
if (GB >= 0) goto L_01AA;
if (COUNT <= N) goto L_01AA;
if (STEP < 9.999999999999971E-7) goto L_01B1;
L_01AA:;
goto U_1016;
L_01B1:;
// 1396 conv = -1; %return
*CONV = -1;
return;
// 1397 nrl7: fb = f; ita = ita-lambda;
U_1016:;
FB = *F;
ITA = ITA - LAMBDA;
// 1398 ->nrl3
goto U_1041;
// 1399 nrl6: %cycle i = 1, 1, n
U_1014:;
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_01B2:
if (I == _final) goto L_01B3;
I += 1;
// 1400 sigma(i) = x(i)-sigma(i)
SIGMA[I - 1] = [/*BUG*/ 0] - SIGMA[I - 1];
// 1401 gamma(i) = g(i)-gamma(i)
GAMMA[I - 1] = G[I - 1] - GAMMA[I - 1];
// 1402 %repeat
goto L_01B2;
}
L_01B3:;
// 1403 sg = dot(sigma, gamma)
SG = DOT(SIGMA, GAMMA);
// 1404 ->nrl13 %unless count>=n
if (COUNT >= N) goto L_01B5;
goto U_1017;
L_01B5:;
// 1405 %return %if sqrt(dot(s, s))<eps %and sqrt(dot(sigma, sigma))<eps
if (_imp_SQRT(DOT(S, S)) >= EPS) goto L_01B6;
if (_imp_SQRT(DOT(SIGMA, SIGMA)) >= EPS) goto L_01B6;
return;
L_01B6:;
// 1406 %comment %end OFTESTFORVANISHINGDERIVATIVE.
// 1407 nrl13:%cycle i = 1, 1, n
U_1017:;
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_01B7:
if (I == _final) goto L_01B8;
I += 1;
// 1408 s(i) = updot(i, h, gamma)
S[I - 1] = UPDOT(I, H, GAMMA);
// 1409 %repeat
goto L_01B7;
}
L_01B8:;
// 1410 ghg = dot(s, gamma)
GHG = DOT(S, GAMMA);
// 1411 k = 1
K = 1;
// 1412 ->nrl8 %if mod(sg)<9.999999999999987@-31 %or %c
if (_imp_MOD(SG) < 9.999999999999987E-31) goto L_01B4;
if (_imp_MOD(GHG) >= 9.999999999999987E-31) goto L_01B9;
L_01B4:;
goto U_1018;
L_01B9:;
// 1413+ mod(ghg)<9.999999999999987@-31
// 1414 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_01BA:
if (I == _final) goto L_01BB;
I += 1;
// 1415 %cycle j = i, 1, n
// %for J = I, 1, N %cycle
{
int _initial = I, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = I, 1, N");
J = _initial - (1);
L_01BD:
if (J == _final) goto L_01BE;
J += 1;
// 1416 h(k) = h(k)+sigma(i)*sigma(j)/sg-s(i)*s(j)/ghg
[/*BUG*/ 0] =
[/*BUG*/ 0] + (double)((SIGMA[I - 1] * SIGMA[J - 1]) / SG) - (double)((S[I - 1] * S[J - 1]) / GHG);
// 1417 k = k+1; %repeat
K = K + 1;
goto L_01BD;
}
L_01BE:;
// 1418 %repeat
goto L_01BA;
}
L_01BB:;
// 1419 ->nrl20 %if iprin=0
if (IPRIN != 0) goto L_01BF;
goto U_1019;
L_01BF:;
// 1420 ->nrl20 %unless count//iprin*iprin=count
if (((long)COUNT / (long)IPRIN) * IPRIN == COUNT) goto L_01C0;
goto U_1019;
L_01C0:;
// 1421 nrl8: newlines(2);
U_1018:;
_imp_NEWLINES(2);
// 1422 printstring("
_imp_PRINTSTRING(_imp_str_literal("\nITERATION "));
// 1423+ ITERATION ")
// 1424 print(count, 2, 0)
_imp_PRINT(COUNT, 2, 0);
// 1425 printstring(":
_imp_PRINTSTRING(_imp_str_literal(":\n NO. OF FN EVALNS = "));
// 1426+ NO. OF FN EVALNS = ")
// 1427 write(fnevalns, 2)
_imp_WRITE(FNEVALNS, 2);
// 1428 newlines(2);
_imp_NEWLINES(2);
// 1429 printstring(" X G
_imp_PRINTSTRING(_imp_str_literal(" X G\n"));
// 1430+ ")
// 1431 %cycle i = 1, 1, n
// %for I = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N");
I = 0;
L_01C1:
if (I == _final) goto L_01C2;
I += 1;
// 1432 printfl(x(i), 6); spaces(4); printfl(g(i), 6); newline
_imp_PRINTFL([/*BUG*/ 0], 6);
_imp_SPACES(4);
_imp_PRINTFL(G[I - 1], 6);
_imp_NEWLINE();
// 1433 %repeat
goto L_01C1;
}
L_01C2:;
// 1434 printstring("
_imp_PRINTSTRING(_imp_str_literal("\n PRESENT MINIMUM VALUE OF FN ="));
// 1435+ PRESENT MINIMUM VALUE OF FN =")
// 1436 printfl(f, 6)
_imp_PRINTFL(*F, 6);
// 1437 newlines(2)
_imp_NEWLINES(2);
// 1438 nrl20:%repeat; ! COUNT LOOP; CONV=-1 ;
U_1019:;
goto L_019E;
}
L_019F:;
// 1439 printstring("
_imp_PRINTSTRING(_imp_str_literal("\nPROBLEM DID NOT CONVERGE IN SPECIFIED NO OF ITERATIONS\nROUTINE TERMINATED."));
// 1440+ PROBLEM DID NOT CONVERGE IN SPECIFIED NO OF ITERATIONS
// 1441+ ROUTINE TERMINATED.")
// 1442 newlines(2)
_imp_NEWLINES(2);
// 1443 %end;
return;
} /* end of IMPFLEPO */
// ?NRL1 unused
// 1444 ; !OF ROUTINE IMP FLEPO
// 1444 !
// 1445 !*END*IMPFLEPO
// 1446
// 1447 !*TOP*IMPGAHER
// 1448 %external %routine imp gaher(%real %function f(%real x), %real eps,
void IMPGAHER(float (*F)(float X), float EPS, float *ANTGRL, float *TEST) {
// 1449+ %real %name antgrl, test)
// 1450 %long %real i4, i6, i9, i12, i15, c
double I4;
double I6;
double I9;
double I12;
double I15;
double C;
// 1451 i4 = .813128354472@-1*(f(1.65068012388)+f(-1.65068012388))+.804914090005* %c
I4 = .813128354472E-1 * (F(1.65068012388) + F(-1.65068012388)) +
.804914090005 * (F(.524647623275) + F(-.524647623275));
// 1452+ (f(.524647623275)+f(-.524647623275))
// 1453 i6 = .453000990551@-2*(f(2.35060497367)+f(-2.35060497367))+.157067320323* %c
I6 = .453000990551E-2 * (F(2.35060497367) + F(-2.35060497367)) +
.157067320323 * (F(1.33584907401) + F(-1.33584907401)) + .724629595224 * (F(.436077411928) + F(-.436077411928));
// 1454+ (f(1.33584907401)+f(-1.33584907401))+.724629595224*(f(.436077411928)+f %c
// 1455+ (-.436077411928))
// 1456 c = mod(mod(i6)-mod(i4))
C = _imp_MOD(_imp_MOD(I6) - _imp_MOD(I4));
// 1457 %if (mod(i6))>1 %then c = c/(mod(i6))
if (_imp_MOD(I6) <= 1) goto L_01C3;
C = (double)(C / _imp_MOD(I6));
L_01C3:;
// 1458 ->nrl1 %unless c<eps
if (C < EPS) goto L_01C4;
goto U_1042;
L_01C4:;
// 1459 antgrl = i6
*ANTGRL = I6;
// 1460 test = c
*TEST = C;
// 1461 %return
return;
// 1462 nrl1:i9 = .396069772632@-4*(f(3.19099320178)+f(-3.19099320178))+ %c
U_1042:;
I9 = .396069772632E-4 * (F(3.19099320178) + F(-3.19099320178)) + .494362427553 -
2 * (F(2.26658058453) + F(-2.26658058453));
// 1463+ .494362427553-2*(f(2.26658058453)+f(-2.26658058453))
// 1464 i9 = i9+.884745273943-1*(f(1.46855328921)+f(-1.46855328921))+ %c
I9 = I9 + .884745273943 - 1 * (F(1.46855328921) + F(-1.46855328921)) +
.432651559002 * (F(.723551018752) + F(-.723551018752)) + .720235215606 * F(0);
// 1465+ .432651559002*(f(.723551018752)+f(-.723551018752))+.720235215606*f(0)
// 1466 c = mod(mod(i9)-mod(i6))
C = _imp_MOD(_imp_MOD(I9) - _imp_MOD(I6));
// 1467 %if (mod(i9))>1 %then c = c/(mod(i9))
if (_imp_MOD(I9) <= 1) goto L_01C5;
C = (double)(C / _imp_MOD(I9));
L_01C5:;
// 1468 ->nrl2 %unless c<eps
if (C < EPS) goto L_01C6;
goto U_1043;
L_01C6:;
// 1469 antgrl = i9
*ANTGRL = I9;
// 1470 test = c
*TEST = C;
// 1471 %return
return;
// 1472 nrl2:i12 = .265855168436@-6*(f(3.88972489787)+f(-3.88972489787))+ %c
U_1043:;
I12 = .265855168436E-6 * (F(3.88972489787) + F(-3.88972489787)) +
.857368704358E-4 * (F(3.02063702512) + F(-3.02063702512)) +
.390539058463E-2 * (F(2.27950708050) + F(-2.27950708050));
// 1473+ .857368704358@-4*(f(3.02063702512)+f(-3.02063702512))+.390539058463@-2* %c
// 1474+ (f(2.27950708050)+f(-2.27950708050))
// 1475 i12 = i12+.516079856158@-1*(f(1.59768263515)+f(-1.59768263515))+ %c
I12 = I12 + .516079856158E-1 * (F(1.59768263515) + F(-1.59768263515)) +
.260492310264 * (F(.947788391240) + F(-.947788391240)) + .570135236262 * (F(.314240376254) + F(-.314240376254));
// 1476+ .260492310264*(f(.947788391240)+f(-.947788391240))+.570135236262*(f( %c
// 1477+ .314240376254)+f(-.314240376254))
// 1478 c = mod(mod(i12)-mod(i9))
C = _imp_MOD(_imp_MOD(I12) - _imp_MOD(I9));
// 1479 %if (mod(i12))>1 %then c = c/(mod(i12))
if (_imp_MOD(I12) <= 1) goto L_01C7;
C = (double)(C / _imp_MOD(I12));
L_01C7:;
// 1480 ->nrl3 %unless c<eps
if (C < EPS) goto L_01C8;
goto U_1011;
L_01C8:;
// 1481 antgrl = i12
*ANTGRL = I12;
// 1482 test = c
*TEST = C;
// 1483 %return
return;
// 1484 nrl3:i15 = .152247580425@-8*(f(4.49999070731)+f(-4.49999070731))+ %c
U_1011:;
I15 = .152247580425E-8 * (F(4.49999070731) + F(-4.49999070731)) +
.105911554771E-5 * (F(3.66995037340) + F(-3.66995037340));
// 1485+ .105911554771@-5*(f(3.66995037340)+f(-3.66995037340))
// 1486 i15 = i15+.100004441232@-3*(f(2.96716692791)+f(-2.96716692791))+ %c
I15 = I15 + .100004441232E-3 * (F(2.96716692791) + F(-2.96716692791)) +
.277806884291E-2 * (F(2.32573248617) + F(-2.32573248617));
// 1487+ .277806884291@-2*(f(2.32573248617)+f(-2.32573248617))
// 1488 i15 = i15+.307800338725@-1*(f(1.71999257519)+f(-1.71999257519))+ %c
I15 = I15 + .307800338725E-1 * (F(1.71999257519) + F(-1.71999257519)) +
.158488915796 * (F(1.13611558521) + F(-1.13611558521));
// 1489+ .158488915796*(f(1.13611558521)+f(-1.13611558521))
// 1490 i15 = i15+.412028687499*(f(.565069583256)+f(-.565069583256))+ %c
I15 = I15 + .412028687499 * (F(.565069583256) + F(-.565069583256)) + .564100308726 * F(0);
// 1491+ .564100308726*f(0)
// 1492 c = mod(mod(i15)-mod(i12))
C = _imp_MOD(_imp_MOD(I15) - _imp_MOD(I12));
// 1493 %if (mod(i15))>1 %then c = c/(mod(i15))
if (_imp_MOD(I15) <= 1) goto L_01C9;
C = (double)(C / _imp_MOD(I15));
L_01C9:;
// 1494 antgrl = i15
*ANTGRL = I15;
// 1495 test = c
*TEST = C;
// 1496 %return %if c<eps
if (C >= EPS) goto L_01CA;
return;
L_01CA:;
// 1497 printstring("
_imp_PRINTSTRING(_imp_str_literal("\n\nMETHOD FAILED TO SATISFY GIVEN ACCURACY. LAST"));
// 1498+
// 1499+ METHOD FAILED TO SATISFY GIVEN ACCURACY. LAST")
// 1500 printstring(" ESTIMATE OF INTEGRAL IS ")
_imp_PRINTSTRING(_imp_str_literal(" ESTIMATE OF INTEGRAL IS "));
// 1501 printfl(antgrl, 10)
_imp_PRINTFL(*ANTGRL, 10);
// 1502 printstring("
_imp_PRINTSTRING(_imp_str_literal("\n\nTHE ERROR OF THIS ANSWER IS "));
// 1503+
// 1504+ THE ERROR OF THIS ANSWER IS ")
// 1505 printfl(test, 10)
_imp_PRINTFL(*TEST, 10);
// 1506 %return
return;
// 1507 %end; !OF IMP GAHER
} /* end of IMPGAHER */
// 1508 !*END*IMPGAHER
// 1509
// 1510 !*TOP*IMPGALAG
// 1511 %external %routine imp galag(%real %function f(%real x), %real eps,
void IMPGALAG(float (*F)(float X), float EPS, float *ANTGRL, float *TEST) {
// 1512+ %real %name antgrl, test)
// 1513 %long %real s4, s6, s9, s12, s15, c
double S4;
double S6;
double S9;
double S12;
double S15;
double C;
// 1514 s4 = 6.03154104342@-1*f(0.322547689619)+3.57418692438@-1*f(1.74576110116) %c
S4 = 6.03154104342E-1 * F(0.322547689619) + 3.57418692438E-1 * F(1.74576110116) +
3.88879085150E-2 * F(4.53662029692) + 5.39294705561E-4 * F(9.395070912301);
// 1515+ +3.88879085150@-2*f(4.53662029692)+5.39294705561@-4*f(9.395070912301)
// 1516 s6 = 4.58964673950@-1*f(.222846604179)+4.17000830772@-1*f(1.18893210167)
S6 = 4.58964673950E-1 * F(.222846604179) + 4.17000830772E-1 * F(1.18893210167);
// 1517 s6 = s6+1.13373382074@-1*f(2.99273632606)+1.03991974531@-2*f(5 %c
S6 = S6 + 1.13373382074E-1 * F(2.99273632606) + 1.03991974531E-2 * F(5.77514356910) +
2.61017202815E-4 * F(9.83746741838) + 8.98547906430E-7 * F(15.9828739806);
// 1518+ .77514356910)+2.61017202815@-4*f(9.83746741838)+8.98547906430@-7*f(15 %c
// 1519+ .9828739806)
// 1520 c = mod(mod(s6)-mod(s4))
C = _imp_MOD(_imp_MOD(S6) - _imp_MOD(S4));
// 1521 %if (mod(s6))>1 %then c = c/(mod(s6))
if (_imp_MOD(S6) <= 1) goto L_01CB;
C = (double)(C / _imp_MOD(S6));
L_01CB:;
// 1522 ->nrl1 %unless c<eps
if (C < EPS) goto L_01CC;
goto U_1043;
L_01CC:;
// 1523 antgrl = s6
*ANTGRL = S6;
// 1524 test = c
*TEST = C;
// 1525 %return
return;
// 1526 nrl1:s9 = 3.36126421798@-1*f(.152322227732)+4.11213980424@-1*f(.807220022742)
U_1043:;
S9 = 3.36126421798E-1 * F(.152322227732) + 4.11213980424E-1 * F(.807220022742);
// 1527 s9 = s9+1.99287525371@-1*f(2.00513515562)+4.74605627657@-2*f(3.78347397333)
S9 = S9 + 1.99287525371E-1 * F(2.00513515562) + 4.74605627657E-2 * F(3.78347397333);
// 1528 s9 = s9+5.59962661079@-3*f(6.20495677788)+3.05249767093@-4*f(9.37298525169)
S9 = S9 + 5.59962661079E-3 * F(6.20495677788) + 3.05249767093E-4 * F(9.37298525169);
// 1529 s9 = s9+6.59212302608@-6*f(13.4662369111)+4.11076933035@-8*f(18.8335977890)
S9 = S9 + 6.59212302608E-6 * F(13.4662369111) + 4.11076933035E-8 * F(18.8335977890);
// 1530 s9 = s9+3.29087403035@-11*f(26.3740718909)
S9 = S9 + 3.29087403035E-11 * F(26.3740718909);
// 1531 c = mod(mod(s9)-mod(s6))
C = _imp_MOD(_imp_MOD(S9) - _imp_MOD(S6));
// 1532 %if (mod(s9))>1 %then c = c/(mod(s9))
if (_imp_MOD(S9) <= 1) goto L_01CD;
C = (double)(C / _imp_MOD(S9));
L_01CD:;
// 1533 ->nrl2 %unless c<eps
if (C < EPS) goto L_01CE;
goto U_1011;
L_01CE:;
// 1534 antgrl = s9
*ANTGRL = S9;
// 1535 test = c
*TEST = C;
// 1536 %return
return;
// 1537 nrl2:s12 = 2.64731371055@-1*f(.115722117358)+3.77759275873@-1*f(.61175748452) %c
U_1011:;
S12 = 2.64731371055E-1 * F(.115722117358) + 3.77759275873E-1 * F(.61175748452) + 9.04492222117E-2 * F(2.83375133774);
// 1538+ +9.04492222117@-2*f(2.83375133774)
// 1539 s12 = s12+2.44082011320@-1*f(1.51261026977)
S12 = S12 + 2.44082011320E-1 * F(1.51261026977);
// 1540 s12 = s12+2.01023811546@-2*f(4.59922763942)+2.66397354187@-3*f(6 %c
S12 = S12 + 2.01023811546E-2 * F(4.59922763942) + 2.66397354187E-3 * F(6.84452545312);
// 1541+ .84452545312)
// 1542 s12 = s12+2.03231592663@-4*f(9.62131684246)+8.36505585682@-6*f(13 %c
S12 = S12 + 2.03231592663E-4 * F(9.62131684246) + 8.36505585682E-6 * F(13.0060549933);
// 1543+ .0060549933)
// 1544 s12 = s12+1.66849387654@-7*f(17.1168551874)+1.34239103051@-9*f(22 %c
S12 = S12 + 1.66849387654E-7 * F(17.1168551874) + 1.34239103051E-9 * F(22.1510903794);
// 1545+ .1510903794)
// 1546 s12 = s12+3.06160163503@-12*f(28.4879672510)+8.14807746743@-16*f(37 %c
S12 = S12 + 3.06160163503E-12 * F(28.4879672510) + 8.14807746743E-16 * F(37.0991210445);
// 1547+ .0991210445)
// 1548 c = mod(mod(s12)-mod(s9))
C = _imp_MOD(_imp_MOD(S12) - _imp_MOD(S9));
// 1549 %if (mod(s12))>1 %then c = c/(mod(s12))
if (_imp_MOD(S12) <= 1) goto L_01CF;
C = (double)(C / _imp_MOD(S12));
L_01CF:;
// 1550 ->nrl3 %unless c<eps
if (C < EPS) goto L_01D0;
goto U_1012;
L_01D0:;
// 1551 antgrl = s12
*ANTGRL = S12;
// 1552 test = c
*TEST = C;
// 1553 %return
return;
// 1554 nrl3:s15 = 2.18234885940@-1*f(.093307812017)+3.42210177923@-1*f(.492691740302)
U_1012:;
S15 = 2.18234885940E-1 * F(.093307812017) + 3.42210177923E-1 * F(.492691740302);
// 1555 s15 = s15+2.63027577942@-1*f(1.21559541207)+1.26425818106@-1*f(2 %c
S15 = S15 + 2.63027577942E-1 * F(1.21559541207) + 1.26425818106E-1 * F(2.26994952620);
// 1556+ .26994952620)
// 1557 s15 = s15+4.02068649210@-2*f(3.66762272175)+8.56387780361@-3*f(5 %c
S15 = S15 + 4.02068649210E-2 * F(3.66762272175) + 8.56387780361E-3 * F(5.42533662741);
// 1558+ .42533662741)
// 1559 s15 = s15+1.2124361472@-3*f(7.56591622661)+1.11674392344@-4*f(10.1202285680)
S15 = S15 + 1.2124361472E-3 * F(7.56591622661) + 1.11674392344E-4 * F(10.1202285680);
// 1560 s15 = s15+6.45992676202@-6*f(13.1302824822)+2.22631690710@-7*f(16 %c
S15 = S15 + 6.45992676202E-6 * F(13.1302824822) + 2.22631690710E-7 * F(16.6544077083);
// 1561+ .6544077083)
// 1562 s15 = s15+4.22743038498@-9*f(20.7764788994)+3.92189726704@-11*f(25 %c
S15 = S15 + 4.22743038498E-9 * F(20.7764788994) + 3.92189726704E-11 * F(25.6238942267);
// 1563+ .6238942267)
// 1564 s15 = s15+1.45651526407@-13*f(31.4075191697)+1.48302705111@-16*f(38 %c
S15 = S15 + 1.45651526407E-13 * F(31.4075191697) + 1.48302705111E-16 * F(38.5306833065);
// 1565+ .5306833065)
// 1566 s15 = s15+1.60059490621@-20*f(48.0260855727)
S15 = S15 + 1.60059490621E-20 * F(48.0260855727);
// 1567 c = mod(mod(s15)-mod(s12))
C = _imp_MOD(_imp_MOD(S15) - _imp_MOD(S12));
// 1568 %if (mod(s15))>1 %then c = c/(mod(s15))
if (_imp_MOD(S15) <= 1) goto L_01D1;
C = (double)(C / _imp_MOD(S15));
L_01D1:;
// 1569 antgrl = s15
*ANTGRL = S15;
// 1570 test = c
*TEST = C;
// 1571 ->nrl4 %unless c<eps
if (C < EPS) goto L_01D2;
goto U_1013;
L_01D2:;
// 1572 %return
return;
// 1573 nrl4:
U_1013:;
// 1574 printstring("
_imp_PRINTSTRING(_imp_str_literal("\n\nMETHOD FAILED. LAST APPROXIMATION TO VALUE OF INTEGRAL IS "));
// 1575+
// 1576+ METHOD FAILED. LAST APPROXIMATION TO VALUE OF INTEGRAL IS ")
// 1577 printfl(antgrl, 10)
_imp_PRINTFL(*ANTGRL, 10);
// 1578 printstring("
_imp_PRINTSTRING(_imp_str_literal("\n\nTHE ERROR OF THIS ANSWER IS "));
// 1579+
// 1580+ THE ERROR OF THIS ANSWER IS ")
// 1581 printfl(test, 10)
_imp_PRINTFL(*TEST, 10);
// 1582 %return
return;
// 1583 %end; !OF IMP GALAG
} /* end of IMPGALAG */
// 1584 !*END*IMPGALAG
// 1585
// 1586 !*TOP*IMPGALEG
// 1587 %external %routine imp galeg(%real %function f(%long %real x), %real xup, xlo,
void IMPGALEG(float (*F)(double X), float XUP, float XLO, float *TIN, float *EPS, int LIM) {
// 1588+ %real %name tin, eps, %integer lim)
// 1589 %long %real i5, i6, int5, int6, h, a, xu, xl
double I5;
double I6;
double INT5;
double INT6;
double H;
double A;
double XU;
double XL;
// 1590 %integer cntr, i
int CNTR;
int I;
// 1591 %routine comp(%long %real a, %long %real %name i5, i6, %real %function f(%long %real x))
void COMP(double A, double *I5, double *I6, float (*F)(double X)) {
// 1592 %long %real c
double C;
// 1593 c = .45308992*h
C = .45308992 * H;
// 1594 i5 = .11846344*(f(a+c)+f(a-c))
*I5 = .11846344 * (F(A + C) + F(A - C));
// 1595 c = .26923466*h
C = .26923466 * H;
// 1596 i5 = i5+0.23931434*(f(a+c)+f(a-c))
*I5 = *I5 + 0.23931434 * (F(A + C) + F(A - C));
// 1597 i5 = h*(i5+.28444444*f(a))
*I5 = H * (*I5 + .28444444 * F(A));
// 1598 c = .46623476*h
C = .46623476 * H;
// 1599 i6 = .85662246@-1*(f(a+c)+f(a-c))
*I6 = .85662246E-1 * (F(A + C) + F(A - C));
// 1600 c = .33060469*h
C = .33060469 * H;
// 1601 i6 = i6+.18038079*(f(a+c)+f(a-c))
*I6 = *I6 + .18038079 * (F(A + C) + F(A - C));
// 1602 c = .11930959*h
C = .11930959 * H;
// 1603 i6 = (i6+0.23395697*(f(a+c)+f(a-c)))*h
*I6 = (*I6 + 0.23395697 * (F(A + C) + F(A - C))) * H;
// 1604 %end; !END OF COMP
return;
} /* end of COMP */
// 1605 cntr = 0
CNTR = 0;
// 1606 a = .5*(xup+xlo)
A = .5 * (XUP + XLO);
// 1607 h = xup-xlo
H = XUP - XLO;
// 1608 comp(a, i5, i6, f)
COMP(A, /*?*/ &I5, /*?*/ &I6, F);
// 1609 int5 = i5
INT5 = I5;
// 1610 int6 = i6
INT6 = I6;
// 1611 nrl1:
U_103B:;
// 1612 ->nrl3 %unless mod(int6-int5)<eps
if (_imp_MOD(INT6 - INT5) < *EPS) goto L_01D3;
goto U_103C;
L_01D3:;
// 1613 tin = int6
*TIN = INT6;
// 1614 eps = mod(int6-int5)
*EPS = _imp_MOD(INT6 - INT5);
// 1615 %return
return;
// 1616 nrl3:cntr = cntr+1
U_103C:;
CNTR = CNTR + 1;
// 1617 ->nrl2 %if cntr<=lim
if (CNTR > LIM) goto L_01D4;
goto U_103D;
L_01D4:;
// 1618 newlines(3)
_imp_NEWLINES(3);
// 1619 printstring("
_imp_PRINTSTRING(_imp_str_literal("\nMETHOD FAILED: REQUIRED ACCURACY NOT ATTAINED AFTER"));
// 1620+ METHOD FAILED: REQUIRED ACCURACY NOT ATTAINED AFTER")
// 1621 write(lim, 5)
_imp_WRITE(LIM, 5);
// 1622 printstring(" BISECTIONS OF THE INTERVAL")
_imp_PRINTSTRING(_imp_str_literal(" BISECTIONS OF THE INTERVAL"));
// 1623 newlines(3)
_imp_NEWLINES(3);
// 1624 %return
return;
// 1625 nrl2:xu = xup
U_103D:;
XU = XUP;
// 1626 int5 = 0
INT5 = 0;
// 1627 int6 = 0
INT6 = 0;
// 1628 h = h/2
H = (double)(H / 2);
// 1629 %cycle i = 1, 1, 2^^cntr
// %for I = 1, 1, _imp_IEXP(2, CNTR) %cycle
{
int _final = _imp_IEXP(2, CNTR);
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, _imp_IEXP(2, CNTR)");
I = 0;
L_01D5:
if (I == _final) goto L_01D6;
I += 1;
// 1630 xl = xu-h
XL = XU - H;
// 1631 a = .5*(xu+xl)
A = .5 * (XU + XL);
// 1632 comp(a, i5, i6, f)
COMP(A, /*?*/ &I5, /*?*/ &I6, F);
// 1633 int5 = int5+i5
INT5 = INT5 + I5;
// 1634 int6 = int6+i6
INT6 = INT6 + I6;
// 1635 xu = xl
XU = XL;
// 1636 %repeat
goto L_01D5;
}
L_01D6:;
// 1637 ->nrl1
goto U_103B;
// 1638 %end; ! OF ROUTINE IMP GALEG
} /* end of IMPGALEG */
// 1639 !*END*IMPGALEG
// 1640
// 1641 !*TOP*IMPLOGFCT
// 1642 %external %real %function imp logfct(%integer i, j)
float IMPLOGFCT(int I, int J) {
// 1643 %integer n, k
int N;
int K;
// 1644 %real x, p
float X;
float P;
// 1645 %switch s(0:20)
int S_index;
int S_line;
char *S_file;
static const void *S[21] = {&&case_S_missing, &&case_S_missing, &&case_S_missing, &&case_S_missing, &&case_S_missing,
&&case_S_missing, &&case_S_missing, &&case_S_missing, &&case_S_missing, &&case_S_missing,
&&case_S_missing, &&case_S_missing, &&case_S_missing, &&case_S_missing, &&case_S_missing,
&&case_S_missing, &&case_S_missing, &&case_S_missing, &&case_S_missing, &&case_S_missing,
&&case_S_missing};
// 1646 %if j=1 %then p = 0.4342944819 %else p = 1.0
if (J != 1) goto L_01D7;
P = 0.4342944819;
goto L_01D8;
L_01D7:;
P = 1.0;
L_01D8:;
// 1647 n = i//50; %if n>20 %then n = 20
N = ((long)I / (long)50);
if (N <= 20) goto L_01D9;
N = 20;
L_01D9:;
// 1648 %if i<0 %then %start
if (I >= 0) goto L_01DA;
// 1649 printstring(" ARGUMENT I IS NEGATIVE. ")
_imp_PRINTSTRING(_imp_str_literal(" ARGUMENT I IS NEGATIVE. "));
// 1650 newline
_imp_NEWLINE();
// 1651 %monitor
/*_imp_monitor(__LINE__, __FILE__, PRETTY_FUNCTION);*/
// 1652 %stop
_imp_stop(0);
// 1653 %finish
L_01DA:;
// 1654 ->s(n)
S_line = __LINE__;
S_file = __FILE__;
S_index = N;
goto case_S_dispatch;
// 1655 s(0):x = 0;
case_S_0:;
X = 0;
// 1656 ->nrl1
goto U_1042;
// 1657 s(1):x = 148.4777669516;
case_S_1:;
X = 148.4777669516;
// 1658 ->nrl1
goto U_1042;
// 1659 s(2):x = 363.7393755532;
case_S_2:;
X = 363.7393755532;
// 1660 ->nrl1
goto U_1042;
// 1661 s(3):x = 605.0201058475;
case_S_3:;
X = 605.0201058475;
// 1662 ->nrl1
goto U_1042;
// 1663 s(4):x = 863.2319871918;
case_S_4:;
X = 863.2319871918;
// 1664 ->nrl1
goto U_1042;
// 1665 s(5):x = 1134.0452317854;
case_S_5:;
X = 1134.0452317854;
// 1666 ->nrl1
goto U_1042;
// 1667 s(6):x = 1414.058499460;
case_S_6:;
X = 1414.058499460;
// 1668 ->nrl1
goto U_1042;
// 1669 s(7):x = 1704.1247472880;
case_S_7:;
X = 1704.1247472880;
// 1670 ->nrl1
goto U_1042;
// 1671 s(8):x = 2000.5006979990;
case_S_8:;
X = 2000.5006979990;
// 1672 ->nrl1
goto U_1042;
// 1673 s(9):x = 2303.1351597674;
case_S_9:;
X = 2303.1351597674;
// 1674 ->nrl1
goto U_1042;
// 1675 s(10):x = 2611.3304584723;
case_S_10:;
X = 2611.3304584723;
// 1676 ->nrl1
goto U_1042;
// 1677 s(11):x = 2924.5291022234;
case_S_11:;
X = 2924.5291022234;
// 1678 ->nrl1
goto U_1042;
// 1679 s(12):x = 3242.2753354039;
case_S_12:;
X = 3242.2753354039;
// 1680 ->nrl1
goto U_1042;
// 1681 s(13):x = 3564.1895888475;
case_S_13:;
X = 3564.1895888475;
// 1682 ->nrl1
goto U_1042;
// 1683 s(14):x = 3889.9508323156;
case_S_14:;
X = 3889.9508323156;
// 1684 ->nrl1
goto U_1042;
// 1685 s(15):x = 4219.2839911876;
case_S_15:;
X = 4219.2839911876;
// 1686 ->nrl1
goto U_1042;
// 1687 s(16):x = 4551.9507307617;
case_S_16:;
X = 4551.9507307617;
// 1688 ->nrl1
goto U_1042;
// 1689 s(17):x = 4887.7425519021;
case_S_17:;
X = 4887.7425519021;
// 1690 ->nrl1
goto U_1042;
// 1691 s(18):x = 5226.4755155920;
case_S_18:;
X = 5226.4755155920;
// 1692 ->nrl1
goto U_1042;
// 1693 s(19):x = 5567.9861427925;
case_S_19:;
X = 5567.9861427925;
// 1694 ->nrl1
goto U_1042;
// 1695 s(20):x = 5912.1281786647
case_S_20:;
X = 5912.1281786647;
// 1696 nrl1:%if i=50*n %then %result = x*p
U_1042:;
if (I != 50 * N) goto L_01DC;
return X * P;
L_01DC:;
// 1697 %cycle k = (50*n)+1, 1, i
// %for K = 50 * N + 1, 1, I %cycle
{
int _initial = 50 * N + 1, _final = I;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 50 * N + 1, 1, I");
K = _initial - (1);
L_01DD:
if (K == _final) goto L_01DE;
K += 1;
// 1698 x = x+log(k)
X = X + _imp_LOG(K);
// 1699 %repeat
goto L_01DD;
}
L_01DE:;
// 1700 %result = x*p
return X * P;
// 1701 %end
case_S_dispatch:
if (S_index < 0 || S_index > 20)
_imp_signalf(6, 3, S_index, "->S(%d) out of bounds at %s line %d", S_index, S_file, S_line);
goto *S[S_index - 0];
case_S_missing:
_imp_signalf(8, 2, S_index, "switch jump to undefined label S(%d): at %s line %d", S_index, S_file, S_line);
} /* end of IMPLOGFCT */
// 1702 !*TOP*IMPDAPRO
// 1703 %external %routine imp dapro(%long %real %function f(%long %real x, %long %real %array %name v), %integer k0, ma, me, n,
void IMPDAPRO(double (*F)(double X, double *V), int K0, int MA, int ME, int N, int *M, double *S, double *X, double *Y,
double *V, double *TOL, double *SF) {
auto void CHOLESKY(int N, double *A, double *B, int *FAULT);
// 1704+ %integer %name m, %long %real %name s, %long %real %array %name x, y, v,
// 1705+ tol, sf)
// 1706 %routine %spec cholesky(%integer n, %long %real %array (2) %name a, %long %real %array %name b,
// 1707+ %integer %name fault)
// 1708 %integer i, k, l, loop, ind
int I;
int K;
int L;
int LOOP;
int IND;
auto double AF(int J, double X);
// 1709 %long %real %function %spec af(%integer j, %long %real x)
// 1710 %long %real %array dy(ma:me), b(1:n), c(1:n, 1:n)
double DY[ME - MA + 1];
double B[N - 1 + 1];
double C[N - 1 + 1][N - 1 + 1];
// 1711 %long %real z, z1
double Z;
double Z1;
// 1712 %if n<2 %then ->nrl100
if (N >= 2) goto L_01DF;
goto U_1004;
L_01DF:;
// 1713 loop = 0
LOOP = 0;
// 1714 nrl1:loop = loop+1
U_1005:;
LOOP = LOOP + 1;
// 1715 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_01E0:
if (K == _final) goto L_01E1;
K += 1;
// 1716 b(k) = 0
B[K - 1] = 0;
// 1717 %cycle l = 1, 1, n
// %for L = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle L = 1, 1, N");
L = 0;
L_01E2:
if (L == _final) goto L_01E3;
L += 1;
// 1718 c(k, l) = 0
C[L - 1] = 0;
// 1719 %repeat
goto L_01E2;
}
L_01E3:;
// 1720 %repeat
goto L_01E0;
}
L_01E1:;
// 1721 %cycle i = ma, 1, me
// %for I = MA, 1, ME %cycle
{
int _initial = MA, _final = ME;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = MA, 1, ME");
I = _initial - (1);
L_01E4:
if (I == _final) goto L_01E5;
I += 1;
// 1722 dy(i) = y(i)-f(x(i), v)
DY[I - MA] = [/*BUG*/ 0] - F([/*BUG*/ 0], V);
// 1723 %if k0=0 %then ->nrl101
if (K0 != 0) goto L_01E6;
goto U_1006;
L_01E6:;
// 1724 dy(i) = dy(i)-v(k0)
DY[I - MA] = DY[I - MA] - [/*BUG*/ 0];
// 1725 nrl101:%repeat
U_1006:;
goto L_01E4;
}
L_01E5:;
// 1726 %cycle i = ma, 1, me
// %for I = MA, 1, ME %cycle
{
int _initial = MA, _final = ME;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = MA, 1, ME");
I = _initial - (1);
L_01E7:
if (I == _final) goto L_01E8;
I += 1;
// 1727 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_01E9:
if (K == _final) goto L_01EA;
K += 1;
// 1728 z = af(k, x(i))
Z = AF(K, [/*BUG*/ 0]);
// 1729 b(k) = b(k)+dy(i)*z
B[K - 1] = B[K - 1] + DY[I - MA] * Z;
// 1730 %cycle l = k, 1, n
// %for L = K, 1, N %cycle
{
int _initial = K, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle L = K, 1, N");
L = _initial - (1);
L_01EB:
if (L == _final) goto L_01EC;
L += 1;
// 1731 z1 = z^2
Z1 = _imp_REXP(Z, 2);
// 1732 %if l=k %then ->nrl102
if (L != K) goto L_01ED;
goto U_1000;
L_01ED:;
// 1733 z1 = af(l, x(i))*z
Z1 = AF(L, [/*BUG*/ 0]) * Z;
// 1734 nrl102: c(l, k) = c(l, k)+z1
U_1000:;
C[K - 1] = C[K - 1] + Z1;
// 1735 %repeat
goto L_01EB;
}
L_01EC:;
// 1736 %repeat
goto L_01E9;
}
L_01EA:;
// 1737 %repeat
goto L_01E7;
}
L_01E8:;
// 1738 cholesky(n, c, b, ind)
CHOLESKY(N, C, B, /*?*/ &IND);
// 1739 %if ind=1 %then ->nrl2
if (IND != 1) goto L_01EE;
goto U_1007;
L_01EE:;
// 1740 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_01EF:
if (K == _final) goto L_01F0;
K += 1;
// 1741 v(k) = v(k)+b(k)
[/*BUG*/ 0] = [/*BUG*/ 0] + B[K - 1];
// 1742 %repeat
goto L_01EF;
}
L_01F0:;
// 1743 %cycle k = 1, 1, n
// %for K = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, N");
K = 0;
L_01F1:
if (K == _final) goto L_01F2;
K += 1;
// 1744 %if (mod(b(k)))>tol(k) %then ->nrl103
if (_imp_MOD(B[K - 1]) <= [/*BUG*/ 0]) goto L_01F3;
goto U_1008;
L_01F3:;
// 1745 %repeat
goto L_01F1;
}
L_01F2:;
// 1746 m = loop
*M = LOOP;
// 1747 ->nrl3
goto U_1009;
// 1748 nrl103:
U_1008:;
// 1749 %if loop<m %then ->nrl1
if (LOOP >= *M) goto L_01F4;
goto U_1005;
L_01F4:;
// 1750 nrl2:m = -loop-1
U_1007:;
*M = -LOOP - 1;
// 1751 nrl3:s = 0
U_1009:;
*S = 0;
// 1752 %cycle i = ma, 1, me
// %for I = MA, 1, ME %cycle
{
int _initial = MA, _final = ME;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = MA, 1, ME");
I = _initial - (1);
L_01F5:
if (I == _final) goto L_01F6;
I += 1;
// 1753 sf(i) = f(x(i), v)
[/*BUG*/ 0] = F([/*BUG*/ 0], V);
// 1754 %if k0=0 %then ->nrl105
if (K0 != 0) goto L_01F7;
goto U_100A;
L_01F7:;
// 1755 sf(i) = sf(i)+v(k0)
[/*BUG*/ 0] = [/*BUG*/ 0] + [/*BUG*/ 0];
// 1756 nrl105:s = s+(y(i)-sf(i))^2
U_100A:;
*S = *S + _imp_REXP([/*BUG*/ 0] - [/*BUG*/ 0], 2);
// 1757 %repeat
goto L_01F5;
}
L_01F6:;
// 1758 s = sqrt(s/(me-ma+1))
*S = _imp_SQRT((double)(*S / (ME - MA + 1)));
// 1759 %return
return;
// 1760 nrl100:newline; spaces(6)
U_1004:;
_imp_NEWLINE();
_imp_SPACES(6);
// 1761 printstring("THIS ERROR OCCURRED IN ROUTINE ===IMP DAPRO===")
_imp_PRINTSTRING(_imp_str_literal("THIS ERROR OCCURRED IN ROUTINE ===IMP DAPRO==="));
// 1762 newline
_imp_NEWLINE();
// 1763 %long %real %function af(%integer j, %long %real x)
double AF(int J, double X) {
// 1764 %long %real zw1, zw2, h
double ZW1;
double ZW2;
double H;
// 1765 %if j=k0 %then %result = 1
if (J != K0) goto L_01F8;
return 1;
L_01F8:;
// 1766 h = tol(j); %if h<5@-2 %then h = 5@-2
H = [/*BUG*/ 0];
if (H >= 5E-2) goto L_01F9;
H = 5E-2;
L_01F9:;
// 1767 %if 0.1<h %then h = 0.1
if (0.1 >= H) goto L_01FA;
H = 0.1;
L_01FA:;
// 1768 v(j) = v(j)+h/20
[/*BUG*/ 0] = [/*BUG*/ 0] + (double)(H / 20);
// 1769 zw1 = f(x, v)
ZW1 = F(X, V);
// 1770 v(j) = v(j)-h/10
[/*BUG*/ 0] = [/*BUG*/ 0] - (double)(H / 10);
// 1771 zw2 = f(x, v)
ZW2 = F(X, V);
// 1772 v(j) = v(j)+h/20
[/*BUG*/ 0] = [/*BUG*/ 0] + (double)(H / 20);
// 1773 %result = 10*(zw1-zw2)/h
return (double)((10 * (ZW1 - ZW2)) / H);
// 1774 %end
} /* end of AF */
// 1775 %routine cholesky(%integer n, %long %real %array (2) %name a, %long %real %array %name b,
void CHOLESKY(int N, double *A, double *B, int *FAULT) {
// 1776+ %integer %name fault)
// 1777 %long %real c, sum
double C;
double SUM;
// 1778 %integer i, j, k
int I;
int J;
int K;
// 1779 fault = 1
*FAULT = 1;
// 1780 %if n<2 %then ->nrl10
if (N >= 2) goto L_01FB;
goto U_101B;
L_01FB:;
// 1781 %cycle i = 1, 1, n-1
// %for I = 1, 1, N - 1 %cycle
{
int _final = N - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, N - 1");
I = 0;
L_01FC:
if (I == _final) goto L_01FD;
I += 1;
// 1782 %cycle j = i+1, 1, n
// %for J = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = I + 1, 1, N");
J = _initial - (1);
L_01FE:
if (J == _final) goto L_01FF;
J += 1;
// 1783 a(i, j) = 0
[/*BUG*/ 0] = 0;
// 1784 %repeat
goto L_01FE;
}
L_01FF:;
// 1785 %repeat
goto L_01FC;
}
L_01FD:;
// 1786 %if a(1, 1)<=0 %then %return
if ([/*BUG*/ 0] > 0) goto L_0200;
return;
L_0200:;
// 1787 c = sqrt(a(1, 1))
C = _imp_SQRT([/*BUG*/ 0]);
// 1788 a(1, 1) = c
[/*BUG*/ 0] = C;
// 1789 %cycle k = 2, 1, n
// %for K = 2, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(2, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 2, 1, N");
K = 1;
L_0201:
if (K == _final) goto L_0202;
K += 1;
// 1790 a(1, k) = a(k, 1)/c
[/*BUG*/ 0] = (double)([/*BUG*/ 0] / C);
// 1791 %repeat
goto L_0201;
}
L_0202:;
// 1792 %cycle i = 2, 1, n
// %for I = 2, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(2, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 2, 1, N");
I = 1;
L_0203:
if (I == _final) goto L_0204;
I += 1;
// 1793 sum = 0
SUM = 0;
// 1794 %cycle j = 1, 1, i-1
// %for J = 1, 1, I - 1 %cycle
{
int _final = I - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, I - 1");
J = 0;
L_0205:
if (J == _final) goto L_0206;
J += 1;
// 1795 sum = sum+a(j, i)*a(j, i)
SUM = SUM + [/*BUG*/ 0] * [/*BUG*/ 0];
// 1796 %repeat
goto L_0205;
}
L_0206:;
// 1797 %if a(i, i)<=sum %then %return
if ([/*BUG*/ 0] > SUM) goto L_0207;
return;
L_0207:;
// 1798 a(i, i) = sqrt(a(i, i)-sum)
[/*BUG*/ 0] = _imp_SQRT([/*BUG*/ 0] - SUM);
// 1799 %if i=n %then ->nrl1
if (I != N) goto L_0208;
goto U_1044;
L_0208:;
// 1800 %cycle k = i+1, 1, n
// %for K = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = I + 1, 1, N");
K = _initial - (1);
L_0209:
if (K == _final) goto L_020A;
K += 1;
// 1801 sum = 0
SUM = 0;
// 1802 %cycle j = 1, 1, i-1
// %for J = 1, 1, I - 1 %cycle
{
int _final = I - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle J = 1, 1, I - 1");
J = 0;
L_020B:
if (J == _final) goto L_020C;
J += 1;
// 1803 sum = sum+a(j, i)*a(j, k)
SUM = SUM + [/*BUG*/ 0] * [/*BUG*/ 0];
// 1804 %repeat
goto L_020B;
}
L_020C:;
// 1805 a(i, k) = (a(k, i)-sum)/a(i, i)
[/*BUG*/ 0] = (double)(([/*BUG*/ 0] - SUM) / [/*BUG*/ 0]);
// 1806 %repeat
goto L_0209;
}
L_020A:;
// 1807 nrl1: %repeat
U_1044:;
goto L_0203;
}
L_0204:;
// 1808 b(1) = b(1)/a(1, 1)
[/*BUG*/ 0] = (double)([/*BUG*/ 0] / [/*BUG*/ 0]);
// 1809 %cycle i = 2, 1, n
// %for I = 2, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(2, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 2, 1, N");
I = 1;
L_020D:
if (I == _final) goto L_020E;
I += 1;
// 1810 sum = 0
SUM = 0;
// 1811 %cycle k = 1, 1, i-1
// %for K = 1, 1, I - 1 %cycle
{
int _final = I - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, I - 1");
K = 0;
L_020F:
if (K == _final) goto L_0210;
K += 1;
// 1812 sum = sum+a(k, i)*b(k)
SUM = SUM + [/*BUG*/ 0] * [/*BUG*/ 0];
// 1813 %repeat
goto L_020F;
}
L_0210:;
// 1814 b(i) = (b(i)-sum)/a(i, i)
[/*BUG*/ 0] = (double)(([/*BUG*/ 0] - SUM) / [/*BUG*/ 0]);
// 1815 %repeat
goto L_020D;
}
L_020E:;
// 1816 b(n) = b(n)/a(n, n)
[/*BUG*/ 0] = (double)([/*BUG*/ 0] / [/*BUG*/ 0]);
// 1817 %cycle i = n-1, -1, 1
// %for I = N - 1, -1, 1 %cycle
{
int _initial = N - 1;
if (invalid_for_loop(_initial, -1, 1)) _imp_signal(5, 1, 0, "Illegal cycle I = N - 1, -1, 1");
I = _initial - (-1);
L_0211:
if (I == 1) goto L_0212;
I -= 1;
// 1818 sum = 0
SUM = 0;
// 1819 %cycle k = i+1, 1, n
// %for K = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = I + 1, 1, N");
K = _initial - (1);
L_0213:
if (K == _final) goto L_0214;
K += 1;
// 1820 sum = sum+a(i, k)*b(k)
SUM = SUM + [/*BUG*/ 0] * [/*BUG*/ 0];
// 1821 %repeat
goto L_0213;
}
L_0214:;
// 1822 b(i) = (b(i)-sum)/a(i, i)
[/*BUG*/ 0] = (double)(([/*BUG*/ 0] - SUM) / [/*BUG*/ 0]);
// 1823 %repeat
goto L_0211;
}
L_0212:;
// 1824 fault = 0
*FAULT = 0;
// 1825 %return
return;
// 1826 nrl10:
U_101B:;
// 1827 printstring("
_imp_PRINTSTRING(_imp_str_literal("\nN TOO SMALL, NAMELY N = "));
// 1828+ N TOO SMALL, NAMELY N = ")
// 1829 print(n, 3, 0)
_imp_PRINT(N, 3, 0);
// 1830 printstring(" THIS ERROR OCCURRED IN ROUTINE ")
_imp_PRINTSTRING(_imp_str_literal(" THIS ERROR OCCURRED IN ROUTINE "));
// 1831 printstring( %c
_imp_PRINTSTRING(_imp_str_literal("== CHOLESKY == CH\n"));
// 1832+ "== CHOLESKY == CH
// 1833+ ")
// 1834 %end
return;
} /* end of CHOLESKY */
// 1835 %end
return;
} /* end of IMPDAPRO */
// 1836 !*END*IMPDAPRO
// 1837
// 1838 !*TOP*IMPDCHLIN
// 1839 %external %routine imp dchlin(%integer m, n, %integer %name ifault,
void IMPDCHLIN(int M, int N, int *IFAULT, double *A, double *B, double *X, double *EPS) {
// 1840+ %long %real %array (2) %name a, %long %real %array %name b, x, %long %real %name eps)
// 1841 %long %real max, rho, sum, rhold
double MAX;
double RHO;
double SUM;
double RHOLD;
// 1842 %integer first, u, v, w, k, np, i
int FIRST;
int U;
int V;
int W;
int K;
int NP;
int I;
// 1843 %long %real %array apt(1:n, 1:n+1), llt(1:n, 1:n), s(1:n+1), bt(1:m)
double APT[N - 1 + 1][N + 1 - 1 + 1];
double LLT[N - 1 + 1][N - 1 + 1];
double S[N + 1 - 1 + 1];
double BT[M - 1 + 1];
// 1844 %integer %array j(1:n+2), ref(1:n+1)
int J[N + 2 - 1 + 1];
int REF[N + 1 - 1 + 1];
auto void SOLVE(int TR, double *IN, double *OUT);
auto void MULT(int TR, int R, int C, double *A, double *IN, double *OUT);
auto void REPLACE(int IN, int OUT);
auto int SIGN(double X);
// 1845 %routine %spec solve(%integer tr, %long %real %array %name in, out)
// 1846 %routine %spec mult(%integer tr, r, c, %long %real %array (2) %name a, %long %real %array %name in, out)
// 1847 %routine %spec replace(%integer in, out)
// 1848 %integer %function %spec sign(%long %real x)
// 1849 ifault = 1
*IFAULT = 1;
// 1850 %if n<2 %or m<=n %then ->er
if (N < 2) goto L_0003;
if (M > N) goto L_0215;
L_0003:;
goto U_1008;
L_0215:;
// 1851 first = 0
FIRST = 0;
// 1852 np = n+1
NP = N + 1;
// 1853 %cycle u = 1, 1, np
// %for U = 1, 1, NP %cycle
{
int _final = NP;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 1, 1, NP");
U = 0;
L_0216:
if (U == _final) goto L_0217;
U += 1;
// 1854 ref(u) = u
REF[U - 1] = U;
// 1855 %repeat
goto L_0216;
}
L_0217:;
// 1856 c1:%cycle v = 1, 1, np
U_1009:;
// %for V = 1, 1, NP %cycle
{
int _final = NP;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle V = 1, 1, NP");
V = 0;
L_0219:
if (V == _final) goto L_021A;
V += 1;
// 1857 w = ref(v)
W = REF[V - 1];
// 1858 %cycle u = 1, 1, n
// %for U = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 1, 1, N");
U = 0;
L_021B:
if (U == _final) goto L_021C;
U += 1;
// 1859 apt(u, v) = a(w, u)
APT[V - 1] = [/*BUG*/ 0];
// 1860 %repeat
goto L_021B;
}
L_021C:;
// 1861 bt(v) = b(w)
BT[V - 1] = [/*BUG*/ 0];
// 1862 %repeat
goto L_0219;
}
L_021A:;
// 1863 %cycle u = 1, 1, n
// %for U = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 1, 1, N");
U = 0;
L_021D:
if (U == _final) goto L_021E;
U += 1;
// 1864 %cycle v = 1, 1, u
// %for V = 1, 1, U %cycle
{
int _final = U;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle V = 1, 1, U");
V = 0;
L_021F:
if (V == _final) goto L_0220;
V += 1;
// 1865 sum = 0
SUM = 0;
// 1866 %cycle w = 1, 1, np
// %for W = 1, 1, NP %cycle
{
int _final = NP;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle W = 1, 1, NP");
W = 0;
L_0221:
if (W == _final) goto L_0222;
W += 1;
// 1867 sum = sum+apt(u, w)*apt(v, w)
SUM = SUM + APT[W - 1] * APT[W - 1];
// 1868 %repeat
goto L_0221;
}
L_0222:;
// 1869 llt(u, v) = sum
LLT[V - 1] = SUM;
// 1870 %repeat
goto L_021F;
}
L_0220:;
// 1871 %repeat
goto L_021D;
}
L_021E:;
// 1872 %if llt(1, 1)<=0 %then ->c2
if (LLT[0] > 0) goto L_0223;
goto U_100A;
L_0223:;
// 1873 rho = sqrt(llt(1, 1))
RHO = _imp_SQRT(LLT[0]);
// 1874 llt(1, 1) = rho
LLT[0] = RHO;
// 1875 %cycle w = 2, 1, n
// %for W = 2, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(2, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle W = 2, 1, N");
W = 1;
L_0224:
if (W == _final) goto L_0225;
W += 1;
// 1876 llt(1, w) = llt(w, 1)/rho
LLT[W - 1] = (double)(LLT[0] / RHO);
// 1877 %repeat
goto L_0224;
}
L_0225:;
// 1878 %cycle u = 2, 1, n
// %for U = 2, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(2, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 2, 1, N");
U = 1;
L_0226:
if (U == _final) goto L_0227;
U += 1;
// 1879 rho = llt(u, u)
RHO = LLT[U - 1];
// 1880 %cycle v = 1, 1, u-1
// %for V = 1, 1, U - 1 %cycle
{
int _final = U - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle V = 1, 1, U - 1");
V = 0;
L_0228:
if (V == _final) goto L_0229;
V += 1;
// 1881 rho = rho-llt(v, u)*llt(v, u)
RHO = RHO - LLT[U - 1] * LLT[U - 1];
// 1882 %repeat
goto L_0228;
}
L_0229:;
// 1883 %if rho<=0 %then ->c2
if (RHO > 0) goto L_022A;
goto U_100A;
L_022A:;
// 1884 rho = sqrt(rho)
RHO = _imp_SQRT(RHO);
// 1885 llt(u, u) = rho
LLT[U - 1] = RHO;
// 1886 %if u=n %then ->c3
if (U != N) goto L_022B;
goto U_100B;
L_022B:;
// 1887 %cycle w = u+1, 1, n
// %for W = U + 1, 1, N %cycle
{
int _initial = U + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle W = U + 1, 1, N");
W = _initial - (1);
L_022C:
if (W == _final) goto L_022D;
W += 1;
// 1888 sum = 0
SUM = 0;
// 1889 %cycle v = 1, 1, u-1
// %for V = 1, 1, U - 1 %cycle
{
int _final = U - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle V = 1, 1, U - 1");
V = 0;
L_022E:
if (V == _final) goto L_022F;
V += 1;
// 1890 sum = sum+llt(v, u)*llt(v, w)
SUM = SUM + LLT[U - 1] * LLT[W - 1];
// 1891 %repeat
goto L_022E;
}
L_022F:;
// 1892 llt(u, w) = (llt(w, u)-sum)/rho
LLT[W - 1] = (double)((LLT[U - 1] - SUM) / RHO);
// 1893 %repeat
goto L_022C;
}
L_022D:;
// 1894 %repeat
goto L_0226;
}
L_0227:;
// 1895 c3:solve(1, bt, x)
U_100B:;
SOLVE(1, BT, X);
// 1896 mult(-1, np, n, apt, x, s)
MULT(-1, NP, N, APT, X, S);
// 1897 sum = 0
SUM = 0;
// 1898 max = 0
MAX = 0;
// 1899 %cycle u = 1, 1, np
// %for U = 1, 1, NP %cycle
{
int _final = NP;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 1, 1, NP");
U = 0;
L_0230:
if (U == _final) goto L_0231;
U += 1;
// 1900 rho = bt(u)-s(u)
RHO = BT[U - 1] - S[U - 1];
// 1901 s(u) = rho
S[U - 1] = RHO;
// 1902 sum = sum+rho*rho
SUM = SUM + RHO * RHO;
// 1903 max = max+mod(rho)
MAX = MAX + _imp_MOD(RHO);
// 1904 %repeat
goto L_0230;
}
L_0231:;
// 1905 %if max<=eps/np %then rho = 0
if (MAX > (double)(*EPS / NP)) goto L_0232;
RHO = 0;
L_0232:;
// 1906 %if max>eps/np %then rho = sum/max
if (MAX <= (double)(*EPS / NP)) goto L_0233;
RHO = (double)(SUM / MAX);
L_0233:;
// 1907 %if first=0 %then ->c4
if (FIRST != 0) goto L_0234;
goto U_1041;
L_0234:;
// 1908 %if rho>rhold %then ->c4
if (RHO <= RHOLD) goto L_0235;
goto U_1041;
L_0235:;
// 1909 c2:%if first=-1 %then %return
U_100A:;
if (FIRST != -1) goto L_0236;
return;
L_0236:;
// 1910 %if first=0 %then ->c5
if (FIRST != 0) goto L_0237;
goto U_1014;
L_0237:;
// 1911 i = i+1
I = I + 1;
// 1912 %if j(i)=0 %then ->c5
if (J[I - 1] != 0) goto L_0238;
goto U_1014;
L_0238:;
// 1913 replace(j(i-1), j(i))
REPLACE(J[I - 1 - 1], J[I - 1]);
// 1914 ->c1
goto U_1009;
// 1915 c5:first = -1
U_1014:;
FIRST = -1;
// 1916 %cycle u = 1, 1, np
// %for U = 1, 1, NP %cycle
{
int _final = NP;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 1, 1, NP");
U = 0;
L_0239:
if (U == _final) goto L_023A;
U += 1;
// 1917 ref(u) = m-n+u-1
REF[U - 1] = M - N + U - 1;
// 1918 %repeat
goto L_0239;
}
L_023A:;
// 1919 ->c1
goto U_1009;
// 1920 c4:rhold = rho
U_1041:;
RHOLD = RHO;
// 1921 %if rho<eps %then ->c6
if (RHO >= *EPS) goto L_023B;
goto U_1015;
L_023B:;
// 1922 %cycle u = 1, 1, np
// %for U = 1, 1, NP %cycle
{
int _final = NP;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 1, 1, NP");
U = 0;
L_023C:
if (U == _final) goto L_023D;
U += 1;
// 1923 bt(u) = bt(u)-rho*sign(s(u))
BT[U - 1] = BT[U - 1] - RHO * SIGN(S[U - 1]);
// 1924 %repeat
goto L_023C;
}
L_023D:;
// 1925 solve(1, bt, x)
SOLVE(1, BT, X);
// 1926 c6:mult(1, m, n, a, x, bt)
U_1015:;
MULT(1, M, N, A, X, BT);
// 1927 max = 0
MAX = 0;
// 1928 k = 0
K = 0;
// 1929 %cycle u = 1, 1, m
// %for U = 1, 1, M %cycle
{
int _final = M;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 1, 1, M");
U = 0;
L_023E:
if (U == _final) goto L_023F;
U += 1;
// 1930 sum = b(u)-bt(u)
SUM = [/*BUG*/ 0] - BT[U - 1];
// 1931 %if max>mod(sum) %then ->c7
if (MAX <= _imp_MOD(SUM)) goto L_0240;
goto U_1016;
L_0240:;
// 1932 max = mod(sum)
MAX = _imp_MOD(SUM);
// 1933 k = u
K = U;
// 1934 w = sign(sum)
W = SIGN(SUM);
// 1935 c7:%repeat
U_1016:;
goto L_023E;
}
L_023F:;
// 1936 %if mod(max-rho)<eps %then ->c8
if (_imp_MOD(MAX - RHO) >= *EPS) goto L_0241;
goto U_1017;
L_0241:;
// 1937 %cycle u = 1, 1, n
// %for U = 1, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 1, 1, N");
U = 0;
L_0242:
if (U == _final) goto L_0243;
U += 1;
// 1938 x(u) = w*a(k, u)
[/*BUG*/ 0] = W * [/*BUG*/ 0];
// 1939 %repeat
goto L_0242;
}
L_0243:;
// 1940 solve(0, x, x)
SOLVE(0, X, X);
// 1941 mult(-1, np, n, apt, x, bt)
MULT(-1, NP, N, APT, X, BT);
// 1942 max = 0
MAX = 0;
// 1943 %cycle u = 1, 1, np
// %for U = 1, 1, NP %cycle
{
int _final = NP;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 1, 1, NP");
U = 0;
L_0244:
if (U == _final) goto L_0245;
U += 1;
// 1944 j(u) = 0
J[U - 1] = 0;
// 1945 %if mod(s(u))<eps*mod(bt(u)) %then ->c9
if (_imp_MOD(S[U - 1]) >= *EPS * _imp_MOD(BT[U - 1])) goto L_0246;
goto U_1018;
L_0246:;
// 1946 sum = bt(u)/s(u)
SUM = (double)(BT[U - 1] / S[U - 1]);
// 1947 %if max>=sum %then ->c9
if (MAX < SUM) goto L_0247;
goto U_1018;
L_0247:;
// 1948 max = sum
MAX = SUM;
// 1949 j(1) = u
J[0] = U;
// 1950 c9:%repeat
U_1018:;
goto L_0244;
}
L_0245:;
// 1951 j(np+1) = 0
J[NP + 1 - 1] = 0;
// 1952 %if j(1)=0 %then ->l1
if (J[0] != 0) goto L_0248;
goto U_1019;
L_0248:;
// 1953 i = 0
I = 0;
// 1954 %cycle u = 1, 1, np
// %for U = 1, 1, NP %cycle
{
int _final = NP;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle U = 1, 1, NP");
U = 0;
L_0249:
if (U == _final) goto L_024A;
U += 1;
// 1955 %if mod(s(u))<eps*mod(bt(u)) %then ->l2
if (_imp_MOD(S[U - 1]) >= *EPS * _imp_MOD(BT[U - 1])) goto L_024B;
goto U_101A;
L_024B:;
// 1956 sum = bt(u)/s(u)
SUM = (double)(BT[U - 1] / S[U - 1]);
// 1957 %if mod(max-sum)>eps %then ->l2
if (_imp_MOD(MAX - SUM) <= *EPS) goto L_024C;
goto U_101A;
L_024C:;
// 1958 i = i+1
I = I + 1;
// 1959 j(i) = ref(u)
J[I - 1] = REF[U - 1];
// 1960 l2:%repeat
U_101A:;
goto L_0249;
}
L_024A:;
// 1961 ->l3
goto U_101B;
// 1962 l1:j(1) = ref(1)
U_1019:;
J[0] = REF[0];
// 1963 l3:replace(k, j(1))
U_101B:;
REPLACE(K, J[0]);
// 1964 i = 1
I = 1;
// 1965 first = 1
FIRST = 1;
// 1966 ->c1
goto U_1009;
// 1967 c8:ifault = 0
U_1017:;
*IFAULT = 0;
// 1968 eps = rho
*EPS = RHO;
// 1969 %return
return;
// 1970 er:newline
U_1008:;
_imp_NEWLINE();
// 1971 printstring("M OR N TOO SMALL, NAMELY M = ")
_imp_PRINTSTRING(_imp_str_literal("M OR N TOO SMALL, NAMELY M = "));
// 1972 print(m, 3, 0)
_imp_PRINT(M, 3, 0);
// 1973 printstring(" N = ")
_imp_PRINTSTRING(_imp_str_literal(" N = "));
// 1974 print(n, 3, 0)
_imp_PRINT(N, 3, 0);
// 1975 newline
_imp_NEWLINE();
// 1976 printstring("THIS ERROR OCCURRED IN ROUTINE ")
_imp_PRINTSTRING(_imp_str_literal("THIS ERROR OCCURRED IN ROUTINE "));
// 1977 printstring("== IMP DCHLIN ==")
_imp_PRINTSTRING(_imp_str_literal("== IMP DCHLIN =="));
// 1978 newline
_imp_NEWLINE();
// 1979 %routine solve(%integer tr, %long %real %array %name in, out)
void SOLVE(int TR, double *IN, double *OUT) {
// 1980 %long %real s
double S;
// 1981 %integer i, k
int I;
int K;
// 1982 %if tr=1 %then mult(1, n, np, apt, in, out)
if (TR != 1) goto L_024D;
MULT(1, N, NP, APT, IN, OUT);
L_024D:;
// 1983 out(1) = out(1)/llt(1, 1)
[/*BUG*/ 0] = (double)([/*BUG*/ 0] / LLT[0]);
// 1984 %cycle i = 2, 1, n
// %for I = 2, 1, N %cycle
{
int _final = N;
if (invalid_for_loop(2, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 2, 1, N");
I = 1;
L_024E:
if (I == _final) goto L_024F;
I += 1;
// 1985 s = 0
S = 0;
// 1986 %cycle k = 1, 1, i-1
// %for K = 1, 1, I - 1 %cycle
{
int _final = I - 1;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, I - 1");
K = 0;
L_0250:
if (K == _final) goto L_0251;
K += 1;
// 1987 s = s+llt(k, i)*out(k)
S = S + LLT[I - 1] * [/*BUG*/ 0];
// 1988 %repeat
goto L_0250;
}
L_0251:;
// 1989 out(i) = (out(i)-s)/llt(i, i)
[/*BUG*/ 0] = (double)(([/*BUG*/ 0] - S) / LLT[I - 1]);
// 1990 %repeat
goto L_024E;
}
L_024F:;
// 1991 out(n) = out(n)/llt(n, n)
[/*BUG*/ 0] = (double)([/*BUG*/ 0] / LLT[N - 1]);
// 1992 %cycle i = n-1, -1, 1
// %for I = N - 1, -1, 1 %cycle
{
int _initial = N - 1;
if (invalid_for_loop(_initial, -1, 1)) _imp_signal(5, 1, 0, "Illegal cycle I = N - 1, -1, 1");
I = _initial - (-1);
L_0252:
if (I == 1) goto L_0253;
I -= 1;
// 1993 s = 0
S = 0;
// 1994 %cycle k = i+1, 1, n
// %for K = I + 1, 1, N %cycle
{
int _initial = I + 1, _final = N;
if (invalid_for_loop(_initial, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = I + 1, 1, N");
K = _initial - (1);
L_0254:
if (K == _final) goto L_0255;
K += 1;
// 1995 s = s+llt(i, k)*out(k)
S = S + LLT[K - 1] * [/*BUG*/ 0];
// 1996 %repeat
goto L_0254;
}
L_0255:;
// 1997 out(i) = (out(i)-s)/llt(i, i)
[/*BUG*/ 0] = (double)(([/*BUG*/ 0] - S) / LLT[I - 1]);
// 1998 %repeat
goto L_0252;
}
L_0253:;
// 1999 %end
return;
} /* end of SOLVE */
// 2000 %routine mult(%integer tr, r, c, %long %real %array (2) %name a, %long %real %array %name in, out)
void MULT(int TR, int R, int C, double *A, double *IN, double *OUT) {
// 2001 %long %real s, m
double S;
double M;
// 2002 %integer i, k
int I;
int K;
// 2003 %cycle i = 1, 1, r
// %for I = 1, 1, R %cycle
{
int _final = R;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, R");
I = 0;
L_0256:
if (I == _final) goto L_0257;
I += 1;
// 2004 s = 0
S = 0;
// 2005 %cycle k = 1, 1, c
// %for K = 1, 1, C %cycle
{
int _final = C;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, C");
K = 0;
L_0258:
if (K == _final) goto L_0259;
K += 1;
// 2006 %if tr=1 %then m = a(i, k)
if (TR != 1) goto L_025A;
M = [/*BUG*/ 0];
L_025A:;
// 2007 %if tr#1 %then m = a(k, i)
if (TR == 1) goto L_025B;
M = [/*BUG*/ 0];
L_025B:;
// 2008 s = s+m*in(k)
S = S + M * [/*BUG*/ 0];
// 2009 %repeat
goto L_0258;
}
L_0259:;
// 2010 out(i) = s
[/*BUG*/ 0] = S;
// 2011 %repeat
goto L_0256;
}
L_0257:;
// 2012 %end
return;
} /* end of MULT */
// 2013 %routine replace(%integer in, out)
void REPLACE(int IN, int OUT) {
// 2014 %integer k
int K;
// 2015 %cycle k = 1, 1, np
// %for K = 1, 1, NP %cycle
{
int _final = NP;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle K = 1, 1, NP");
K = 0;
L_025C:
if (K == _final) goto L_025D;
K += 1;
// 2016 %if ref(k)=out %then ->u1
if (REF[K - 1] != OUT) goto L_025E;
goto U_1045;
L_025E:;
// 2017 %repeat
goto L_025C;
}
L_025D:;
// 2018 u1: ref(k) = in
U_1045:;
REF[K - 1] = IN;
// 2019 %end
return;
} /* end of REPLACE */
// 2020 %integer %function sign(%long %real x)
int SIGN(double X) {
// 2021 %if x>0 %then %result = 1
if (X <= 0) goto L_025F;
return 1;
L_025F:;
// 2022 %if x<0 %then %result = -1
if (X >= 0) goto L_0260;
return -1;
L_0260:;
// 2023 %result = 0
return 0;
// 2024 %end
} /* end of SIGN */
// 2025 %end
return;
} /* end of IMPDCHLIN */
// 2026 !*END*IMPDCHLIN
// 2027
// 2028 !*TOP*IMPDELI3
// 2029 %external %long %real %function imp deli3(%long %real x, a, b)
double IMPDELI3(double X, double A, double B) {
// 2030 %long %real am, ap, c, ca, cb, d, de, e, f, fa, g, h, hh, ln2, p1, pi, pm,
double AM;
double AP;
double C;
double CA;
double CB;
double D;
double DE;
double E;
double F;
double FA;
double G;
double H;
double HH;
double LN2;
double P1;
double PI;
double PM;
double PZ;
double Q;
double R;
double S;
double T;
double U;
double V;
double W;
double Y;
double YE;
double Z;
double ZD;
// 2031+ pz, q, r, s, t, u, v, w, y, ye, z, zd
// 2032 %integer k, l, m, n, bo, bk
int K;
int L;
int M;
int N;
int BO;
int BK;
// 2033 %long %real %array ra, rb, rr(2:14)
double RA[13];
double RB[13];
double RR[13];
auto double ARSINH(double X);
auto int SIGN(double X);
// 2034 %long %real %function %spec arsinh(%long %real x)
// 2035 %integer %function %spec sign(%long %real x)
// 2036 %if x=0 %then %result = 0
if (X != 0) goto L_0261;
return 0;
L_0261:;
// 2037 ln2 = 0.6931471805599
LN2 = 0.6931471805599;
// 2038 pi = 3.141592653590
PI = 3.141592653590;
// 2039 ca = 1@-6; cb = 1@-14
CA = 1E-6;
CB = 1E-14;
// 2040 hh = x*x; f = b*hh
HH = X * X;
F = B * HH;
// 2041 %if a#0 %then s = a %else s = ca/(1+mod(x))
if (A == 0) goto L_0262;
S = A;
goto L_0263;
L_0262:;
S = (double)(CA / (1 + _imp_MOD(X)));
L_0263:;
// 2042 t = s*s; pm = t/2
T = S * S;
PM = (double)(T / 2);
// 2043 e = hh*t;
E = HH * T;
// 2044 z = mod(f)
Z = _imp_MOD(F);
// 2045 r = mod(b)
R = _imp_MOD(B);
// 2046 h = 1+hh
H = 1 + HH;
// 2047 %if e<0.1 %and z<0.1 %and t<1 %and r<1 %then ->sl
if (E >= 0.1) goto L_0264;
if (Z >= 0.1) goto L_0264;
if (T >= 1) goto L_0264;
if (R >= 1) goto L_0264;
goto U_1045;
L_0264:;
// 2048 w = 1+f
W = 1 + F;
// 2049 %if w=0 %then ->er
if (W != 0) goto L_0265;
goto U_1046;
L_0265:;
// 2050 %if b#0 %then p1 = b %else p1 = cb/hh
if (B == 0) goto L_0266;
P1 = B;
goto L_0267;
L_0266:;
P1 = (double)(CB / HH);
L_0267:;
// 2051 s = mod(s)
S = _imp_MOD(S);
// 2052 y = mod(x)
Y = _imp_MOD(X);
// 2053 g = p1-1; %if g=0 %then g = cb
G = P1 - 1;
if (G != 0) goto L_0268;
G = CB;
L_0268:;
// 2054 f = p1-t; %if f=0 %then f = cb*t
F = P1 - T;
if (F != 0) goto L_0269;
F = CB * T;
L_0269:;
// 2055 am = 1-t; ap = 1+e
AM = 1 - T;
AP = 1 + E;
// 2056 r = p1*h; fa = g/(f*p1)
R = P1 * H;
FA = (double)(G / (F * P1));
// 2057 %if fa>0 %then bo = 1 %else bo = -1
if (FA <= 0) goto L_026A;
BO = 1;
goto L_026B;
L_026A:;
BO = -1;
L_026B:;
// 2058 fa = mod(fa)
FA = _imp_MOD(FA);
// 2059 pz = mod(g*f)
PZ = _imp_MOD(G * F);
// 2060 de = sqrt(pz)
DE = _imp_SQRT(PZ);
// 2061 q = sqrt(mod(p1))
Q = _imp_SQRT(_imp_MOD(P1));
// 2062 %if pm>0.5 %then pm = 0.5
if (PM <= 0.5) goto L_026C;
PM = 0.5;
L_026C:;
// 2063 pm = p1-pm
PM = P1 - PM;
// 2064 %if pm>=0 %then %start
if (PM < 0) goto L_026D;
// 2065 v = y*de; u = sqrt(r*ap)
V = Y * DE;
U = _imp_SQRT(R * AP);
// 2066 %if g<0 %then v = -v
if (G >= 0) goto L_026E;
V = -V;
L_026E:;
// 2067 d = 1/q; c = 1
D = (double)(1 / Q);
C = 1;
// 2068 %finish %else %start
goto L_026F;
L_026D:;
// 2069 u = sqrt(h*ap*pz)
U = _imp_SQRT(H * AP * PZ);
// 2070 ye = y*q; v = am*ye
YE = Y * Q;
V = AM * YE;
// 2071 q = -de/g; d = -am/de
Q = -((double)(DE / G));
D = -((double)(AM / DE));
// 2072 c = 0; pz = ap-r
C = 0;
PZ = AP - R;
// 2073 %finish
L_026F:;
// 2074 %if bo=1 %then %start
if (BO != 1) goto L_0270;
// 2075 r = v/u
R = (double)(V / U);
// 2076 z = 1; k = 1
Z = 1;
K = 1;
// 2077 %if pm<0 %then %start
if (PM >= 0) goto L_0271;
// 2078 h = y*sqrt(h/(ap*fa))
H = Y * _imp_SQRT((double)(H / (AP * FA)));
// 2079 h = 1/h-h
H = (double)(1 / H) - H;
// 2080 z = h-2*r; r = 2+r*h
Z = H - 2 * R;
R = 2 + R * H;
// 2081 %if r=0 %then r = cb
if (R != 0) goto L_0272;
R = CB;
L_0272:;
// 2082 %if z=0 %then z = h*cb
if (Z != 0) goto L_0273;
Z = H * CB;
L_0273:;
// 2083 w = r/z; z = w
W = (double)(R / Z);
Z = W;
// 2084 r = w; w = pz
R = W;
W = PZ;
// 2085 %finish
L_0271:;
// 2086 u = u/w; v = v/w
U = (double)(U / W);
V = (double)(V / W);
// 2087 %finish %else %start
goto L_0274;
L_0270:;
// 2088 t = u+mod(v)
T = U + _imp_MOD(V);
// 2089 bk = 1
BK = 1;
// 2090 %if p1<0 %then %start
if (P1 >= 0) goto L_0275;
// 2091 de = v/pz; ye = 2*u*ye
DE = (double)(V / PZ);
YE = 2 * U * YE;
// 2092 u = t/pz; v = (-f-g*e)/t
U = (double)(T / PZ);
V = (double)((-F - G * E) / T);
// 2093 t = pz*mod(w)
T = PZ * _imp_MOD(W);
// 2094 z = (hh*r*f-g*ap+ye)/t
Z = (double)((HH * R * F - G * AP + YE) / T);
// 2095 ye = ye/t
YE = (double)(YE / T);
// 2096 %finish %else %start
goto L_0276;
L_0275:;
// 2097 de = v/w; ye = 0
DE = (double)(V / W);
YE = 0;
// 2098 u = (e+p1)/t
U = (double)((E + P1) / T);
// 2099 v = t/w; z = 1
V = (double)(T / W);
Z = 1;
// 2100 %finish
L_0276:;
// 2101 %if s>1 %then %start
if (S <= 1) goto L_0277;
// 2102 h = u; u = v
H = U;
U = V;
// 2103 v = h
V = H;
// 2104 %finish
L_0277:;
// 2105 %finish
L_0274:;
// 2106 y = 1/y; e = s
Y = (double)(1 / Y);
E = S;
// 2107 n = 1; t = 1
N = 1;
T = 1;
// 2108 l = 0; m = 0
L = 0;
M = 0;
// 2109 mr:y = y-e/y
U_1047:;
Y = Y - (double)(E / Y);
// 2110 %if y=0 %then y = cb*sqrt(e)
if (Y != 0) goto L_0278;
Y = CB * _imp_SQRT(E);
L_0278:;
// 2111 f = c; c = c+d/q
F = C;
C = C + (double)(D / Q);
// 2112 g = e/q; d = 2*(f*g+d)
G = (double)(E / Q);
D = 2 * (F * G + D);
// 2113 q = q+g; g = t
Q = Q + G;
G = T;
// 2114 t = s+t; n = 2*n
T = S + T;
N = 2 * N;
// 2115 m = 2*m
M = 2 * M;
// 2116 %if bo=1 %then %start
if (BO != 1) goto L_0279;
// 2117 %if z<0 %then m = m+k
if (Z >= 0) goto L_027A;
M = M + K;
L_027A:;
// 2118 k = sign(r); h = e/(u*u+v*v)
K = SIGN(R);
H = (double)(E / (U * U + V * V));
// 2119 u = u*(1+h); v = v*(1-h)
U = U * (1 + H);
V = V * (1 - H);
// 2120 %finish %else %start
goto L_027B;
L_0279:;
// 2121 r = u/v; h = z*r
R = (double)(U / V);
H = Z * R;
// 2122 z = h*z; hh = e/v
Z = H * Z;
HH = (double)(E / V);
// 2123 %if bk=1 %then %start
if (BK != 1) goto L_027C;
// 2124 de = de/u
DE = (double)(DE / U);
// 2125 ye = ye*(h+1/h)+de*(1+r)
YE = YE * (H + (double)(1 / H)) + DE * (1 + R);
// 2126 de = de*(u-hh)
DE = DE * (U - HH);
// 2127 %if mod(ye)<1 %then bk = 1 %else bk = -1
if (_imp_MOD(YE) >= 1) goto L_027D;
BK = 1;
goto L_027E;
L_027D:;
BK = -1;
L_027E:;
// 2128 %finish %else %start
goto L_027F;
L_027C:;
// 2129 z = log(z)
Z = _imp_LOG(Z);
// 2130 k = intpt(z/ln2)+1
K = _imp_INTPT((double)(Z / LN2)) + 1;
// 2131 z = exp(z-k*ln2)
Z = _imp_EXP(Z - K * LN2);
// 2132 m = m+k
M = M + K;
// 2133 %finish
L_027F:;
// 2134 %finish
L_027B:;
// 2135 %if mod(g-s)>ca*g %then %start
if (_imp_MOD(G - S) <= CA * G) goto L_0280;
// 2136 %if bo=1 %then %start
if (BO != 1) goto L_0281;
// 2137 g = (1/r-r)/2
G = (double)(((double)(1 / R) - R) / 2);
// 2138 hh = u+v*g
HH = U + V * G;
// 2139 h = g*u-v
H = G * U - V;
// 2140 %if hh=0 %then hh = u*cb
if (HH != 0) goto L_0282;
HH = U * CB;
L_0282:;
// 2141 %if h=0 %then h = v*cb
if (H != 0) goto L_0283;
H = V * CB;
L_0283:;
// 2142 z = r*h
Z = R * H;
// 2143 r = hh/h
R = (double)(HH / H);
// 2144 %finish %else %start
goto L_0284;
L_0281:;
// 2145 u = u+e/u
U = U + (double)(E / U);
// 2146 v = v+hh
V = V + HH;
// 2147 %finish
L_0284:;
// 2148 s = 2*sqrt(e)
S = 2 * _imp_SQRT(E);
// 2149 e = s*t
E = S * T;
// 2150 l = 2*l
L = 2 * L;
// 2151 %if y<0 %then l = l+1
if (Y >= 0) goto L_0285;
L = L + 1;
L_0285:;
// 2152 ->mr
goto U_1047;
// 2153 %finish
L_0280:;
// 2154 %if y<0 %then l = l+1
if (Y >= 0) goto L_0286;
L = L + 1;
L_0286:;
// 2155 e = arctan(1, t/y)+pi*l
E = _imp_ARCTAN(1, (double)(T / Y)) + PI * L;
// 2156 e = e*(c*t+d)/(t*(t+q))
E = (double)((E * (C * T + D)) / (T * (T + Q)));
// 2157 %if bo=1 %then %start
if (BO != 1) goto L_0287;
// 2158 h = v/(t+u)
H = (double)(V / (T + U));
// 2159 z = 1-r*h
Z = 1 - R * H;
// 2160 h = r+h
H = R + H;
// 2161 %if z=0 %then z = cb
if (Z != 0) goto L_0288;
Z = CB;
L_0288:;
// 2162 %if z<0 %then m = m+sign(h)
if (Z >= 0) goto L_0289;
M = M + SIGN(H);
L_0289:;
// 2163 s = arctan(1, h/z)+m*pi
S = _imp_ARCTAN(1, (double)(H / Z)) + M * PI;
// 2164 %finish %else %start
goto L_028A;
L_0287:;
// 2165 %if bk=1 %then s = arsinh(ye) %else s = log(z)+m*ln2
if (BK != 1) goto L_028B;
S = ARSINH(YE);
goto L_028C;
L_028B:;
S = _imp_LOG(Z) + M * LN2;
L_028C:;
// 2166 s = s/2
S = (double)(S / 2);
// 2167 %finish
L_028A:;
// 2168 e = (e+sqrt(fa)*s)/n
E = (double)((E + _imp_SQRT(FA) * S) / N);
// 2169 %result = e*sign(x)
return E * SIGN(X);
// 2170 sl:%cycle k = 2, 1, 10
U_1045:;
// %for K = 2, 1, 10 %cycle
{
K = 1;
L_028D:
if (K == 10) goto L_028E;
K += 1;
// 2171 rb(k) = 0.5/k
RB[K - 2] = (double)(0.5 / K);
// 2172 ra(k) = 1-rb(k)
RA[K - 2] = 1 - RB[K - 2];
// 2173 %repeat
goto L_028D;
}
L_028E:;
// 2174 zd = 0.0454545454545
ZD = 0.0454545454545;
// 2175 s = b+pm
S = B + PM;
// 2176 %cycle k = 2, 1, 10
// %for K = 2, 1, 10 %cycle
{
K = 1;
L_028F:
if (K == 10) goto L_0290;
K += 1;
// 2177 rr(k) = s
RR[K - 2] = S;
// 2178 pm = pm*t*ra(k)
PM = PM * T * RA[K - 2];
// 2179 s = s*b+pm
S = S * B + PM;
// 2180 %repeat
goto L_028F;
}
L_0290:;
// 2181 u = s*zd
U = S * ZD;
// 2182 s = u
S = U;
// 2183 bo = -1
BO = -1;
// 2184 %cycle k = 10, -1, 2
// %for K = 10, -1, 2 %cycle
{
K = 11;
L_0291:
if (K == 2) goto L_0292;
K -= 1;
// 2185 u = u+(rr(k)-u)*rb(k)
U = U + (RR[K - 2] - U) * RB[K - 2];
// 2186 bo = -bo
BO = -BO;
// 2187 v = -u*bo
V = -(U * BO);
// 2188 s = s*hh+v
S = S * HH + V;
// 2189 %repeat
goto L_0291;
}
L_0292:;
// 2190 %if bo=1 %then s = -s
if (BO != 1) goto L_0293;
S = -S;
L_0293:;
// 2191 u = (u+1)/2
U = (double)((U + 1) / 2);
// 2192 %result = (u-s*h)*sqrt(h)*x+u*arsinh(x)
return (U - S * H) * _imp_SQRT(H) * X + U * ARSINH(X);
// 2193 er:newline
U_1046:;
_imp_NEWLINE();
// 2194 printstring("1+B*X*X = 0, NAMELY B = ")
_imp_PRINTSTRING(_imp_str_literal("1+B*X*X = 0, NAMELY B = "));
// 2195 printfl(b, 10)
_imp_PRINTFL(B, 10);
// 2196 printstring(" X = ")
_imp_PRINTSTRING(_imp_str_literal(" X = "));
// 2197 printfl(x, 10)
_imp_PRINTFL(X, 10);
// 2198 newline
_imp_NEWLINE();
// 2199 printstring("THIS ERROR OCCURRED IN ROUTINE == IMP DELI3 ==")
_imp_PRINTSTRING(_imp_str_literal("THIS ERROR OCCURRED IN ROUTINE == IMP DELI3 =="));
// 2200 newline
_imp_NEWLINE();
// 2201 %result = 1@17
return 1E17;
// 2202 %integer %function sign(%long %real x)
int SIGN(double X) {
// 2203 %if x>0 %then %result = 1
if (X <= 0) goto L_0294;
return 1;
L_0294:;
// 2204 %if x=0 %then %result = 0
if (X != 0) goto L_0295;
return 0;
L_0295:;
// 2205 %result = -1
return -1;
// 2206 %end
} /* end of SIGN */
// 2207 %long %real %function arsinh(%long %real x)
double ARSINH(double X) {
// 2208 %long %real a, y, b0, b1, b2
double A;
double Y;
double B0;
double B1;
double B2;
// 2209 %integer i
int I;
// 2210 %long %real %array c(0:13)
double C[14];
// 2211 a = x*x
A = X * X;
// 2212 %if a<1 %then %start
if (A >= 1) goto L_0296;
// 2213 c(0) = 1.8717996008826; c(1) = -0.0588117611900
C[0] = 1.8717996008826;
C[1] = -0.0588117611900;
// 2214 c(2) = 0.0047274654322; c(3) = -0.0004938363163
C[2] = 0.0047274654322;
C[3] = -0.0004938363163;
// 2215 c(4) = 0.0000585062071; c(5) = -0.0000074669983
C[4] = 0.0000585062071;
C[5] = -0.0000074669983;
// 2216 c(6) = 0.0000010011694; c(7) = -0.0000001390354
C[6] = 0.0000010011694;
C[7] = -0.0000001390354;
// 2217 c(8) = 0.0000000198232; c(9) = -0.0000000028847
C[8] = 0.0000000198232;
C[9] = -0.0000000028847;
// 2218 c(10) = 0.0000000004267; c(11) = -0.0000000000640
C[10] = 0.0000000004267;
C[11] = -0.0000000000640;
// 2219 c(12) = 0.0000000000097; c(13) = -0.0000000000015
C[12] = 0.0000000000097;
C[13] = -0.0000000000015;
// 2220 b1 = 0.0000000000002
B1 = 0.0000000000002;
// 2221 y = 2*a-1
Y = 2 * A - 1;
// 2222 b0 = 2*y*b1+c(13)
B0 = 2 * Y * B1 + C[13];
// 2223 %cycle i = 12, -1, 0
// %for I = 12, -1, 0 %cycle
{
I = 13;
L_0297:
if (I == 0) goto L_0298;
I -= 1;
// 2224 b2 = b1
B2 = B1;
// 2225 b1 = b0
B1 = B0;
// 2226 b0 = c(i)-b2+2*y*b1
B0 = C[I] - B2 + 2 * Y * B1;
// 2227 %repeat
goto L_0297;
}
L_0298:;
// 2228 %result = x*(b0-b2)/2
return (double)((X * (B0 - B2)) / 2);
// 2229 %finish
L_0296:;
// 2230 y = mod(x)
Y = _imp_MOD(X);
// 2231 %if a>1@16 %then b0 = log(y)+0.6931471805599 %else b0 = log(y+sqrt(1+a))
if (A <= 1E16) goto L_0299;
B0 = _imp_LOG(Y) + 0.6931471805599;
goto L_029A;
L_0299:;
B0 = _imp_LOG(Y + _imp_SQRT(1 + A));
L_029A:;
// 2232 %if x<0 %then b0 = -b0
if (X >= 0) goto L_029B;
B0 = -B0;
L_029B:;
// 2233 %result = b0
return B0;
// 2234 %end
} /* end of ARSINH */
// 2235 %end
_imp_signal(8, 4, 0, "result missing");
} /* end of IMPDELI3 */
// 2236 !*END*IMPDELI3
// 2237
// 2238 %end %of %file