Золотой треугольник

PROCEDURE InitPrime*;
VAR k,j,p,q: LONGINT; step: SHORTINT; init: BOOLEAN; set: LONGINT;
BEGIN (* PRE: r.int - верхняя граница диапазона чисел *)
ASSERT (r.int > 3);
primeDesc.this := 1; primeDesc.count := 0;
primeDesc.max := (r.int DIV 2) - 1;
(* только среди нечетных *)
set := (primeDesc.max DIV MAX(SET)) + 1;
(* кол-во наборов *)
NEW(isPrime,set);
(* формируем динамический массив множеств *)
FOR k := 1 TO primeDesc.max DO
(* выставляем бит; isPrime<0> не используется *)
INCL(isPrime[k DIV MAX(SET)], k MOD MAX(SET));
END;
k := 1; j := 1; p := 3; q := 4;
(* p = 2*j + 1; q = 2*j + 2*j^2 *)
step := 1; (* шаг алгоритма *)
LOOP (* цикл фиктивный: это переключатель шагов *)
CASE step OF
1: IF ~((j MOD MAX(SET)) IN
isPrime[j DIV MAX(SET)]) THEN
step := 3;
ELSE k := q; step := 2;
END;
| 2: IF (k <= primeDesc.max) THEN
EXCL(isPrime[k DIV MAX(SET)],
k MOD MAX(SET));
INC(k,p); (* еще раз *)
ELSE step := 3;
END;
| 3: INC(j); INC(p,2); q := q + 2*p - 2;
IF (j <= primeDesc.max) THEN step := 1;
ELSE EXIT (* выходим из цикла LOOP *)
END;
ELSE (* других вариантов нет *)
END;
END;
init := FALSE; sPrimeIsEmpty := TRUE;
FOR k := 1 TO primeDesc.max DO
IF ((k MOD MAX(SET)) IN isPrime[k DIV MAX(SET)])
THEN
IF ~init THEN
primeDesc.this := k; r.result := 2*k + 1;
init := TRUE;
sPrimeIsEmpty := FALSE;
END;
INC(primeDesc.count);
END;
END; (* POST: r.result содержит первое простое число *)
END InitPrime;