code 41569;
real procedure BINORCOR(R,RHO,N);
value R, RHO, N; real R, RHO, N;
begin
real procedure SAMCORBIVNORDEN(R,RHO,N);
value R, RHO, N; real R, RHO, N;
begin
if ABS(R) ≥ 1 then SAMCORBIVNORDEN := 0
else
begin
real W1, W3, Y1, Y2, Y3, Y4, N1, R2, RRHO,
R2RHO2, W2, PI, RHO2;
R2 := R × R;
RHO2 := RHO × RHO;
RRHO := R × RHO;
R2RHO2 := R2 × RHO2;
W1 := SQRT(1 - R2);
W2 := SQRT(1 - RHO2);
W3 := SQRT(1 - RHO2 × R2);
PI := ARCTAN(1) × 4;
N1 := N - 1;
if N < 15 then
begin
real array SB[3:N]; integer I;
SB[3] :=
(1-RHO2) / PI / W1 × (1 + RRHO ×
ARCCOS(-RRHO) / W3) / (1 - R2RHO2);
if N =3 then SAMCORBIVNORDEN := SB[3]
else
begin
SB[4] :=
(1 - RHO2) × W2 / PI × (3 × RRHO +
(1 + 2 × R2RHO2) × ARCCOS(-RRHO) / W3)
/ (1 - R2RHO2) / (1 - R2RHO2);
if N = 4 then SAMCORBIVNORDEN := SB[4]
else
begin for I := 5 step 1 until N do
SB[I] :=
(2 × I - 5) / (I - 3) × RRHO × W1 × W2
/ (1 - R2RHO2) × SB[I-1] +
(I - 3) / (I - 4) × (1 - RHO2) ×
(1 - R2) / (1 - R2RHO2) × SB[I-2];
SAMCORBIVNORDEN := SB[N];
end;
end;
end
else
begin
Y1 := (RRHO + 2) / 8;
Y2 := (3 × RRHO + 2) × (3 × RRHO + 2) / 128;
Y3 := (((15 × RRHO + 18) × RRHO - 4)
× RRHO - 8) × 5 / 1024;
Y4 := ((((3675 × RRHO + 4200) × RRHO - 2520)
× RRHO - 3360) × RRHO - 336) / 32768;
SAMCORBIVNORDEN :=
(N - 2) / SQRT(N - 1) × (1 - RHO2) × W2
× (W1 × W2 / (1 - RRHO)) ⭡ N1 ×
SQRT((1 - RRHO) / 2 / PI)
/ ((1 - R2) × W1 × (1 - RHO2) × W2)
× ((((Y4 / N1 + Y3) / N1 + Y2) / N1 + Y1)
/ N1 + 1);
end;
end;
end SAMCORBIVNORDEN;
if ABS(RHO) ≥ 1 then
STATAL3 ERROR(“BINORCOR”, 2, RHO)
else if N > ENTIER(N) ∨ N < 3 then
STATAL3 ERROR(“BINORCOR”, 3, N)
else if R ≤ -1 then BINORCOR := 0
else if R ≥ 1 then BINORCOR := 1
else if RHO = 0 then
BINORCOR :=
STUDENT(R × SQRT((N - 2) / (1 - R × R)), N - 2)
else
if N ≤ 500 then
begin real W1, W2, W3, PI, R2, RHO2, RHO3,
RHO4, RRHO, R2RHO2;
R2 := R × R;
RHO2 := RHO × RHO;
RHO3 := RHO2 × RHO;
RHO4 := RHO2 × RHO2;
RRHO := R × RHO;
R2RHO2 := R2 × RHO2;
W1 := SQRT(1 - R2);
W2 := SQRT(1 - RHO2);
W3 := SQRT(1 - RHO2 × R2);
PI := ARCTAN(1) × 4;
if N = 3 then
BINORCOR :=
(ARCCOS(-R) - RHO × W1 / W3 × ARCCOS(-RRHO)) / PI
else if N = 4 then
BINORCOR :=
W1 × W2 × SAMCORBIVNORDEN(R, RHO, 3) / RHO -
(W2 / RHO - ARCCOS(RHO)) / PI
else if N = 5 then
BINORCOR :=
W1 × W2 × SAMCORBIVNORDEN(R, RHO, 4) / 2 / RHO
- R × (1 - R2) / 2 × SAMCORBIVNORDENC(R, RHO, 3)
- W1 × (1 + RHO2) / (2 × PI × RHO) × ARCCOS(-RRHO)
/ W3 + ARCCOS(-R) / PI
else
if N = 6 then
BINORCOR :=
W2 × (1 - 4 × RHO2) / (3 × PI × RHO3)
+ ARCCOS(RHO) / PI
- (1 - RHO2) × W1 × W2 / 3 / RHO3 ×
SAMCORBIVNORDEN(R, RHO, 3)
+ (1 - RHO2) × R / 3 / RHO2 ×
SAMCORBIVNORDEN(R, RHO, 4)
+ W1 × W2 / 3 / RHO × SAMCORBIVNORDEN(R, RHO, 5)
else if N = 7 then
BINORCOR :=
ARCCOS(-R) / PI - (3 + 6 × RHO2 - RHO4) ×
ARCCOS(-RRHO) / W3 × W1 / (8 × PI × RHO)
-R × (1 - R2) × (4 - 3 × RHO2 + 3 × RHO4) / 8 /
RHO2 × SAMCORBIVNORDEN(R, RHO, 3)
- R2 × W1 × W2 × (2 - RHO2) / 8 / RHO ×
SAMCORBIVNORDEN(R, RHO, 4)
+ (1 - RHO2) × R / 4 / RHO2 ×
SAMCORBIVNORDEN(R, RHO, 5)
+ W1 × W2 / 4 / RHO × SAMCORBIVNORDEN(R, RHO, 6)
else if N = 8 then
BINORCOR := ARCCOS(RHO)
/ PI - W2 × (3 - 11 × RHO2 + 23 × RHO4)
/ 15 / PI / RHO4 / RHO
+ (W2 / RHO) ⭡ 5 × W1 / 5 ×
SAMCORBIVNORDEN(R, RHO, 3)
- R × (1-RHO2) × (1-RHO2) / 5 / RHO4 ×
SAMCORBIVNORDEN(R, RHO, 4)
+ (3 × R2 - 1) × (1 - RHO2) × W2 / W1 / 15 / RHO3
× SAMCORBIVNORDEN(R, RHO, 5)
+ (1 - RHO2) × R / 5 / RHO2 ×
SAMCORBIVNORDEN(R, RHO, 6)
+ W1 × W2 / 5 / RHO × SAMCORBIVNORDENC(R, RHO, 7)
else
begin real array E[1:3]; real X;
E[1] := E[2] := 10-6;
BINORCOR :=
STUDENT(-RHO × SQRT(N - 1) / W2, N - 1)
+ QADRAT(X, 0, R, SAMCORBIVNORDEN(X, RHO, N), E)
end;
end
else
begin real R2, RHO2, RHO3; integer N1, N2, N3;
N1 := N - 1; N2 := N1 × N1; N3 := N1 × N2;
R2 := R × R;
RHO2 := RHO × RHO;
RHO3 := RHO2 × RHO;
BINORCOR :=
if ABS(RHO) ≤ .7 then
PHI((R × SQRT((N - 2.5) / (1 - R2)) -
RHO × SQRT((N - 1.5) / (1 - RHO2)))
/ SQRT(1 + RHO2 / 2 / (1 - RHO2) + R2 /
2 / (1 - R2)))
else
PHI((.5 × LN((1 + R) × (1 - RHO) / (1 - R) /
(1 + RHO)) - RHO / 2/ N1 -
(5 × RHO + RHO3) / 8 / N2 - (11 × RHO +
(2 + RHO2) × 3 × RHO3) / 16 / N3)
/ SQRT(1 / N1 + (4 - RHO2) / 2 / N2 +
(22 - (6 - 3 × RHO2) × RHO2) / 6 / N3));
end;
end BINORCOR;
eop