donuts.cc

#include <iostream>
#include <iomanip>
#include <type_traits>

constexpr int N = 10;
constexpr int max_k = 100;
double dp[N][max_k+1];

#ifdef DEBUG
template<class T>
constexpr T abs2(T val){
  static_assert( std::is_arithmetic<T>::value, "argument must be arithmetic type" );
  return val < 0 ? -val : val;
}

template<class T>
constexpr int _num_digits(T n){
  return n < 10 ? 1 : _num_digits(n/10) + 1;
}

template<class T> 
constexpr int num_digits(T n){
  static_assert( std::is_integral<T>::value, "argument must be integral" );
  return _num_digits(abs2(n));
}

constexpr int N_digit = num_digits(N);
constexpr int max_k_digit = num_digits(max_k);
#endif

int main(int argc, char const* argv[])
{
  dp[0][0] = 1.0;
  for(int n = 1; n < N; n++){
    for(int k = n; k <= n + (max_k - N); k++){
      dp[n][k] = dp[n-1][k-1] * (double)(N-n+1)/N + dp[n][k-1] * (double)n/N;
#ifdef DEBUG
      std::cerr << "n:" << std::setw(N_digit) << n 
          << ", k:" << std::setw(max_k_digit) << k << ", dp[n][k]:" << dp[n][k] << std::endl;
#endif
    }
  }

  double acc = 0.0;
  for(int k = 1; k <= max_k; k++){
    double prob = dp[N-1][k-1] * 1.0/N;
    acc += prob;
    std::cout << k << "," << prob << "," << acc << std::endl;
  }

  return 0;
}