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 (>).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 = ;
for ( ; q > ; c += (q /= p) );
return c;
}
int main() {
vector<bool> isPrime( m + 1);
FOR (q, isPrime) *q = true;
isPrime[] = 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 = ;
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;
}