code 41020;
real procedure WILCOX(X,M,N);
value X,M,N; real X,M,N;
begin
    integer procedure MIN(A,B);
    value A,B; integer A,B;
        MIN := if A > B then B else A;
    integer procedure MAX(A,B);
    value A,B; integer A,B;
        MAX := if A > B then A else B;
    real WP1;    Boolean X EVEN, RIGHT;
    integer MN;
    if M < 0 ∨ ENTIER(M) < M then
    STATAL3 ERROR(“WILCOX”, 2, M);
    if N < 0 ∨ ENTIER(N) < N then
    STATAL3 ERROR(“WILCOX”, 3, N);
    MN := M × N;
    X := ENTIER(X/2);
    if X < MN/2 then RIGHT := false else
    begin RIGHT := true; X := MN-X-1 end;
    X EVEN := ENTIER(X/2) × 2 = X;
    M := MIN(M,N); N := MN/M;
    if X < 0 then WP1 := 0 else
    if M = 1 then WP1 := (X+1) / (N+1) else
    if M = 2 then WP1 := (if X EVEN then
        (X+2)×(X+2) / (2×(N+1)×(N+2)) else
        (X+1)×(X+3) / (2×(N+1)×(N+2)) )     else
    if 2×X = MN - 1 then WP1 := .5 else
    if MN > 400 then
    begin integer NOEM, N2, N3, N4, M2, M3, M4;
        real F0Y, F3Y, F5Y, F7Y, T3, T5, T7, Y, Y2;
        Y := (2×X + 1 - MN) / SQRT(MN × (M + N + 1) / 3);
        Y2 := Y × Y;
        NOEF := 10 × MN × (M + N + 1); F0Y := PHIDENS(Y);
        N2 := N × N; N3 := N2 × N; N4 := N2 × N2;
        M2 := M × M; M3 := M2 × M; M4 := M2 × M2;
        F3Y := Y × (3 - Y2);
        F5Y := Y × (-15 + Y2 × (10 - Y2));
        F7Y := Y × (105 - Y2 × (105 - Y2 × (21 - Y2)));
        T3 := (M2 + N2 + MN + M + N) / NOEM / 2;
        T5 := ( 2 × (M4 + N4) + 4 +
                  (M3 × N + N3 × M + M3 + N3) +
                  6 × M2 × N2 + 7 × MN × (M + N) + M2 + N2 +
                  2 × MN - M - N) / (NOEM × NOEM × 2.1);
        T7 := (M2 + N2 + MN + M + N) ⭡ 2
              / (NOEM × NOEM × 8);
        WP1 := MAX(PHI(Y) - F0Y ×
               (T3 × F3Y - T5 × F5Y - T7 × F7Y), 0);
    end else
    begin integer I,J,W,UP,UP1,UP2; real H1,H2;
        if N × (X+1) ≤ 12000 then
        begin M := N; N := MN / M end;
        begin real array WP[0:X, 1:M];
            UP2 := MIN(M, ENTIER((MN-X-1)/(N-1)));
            for I := MAX(2, -ENTIER(X/2-M)) step 1
            until UP2 do
            begin UP:= X-(M-I)×2; UP1 := MIN(UP,I-1);
                   H1 := 1/(I+1);
                   for W:= MAX(0, X-(M-I)×N) step 1
                   until UP1 do
                   WP[W,I] := H1 × (W+1);
            end;
            for J := 2 step 1 until N do
            begin UP:= MIN(X-(M-1)×J , J-1);
                H2 := 1/(J+1);
                for W:= MAX(0, X-(M-2)×N-J) step 1
                until UP do
                WP[W,1] := H2 × (W+1);
                UP2 := (if J×M < X+1 then
                  ENTIER((MN-X-1)/(N-J)) else M);
                for I := MAX(2,-ENTIER(X/J-M)) step 1
               until UP2 do
               begin UP:= X - (M-I)×J;
                  H1:= J/(I+J); H2 := I/(I+J);
                  UP1 := MIN(UP,J-1);
                  for  W := MAX(0, X-(M-I)×N) step 1
                  until UP1 do
                    WP[W,I] := WP[W,I]×H1;
                  UP1 := MIN(UP,I×J-I-1);
                  for  W := MAX(J, X-(M-I)×N) step 1
                  until UP1 do
                       WP[W,I] := WP[W,I]×H1 + WP[W-J,I-1]×H2;
                    UP1 := MIN(UP,I×J-1);
                    for W:= MAX(I×J-I, X-(M-I)×N) step 1
                    until UP1 do
                      WP[W,I] := H1 + WP[W-J,I-1]×H2;
                   end
              end;
              WP1 := WP[X,M];
          end;
     end;
     WILCOX := if RIGHT then 1-WP1 else WP1;
end WILCOXCDF;
eop