結果

問題 No.248 ミラー君の宿題
コンテスト
ユーザー Min_25
提出日時 2015-05-18 01:38:05
言語 C++11（廃止可能性あり）
(gcc 13.3.0)
結果
AC  
実行時間 156 ms / 5,000 ms
コード長 4,375 bytes
コンパイル時間 741 ms
コンパイル使用メモリ 86,380 KB
実行使用メモリ 16,328 KB
最終ジャッジ日時 2024-07-06 05:15:14
合計ジャッジ時間 2,079 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge3
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 18
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <ctime>
#include <cassert>

#include <iostream>
#include <utility>
#include <algorithm>
#include <queue>
#include <functional>
#include <vector>
#include <map>
#include <set>
#include <complex>

#define getchar getchar_unlocked
#define putchar putchar_unlocked

using namespace std;

typedef long long int64;
typedef long long unsigned uint64;
typedef long double float80;
typedef unsigned short uint16;
typedef unsigned uint;
typedef unsigned char uint8;

const uint N_PRIMES = 15;
const uint STATE_MAX = 1 << N_PRIMES;
const uint BITS = 65;

typedef double prob_t;

uint64 primes[N_PRIMES];
uint prime_twos[STATE_MAX];
uint state_max_twos[STATE_MAX];

prob_t probs_cumu[STATE_MAX][BITS];
prob_t probs_same[STATE_MAX][BITS];
prob_t probs_pos[STATE_MAX];
prob_t probs_neg[STATE_MAX];

prob_t state_prods[STATE_MAX];
prob_t state_phis[STATE_MAX];

prob_t pows[BITS];
prob_t dp[STATE_MAX];

uint conv[STATE_MAX];

inline uint ilog2(uint64 x) {
  union Data {
    uint64 u64;
    double d;
  } n;
  n.d = double(x) + 0.5;
  return (n.u64 >> 52) - 1023;
}

void init(uint N) {
  const uint state_max = 1 << N;

  pows[0] = 1.0;
  for (uint i = 1; i < BITS; ++i) {
    pows[i] = pows[i-1] * 0.5;
  }

  uint max_two = 0;
  for (uint i = 0; i < N; ++i) {
    uint64 p = primes[i] - 1;
    prime_twos[i] = ilog2(p & -p);
    max_two = max(max_two, prime_twos[i]);
  }

  for (uint i = 0; i < state_max; ++i) {
    for (uint j = 0; j <= max_two; ++j) {
      probs_cumu[i][j] = 1.0;
      probs_same[i][j] = 1.0;
    }
    probs_pos[i] = 1.0;
    probs_neg[i] = 1.0;
    state_prods[i] = 1.0;
    state_phis[i] = 1.0;
  }

  for (uint state = 1; state < state_max; ++state) {
    uint f = state & -state;
    uint pstate = state ^ f;

    uint idx = ilog2(f);
    uint two = prime_twos[idx];

    state_max_twos[state] = max(state_max_twos[pstate], two);

    prob_t p = primes[idx];

    probs_pos[state] = probs_pos[pstate] / p;
    probs_neg[state] = probs_neg[pstate] * (1. - 1. / p);

    state_prods[state] = state_prods[pstate] * p;
    state_phis[state]  = state_phis[pstate] * (p - 1.0);

    prob_t pw = pows[two];
    prob_t q = pw;

    for (uint i = 0; i <= two; ++i) {
      probs_cumu[state][i] = probs_cumu[pstate][i] * q;
      probs_same[state][i] = probs_same[pstate][i] * pw;
      if (i > 0) {
        pw *= 2.0;
      }
      q += pw;
    }

    for (uint i = two + 1; i <= max_two; ++i) {
      probs_cumu[state][i] = probs_cumu[pstate][i];
      probs_same[state][i] = 0.0;
    }
  }
}

uint pop_count(uint n) {
  return __builtin_popcount(n);
}

prob_t calc(const uint N) {
  const uint total = 1 << N;
  for (uint state = 1; state < total; ++state) {
    uint bits = pop_count(state);
    if (bits == 1) {
      dp[state] = 1.0;
      continue;
    }

    uint rest = state;
    uint conv_size = 1;
    conv[0] = 0;
    while (rest) {
      uint mask = rest & -rest;
      for (uint i = conv_size; i < 2 * conv_size; ++i) {
        conv[i] = conv[i - conv_size] | mask;
      }
      conv_size <<= 1;
      rest ^= mask;
    }

    prob_t Q = state_phis[state];
    prob_t P = state_prods[state];

    const uint state_total = 1 << bits;

    prob_t expected_count = 1.0;
    prob_t prob_not_repeat = 0.0;

    for (uint i = 1; i < state_total - 1; ++i) {
      uint left_state = conv[i];
      uint right_state = conv[state_total - 1 - i];

      prob_t p1 = 0.0;
      for (uint j = 0; j < state_max_twos[right_state]; ++j) {
        p1 += probs_cumu[left_state][j] * probs_same[right_state][j + 1];
      }
      p1 *= Q / P;
      
      prob_t p2 = probs_pos[left_state] * probs_neg[right_state];
      prob_t p3 = p1 + p2;

      expected_count += p3 * (dp[left_state] + dp[right_state]);
      prob_not_repeat += p3;
    }
    expected_count /= prob_not_repeat;
    dp[state] = expected_count;
  }
  return dp[total - 1];
}

void solve() {
  uint T; assert(1 == scanf("%u", &T));
  assert(1 <= T && T <= 1000);
  for (; T; --T) {
    uint N; assert(1 == scanf("%u", &N));
    assert(1 <= N && N <= 13);
    for (uint i = 0; i < N; ++i) {
      assert(1 == scanf("%llu", &primes[i]));
      assert(3 <= primes[i] && primes[i] < 1ull << 63);
    }
    init(N);
    prob_t ans = calc(N);
    // WA
    printf("%.8lf\n", ans);
  }
}

int main() {
  solve();
  return 0;
}
0