Problem 231

231 The prime factorisation of binomial coefficients
Problem 231 - Project Euler
二項係数の約数和を求める問題．
結局は含まれる各素数の数が分かれば良いので，ふるいで素数を求めれば，あとは簡単．

import Number (primesUpTo)

count :: (Integral a) => a -> a -> a
count n p= sum.takeWhile (>0).tail.iterate (`div` p) $ n

p231 m n = sum.map f.primesUpTo $ m
    where f p = p * ( count m p - count n p - count (m - n) p )

main :: IO ()
main = print $ p231 (2*10^7) (15*10^6)

やはり，Haskellではあまり速くない．
C++とかで実装したら速いんだろうな．

追記

C++で実装してみた．
一瞬すぎる…

#include <iostream>
#include <vector>
#include <cmath>
#define FOR(x, xs) for (__typeof( (xs).begin() ) x = (xs).begin(); x != (xs).end(); x++)

using namespace std;

static const int m = 20000000;
static const int n = 15000000;
static const int u = sqrt(m);

int count (int q, int p) {
  int c = 0;
  for ( ; q > 0; c += (q /= p) );
  return c;
}

int main() {
  vector<bool> isPrime( m + 1);
  FOR (q, isPrime) *q = true;
  isPrime[0] = false;
  isPrime[1] = false;

  for (int q = 2; q <= u; q++)
    if ( isPrime[q] )
      for (int s = q * q; s <= m; s += q)
	isPrime[s] = false;

  long long s = 0;
  for (int q = 2; q <= m; q++)
    if ( isPrime[q] ) s += q * ( count(m , q) - count(n, q) - count(m - n, q));

  cout << s <<endl;
}