100만까지의 소수 나열

'휠 체(Wheel Sieve)'라고 하는 고속화 기법을 사용하고 있습니다. 아마도요.
영어 해설은 잘 이해하지 못했지만, 위키백과에 있던 이미지를 보고 추측하여 대충 구성해 보았습니다.

/* 먼저 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};
/* 체의 인덱스에서 휠상의 위치를 뺀 테이블도 만들고 */
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의 인덱스를 실제 수치로 변환합니다 */
　　　　　　j = ((num * num) - 3) / 2; 　　　　/* 이는 제곱을 역변환하여 인덱스로 되돌리는 것입니다 */
　　　　　　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];