素数(4題)

プログラミングコンテストチャレンジブック [第2版] ~問題解決のアルゴリズム活用力とコーディングテクニックを鍛える~
この本(通称:蟻本)の章末で紹介されている練習問題をひたすらといている。

素数の問題。エラトステネスのふるい。有名だし、これといってコメントがないんだけど。

あ、エラトステネスのふるいの計算量は O(n log log n) 程度らしい。蟻本によると。

それと、数学系のこういう問題って、動的計画法とかと違って、「どうやってアルゴリズムを組み立てるか」をあまり考えなくて済むので、一度実装できると類題は瞬殺で終わることが多いね。
Spaghetti Source - 各種アルゴリズムの C++ による実装 みたいなサイトから基本的な関数を拝借すれば、解ける問題の幅が広がる。

AOJ0009: Prime Number

解法

エラトステネスのふるい。

#include <fstream>
#include <iostream>
#include <cmath>
#include <vector>
#include <queue>
#include <utility>
#include <algorithm>
using namespace std;
const int N = 1000000;

int main() {
  // エラトステネスのふるいを作成
  bool isPrime[N]; // ふるい
  int cntPrime[N]; // n個以下の素数の個数
  for (int i = 0; i < N; i++)
    isPrime[i] = true;
  isPrime[0] = false; isPrime[1] = false;
  cntPrime[0] = 0;
  for (int i = 1; i < N; i++) {
    if (isPrime[i]) {
      cntPrime[i] = cntPrime[i-1] + 1;
      for (int j = 2 * i; j < N; j += i)
        isPrime[j] = false;
    } else {
      cntPrime[i] = cntPrime[i - 1];
    }
  }

  // 入力と出力
  int n;
  while (cin >> n) {
    cout << cntPrime[n] << endl;
  }
}

POJ3126: Prime Path

3126 -- Prime Path

解法

4桁の素数の表を作るのはエラトステネスのふるいでよい。
最短パスを見つけるためには幅優先探索で。

#include <fstream>
#include <iostream>
#include <vector>
#include <queue>
#include <utility>
#include <algorithm>
using namespace std;
const int N = 10000;
typedef pair<int, int> P;

// ---- n ** m -----
int pow(int n, int m) {
  int res = 1;
  for (int i = 0; i < m; i++)
    res *= n;
  return res;
}

int main() {
  // エラトステネスのふるいで素数の表を作る
  bool isPrime[N];
  for (int i = 0; i < N; i++) {
    isPrime[i] = true;
  }
  isPrime[0] = false; isPrime[1] = false;
  for (int i = 0; i < N; i++) {
    if (isPrime[i]) {
      for (int j = i * 2; j < N; j += i) {
        isPrime[j] = false;
      }
    }
  }

  // 入力
  bool visited[N]; // すでに訪れた素数
  int n, p, q;
  cin >> n;
  while (n--) {
    // visited の初期化
    for (int i = 0; i < N; i++) {
      visited[i] = false;
    }
    cin >> p >> q;
    // 幅優先探索
    // pair (素数, 距離)
    queue<P> que;
    que.push(P(p, 0));
    visited[p] = true;
    int res_d = -1;
    while(!que.empty()) {
      P x = que.front();
      que.pop();
      int _p = x.first;
      int d = x.second;
      if (_p == q) {
        res_d = d;
        break;
      }
      int a, b, c, np;
      // p の一桁を違う数に置換した数を作る
      for (int i = 0; i < 4; i++) {
        a = pow(10, i);
        b = (_p / a) % 10;
        c = _p - b * a;
        for (int j = 0; j < 10; j++) {
          np = c + a * j;
          if (np >= 1000 && np != p && isPrime[np] && !visited[np]) {
            visited[np] = true;
            que.push(P(np, d + 1));
          }
        }
      }
    }
    // 出力
    if (res_d >= 0) {
      cout << res_d << endl;
    } else {
      cout << "Impossible" << endl;
    }
  }
}

POJ3421: X-factor chains

問題概要

3421 -- X-factor Chains
正の整数 X に対して長さ m の X-factor chain を次のような数列と定義する。
1 = X0, X1, X2, …, Xm = X
数列は次の性質を満たす。
Xi < Xi+1 かつ Xi は Xi+1 を割り切る。
X-factor chain の最大の長さと、そのような長さの chain がいくつあるかを求める。
x <= 2 ^ 20

解法

X-factor chain は、1から始めて少しずつ素数をかけて X にすれば、最も長いものが作れる。
よって、最大の長さは、X を素因数分解したときの素数の個数に等しい。
また、長さが最大となる chain の個数は、素数をかける順番に依存するので、同じものを含む順列で計算する。
必要な素数は 2 ^ 10 以下なので、エラトステネスのふるいが速い。

#include <fstream>
#include <iostream>
#include <cmath>
#include <vector>
#include <queue>
#include <utility>
#include <algorithm>
using namespace std;
const int N = 1024;

// n!
long long factorial(int n) {
  long long res = 1;
  while (n > 1) {
    res *= n--;
  }
  return res;
}

int main() {
  ifstream cin("../test.txt");
  // エラトステネスのふるい
  bool isPrime[N];
  vector<int> prime; // i番目の素数
  for (int i = 0; i < N; i++)
    isPrime[i] = true;
  for (int i = 2; i < N; i++) {
    if (isPrime[i]) {
      prime.push_back(i);
      for (int j = 2 * i; j < N; j += i) {
        isPrime[j] = false;
      }
    }
  }
  int M = prime.size();

  int n;
  while (cin >> n) {
    // n を素因数分解する
    vector<int> factors(M);
    fill(factors.begin(), factors.end(), 0);
    int cnt_len = 0;
    for (int i = 0; i < M; i++) {
      while (n % prime[i] == 0) {
        n /= prime[i];
        factors[i]++;
        cnt_len++;
      }
      if (n == 1) break;
    }
    // 残った n が素数の場合
    if (n > 1) {
      cnt_len++;
    }
    // 最長 chain の個数を数える
    long long cnt_chains = factorial(cnt_len);
    for (int i = 0; i < M; i++) {
      if (factors[i] > 0) {
        cnt_chains /= factorial(factors[i]);
      }
    }
    // 出力
    cout << cnt_len << " " << cnt_chains << endl;
  }
}

POJ3292: Semi-prime H-numbers

3292 -- Semi-prime H-numbers

問題概要

H-数を 4 * n + 1 で表せる正の整数とする。つまり、H-数とは 1, 5, 9, 13, ... である。
H-数は乗法について閉じている(2つの H-数の積はまた H-数になる)ので、ふつうの整数と同じように1や素数や合成数を定義できる。
・1 はそのまま 1 とする。
・H-数 h が "H-素数"であるとは、H-数の積に分解すると 1 * h の形にしかならないことをいう。ただし 1 は H-素数ではない。
・H-数 h が"H-合成数"であるとは、h を H-素数の積で表せることをいう。
さて、さらに H-擬素数というものを定義する。
・H-数 h が"H-擬素数"であろとは、ちょうど2つの H-素数の積で表せることをいう。
たとえば、 25 = 5 * 5 は H-擬素数であるが、 125 = 5 * 5 * 5 は H-擬素数でない。
与えられた H-数 h に対して、h 以下の H-擬素数の個数を求めるプログラムを書く。

解法

変わった問題だが、ふつうにエラトステネスのふるいが通用する。
1000001 以下の素数を記録しておき、h 以下の擬素数の個数も先に計算しておく。

#include <fstream>
#include <iostream>
#include <cmath>
#include <vector>
#include <queue>
#include <utility>
#include <algorithm>
using namespace std;
const int N = 250001;

// m = 4 * i + 1 から i を求める
int to_i(int m) {
  return (m - 1) / 4;
}

int main() {
  #ifndef ONLINE_JUDGE
  ifstream cin("../test.txt");
  #endif
  // エラトステネスのふるい
  bool isPrime[N]; // isPrime[i] = true は 4 * i + 1 が H-素数ということ
  vector<int> prime; // i 番目の H-素数
  for (int i = 0; i < N; i++)
      isPrime[i] = true;
  isPrime[0] = false;
  int m;
  for (int i = 0; i < N; i++) {
    if (isPrime[i]) {
      m = 4 * i + 1;
      prime.push_back(m);
      for (int j = 5; to_i(m * j) < N; j += 4) {
        isPrime[to_i(m * j)] = false;
      }
    }
  }

  int n = prime.size(); // 1000001以下の H-素数の個数
  // 4 * i + 1 が擬素数かどうか
  bool is_semi_prime[N];
  for (int i = 0; i < N; i++)
    is_semi_prime[i] = false;
  long long k;
  for (int i = 0; i < n; i++) {
    for (int j = i; j < n; j++) {
      k = (long long)prime[i] * prime[j];
      if (k <= 1000001 && k >= 0) {
        is_semi_prime[to_i(k)] = true;
      } else {
        break;
      }
    }
  }

  // 擬素数を数えておく
  int cnt_semi_primes[N];
  cnt_semi_primes[0] = 0;
  for (int i = 1; i < N; i++) {
    cnt_semi_primes[i] = cnt_semi_primes[i-1];
    if (is_semi_prime[i]) cnt_semi_primes[i]++;
  }

  // 入力と出力
  while (true) {
    int h;
    cin >> h;
    if (h == 0) break;
    cout << h << " " << cnt_semi_primes[to_i(h)] << endl;
  }
}