ログイン
編集不可のページディスカッション情報添付ファイル
iz/競技プログラミング/Project Euler 023 Non-abundant sums

MMA

http://projecteuler.net/problem=23

2つの過剰数の和で表せない整数の総和を求めよ。ただし28123より大きい整数についてそのような整数は存在しないことが保証されている。

Rubyで適当に書いたらかなり時間がかかったので最適化。最後の二重ループとfactorでmapを使っているところに高速化の余地あり。

   1 #include <iostream>
   2 #include <vector>
   3 #include <map>
   4 #include <cstring>
   5 using namespace std;
   6 
   7 #define LIMIT 28123
   8 
   9 bool isPrime[LIMIT+1];
  10 vector<int> primes;
  11 vector<int> abundants;
  12 
  13 void factor(int n, map<int, int>& fs) {
  14     for (int i = 0; primes[i] * primes[i] <= n; i++) {
  15         while (n % primes[i] == 0) {
  16             n /= primes[i];
  17             fs[primes[i]]++;
  18         }
  19     }
  20     if (n != 1) fs[n]++;
  21 }
  22 
  23 int power(int x, int n) { // x ^ n
  24     int p = 1;
  25     for (int i = 0; i < n; i++) {
  26         p *= x;
  27     }
  28     return p;
  29 }
  30 
  31 int divisorSum(int n) {
  32     int sum = 1;
  33     map<int, int> fs;
  34     factor(n, fs);
  35     for (map<int, int>::iterator it = fs.begin(); it != fs.end(); it++) {
  36         // 約数の和: p ^ (o + 1) / (p - 1)
  37         //     ただし p: 素因数, o: 出現数
  38         sum *= (power(it->first, it->second+1) - 1) / (it->first - 1); 
  39     }
  40     return sum - n; // その数自身は除く
  41 }
  42 
  43 inline bool isAbundant(int n) {
  44     return divisorSum(n) > n;
  45 }
  46 
  47 void init() {
  48     memset(isPrime, 1, sizeof(isPrime));
  49     for (int i = 2; i <= LIMIT; i++) {
  50         if (isPrime[i]) {
  51             primes.push_back(i);
  52             for (int j = i + i; j <= LIMIT; j += i) {
  53                 isPrime[j] = false;
  54             }
  55         }
  56     }
  57     for (int i = 1; i <= LIMIT; i++) {
  58         if (isAbundant(i)) {
  59             abundants.push_back(i);
  60         }
  61     }
  62 }
  63 
  64 int main() {
  65     init();
  66     bool possible[LIMIT+1] = {};
  67     for (int i = 0; i < abundants.size(); i++) {
  68         for (int j = i; j < abundants.size(); j++) {
  69             if (abundants[i] + abundants[j] > LIMIT) break;
  70             possible[abundants[i]+abundants[j]] = true;
  71         }
  72     }
  73     int sum = 0;
  74     for (int i = 1; i <= LIMIT; i++) {
  75         if (!possible[i]) {
  76             sum += i;
  77         }
  78     }
  79     cout << sum << endl;
  80     return 0;
  81 }

iz/競技プログラミング/Project Euler 023 Non-abundant sums (最終更新日時 2012-12-22 01:38:21 更新者 iz)