100万までの素数を列挙

Wheel Sieveとかいう高速化テクニックを使っています。多分。
英語の解説はよく理解できなかったのですが、Wikipediaにあった画像から類推して適当に組立ててみました。

/* まず2*3*5*7=210までの自然数から2,3,5,7の倍数を排除した配列Wheelを用意。処理の関係で211までありますが、これは「210+1」であって、次の素数なのは偶々です。 */
numeric wheel[] = {1,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,121,127,131,137,139,143,149,151,157,163,167,169,173,179,181,187,191,193,197,199,209,211};
/* Wheelの要素数は */
whlen = 48;
/* 差分(の1/2)の配列を用意して（都合で開始が1つずれてます。回転処理が面倒なので2巡分あります）*/
numeric ntb105[] = {1,2,1,2,3,1,3,2,1,2,3,3,1,3,2,1,3,2,3,4,2,1,2,1,2,4,3,2,3,1,2,3,1,3,3,2,1,2,3,1,3,2,1,2,1,5,1,5,1,2,1,2,3,1,3,2,1,2,3,3,1,3,2,1,3,2,3,4,2,1,2,1,2,4,3,2,3,1,2,3,1,3,3,2,1,2,3,1,3,2,1,2,1,5,1,5};
/* 篩のindexからWheel上の位置を引くテーブルも作って */
numeric nrev105[105];
for (k = 1; k ＜= whlen; k = k + 1) {
　　nrev105[(wheel[k] - 3) / 2] = k - 1;
}
/* あとは適当に回す、みたいな */
numeric numx[48];
idx = 4;
while (idx ＜ enphsiev) {
　　for (z = 0; z ＜ whlen; z = z + 1) {
　　　　if (sieve[idx] == 0) {
　　　　　　m = nrev105[mod(idx, 105)];　　　　/* wheel開始位置 */
　　　　　　num = idx * 2 + 3; 　　　　　　　　/* この式はsieveのindexから実際の数値に変換します */
　　　　　　j = ((num * num) - 3) / 2; 　　　　/* これは2乗を逆変換してindexに戻してます */
　　　　　　for (k = 0; k ＜ whlen; k = k + 1) {
　　　　　　　　numx[k] = ntb105[m + k] * num; /* 掛け算は重いのでテーブルにしておく */
　　　　　　}
　　　　　　while (j ＜ sievemax) {　　　　　　/* 篩のメインループ */
　　　　　　　　for (k = 0; k ＜ whlen; k = k + 1) {
　　　　　　　　　　sieve[j] = 1;
　　　　　　　　　　j = j + numx[k];
　　　　　　　　}
　　　　　　}
　　　　}
　　　　idx = idx + ntb105[z];
　　　　if (idx >= enphsiev) {break;}
　　}
}
※篩の配列sieve[idx]は、idx=0 → 3，idx=1 → 5，idx=2 → 7 のように対応してます。
※ちなみにメインループで sieve[j] の引数 j が最大で whlen-1 だけ maxsieve を超過するので何らかの対策が必要です。
※完成するには以下のこれと、結果表示する処理が必要です。
/* 入力変数をUpperLimitとすると */
sievemax = ceiling(UpperLimit / 2);
enphsiev = int(sqrt(UpperLimit) / 2);
/* 入力数値は変更するまで定数扱いになるらしく配列の引数に使える */
numeric sieve[UpperLimit];