code 41022;
real procedure WILCOXPROB(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;
    integer MN;
    if M < 0 ∨ ENTIER(M) < M then
    STATAL3 ERROR(“WILCOXPROB”, 2, M);
    if N < 0 ∨ ENTIER(N) < N then
    STATAL3 ERROR(“WILCOXPROB”, 3, N);
    MN := M × N;
    X := MIN(X/2, MN - X/2);
    M := MIN(M,N); N := MN/M;
    if ENTIER(X) < X then WP1 := 0 else
    if X < 0 then WP1 := 0 else
    if MN = 0 then WP1 := 1 else
    if M = 1 then WP1 := 1/(N+1) else
    if M = 2 then
             WP1 := ENTIER(X/2+1) / ((N+1)×(N+2)/2) else
     if MN > 400 then
         WP1 := WILCOX(2×X,M,N) - WILCOX(2×X-2,M,N) 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)/(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);
                      H1 := 1/(I+1);
                      for W:= MAX(0, X-(M-I)×N) step 1
                      until UP1 do
                      WP[W,I] := H1;
                  end;
                  for J := 2 step 1 until N do
                  begin UP:= MIN(X-(M-1)×J , J); H2 := 1/(J+1);
                    for W:= MAX(0, X-(M-2)×N-J) step 1
                    until UP do
                      WP[W,1] := H2;
                      UP2 := if J×M<X+1 then
                          ENTIER((MN-X)/(N-J)) else M;
                      for I := MAX(2,-ENTIER(X/2-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,1]×H1;
                        UP1 := MIN(UP,I×J-I);
                        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);
                        for W:=MAX(I×J-I+1,X-(M-I)×N) step 1
                        until UP1 do
                          WP[W,I] := WP[W-J,I-1]×H2;
                      end
                  end;
                  WP1   := WP[X,M];
             end;
     end;
     WILCOXPROB := WP1
end WILCOXPROB;
eop