C言語でアトキンの篩を実装したい〈初心者〉

Question

この春から大学でプログラムを学び始めた者です。
講義にて、素数の演算と表示に関する話がありました。
その中で、エラトステネスの篩とアトキンの篩について紹介があり、どちらの方がどのくらい早いのかと、入力した数以下の素数と処理時間を表示するものを組もうとしています。
しかし、エラトステネスの篩は出来たのですが、アトキンの篩が上手いこと組めずにいます。
そこで、ネットから参考となるコードを探したのですが、C言語で記述されたものを見つけることができませんでした。
次に、C++のコードをC言語に出来ないかと試み、エラーは消えたのですが、知識が足りないためギャップを埋めきれなかったようで、前述した素数と処理時間の表示という動作をしませんでした。（素数が表示されず処理時間のみ表示される）
どう違っているのか教えて下さい！

該当のソースコード

#include <stdio.h>
#include <math.h>
#include <sys/time.h>

#define M 1000
#define true 1
#define false 0

int main(){
    struct timeval s, e;
    gettimeofday(&s, NULL);

  /*process*/
    void atkin() {
      int  n,z,y,x,k;
      double N = 1000;
      int isprime[200];
  for (z = 1; z <= 5; z += 4) {
    for (y = z; y <= sqrt(N); y += 6) {
      for (x = 1; x <= sqrt(N) && (n = 4*x*x+y*y) <= M; ++x)
        isprime[n] = !isprime[n];
      for (x = y+1; x <= sqrt(N) && (n = 3*x*x-y*y) <= M; x += 2)
        isprime[n] = !isprime[n];
    }
  }
  for (z = 2; z <= 4; z += 2) {
    for (y = z; y <= sqrt(N); y += 6) {
      for (x = 1; x <= sqrt(N) && (n = 3*x*x+y*y) <= M; x += 2)
        isprime[n] = !isprime[n];
      for (x = y+1; x <= sqrt(N) && (n = 3*x*x-y*y) <= M; x += 2)
        isprime[n] = !isprime[n];
    }
  }
  for (y = 3; y <= sqrt(N); y += 6) {
    for (z = 1; z <= 2; ++z) {
      for (x = z; x <= sqrt(N) && (n = 4*x*x+y*y) <= M; x += 3)
        isprime[n] = !isprime[n];
    }
  }
  for (n = 5; n <= sqrt(N); ++n)
    if (isprime[n])
      for (k = n*n; k <= M; k+=n*n)
        isprime[k] = false;
  isprime[2] = isprime[3] = true;
}
  
 gettimeofday(&e, NULL);
    printf("time = %lf\n", (e.tv_sec - s.tv_sec) + (e.tv_usec - s.tv_usec)*1.0E-6);
    return 0;
}

参考にしたサイト

時間の計測
アトキンの篩のやり方
アトキンの篩
アトキンの篩
C言語のtrueとfalseについて
C++の記号一覧

補足情報（FW/ツールのバージョンなど）

Oracle VM VirtualBox
Oracle VM VirtualBox Extension Pack
Cent OS 7
C99以前

Comments

[jimbe]

URL を「リンクの挿入」での記述に編集して頂けると, クリックで飛べるようになります.

> 思った通りの動作をしませんでした
どう思ったのにどう動作したのでしょうか. その原因はどこまで追跡されていますでしょうか.

[cateye]

最適化されれば消えるかもですが、sqrt(N)は先に求めておいた方がいいのでは？

[moredeep]

配列を初期化していないせいな気がします。(自信０)

[cateye]

>思った通りの動作・・・とは？
素数の数が合わないとかでしょうか？

Answer 1

まずは関数の定義など基本的な文法から学び直してください。
また、定義に立ち返ってアルゴリズムをしっかり確認してみてください。

処理が複雑なら分解して関数化する、スコープを小さくして間違いにくくする、というのも手です。

具体的に何が問題か指摘するのが大変なので、以下にアトキンの篩のコードを書いておきます。

#include <stdio.h>

#define N 1000000
#define TRUE 1
#define FALSE 0

int is_in(int num, int set[], int length) {
  for(int i = 0; i < length; i++) {
    if(num == set[i]) return TRUE;
  }
  return FALSE;
}

void toggle(int* target) {
  if(*target == TRUE) {
    *target = FALSE;
  }
  *target = TRUE;
}

void step1(int is_prime[]) {
  int n, n60;
  int s[] = {1,13,17,29,37,41,49,53};
  int s_len = sizeof(s) / sizeof(s[0]);

  for(int x = 1; 4 * x * x < N; x++) {
    for(int y = 1; y < N; y = y + 2) {
      n = 4 * x * x + y * y;
      if (n > N) break;
      n60 = n % 60;
      if(is_in(n60, s, s_len) == TRUE) {
        toggle(is_prime + n);
      }
    }
  }
}

void step2(int is_prime[]) {
  int n, n60;
  int s[] = {7,19,31,43};
  int s_len = sizeof(s) / sizeof(s[0]);

  for(int x = 1; 3 * x * x < N; x = x + 2) {
    for(int y = 2; y < N; y = y + 2) {
      n = 3 * x * x + y * y;
      if (n > N) break;
      n60 = n % 60;
      if(is_in(n60, s, s_len) == TRUE) {
        toggle(is_prime + n);
      }
    }
  }
}

void step3(int is_prime[]) {
  int n, n60;
  int s[] = {11,23,47,59};
  int s_len = sizeof(s) / sizeof(s[0]);

  for(int x = 2; 3 * x * x < N; x++) {
    for(int y = x - 1; y > 0; y = y - 2) {
      n = 3 * x * x - y * y;
      if (n > N) break;
      n60 = n % 60;
      if(is_in(n60, s, s_len) == TRUE) {
        toggle(is_prime + n);
      }
    }
  }
}

void eliminate_composites(int is_prime[]) {
  int s[] = {1,7,11,13,17,19,23,29,31,37,41,43,47,49,53,59};
  int s_len = sizeof(s) / sizeof(s[0]);
  int n;

  for(n = 7; n*n < N; n++) {
    if( is_in(n % 60, s, s_len) == TRUE && is_prime[n] == TRUE) {
      for(int i = 1; i < N/(n*n); i++)
        is_prime[n*n*i] = FALSE;
    }
  }
}


int main(int argc, char* argv[]) {
  int is_prime[N+1] = {FALSE};

  step1(is_prime);
  step2(is_prime);
  step3(is_prime);
  eliminate_composites(is_prime);

  printf("2\n3\n5\n");
  for(int n = 7; n <= N; n++) {
    if(is_prime[n] == TRUE)
      printf("%d\n", n);
  }

  return 0;
}

Answer 2

ざっとみて（処理の正当性は別として）

１．main() の中に void atkin() の定義を内包しているのは誤りなので main() の外に atkin() の定義内容を全部追い出してください。

２．atkin() の中で int isprime[200]; と宣言してありますが、atkin() の中で isprime[] にアクセスするnやkの値が200以上になる事がありメモリを破壊しそうです。
(これくらいの大きさならほとんどの環境で大丈夫だとは思いますが関数内であまりにも巨大すぎる配列をオート変数として確保するとスタックを破壊するリスクもあるので大きめな配列を確保するのならグローバル変数にしてしまう方がいいかもしれません)

３．isprime配列を初期化せずに参照しているので先に初期化しておくと良いかもしれません(ローカル変数は初期値不定です)