#include <bits/stdc++.h>
using namespace std;
template <class T, class U>
ostream& operator<<(ostream& o, const pair<T, U>& p) {
  o << "(" << p.first << "," << p.second << ")";
  return o;
}
template <class T>
ostream& operator<<(ostream& o, const vector<T>& v) {
  o << "[";
  for (T t : v) {
    o << t << ",";
  }
  o << "]";
  return o;
}
long double fact[2010];
long double pow(long double p, int x) {
  long double res = 1;
  while (x > 0) {
    if (x % 2) res *= p;
    p *= p;
    x /= 2;
  }
  return res;
}
long double fac(int x) {
  if (x <= 2000)
    return fact[x];
  else
    return fact[2000];
}
int solve() {
  int n = 250;
  vector<int> vec(n);
  for (int i = (int)(0); i < (int)(n); i++) cin >> vec[i];
  int maxid = max_element((vec).begin(), (vec).end()) - vec.begin();
  int minid = min_element((vec).begin(), (vec).end()) - vec.begin();
  double p =
      (accumulate((vec).begin(), (vec).end(), 0) - vec[maxid] - vec[minid]) /
      (double)(n - 2);
  double var = 0;
  for (int i = (int)(0); i < (int)(n); i++)
    if (i != maxid && i != minid) var += (vec[i] - p) * (vec[i] - p);
  var /= (double)(n - 2);
  double diffp = abs(p - var);
  double diffu = abs(p * p / 3.0 - var);
  if (diffp < diffu) {
    long double bp = 0;
    int pp = -1;
    for (int i = (int)(10); i < (int)(1001); i++) {
      long double tp = 1.0;
      long double ee = exp((long double)-i);
      for (int j = (int)(0); j < (int)(n); j++) {
        tp *= pow(i, vec[j]) * ee / fac(vec[j]);
        if (tp < bp) break;
      }
      if (tp > bp) bp = tp, pp = i;
    }
    return pp;
  } else {
    long double bp = 0;
    int pp = -1;
    for (int i = (int)(10); i < (int)(1001); i++) {
      long double tp = 1;
      long double tmp = 1.0 / ((long double)2.0 * i + 1.0);
      for (int j = (int)(0); j < (int)(n); j++) {
        if (vec[j] > 2 * i)
          tp = 0;
        else
          tp *= tmp;
        if (tp < bp) break;
      }
      if (tp > bp) bp = tp, pp = i;
    }
    return pp;
  }
}
int main() {
  int t;
  cin >> t;
  fact[0] = 1;
  for (int i = (int)(0); i < (int)(2000); i++)
    fact[i + 1] = fact[i] * (long double)(i + 1);
  for (int i = (int)(0); i < (int)(t); i++) {
    printf("%d\n", 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
#include <bits/stdc++.h>
using namespace std;
template <class T, class U>
ostream& operator<<(ostream& o, const pair<T, U>& p) {
  o << "(" << p.first << "," << p.second << ")";
  return o;
}
template <class T>
ostream& operator<<(ostream& o, const vector<T>& v) {
  o << "[";
  for (T t : v) {
    o << t << ",";
  }
  o << "]";
  return o;
}
long double fact[2010];
long double pow(long double p, int x) {
  long double res = 1;
  while (x > 0) {
    if (x % 2) res *= p;
    p *= p;
    x /= 2;
  }
  return res;
}
long double fac(int x) {
  if (x <= 2000)
    return fact[x];
  else
    return fact[2000];
}
int solve() {
  int n = 250;
  vector<int> vec(n);
  for (int i = (int)(0); i < (int)(n); i++) cin >> vec[i];
  int maxid = max_element((vec).begin(), (vec).end()) - vec.begin();
  int minid = min_element((vec).begin(), (vec).end()) - vec.begin();
  double p =
      (accumulate((vec).begin(), (vec).end(), 0) - vec[maxid] - vec[minid]) /
      (double)(n - 2);
  double var = 0;
  for (int i = (int)(0); i < (int)(n); i++)
    if (i != maxid && i != minid) var += (vec[i] - p) * (vec[i] - p);
  var /= (double)(n - 2);
  double diffp = abs(p - var);
  double diffu = abs(p * p / 3.0 - var);
  if (diffp < diffu) {
    long double bp = 0;
    int pp = -1;
    for (int i = (int)(10); i < (int)(1001); i++) {
      long double tp = 1.0;
      long double ee = exp((long double)-i);
      for (int j = (int)(0); j < (int)(n); j++) {
        tp *= pow(i, vec[j]) * ee / fac(vec[j]);
        if (tp < bp) break;
      }
      if (tp > bp) bp = tp, pp = i;
    }
    return pp;
  } else {
    long double bp = 0;
    int pp = -1;
    for (int i = (int)(10); i < (int)(1001); i++) {
      long double tp = 1;
      long double tmp = 1.0 / ((long double)2.0 * i + 1.0);
      for (int j = (int)(0); j < (int)(n); j++) {
        if (vec[j] > 2 * i)
          tp = 0;
        else
          tp *= tmp;
        if (tp < bp) break;
      }
      if (tp > bp) bp = tp, pp = i;
    }
    return pp;
  }
}
int main() {
  int t;
  cin >> t;
  fact[0] = 1;
  for (int i = (int)(0); i < (int)(2000); i++)
    fact[i + 1] = fact[i] * (long double)(i + 1);
  for (int i = (int)(0); i < (int)(t); i++) {
    printf("%d\n", solve());
  }
  return 0;
}

```