#include <bits/stdc++.h>
using namespace std;
static bool debug = false;
template <typename T, typename Q, typename S>
bool Bounded(const T &x, const Q &a, const S &b) {
  return a <= x && x <= b;
}
const long double pi = 4 * atan(long double(1.0));
const long double Inf = numeric_limits<long double>::infinity();
long double isapprox(
    long double x, long double y,
    long double rtol = sqrt(numeric_limits<long double>::epsilon()),
    long double atol = long double(0.0), bool nans = false) {
  return x == y || (nans && isnan(x) && isnan(y)) ||
         (isfinite(x) && isfinite(y) &&
          abs(x - y) <= atol + rtol * max(abs(x), abs(y)));
}
template <typename T>
int sz(const T &x) {
  return (int)x.size();
}
template <typename T, typename Q>
bool mem(const T &s, const Q &x) {
  return s.find(x) != s.end();
}
template <typename T>
void push(vector<T> &c, const T &e) {
  c.push_back(e);
}
template <typename T>
void push(queue<T> &c, const T &e) {
  c.push(e);
}
template <typename T>
void push(set<T> &c, const T &e) {
  c.insert(e);
}
template <typename T>
T pop(vector<T> &v) {
  assert(!v.empty());
  T u = v.back();
  v.pop_back();
  return u;
}
template <typename T>
T pop(queue<T> &v) {
  assert(!v.empty());
  T u = v.front();
  v.pop();
  return u;
}
template <typename T, typename... Args>
auto emplace(vector<T> &c, Args &&...args)
    -> decltype(c.emplace_back(args...)) {
  c.emplace_back(args...);
}
template <typename T, typename... Args>
auto emplace(queue<T> &c, Args &&...args) -> decltype(c.emplace(args...)) {
  c.emplace(args...);
}
template <typename T, typename... Args>
auto emplace(set<T> &c, Args &&...args) -> decltype(c.emplace(args...)) {
  c.emplace(args...);
}
template <typename T>
void umin(T &x, const T &y) {
  x = min(x, y);
}
template <typename T>
void umax(T &x, const T &y) {
  x = max(x, y);
}
template <typename T>
int sign(const T &x) {
  return x == 0 ? 0 : x > 0 ? 1 : -1;
}
template <typename T>
T square(const T &x) {
  return x * x;
}
template <typename S, typename T>
istream &operator>>(istream &in, pair<S, T> &p) {
  return in >> p.first >> p.second;
}
template <typename T>
T read() {
  T x;
  cin >> x;
  return x;
}
template <typename T, typename I = typename T::iterator>
void read(T &cnt) {
  for (auto &x : cnt) cin >> x;
}
template <typename S, typename T>
ostream &operator<<(ostream &out, const pair<S, T> &p) {
  return out << "{" << p.first << ", " << p.second << "}";
}
template <typename T>
struct show_container {
  const T &container;
  show_container(const T &container_) : container(container_) {}
};
template <typename T>
ostream &operator<<(ostream &o, const show_container<T> &thing) {
  bool first = true;
  o << "{";
  for (const auto &x : thing.container) {
    if (!first) o << ", ";
    first = false;
    o << x;
  }
  o << "}";
  return o;
}
template <typename T>
show_container<T> show(const T &container) {
  return show_container<T>(container);
}
struct dout_ {
} dout;
template <typename T>
dout_ &operator<<(dout_ &out, const T &thing) {
  return out;
}
dout_ &operator<<(dout_ &out, ostream &(*thing)(ostream &)) { return out; }
struct outboth {
  ostream &a;
  ostream &b;
  outboth(ostream &a_, ostream &b_) : a(a_), b(b_) {}
};
template <typename T>
outboth &operator<<(outboth &both, const T &thing) {
  both.a << thing;
  both.b << thing;
  return both;
}
outboth &operator<<(outboth &both, ostream &(*thing)(ostream &)) {
  both.a << thing;
  both.b << thing;
  return both;
}
long double logFactorial(int x) {
  static map<int, long double> memo;
  if (mem(memo, x)) return memo[x];
  return memo[x] = lgamma(long double(x + 1));
}
long double logpPoisson(const vector<int> &X, int P) {
  long double ans = 0;
  for (int x : X) ans += x * log(long double(P)) - P - logFactorial(x);
  return ans;
}
long double logpUniform(const vector<int> &X, int P) {
  for (int x : X)
    if (!Bounded(x, 0, 2 * P)) return -Inf;
  return -sz(X) * log(2 * P + long double(1));
}
void solve(ostream &rout) {
  int V;
  cin >> V;
  if (!cin) return;
  const int N = 250;
  for (int v = 0; v < int(V); ++v) {
    vector<int> X(N);
    read(X);
    int xmin = *min_element(X.begin(), X.end()),
        xmax = *max_element(X.begin(), X.end());
    long double sumX = 0, sumLogFactX = 0;
    for (int x : X) {
      sumX += x;
      sumLogFactX += logFactorial(x);
    }
    long double pP = -Inf, pU = -Inf;
    int whichP = -1, whichU = -1;
    for (int p = 10; p <= 1000; ++p) {
      long double logprobPoisson =
          sumX * log(long double(p)) - N * p - sumLogFactX;
      if (logprobPoisson > pP) {
        pP = logprobPoisson;
        whichP = p;
      }
      if (xmin >= 0 && xmax <= 2 * p + 1) {
        long double logprobUniform = -N * log(long double(2 * p + 1));
        if (logprobUniform > pU) {
          pU = logprobUniform;
          whichU = p;
        }
      }
    }
    rout << (pP > pU ? whichP : whichU) << endl;
  }
}
int main() {
  std::ios::sync_with_stdio(false);
  cin.tie(nullptr);
  solve(cout);
  return 0;
}

### Prompt

Create a solution in cpp for 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
#include <bits/stdc++.h>
using namespace std;
static bool debug = false;
template <typename T, typename Q, typename S>
bool Bounded(const T &x, const Q &a, const S &b) {
  return a <= x && x <= b;
}
const long double pi = 4 * atan(long double(1.0));
const long double Inf = numeric_limits<long double>::infinity();
long double isapprox(
    long double x, long double y,
    long double rtol = sqrt(numeric_limits<long double>::epsilon()),
    long double atol = long double(0.0), bool nans = false) {
  return x == y || (nans && isnan(x) && isnan(y)) ||
         (isfinite(x) && isfinite(y) &&
          abs(x - y) <= atol + rtol * max(abs(x), abs(y)));
}
template <typename T>
int sz(const T &x) {
  return (int)x.size();
}
template <typename T, typename Q>
bool mem(const T &s, const Q &x) {
  return s.find(x) != s.end();
}
template <typename T>
void push(vector<T> &c, const T &e) {
  c.push_back(e);
}
template <typename T>
void push(queue<T> &c, const T &e) {
  c.push(e);
}
template <typename T>
void push(set<T> &c, const T &e) {
  c.insert(e);
}
template <typename T>
T pop(vector<T> &v) {
  assert(!v.empty());
  T u = v.back();
  v.pop_back();
  return u;
}
template <typename T>
T pop(queue<T> &v) {
  assert(!v.empty());
  T u = v.front();
  v.pop();
  return u;
}
template <typename T, typename... Args>
auto emplace(vector<T> &c, Args &&...args)
    -> decltype(c.emplace_back(args...)) {
  c.emplace_back(args...);
}
template <typename T, typename... Args>
auto emplace(queue<T> &c, Args &&...args) -> decltype(c.emplace(args...)) {
  c.emplace(args...);
}
template <typename T, typename... Args>
auto emplace(set<T> &c, Args &&...args) -> decltype(c.emplace(args...)) {
  c.emplace(args...);
}
template <typename T>
void umin(T &x, const T &y) {
  x = min(x, y);
}
template <typename T>
void umax(T &x, const T &y) {
  x = max(x, y);
}
template <typename T>
int sign(const T &x) {
  return x == 0 ? 0 : x > 0 ? 1 : -1;
}
template <typename T>
T square(const T &x) {
  return x * x;
}
template <typename S, typename T>
istream &operator>>(istream &in, pair<S, T> &p) {
  return in >> p.first >> p.second;
}
template <typename T>
T read() {
  T x;
  cin >> x;
  return x;
}
template <typename T, typename I = typename T::iterator>
void read(T &cnt) {
  for (auto &x : cnt) cin >> x;
}
template <typename S, typename T>
ostream &operator<<(ostream &out, const pair<S, T> &p) {
  return out << "{" << p.first << ", " << p.second << "}";
}
template <typename T>
struct show_container {
  const T &container;
  show_container(const T &container_) : container(container_) {}
};
template <typename T>
ostream &operator<<(ostream &o, const show_container<T> &thing) {
  bool first = true;
  o << "{";
  for (const auto &x : thing.container) {
    if (!first) o << ", ";
    first = false;
    o << x;
  }
  o << "}";
  return o;
}
template <typename T>
show_container<T> show(const T &container) {
  return show_container<T>(container);
}
struct dout_ {
} dout;
template <typename T>
dout_ &operator<<(dout_ &out, const T &thing) {
  return out;
}
dout_ &operator<<(dout_ &out, ostream &(*thing)(ostream &)) { return out; }
struct outboth {
  ostream &a;
  ostream &b;
  outboth(ostream &a_, ostream &b_) : a(a_), b(b_) {}
};
template <typename T>
outboth &operator<<(outboth &both, const T &thing) {
  both.a << thing;
  both.b << thing;
  return both;
}
outboth &operator<<(outboth &both, ostream &(*thing)(ostream &)) {
  both.a << thing;
  both.b << thing;
  return both;
}
long double logFactorial(int x) {
  static map<int, long double> memo;
  if (mem(memo, x)) return memo[x];
  return memo[x] = lgamma(long double(x + 1));
}
long double logpPoisson(const vector<int> &X, int P) {
  long double ans = 0;
  for (int x : X) ans += x * log(long double(P)) - P - logFactorial(x);
  return ans;
}
long double logpUniform(const vector<int> &X, int P) {
  for (int x : X)
    if (!Bounded(x, 0, 2 * P)) return -Inf;
  return -sz(X) * log(2 * P + long double(1));
}
void solve(ostream &rout) {
  int V;
  cin >> V;
  if (!cin) return;
  const int N = 250;
  for (int v = 0; v < int(V); ++v) {
    vector<int> X(N);
    read(X);
    int xmin = *min_element(X.begin(), X.end()),
        xmax = *max_element(X.begin(), X.end());
    long double sumX = 0, sumLogFactX = 0;
    for (int x : X) {
      sumX += x;
      sumLogFactX += logFactorial(x);
    }
    long double pP = -Inf, pU = -Inf;
    int whichP = -1, whichU = -1;
    for (int p = 10; p <= 1000; ++p) {
      long double logprobPoisson =
          sumX * log(long double(p)) - N * p - sumLogFactX;
      if (logprobPoisson > pP) {
        pP = logprobPoisson;
        whichP = p;
      }
      if (xmin >= 0 && xmax <= 2 * p + 1) {
        long double logprobUniform = -N * log(long double(2 * p + 1));
        if (logprobUniform > pU) {
          pU = logprobUniform;
          whichU = p;
        }
      }
    }
    rout << (pP > pU ? whichP : whichU) << endl;
  }
}
int main() {
  std::ios::sync_with_stdio(false);
  cin.tie(nullptr);
  solve(cout);
  return 0;
}

```