code 41556;
real procedure KOLSMIR(D, XSIZE, YSIZE, EPS);
value D, XSIZE, YSIZE, EPS; real D, XSIZE, YSIZE, EPS;
begin integer I, KGV, M, N;
      integer procedure GGD(M, N);
      value M, N; integer M, N;
      GGD:= if N = 0 then M else GGD(N, M - M ÷ N × N);
      
      procedure APPROX;
      begin integer K; real SUM, X, TERM, THETA;
          SUM:= .5; THETA:= (1 + (M / KGV) ⭡ 1.2) / (M + N);
          X:= (I / KGV + THETA) ⭡ 2 × 2 / (1 / M + 1 / N);
          for K:= 1, K + 2 while TERM > EPS do
          begin TERM:= EXP(-X × K × K);
              SUM:= SUM - TERM × (1 - EXP(-X × (2 × K + 1)))
          end;
          KOLSMIR:= 2 × SUM
      end;
      
      procedure EXACT;
      begin integer array LOW[0:N]; array H[0:M];
          integer DMN, MN1, X, Y, UPP; real SUM, BINCOEF;
          BINCOEF:= 1; LOW[0]:= 0; MN1:= M + N + 1;
          DMN:= I × M × N / KGV;
          for X:= 1 step 1 until N do
          begin integer T, TN;
              T:= M × X - DMN; TN:= T ÷ N;
              LOW[X]:= if T ≤ 0 then 0 else
                  if TN = T / N then TN else TN + 1;
              BINCOEF:= BINCOEF × (MN1 - X) / X;
              if BINCOEF > 10318 then
              begin EPS:= 10-4;goto L end
          end;
          H[0]:= 1;
          for Y:= 1 step 1 until M do H[Y]:= 0;
          for X:= 0 step 1 until N do
          begin Y:= LOW[X]; SUM:= H[Y];
              UPP:= M - LOW[N - X];
              for  Y:= Y + 1 step 1 until UPP do
              H[Y]:= SUM:= SUM + H[Y]
          end;
          KOLSMIR:= SUM / BINCOEF
      end;
    if XSIZE ≤ 0 ∨ XSIZE > ENTIER(XSIZE) then
        STATAL3 ERROR(“KOLSMIR”, 2, XSIZE);
    if YSIZE ≤ 0 ∨ YSIZE > ENTIER(YSIZE) then
        STATAL3 ERROR(“KOLSMIR”, 3, YSIZE);
    if XSIZE < YSIZE then
    begin N:= XSIZE; M:= YSIZE end
    else begin M:= XSIZE; N:= YSIZE end;
    if EPS < 0 ∨ EPS > 10-2 then EPS:= 10-3;
    if XSIZE < YSIZE then
    begin N:= XSIZE; M:= YSIZE end
    else begin M:= XSIZE; N:= YSIZE end;
    KGV:= M × N / GGD(M, N);
    I:= ENTIER((1 + 10-14) × D × KGV);
    if EPS ≥ 10-3 then L: APPROX else EXACT;
end KOLSMIR;
eop