code 41021;
real procedure WILCOXINV(PROB, M, N, LEFT);
value PROB, M, N, LEFT; real PROB, M, N; Boolean LEFT;
begin
    integer X, W, WI, MN; real Z; Boolean MN EVEN;
    MN := M × N;
    MN EVEN := ENTIER(MN/2) × 2 = MN;
    X := if MN EVEN then MN/2 - 1 else MN/2 - 1.5;
    PROB := PROB + 10-13×(1-PROB);
    if PROB ≤ 0 ∨ PROB ≥ 1 then
    STATAL3 ERROR(“WILCOXINV”, 1, PROB) else
    if M < 0 ∨ ENTIER(M) < M then
    STATAL3 ERROR(“WILCOXINV”, 2, M) else
    if N < 0 ∨ ENTIER(N) < N then
    STATAL3 ERROR(“WILCOXINV”, 2, N) else
    if MN = 0 then WI := -2 else
    if PROB = .5 then WI := ENTIER((MN-1)/2) × 2 else
    if M = 1 then WI := ENTIER(PROB×(N+1)) × 2 - 2
    else
    if N = 1 then WI := ENTIER(PROB×(M+1)) × 2 - 2
    else if MN > 400 ∨ M=2 ∨ N=2 then
    begin Z := PHINV(PROB) × SQRT(MN×(M+N+1)/3) + MN;
        WI := W := ENTIER(Z/2) × 2;
        if WI < 0 then WI := W := 0;
        if WI > 2×MN then WI := W := 2×MN;
        if WILCOX(W, M, N) ≤ PROB then
        begin for W := W + 2
          while WILCOX(W, M, N) ≤ PROB do WI := W
        end else
        begin for W:= W - 2
           while WILCOX(W, M, N) > PROB do WI := W;
             WI := WI - 2;
        end;
    end else
    begin integer I,J,UP,UP1; real H1,H2;
        Boolean RIGHT, EQUAL; real array WCMN[-1:X+2];
        integer procedure MAX(A,B);
        value A,B; integer A,B;
             MAX := if A > B then A else B;
         integer procedure MIN(A,B);
         value A,B; integer A,B;
             MIN := if A > B then B else A;
         RIGHT := PROB > .5;
            if RIGHT then PROB := 1 - PROB;
            I := MAX(M,N); J := MIN(M,N);
            if I × (X+1) > 12000 then
         begin M := J; N := I end else
         begin M := I; N := J end;
         begin real array WP[0:X, 1:M];
             for I := MAX(2,-ENTIER(X/2-M)) step 1
                until M do
                begin UP:= X-(M-I)×2; UP1 := MIN(UP,I);
                  H1 := 1/(I+1);
                  for W:= 0 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:= 0 step 1 until UP do
                  WP[W,1] := H2;
                  for I := MAX(2,-ENTIER(X/J-M)) step 1
                  until M do
                  begin UP:= X - (M-I)×J;
                    H1:= J/(I+J); H2 := I/(I+J);
                    UP1 := MIN(UP,J-1);
                    for W := 0 step 1 until UP1 do
                      WP[W,I] := WP[W,I]×H1;
                    UP1 := MIN(UP,I×J-I);
                    for W := J 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 := I×J-I+1 step 1 until UP1
                    do  WP[W,I] := WP[W-J,I-1]×H2;
                   end
               end;
               WCMN[-1] := 0;
               WCMN[0] := WP[0,M];
               for W := 1 step 1 until X do
                     WCMN[W] := WCMN[W-1] + WP[W,M];
               if MN EVEN then
                 begin WCMN[X+1] := 1 - WCMN[X];
                         WCMN[X+2] := 1 - WCMN[X-1];
                   end
               else
                 begin WCMN[X+1] := .5;
                         WCMN[X+2] := 1 - WCMN[X];
                   end;
          end;
          WI := PHINV(PROB) × SQRT(MN×(M+N+1)/3) + MN;
          W := ENTIER(WI/2);
          if W < 0 then WI := W := 0 else WI:= W;
          if WCMN[W] ≤ PROB then
          begin for W := W + 1
                while WCMN[W] ≤ PROB do WI := W
          end else
          begin for W := W - 1
              while WCMN[W] > PROB do WI := W;
              WI := WI - 1;
          end;
          EQUAL := WCMN[WI] = PROB;
          if RIGHT then
          begin if EQUAL then WI := 2 × (MN - WI - 1)
                             else WI := 2 × (MN - WI - 2)
          end
          else WI := 2 × WI;
      end;
      WILCOXINV := if LEFT then WI else 2×MN - WI;
end WILCOXINV;
eop