#ifdef DBG1
  #define _GLIBCXX_DEBUG
#endif

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cassert>
#include <cstring>
#include <string>
#include <vector>
#include <algorithm>
#include <map>
#include <set>

#ifdef DBG1
    #define dbg(...) fprintf(stderr, __VA_ARGS__)
#else
    #define dbg(...)
#endif

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef pair <int, int> pii;

const double E = exp(1.0);
const int MINP = 10;
const int MAXP = 1000;

double sqr(double x) {
  return x * x;
}

double log_fact[10000];

void find_mean_uniform(const vector <int> & a) {
  int n = int(a.size());
  int maxA = 0;
  for (int i = 0; i < n; ++i) {
    maxA = max(maxA, a[i]);
  }
  for (int P = MINP; P <= MAXP; ++P) {
    if (maxA > 2 * P) { continue; }
    printf("%d\n", P);
    return;
  }
}

void find_mean_poisson(const vector <int> & a) {
  int n = int(a.size());
  int best = 0;
  double bestProb = -1;
  for (int P = MINP; P <= MAXP; ++P) {
    double logP = log(P + 0.0);
    double curProb = 0.0;
    for (int i = 0; i < n; ++i) {
      curProb += -P + logP * a[i] - log_fact[a[i]];
    }
    if (best == 0 || curProb > bestProb) {
      bestProb = curProb;
      best = P;
    }
  }
  printf("%d\n", best);
}

double PoissonKnuth(double rate) {
  double expRateInv = exp(-rate);
  double k = 0, prod = rand() / double(RAND_MAX);
  while (prod > expRateInv) {
    prod *= rand() / double(RAND_MAX);
    ++k;
  }
  return k;
}

int genRand(int type, int PP) {
  if (type == 0) {
    return rand() % (2 * PP + 1);
  } else {
    return PoissonKnuth(PP);
  }
}

void solve() {
  const int n = 250;
  vector <int> a(n);
  //int PP = rand() % 900 + 10;
  //int type = rand() & 1;
  //PP = 10;
  //dbg("%s %d\n", type == 0 ? "uniform" : "poisson", PP);
  int sum = 0;
  for (int i = 0; i < n; ++i) {
    scanf("%d", &a[i]);
    //a[i] = genRand(type, PP);
    //printf("%d ", a[i]);
    sum += a[i];
  }
  //printf("\n");

  double mean = sum / double(n);
  int P = int(mean + 0.5);
  dbg("mean %.10lf\n", mean);

  int ans = -1;
  double bestDiff = 0.0;
  int ans1 = 0, ans2 = 0;
  for (int lambda = max(10, P - 2); lambda <= P + 2; ++lambda) {
    double log_lambda = log(lambda + 0.0);
    double factor = 1.0;
    for (int i = 1; i <= lambda; ++i) {
      factor *= (i / E);
    }
    double p1 = - n * log(2 * lambda + 1.0);
    double p2 = 0.0;
    for (int i = 0; i < n; ++i) {
      p2 += -lambda + log_lambda * a[i] - log_fact[a[i]];
      /*for (int j = 1; j <= a[i]; ++j) {
        p2 += log(double(lambda) / j);
      }*/
    }
    dbg("p1 %.10lf p2 %.10lf\n", p1, p2);
/*    double curDiff = (p1 > p2 ? p1 / p2 : p2 / p1);
    int curAns = (p1 > p2 ? 0 : 1);
    if (ans == -1 || bestDiff < curDiff) {
      ans = curAns;
      bestDiff = curDiff;
    }*/
    if (p1 > p2) {
      ans1++;
    } else {
      ans2++;
    }
  }
  ans = (ans1 > ans2 ? 0 : 1);
  //printf("%s\n", ans == 0 ? "uniform" : "poisson");
  if (ans == 0) {
    find_mean_uniform(a);
  } else {
    find_mean_poisson(a);
  }
}

int main() {
  for (int i = 1; i < 10000; ++i) {
    log_fact[i] = log_fact[i - 1] + log(i + 0.0);
  }
  int v;
  scanf("%d", &v);
  for (int i = 0; i < v; ++i) {
    solve();
  }
  return 0;
}

### Prompt

In Cpp, your task is to solve the following problem:

Good job! Now that Heidi is able to distinguish between Poisson and uniform distributions, she is in a good position to actually estimate the populations.

Can you help Heidi estimate each village's population?

Input

Same as the easy version.

Output

Output one line per village, in the same order as provided in the input, containing your (integer) population estimate.

Your answer is considered correct if it is an integer that falls into the interval <image>, where P is the real population of the village, used to create the distribution (either Poisson or uniform) from which the marmots drew their answers.

### Response

```cpp
#ifdef DBG1
  #define _GLIBCXX_DEBUG
#endif

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cassert>
#include <cstring>
#include <string>
#include <vector>
#include <algorithm>
#include <map>
#include <set>

#ifdef DBG1
    #define dbg(...) fprintf(stderr, __VA_ARGS__)
#else
    #define dbg(...)
#endif

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef pair <int, int> pii;

const double E = exp(1.0);
const int MINP = 10;
const int MAXP = 1000;

double sqr(double x) {
  return x * x;
}

double log_fact[10000];

void find_mean_uniform(const vector <int> & a) {
  int n = int(a.size());
  int maxA = 0;
  for (int i = 0; i < n; ++i) {
    maxA = max(maxA, a[i]);
  }
  for (int P = MINP; P <= MAXP; ++P) {
    if (maxA > 2 * P) { continue; }
    printf("%d\n", P);
    return;
  }
}

void find_mean_poisson(const vector <int> & a) {
  int n = int(a.size());
  int best = 0;
  double bestProb = -1;
  for (int P = MINP; P <= MAXP; ++P) {
    double logP = log(P + 0.0);
    double curProb = 0.0;
    for (int i = 0; i < n; ++i) {
      curProb += -P + logP * a[i] - log_fact[a[i]];
    }
    if (best == 0 || curProb > bestProb) {
      bestProb = curProb;
      best = P;
    }
  }
  printf("%d\n", best);
}

double PoissonKnuth(double rate) {
  double expRateInv = exp(-rate);
  double k = 0, prod = rand() / double(RAND_MAX);
  while (prod > expRateInv) {
    prod *= rand() / double(RAND_MAX);
    ++k;
  }
  return k;
}

int genRand(int type, int PP) {
  if (type == 0) {
    return rand() % (2 * PP + 1);
  } else {
    return PoissonKnuth(PP);
  }
}

void solve() {
  const int n = 250;
  vector <int> a(n);
  //int PP = rand() % 900 + 10;
  //int type = rand() & 1;
  //PP = 10;
  //dbg("%s %d\n", type == 0 ? "uniform" : "poisson", PP);
  int sum = 0;
  for (int i = 0; i < n; ++i) {
    scanf("%d", &a[i]);
    //a[i] = genRand(type, PP);
    //printf("%d ", a[i]);
    sum += a[i];
  }
  //printf("\n");

  double mean = sum / double(n);
  int P = int(mean + 0.5);
  dbg("mean %.10lf\n", mean);

  int ans = -1;
  double bestDiff = 0.0;
  int ans1 = 0, ans2 = 0;
  for (int lambda = max(10, P - 2); lambda <= P + 2; ++lambda) {
    double log_lambda = log(lambda + 0.0);
    double factor = 1.0;
    for (int i = 1; i <= lambda; ++i) {
      factor *= (i / E);
    }
    double p1 = - n * log(2 * lambda + 1.0);
    double p2 = 0.0;
    for (int i = 0; i < n; ++i) {
      p2 += -lambda + log_lambda * a[i] - log_fact[a[i]];
      /*for (int j = 1; j <= a[i]; ++j) {
        p2 += log(double(lambda) / j);
      }*/
    }
    dbg("p1 %.10lf p2 %.10lf\n", p1, p2);
/*    double curDiff = (p1 > p2 ? p1 / p2 : p2 / p1);
    int curAns = (p1 > p2 ? 0 : 1);
    if (ans == -1 || bestDiff < curDiff) {
      ans = curAns;
      bestDiff = curDiff;
    }*/
    if (p1 > p2) {
      ans1++;
    } else {
      ans2++;
    }
  }
  ans = (ans1 > ans2 ? 0 : 1);
  //printf("%s\n", ans == 0 ? "uniform" : "poisson");
  if (ans == 0) {
    find_mean_uniform(a);
  } else {
    find_mean_poisson(a);
  }
}

int main() {
  for (int i = 1; i < 10000; ++i) {
    log_fact[i] = log_fact[i - 1] + log(i + 0.0);
  }
  int v;
  scanf("%d", &v);
  for (int i = 0; i < v; ++i) {
    solve();
  }
  return 0;
}

```