Language: C++

#include <bits/stdc++.h>
#define BIT(n) (1LL << (n))
#define BITF(n, i) (((n) >> i) & 1)
#define REP(i, n) for (int i = 0; i < n; i++)
#define REPB(i, n) for (int i = 0; i < BIT(n); i++)
#define REPS(i, x) for (int i = 1; i <= x; i++)
#define REPR(i, n) for (int i = n; i >= 0; i--)
#define REPZ(i, x) for (int i = 0; i <= x; i++)
#define FOR(i, m, n) for (int i = m; i < n; i++)
#define FORA(i, n) for (auto &&i : n)
#define FORS(i, m, n) for (int i = m; i <= n; i++)
using namespace std;
#define PRINTARR(x, y)                                     \
  cerr << #x << "=\n";                                     \
  for (auto itr = x; itr != y; itr++) cerr << *itr << " "; \
  cerr << endl;
#define PRINTARR2(x, i0, i1)                                       \
  cerr << #x << "=\n";                                             \
  for (int ii0 = 0; ii0 < i0; ii0++) {                             \
    for (int ii1 = 0; ii1 < i1; ii1++) cerr << x[ii0][ii1] << " "; \
    cerr << endl;                                                  \
  }
#define DUMPOUT cerr
// vector
template <typename T> istream &operator>>(istream &is, vector<T> &vec) {
  for (T &x : vec) is >> x;
  return is;
}
// pair
template <typename T, typename U> ostream &operator<<(ostream &os, pair<T, U> &pair_var) {
  os << "(" << pair_var.first << ", " << pair_var.second << ")";
  return os;
}
// vector
template <typename T> ostream &operator<<(ostream &os, const vector<T> &vec) {
  os << "{";
  REP(i, (int)vec.size())
  os << vec[i] << (i + 1 == (int)vec.size() ? "" : ", ");
  os << "}";
  return os;
}
// map
template <typename T, typename U> ostream &operator<<(ostream &os, map<T, U> &map_var) {
  os << "{";
  FORA(itr, map_var) {
    os << *itr;
    itr++;
    if (itr != map_var.end()) os << ", ";
    itr--;
  }
  os << "}";
  return os;
}
// set
template <typename T> ostream &operator<<(ostream &os, set<T> &set_var) {
  os << "{";
  FORA(itr, set_var) {
    os << *itr;
    itr++;
    if (itr != set_var.end()) os << ", ";
    itr--;
  }
  os << "}";
  return os;
}
void dump_func() { DUMPOUT << endl; }
template <class Head, class... Tail> void dump_func(Head &&head, Tail &&... tail) {
  DUMPOUT << head;
  if (sizeof...(Tail) > 0) DUMPOUT << ", ";
  dump_func(std::move(tail)...);
}
#ifdef DEBUG_
#define DEB
#define DUMP(...)                                                             \
  DUMPOUT << "  " << string(#__VA_ARGS__) << ": "                             \
          << "[" << to_string(__LINE__) << ":" << __FUNCTION__ << "]" << endl \
          << "    ",                                                          \
      dump_func(__VA_ARGS__)
#else
#define DEB if (false)
#define DUMP(...)
#endif

#define ALL(v) v.begin(), v.end()
#define fst first
#define snd second
#define mp make_pair
#define pb push_back
#define epb emplace_back
#define int long long
#define pint pair<int, int>
#define ld long double
using namespace std;
template <class T> bool chmax(T &a, const T &b) {
  if (a < b) {
    a = b;
    return 1;
  }
  return 0;
}
template <class T> bool chmin(T &a, const T &b) {
  if (a > b) {
    a = b;
    return 1;
  }
  return 0;
}
template <class T> using vec = std::vector<T>;
template <class T> void print(const T &x) { cout << x << "\n"; }
const int MOD = 1000000007, INF0 = 1061109567, INF = INF0 * INF0;
const double EPS = 1e-10, PI = acos(-1.0);
const int dx[4] = {1, 0, -1, 0}, dy[4] = {0, 1, 0, -1};
#define MAXN 200200
#define MAXSIZE (1 << 25)

/*
Fast Fourier/Cosine/Sine Transform
    dimension   :one
    data length :power of 2
    decimation  :frequency
    radix       :4, 2
    data        :inplace
    table       :use
functions
    cdft: Complex Discrete Fourier Transform
    rdft: Real Discrete Fourier Transform
    ddct: Discrete Cosine Transform
    ddst: Discrete Sine Transform
    dfct: Cosine Transform of RDFT (Real Symmetric DFT)
    dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
function prototypes
    void cdft(signed, signed, double *, signed *, double *);
    void rdft(signed, signed, double *, signed *, double *);
    void ddct(signed, signed, double *, signed *, double *);
    void ddst(signed, signed, double *, signed *, double *);
    void dfct(signed, double *, double *, signed *, double *);
    void dfst(signed, double *, double *, signed *, double *);


-------- Complex DFT (Discrete Fourier Transform) --------
    [definition]
        <case1>
            X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
        <case2>
            X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
        (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
    [usage]
        <case1>
            ip[0] = 0; // first time only
            cdft(2*n, 1, a, ip, w);
        <case2>
            ip[0] = 0; // first time only
            cdft(2*n, -1, a, ip, w);
    [parameters]
        2*n            :data length (signed)
                        n >= 1, n = power of 2
        a[0...2*n-1]   :input/output data (double *)
                        input data
                            a[2*j] = Re(x[j]),
                            a[2*j+1] = Im(x[j]), 0<=j<n
                        output data
                            a[2*k] = Re(X[k]),
                            a[2*k+1] = Im(X[k]), 0<=k<n
        ip[0...*]      :work area for bit reversal (signed *)
                        length of ip >= 2+sqrt(n)
                        strictly,
                        length of ip >=
                            2+(1<<(signed)(log(n+0.5)/log(2))/2).
                        ip[0],ip[1] are posigneders of the cos/sin table.
        w[0...n/2-1]   :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of
            cdft(2*n, -1, a, ip, w);
        is
            cdft(2*n, 1, a, ip, w);
            for (j = 0; j <= 2 * n - 1; j++) {
                a[j] *= 1.0 / n;
            }
        .


-------- Real DFT / Inverse of Real DFT --------
    [definition]
        <case1> RDFT
            R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
            I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
        <case2> IRDFT (excluding scale)
            a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
                   sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
                   sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
    [usage]
        <case1>
            ip[0] = 0; // first time only
            rdft(n, 1, a, ip, w);
        <case2>
            ip[0] = 0; // first time only
            rdft(n, -1, a, ip, w);
    [parameters]
        n              :data length (signed)
                        n >= 2, n = power of 2
        a[0...n-1]     :input/output data (double *)
                        <case1>
                            output data
                                a[2*k] = R[k], 0<=k<n/2
                                a[2*k+1] = I[k], 0<k<n/2
                                a[1] = R[n/2]
                        <case2>
                            input data
                                a[2*j] = R[j], 0<=j<n/2
                                a[2*j+1] = I[j], 0<j<n/2
                                a[1] = R[n/2]
        ip[0...*]      :work area for bit reversal (signed *)
                        length of ip >= 2+sqrt(n/2)
                        strictly,
                        length of ip >=
                            2+(1<<(signed)(log(n/2+0.5)/log(2))/2).
                        ip[0],ip[1] are posigneders of the cos/sin table.
        w[0...n/2-1]   :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of
            rdft(n, 1, a, ip, w);
        is
            rdft(n, -1, a, ip, w);
            for (j = 0; j <= n - 1; j++) {
                a[j] *= 2.0 / n;
            }
        .


-------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
    [definition]
        <case1> IDCT (excluding scale)
            C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
        <case2> DCT
            C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
    [usage]
        <case1>
            ip[0] = 0; // first time only
            ddct(n, 1, a, ip, w);
        <case2>
            ip[0] = 0; // first time only
            ddct(n, -1, a, ip, w);
    [parameters]
        n              :data length (signed)
                        n >= 2, n = power of 2
        a[0...n-1]     :input/output data (double *)
                        output data
                            a[k] = C[k], 0<=k<n
        ip[0...*]      :work area for bit reversal (signed *)
                        length of ip >= 2+sqrt(n/2)
                        strictly,
                        length of ip >=
                            2+(1<<(signed)(log(n/2+0.5)/log(2))/2).
                        ip[0],ip[1] are posigneders of the cos/sin table.
        w[0...n*5/4-1] :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of
            ddct(n, -1, a, ip, w);
        is
            a[0] *= 0.5;
            ddct(n, 1, a, ip, w);
            for (j = 0; j <= n - 1; j++) {
                a[j] *= 2.0 / n;
            }
        .


-------- DST (Discrete Sine Transform) / Inverse of DST --------
    [definition]
        <case1> IDST (excluding scale)
            S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
        <case2> DST
            S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
    [usage]
        <case1>
            ip[0] = 0; // first time only
            ddst(n, 1, a, ip, w);
        <case2>
            ip[0] = 0; // first time only
            ddst(n, -1, a, ip, w);
    [parameters]
        n              :data length (signed)
                        n >= 2, n = power of 2
        a[0...n-1]     :input/output data (double *)
                        <case1>
                            input data
                                a[j] = A[j], 0<j<n
                                a[0] = A[n]
                            output data
                                a[k] = S[k], 0<=k<n
                        <case2>
                            output data
                                a[k] = S[k], 0<k<n
                                a[0] = S[n]
        ip[0...*]      :work area for bit reversal (signed *)
                        length of ip >= 2+sqrt(n/2)
                        strictly,
                        length of ip >=
                            2+(1<<(signed)(log(n/2+0.5)/log(2))/2).
                        ip[0],ip[1] are posigneders of the cos/sin table.
        w[0...n*5/4-1] :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of
            ddst(n, -1, a, ip, w);
        is
            a[0] *= 0.5;
            ddst(n, 1, a, ip, w);
            for (j = 0; j <= n - 1; j++) {
                a[j] *= 2.0 / n;
            }
        .


-------- Cosine Transform of RDFT (Real Symmetric DFT) --------
    [definition]
        C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
    [usage]
        ip[0] = 0; // first time only
        dfct(n, a, t, ip, w);
    [parameters]
        n              :data length - 1 (signed)
                        n >= 2, n = power of 2
        a[0...n]       :input/output data (double *)
                        output data
                            a[k] = C[k], 0<=k<=n
        t[0...n/2]     :work area (double *)
        ip[0...*]      :work area for bit reversal (signed *)
                        length of ip >= 2+sqrt(n/4)
                        strictly,
                        length of ip >=
                            2+(1<<(signed)(log(n/4+0.5)/log(2))/2).
                        ip[0],ip[1] are posigneders of the cos/sin table.
        w[0...n*5/8-1] :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of
            a[0] *= 0.5;
            a[n] *= 0.5;
            dfct(n, a, t, ip, w);
        is
            a[0] *= 0.5;
            a[n] *= 0.5;
            dfct(n, a, t, ip, w);
            for (j = 0; j <= n; j++) {
                a[j] *= 2.0 / n;
            }
        .


-------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
    [definition]
        S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
    [usage]
        ip[0] = 0; // first time only
        dfst(n, a, t, ip, w);
    [parameters]
        n              :data length + 1 (signed)
                        n >= 2, n = power of 2
        a[0...n-1]     :input/output data (double *)
                        output data
                            a[k] = S[k], 0<k<n
                        (a[0] is used for work area)
        t[0...n/2-1]   :work area (double *)
        ip[0...*]      :work area for bit reversal (signed *)
                        length of ip >= 2+sqrt(n/4)
                        strictly,
                        length of ip >=
                            2+(1<<(signed)(log(n/4+0.5)/log(2))/2).
                        ip[0],ip[1] are posigneders of the cos/sin table.
        w[0...n*5/8-1] :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of
            dfst(n, a, t, ip, w);
        is
            dfst(n, a, t, ip, w);
            for (j = 1; j <= n - 1; j++) {
                a[j] *= 2.0 / n;
            }
        .


Appendix :
    The cos/sin table is recalculated when the larger table required.
    w[] and ip[] are compatible with all routines.
*/

void cdft(signed n, signed isgn, double *a, signed *ip, double *w) {
  void makewt(signed nw, signed *ip, double *w);
  void bitrv2(signed n, signed *ip, double *a);
  void bitrv2conj(signed n, signed *ip, double *a);
  void cftfsub(signed n, double *a, double *w);
  void cftbsub(signed n, double *a, double *w);

  if (n > (ip[0] << 2)) {
    makewt(n >> 2, ip, w);
  }
  if (n > 4) {
    if (isgn >= 0) {
      bitrv2(n, ip + 2, a);
      cftfsub(n, a, w);
    } else {
      bitrv2conj(n, ip + 2, a);
      cftbsub(n, a, w);
    }
  } else if (n == 4) {
    cftfsub(n, a, w);
  }
}

void rdft(signed n, signed isgn, double *a, signed *ip, double *w) {
  void makewt(signed nw, signed *ip, double *w);
  void makect(signed nc, signed *ip, double *c);
  void bitrv2(signed n, signed *ip, double *a);
  void cftfsub(signed n, double *a, double *w);
  void cftbsub(signed n, double *a, double *w);
  void rftfsub(signed n, double *a, signed nc, double *c);
  void rftbsub(signed n, double *a, signed nc, double *c);
  signed nw, nc;
  double xi;

  nw = ip[0];
  if (n > (nw << 2)) {
    nw = n >> 2;
    makewt(nw, ip, w);
  }
  nc = ip[1];
  if (n > (nc << 2)) {
    nc = n >> 2;
    makect(nc, ip, w + nw);
  }
  if (isgn >= 0) {
    if (n > 4) {
      bitrv2(n, ip + 2, a);
      cftfsub(n, a, w);
      rftfsub(n, a, nc, w + nw);
    } else if (n == 4) {
      cftfsub(n, a, w);
    }
    xi = a[0] - a[1];
    a[0] += a[1];
    a[1] = xi;
  } else {
    a[1] = 0.5 * (a[0] - a[1]);
    a[0] -= a[1];
    if (n > 4) {
      rftbsub(n, a, nc, w + nw);
      bitrv2(n, ip + 2, a);
      cftbsub(n, a, w);
    } else if (n == 4) {
      cftfsub(n, a, w);
    }
  }
}

void ddct(signed n, signed isgn, double *a, signed *ip, double *w) {
  void makewt(signed nw, signed *ip, double *w);
  void makect(signed nc, signed *ip, double *c);
  void bitrv2(signed n, signed *ip, double *a);
  void cftfsub(signed n, double *a, double *w);
  void cftbsub(signed n, double *a, double *w);
  void rftfsub(signed n, double *a, signed nc, double *c);
  void rftbsub(signed n, double *a, signed nc, double *c);
  void dctsub(signed n, double *a, signed nc, double *c);
  signed j, nw, nc;
  double xr;

  nw = ip[0];
  if (n > (nw << 2)) {
    nw = n >> 2;
    makewt(nw, ip, w);
  }
  nc = ip[1];
  if (n > nc) {
    nc = n;
    makect(nc, ip, w + nw);
  }
  if (isgn < 0) {
    xr = a[n - 1];
    for (j = n - 2; j >= 2; j -= 2) {
      a[j + 1] = a[j] - a[j - 1];
      a[j] += a[j - 1];
    }
    a[1] = a[0] - xr;
    a[0] += xr;
    if (n > 4) {
      rftbsub(n, a, nc, w + nw);
      bitrv2(n, ip + 2, a);
      cftbsub(n, a, w);
    } else if (n == 4) {
      cftfsub(n, a, w);
    }
  }
  dctsub(n, a, nc, w + nw);
  if (isgn >= 0) {
    if (n > 4) {
      bitrv2(n, ip + 2, a);
      cftfsub(n, a, w);
      rftfsub(n, a, nc, w + nw);
    } else if (n == 4) {
      cftfsub(n, a, w);
    }
    xr = a[0] - a[1];
    a[0] += a[1];
    for (j = 2; j < n; j += 2) {
      a[j - 1] = a[j] - a[j + 1];
      a[j] += a[j + 1];
    }
    a[n - 1] = xr;
  }
}

void ddst(signed n, signed isgn, double *a, signed *ip, double *w) {
  void makewt(signed nw, signed *ip, double *w);
  void makect(signed nc, signed *ip, double *c);
  void bitrv2(signed n, signed *ip, double *a);
  void cftfsub(signed n, double *a, double *w);
  void cftbsub(signed n, double *a, double *w);
  void rftfsub(signed n, double *a, signed nc, double *c);
  void rftbsub(signed n, double *a, signed nc, double *c);
  void dstsub(signed n, double *a, signed nc, double *c);
  signed j, nw, nc;
  double xr;

  nw = ip[0];
  if (n > (nw << 2)) {
    nw = n >> 2;
    makewt(nw, ip, w);
  }
  nc = ip[1];
  if (n > nc) {
    nc = n;
    makect(nc, ip, w + nw);
  }
  if (isgn < 0) {
    xr = a[n - 1];
    for (j = n - 2; j >= 2; j -= 2) {
      a[j + 1] = -a[j] - a[j - 1];
      a[j] -= a[j - 1];
    }
    a[1] = a[0] + xr;
    a[0] -= xr;
    if (n > 4) {
      rftbsub(n, a, nc, w + nw);
      bitrv2(n, ip + 2, a);
      cftbsub(n, a, w);
    } else if (n == 4) {
      cftfsub(n, a, w);
    }
  }
  dstsub(n, a, nc, w + nw);
  if (isgn >= 0) {
    if (n > 4) {
      bitrv2(n, ip + 2, a);
      cftfsub(n, a, w);
      rftfsub(n, a, nc, w + nw);
    } else if (n == 4) {
      cftfsub(n, a, w);
    }
    xr = a[0] - a[1];
    a[0] += a[1];
    for (j = 2; j < n; j += 2) {
      a[j - 1] = -a[j] - a[j + 1];
      a[j] -= a[j + 1];
    }
    a[n - 1] = -xr;
  }
}

void dfct(signed n, double *a, double *t, signed *ip, double *w) {
  void makewt(signed nw, signed *ip, double *w);
  void makect(signed nc, signed *ip, double *c);
  void bitrv2(signed n, signed *ip, double *a);
  void cftfsub(signed n, double *a, double *w);
  void rftfsub(signed n, double *a, signed nc, double *c);
  void dctsub(signed n, double *a, signed nc, double *c);
  signed j, k, l, m, mh, nw, nc;
  double xr, xi, yr, yi;

  nw = ip[0];
  if (n > (nw << 3)) {
    nw = n >> 3;
    makewt(nw, ip, w);
  }
  nc = ip[1];
  if (n > (nc << 1)) {
    nc = n >> 1;
    makect(nc, ip, w + nw);
  }
  m = n >> 1;
  yi = a[m];
  xi = a[0] + a[n];
  a[0] -= a[n];
  t[0] = xi - yi;
  t[m] = xi + yi;
  if (n > 2) {
    mh = m >> 1;
    for (j = 1; j < mh; j++) {
      k = m - j;
      xr = a[j] - a[n - j];
      xi = a[j] + a[n - j];
      yr = a[k] - a[n - k];
      yi = a[k] + a[n - k];
      a[j] = xr;
      a[k] = yr;
      t[j] = xi - yi;
      t[k] = xi + yi;
    }
    t[mh] = a[mh] + a[n - mh];
    a[mh] -= a[n - mh];
    dctsub(m, a, nc, w + nw);
    if (m > 4) {
      bitrv2(m, ip + 2, a);
      cftfsub(m, a, w);
      rftfsub(m, a, nc, w + nw);
    } else if (m == 4) {
      cftfsub(m, a, w);
    }
    a[n - 1] = a[0] - a[1];
    a[1] = a[0] + a[1];
    for (j = m - 2; j >= 2; j -= 2) {
      a[2 * j + 1] = a[j] + a[j + 1];
      a[2 * j - 1] = a[j] - a[j + 1];
    }
    l = 2;
    m = mh;
    while (m >= 2) {
      dctsub(m, t, nc, w + nw);
      if (m > 4) {
        bitrv2(m, ip + 2, t);
        cftfsub(m, t, w);
        rftfsub(m, t, nc, w + nw);
      } else if (m == 4) {
        cftfsub(m, t, w);
      }
      a[n - l] = t[0] - t[1];
      a[l] = t[0] + t[1];
      k = 0;
      for (j = 2; j < m; j += 2) {
        k += l << 2;
        a[k - l] = t[j] - t[j + 1];
        a[k + l] = t[j] + t[j + 1];
      }
      l <<= 1;
      mh = m >> 1;
      for (j = 0; j < mh; j++) {
        k = m - j;
        t[j] = t[m + k] - t[m + j];
        t[k] = t[m + k] + t[m + j];
      }
      t[mh] = t[m + mh];
      m = mh;
    }
    a[l] = t[0];
    a[n] = t[2] - t[1];
    a[0] = t[2] + t[1];
  } else {
    a[1] = a[0];
    a[2] = t[0];
    a[0] = t[1];
  }
}

void dfst(signed n, double *a, double *t, signed *ip, double *w) {
  void makewt(signed nw, signed *ip, double *w);
  void makect(signed nc, signed *ip, double *c);
  void bitrv2(signed n, signed *ip, double *a);
  void cftfsub(signed n, double *a, double *w);
  void rftfsub(signed n, double *a, signed nc, double *c);
  void dstsub(signed n, double *a, signed nc, double *c);
  signed j, k, l, m, mh, nw, nc;
  double xr, xi, yr, yi;

  nw = ip[0];
  if (n > (nw << 3)) {
    nw = n >> 3;
    makewt(nw, ip, w);
  }
  nc = ip[1];
  if (n > (nc << 1)) {
    nc = n >> 1;
    makect(nc, ip, w + nw);
  }
  if (n > 2) {
    m = n >> 1;
    mh = m >> 1;
    for (j = 1; j < mh; j++) {
      k = m - j;
      xr = a[j] + a[n - j];
      xi = a[j] - a[n - j];
      yr = a[k] + a[n - k];
      yi = a[k] - a[n - k];
      a[j] = xr;
      a[k] = yr;
      t[j] = xi + yi;
      t[k] = xi - yi;
    }
    t[0] = a[mh] - a[n - mh];
    a[mh] += a[n - mh];
    a[0] = a[m];
    dstsub(m, a, nc, w + nw);
    if (m > 4) {
      bitrv2(m, ip + 2, a);
      cftfsub(m, a, w);
      rftfsub(m, a, nc, w + nw);
    } else if (m == 4) {
      cftfsub(m, a, w);
    }
    a[n - 1] = a[1] - a[0];
    a[1] = a[0] + a[1];
    for (j = m - 2; j >= 2; j -= 2) {
      a[2 * j + 1] = a[j] - a[j + 1];
      a[2 * j - 1] = -a[j] - a[j + 1];
    }
    l = 2;
    m = mh;
    while (m >= 2) {
      dstsub(m, t, nc, w + nw);
      if (m > 4) {
        bitrv2(m, ip + 2, t);
        cftfsub(m, t, w);
        rftfsub(m, t, nc, w + nw);
      } else if (m == 4) {
        cftfsub(m, t, w);
      }
      a[n - l] = t[1] - t[0];
      a[l] = t[0] + t[1];
      k = 0;
      for (j = 2; j < m; j += 2) {
        k += l << 2;
        a[k - l] = -t[j] - t[j + 1];
        a[k + l] = t[j] - t[j + 1];
      }
      l <<= 1;
      mh = m >> 1;
      for (j = 1; j < mh; j++) {
        k = m - j;
        t[j] = t[m + k] + t[m + j];
        t[k] = t[m + k] - t[m + j];
      }
      t[0] = t[m + mh];
      m = mh;
    }
    a[l] = t[0];
  }
  a[0] = 0;
}

/* -------- initializing routines -------- */

#include <math.h>

void makewt(signed nw, signed *ip, double *w) {
  void bitrv2(signed n, signed *ip, double *a);
  signed j, nwh;
  double delta, x, y;

  ip[0] = nw;
  ip[1] = 1;
  if (nw > 2) {
    nwh = nw >> 1;
    delta = atan(1.0) / nwh;
    w[0] = 1;
    w[1] = 0;
    w[nwh] = cos(delta * nwh);
    w[nwh + 1] = w[nwh];
    if (nwh > 2) {
      for (j = 2; j < nwh; j += 2) {
        x = cos(delta * j);
        y = sin(delta * j);
        w[j] = x;
        w[j + 1] = y;
        w[nw - j] = y;
        w[nw - j + 1] = x;
      }
      bitrv2(nw, ip + 2, w);
    }
  }
}

void makect(signed nc, signed *ip, double *c) {
  signed j, nch;
  double delta;

  ip[1] = nc;
  if (nc > 1) {
    nch = nc >> 1;
    delta = atan(1.0) / nch;
    c[0] = cos(delta * nch);
    c[nch] = 0.5 * c[0];
    for (j = 1; j < nch; j++) {
      c[j] = 0.5 * cos(delta * j);
      c[nc - j] = 0.5 * sin(delta * j);
    }
  }
}

/* -------- child routines -------- */

void bitrv2(signed n, signed *ip, double *a) {
  signed j, j1, k, k1, l, m, m2;
  double xr, xi, yr, yi;

  ip[0] = 0;
  l = n;
  m = 1;
  while ((m << 3) < l) {
    l >>= 1;
    for (j = 0; j < m; j++) {
      ip[m + j] = ip[j] + l;
    }
    m <<= 1;
  }
  m2 = 2 * m;
  if ((m << 3) == l) {
    for (k = 0; k < m; k++) {
      for (j = 0; j < k; j++) {
        j1 = 2 * j + ip[k];
        k1 = 2 * k + ip[j];
        xr = a[j1];
        xi = a[j1 + 1];
        yr = a[k1];
        yi = a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
        j1 += m2;
        k1 += 2 * m2;
        xr = a[j1];
        xi = a[j1 + 1];
        yr = a[k1];
        yi = a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
        j1 += m2;
        k1 -= m2;
        xr = a[j1];
        xi = a[j1 + 1];
        yr = a[k1];
        yi = a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
        j1 += m2;
        k1 += 2 * m2;
        xr = a[j1];
        xi = a[j1 + 1];
        yr = a[k1];
        yi = a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
      }
      j1 = 2 * k + m2 + ip[k];
      k1 = j1 + m2;
      xr = a[j1];
      xi = a[j1 + 1];
      yr = a[k1];
      yi = a[k1 + 1];
      a[j1] = yr;
      a[j1 + 1] = yi;
      a[k1] = xr;
      a[k1 + 1] = xi;
    }
  } else {
    for (k = 1; k < m; k++) {
      for (j = 0; j < k; j++) {
        j1 = 2 * j + ip[k];
        k1 = 2 * k + ip[j];
        xr = a[j1];
        xi = a[j1 + 1];
        yr = a[k1];
        yi = a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
        j1 += m2;
        k1 += m2;
        xr = a[j1];
        xi = a[j1 + 1];
        yr = a[k1];
        yi = a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
      }
    }
  }
}

void bitrv2conj(signed n, signed *ip, double *a) {
  signed j, j1, k, k1, l, m, m2;
  double xr, xi, yr, yi;

  ip[0] = 0;
  l = n;
  m = 1;
  while ((m << 3) < l) {
    l >>= 1;
    for (j = 0; j < m; j++) {
      ip[m + j] = ip[j] + l;
    }
    m <<= 1;
  }
  m2 = 2 * m;
  if ((m << 3) == l) {
    for (k = 0; k < m; k++) {
      for (j = 0; j < k; j++) {
        j1 = 2 * j + ip[k];
        k1 = 2 * k + ip[j];
        xr = a[j1];
        xi = -a[j1 + 1];
        yr = a[k1];
        yi = -a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
        j1 += m2;
        k1 += 2 * m2;
        xr = a[j1];
        xi = -a[j1 + 1];
        yr = a[k1];
        yi = -a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
        j1 += m2;
        k1 -= m2;
        xr = a[j1];
        xi = -a[j1 + 1];
        yr = a[k1];
        yi = -a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
        j1 += m2;
        k1 += 2 * m2;
        xr = a[j1];
        xi = -a[j1 + 1];
        yr = a[k1];
        yi = -a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
      }
      k1 = 2 * k + ip[k];
      a[k1 + 1] = -a[k1 + 1];
      j1 = k1 + m2;
      k1 = j1 + m2;
      xr = a[j1];
      xi = -a[j1 + 1];
      yr = a[k1];
      yi = -a[k1 + 1];
      a[j1] = yr;
      a[j1 + 1] = yi;
      a[k1] = xr;
      a[k1 + 1] = xi;
      k1 += m2;
      a[k1 + 1] = -a[k1 + 1];
    }
  } else {
    a[1] = -a[1];
    a[m2 + 1] = -a[m2 + 1];
    for (k = 1; k < m; k++) {
      for (j = 0; j < k; j++) {
        j1 = 2 * j + ip[k];
        k1 = 2 * k + ip[j];
        xr = a[j1];
        xi = -a[j1 + 1];
        yr = a[k1];
        yi = -a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
        j1 += m2;
        k1 += m2;
        xr = a[j1];
        xi = -a[j1 + 1];
        yr = a[k1];
        yi = -a[k1 + 1];
        a[j1] = yr;
        a[j1 + 1] = yi;
        a[k1] = xr;
        a[k1 + 1] = xi;
      }
      k1 = 2 * k + ip[k];
      a[k1 + 1] = -a[k1 + 1];
      a[k1 + m2 + 1] = -a[k1 + m2 + 1];
    }
  }
}

void cftfsub(signed n, double *a, double *w) {
  void cft1st(signed n, double *a, double *w);
  void cftmdl(signed n, signed l, double *a, double *w);
  signed j, j1, j2, j3, l;
  double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;

  l = 2;
  if (n > 8) {
    cft1st(n, a, w);
    l = 8;
    while ((l << 2) < n) {
      cftmdl(n, l, a, w);
      l <<= 2;
    }
  }
  if ((l << 2) == n) {
    for (j = 0; j < l; j += 2) {
      j1 = j + l;
      j2 = j1 + l;
      j3 = j2 + l;
      x0r = a[j] + a[j1];
      x0i = a[j + 1] + a[j1 + 1];
      x1r = a[j] - a[j1];
      x1i = a[j + 1] - a[j1 + 1];
      x2r = a[j2] + a[j3];
      x2i = a[j2 + 1] + a[j3 + 1];
      x3r = a[j2] - a[j3];
      x3i = a[j2 + 1] - a[j3 + 1];
      a[j] = x0r + x2r;
      a[j + 1] = x0i + x2i;
      a[j2] = x0r - x2r;
      a[j2 + 1] = x0i - x2i;
      a[j1] = x1r - x3i;
      a[j1 + 1] = x1i + x3r;
      a[j3] = x1r + x3i;
      a[j3 + 1] = x1i - x3r;
    }
  } else {
    for (j = 0; j < l; j += 2) {
      j1 = j + l;
      x0r = a[j] - a[j1];
      x0i = a[j + 1] - a[j1 + 1];
      a[j] += a[j1];
      a[j + 1] += a[j1 + 1];
      a[j1] = x0r;
      a[j1 + 1] = x0i;
    }
  }
}

void cftbsub(signed n, double *a, double *w) {
  void cft1st(signed n, double *a, double *w);
  void cftmdl(signed n, signed l, double *a, double *w);
  signed j, j1, j2, j3, l;
  double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;

  l = 2;
  if (n > 8) {
    cft1st(n, a, w);
    l = 8;
    while ((l << 2) < n) {
      cftmdl(n, l, a, w);
      l <<= 2;
    }
  }
  if ((l << 2) == n) {
    for (j = 0; j < l; j += 2) {
      j1 = j + l;
      j2 = j1 + l;
      j3 = j2 + l;
      x0r = a[j] + a[j1];
      x0i = -a[j + 1] - a[j1 + 1];
      x1r = a[j] - a[j1];
      x1i = -a[j + 1] + a[j1 + 1];
      x2r = a[j2] + a[j3];
      x2i = a[j2 + 1] + a[j3 + 1];
      x3r = a[j2] - a[j3];
      x3i = a[j2 + 1] - a[j3 + 1];
      a[j] = x0r + x2r;
      a[j + 1] = x0i - x2i;
      a[j2] = x0r - x2r;
      a[j2 + 1] = x0i + x2i;
      a[j1] = x1r - x3i;
      a[j1 + 1] = x1i - x3r;
      a[j3] = x1r + x3i;
      a[j3 + 1] = x1i + x3r;
    }
  } else {
    for (j = 0; j < l; j += 2) {
      j1 = j + l;
      x0r = a[j] - a[j1];
      x0i = -a[j + 1] + a[j1 + 1];
      a[j] += a[j1];
      a[j + 1] = -a[j + 1] - a[j1 + 1];
      a[j1] = x0r;
      a[j1 + 1] = x0i;
    }
  }
}

void cft1st(signed n, double *a, double *w) {
  signed j, k1, k2;
  double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
  double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;

  x0r = a[0] + a[2];
  x0i = a[1] + a[3];
  x1r = a[0] - a[2];
  x1i = a[1] - a[3];
  x2r = a[4] + a[6];
  x2i = a[5] + a[7];
  x3r = a[4] - a[6];
  x3i = a[5] - a[7];
  a[0] = x0r + x2r;
  a[1] = x0i + x2i;
  a[4] = x0r - x2r;
  a[5] = x0i - x2i;
  a[2] = x1r - x3i;
  a[3] = x1i + x3r;
  a[6] = x1r + x3i;
  a[7] = x1i - x3r;
  wk1r = w[2];
  x0r = a[8] + a[10];
  x0i = a[9] + a[11];
  x1r = a[8] - a[10];
  x1i = a[9] - a[11];
  x2r = a[12] + a[14];
  x2i = a[13] + a[15];
  x3r = a[12] - a[14];
  x3i = a[13] - a[15];
  a[8] = x0r + x2r;
  a[9] = x0i + x2i;
  a[12] = x2i - x0i;
  a[13] = x0r - x2r;
  x0r = x1r - x3i;
  x0i = x1i + x3r;
  a[10] = wk1r * (x0r - x0i);
  a[11] = wk1r * (x0r + x0i);
  x0r = x3i + x1r;
  x0i = x3r - x1i;
  a[14] = wk1r * (x0i - x0r);
  a[15] = wk1r * (x0i + x0r);
  k1 = 0;
  for (j = 16; j < n; j += 16) {
    k1 += 2;
    k2 = 2 * k1;
    wk2r = w[k1];
    wk2i = w[k1 + 1];
    wk1r = w[k2];
    wk1i = w[k2 + 1];
    wk3r = wk1r - 2 * wk2i * wk1i;
    wk3i = 2 * wk2i * wk1r - wk1i;
    x0r = a[j] + a[j + 2];
    x0i = a[j + 1] + a[j + 3];
    x1r = a[j] - a[j + 2];
    x1i = a[j + 1] - a[j + 3];
    x2r = a[j + 4] + a[j + 6];
    x2i = a[j + 5] + a[j + 7];
    x3r = a[j + 4] - a[j + 6];
    x3i = a[j + 5] - a[j + 7];
    a[j] = x0r + x2r;
    a[j + 1] = x0i + x2i;
    x0r -= x2r;
    x0i -= x2i;
    a[j + 4] = wk2r * x0r - wk2i * x0i;
    a[j + 5] = wk2r * x0i + wk2i * x0r;
    x0r = x1r - x3i;
    x0i = x1i + x3r;
    a[j + 2] = wk1r * x0r - wk1i * x0i;
    a[j + 3] = wk1r * x0i + wk1i * x0r;
    x0r = x1r + x3i;
    x0i = x1i - x3r;
    a[j + 6] = wk3r * x0r - wk3i * x0i;
    a[j + 7] = wk3r * x0i + wk3i * x0r;
    wk1r = w[k2 + 2];
    wk1i = w[k2 + 3];
    wk3r = wk1r - 2 * wk2r * wk1i;
    wk3i = 2 * wk2r * wk1r - wk1i;
    x0r = a[j + 8] + a[j + 10];
    x0i = a[j + 9] + a[j + 11];
    x1r = a[j + 8] - a[j + 10];
    x1i = a[j + 9] - a[j + 11];
    x2r = a[j + 12] + a[j + 14];
    x2i = a[j + 13] + a[j + 15];
    x3r = a[j + 12] - a[j + 14];
    x3i = a[j + 13] - a[j + 15];
    a[j + 8] = x0r + x2r;
    a[j + 9] = x0i + x2i;
    x0r -= x2r;
    x0i -= x2i;
    a[j + 12] = -wk2i * x0r - wk2r * x0i;
    a[j + 13] = -wk2i * x0i + wk2r * x0r;
    x0r = x1r - x3i;
    x0i = x1i + x3r;
    a[j + 10] = wk1r * x0r - wk1i * x0i;
    a[j + 11] = wk1r * x0i + wk1i * x0r;
    x0r = x1r + x3i;
    x0i = x1i - x3r;
    a[j + 14] = wk3r * x0r - wk3i * x0i;
    a[j + 15] = wk3r * x0i + wk3i * x0r;
  }
}

void cftmdl(signed n, signed l, double *a, double *w) {
  signed j, j1, j2, j3, k, k1, k2, m, m2;
  double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
  double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;

  m = l << 2;
  for (j = 0; j < l; j += 2) {
    j1 = j + l;
    j2 = j1 + l;
    j3 = j2 + l;
    x0r = a[j] + a[j1];
    x0i = a[j + 1] + a[j1 + 1];
    x1r = a[j] - a[j1];
    x1i = a[j + 1] - a[j1 + 1];
    x2r = a[j2] + a[j3];
    x2i = a[j2 + 1] + a[j3 + 1];
    x3r = a[j2] - a[j3];
    x3i = a[j2 + 1] - a[j3 + 1];
    a[j] = x0r + x2r;
    a[j + 1] = x0i + x2i;
    a[j2] = x0r - x2r;
    a[j2 + 1] = x0i - x2i;
    a[j1] = x1r - x3i;
    a[j1 + 1] = x1i + x3r;
    a[j3] = x1r + x3i;
    a[j3 + 1] = x1i - x3r;
  }
  wk1r = w[2];
  for (j = m; j < l + m; j += 2) {
    j1 = j + l;
    j2 = j1 + l;
    j3 = j2 + l;
    x0r = a[j] + a[j1];
    x0i = a[j + 1] + a[j1 + 1];
    x1r = a[j] - a[j1];
    x1i = a[j + 1] - a[j1 + 1];
    x2r = a[j2] + a[j3];
    x2i = a[j2 + 1] + a[j3 + 1];
    x3r = a[j2] - a[j3];
    x3i = a[j2 + 1] - a[j3 + 1];
    a[j] = x0r + x2r;
    a[j + 1] = x0i + x2i;
    a[j2] = x2i - x0i;
    a[j2 + 1] = x0r - x2r;
    x0r = x1r - x3i;
    x0i = x1i + x3r;
    a[j1] = wk1r * (x0r - x0i);
    a[j1 + 1] = wk1r * (x0r + x0i);
    x0r = x3i + x1r;
    x0i = x3r - x1i;
    a[j3] = wk1r * (x0i - x0r);
    a[j3 + 1] = wk1r * (x0i + x0r);
  }
  k1 = 0;
  m2 = 2 * m;
  for (k = m2; k < n; k += m2) {
    k1 += 2;
    k2 = 2 * k1;
    wk2r = w[k1];
    wk2i = w[k1 + 1];
    wk1r = w[k2];
    wk1i = w[k2 + 1];
    wk3r = wk1r - 2 * wk2i * wk1i;
    wk3i = 2 * wk2i * wk1r - wk1i;
    for (j = k; j < l + k; j += 2) {
      j1 = j + l;
      j2 = j1 + l;
      j3 = j2 + l;
      x0r = a[j] + a[j1];
      x0i = a[j + 1] + a[j1 + 1];
      x1r = a[j] - a[j1];
      x1i = a[j + 1] - a[j1 + 1];
      x2r = a[j2] + a[j3];
      x2i = a[j2 + 1] + a[j3 + 1];
      x3r = a[j2] - a[j3];
      x3i = a[j2 + 1] - a[j3 + 1];
      a[j] = x0r + x2r;
      a[j + 1] = x0i + x2i;
      x0r -= x2r;
      x0i -= x2i;
      a[j2] = wk2r * x0r - wk2i * x0i;
      a[j2 + 1] = wk2r * x0i + wk2i * x0r;
      x0r = x1r - x3i;
      x0i = x1i + x3r;
      a[j1] = wk1r * x0r - wk1i * x0i;
      a[j1 + 1] = wk1r * x0i + wk1i * x0r;
      x0r = x1r + x3i;
      x0i = x1i - x3r;
      a[j3] = wk3r * x0r - wk3i * x0i;
      a[j3 + 1] = wk3r * x0i + wk3i * x0r;
    }
    wk1r = w[k2 + 2];
    wk1i = w[k2 + 3];
    wk3r = wk1r - 2 * wk2r * wk1i;
    wk3i = 2 * wk2r * wk1r - wk1i;
    for (j = k + m; j < l + (k + m); j += 2) {
      j1 = j + l;
      j2 = j1 + l;
      j3 = j2 + l;
      x0r = a[j] + a[j1];
      x0i = a[j + 1] + a[j1 + 1];
      x1r = a[j] - a[j1];
      x1i = a[j + 1] - a[j1 + 1];
      x2r = a[j2] + a[j3];
      x2i = a[j2 + 1] + a[j3 + 1];
      x3r = a[j2] - a[j3];
      x3i = a[j2 + 1] - a[j3 + 1];
      a[j] = x0r + x2r;
      a[j + 1] = x0i + x2i;
      x0r -= x2r;
      x0i -= x2i;
      a[j2] = -wk2i * x0r - wk2r * x0i;
      a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
      x0r = x1r - x3i;
      x0i = x1i + x3r;
      a[j1] = wk1r * x0r - wk1i * x0i;
      a[j1 + 1] = wk1r * x0i + wk1i * x0r;
      x0r = x1r + x3i;
      x0i = x1i - x3r;
      a[j3] = wk3r * x0r - wk3i * x0i;
      a[j3 + 1] = wk3r * x0i + wk3i * x0r;
    }
  }
}

void rftfsub(signed n, double *a, signed nc, double *c) {
  signed j, k, kk, ks, m;
  double wkr, wki, xr, xi, yr, yi;

  m = n >> 1;
  ks = 2 * nc / m;
  kk = 0;
  for (j = 2; j < m; j += 2) {
    k = n - j;
    kk += ks;
    wkr = 0.5 - c[nc - kk];
    wki = c[kk];
    xr = a[j] - a[k];
    xi = a[j + 1] + a[k + 1];
    yr = wkr * xr - wki * xi;
    yi = wkr * xi + wki * xr;
    a[j] -= yr;
    a[j + 1] -= yi;
    a[k] += yr;
    a[k + 1] -= yi;
  }
}

void rftbsub(signed n, double *a, signed nc, double *c) {
  signed j, k, kk, ks, m;
  double wkr, wki, xr, xi, yr, yi;

  a[1] = -a[1];
  m = n >> 1;
  ks = 2 * nc / m;
  kk = 0;
  for (j = 2; j < m; j += 2) {
    k = n - j;
    kk += ks;
    wkr = 0.5 - c[nc - kk];
    wki = c[kk];
    xr = a[j] - a[k];
    xi = a[j + 1] + a[k + 1];
    yr = wkr * xr + wki * xi;
    yi = wkr * xi - wki * xr;
    a[j] -= yr;
    a[j + 1] = yi - a[j + 1];
    a[k] += yr;
    a[k + 1] = yi - a[k + 1];
  }
  a[m + 1] = -a[m + 1];
}

void dctsub(signed n, double *a, signed nc, double *c) {
  signed j, k, kk, ks, m;
  double wkr, wki, xr;

  m = n >> 1;
  ks = nc / n;
  kk = 0;
  for (j = 1; j < m; j++) {
    k = n - j;
    kk += ks;
    wkr = c[kk] - c[nc - kk];
    wki = c[kk] + c[nc - kk];
    xr = wki * a[j] - wkr * a[k];
    a[j] = wkr * a[j] + wki * a[k];
    a[k] = xr;
  }
  a[m] *= c[0];
}

void dstsub(signed n, double *a, signed nc, double *c) {
  signed j, k, kk, ks, m;
  double wkr, wki, xr;

  m = n >> 1;
  ks = nc / n;
  kk = 0;
  for (j = 1; j < m; j++) {
    k = n - j;
    kk += ks;
    wkr = c[kk] - c[nc - kk];
    wki = c[kk] + c[nc - kk];
    xr = wki * a[k] - wkr * a[j];
    a[k] = wkr * a[k] + wki * a[j];
    a[j] = xr;
  }
  a[m] *= c[0];
}

//----//

signed a[MAXN];
signed b[MAXN];
double ab[MAXSIZE];

int binpow(int a, int b, int mo) {
  int ans = 1;
  do {
    if (b & 1) ans = 1ll * ans * a % mo;
    a = 1ll * a * a % mo;
  } while (b >>= 1);
  return ans;
}

const int FACTSIZE = 10100;
int fact[FACTSIZE];
int revfact[FACTSIZE];

int ncr(int n, int r) {
  if (n < r) return 0;
  int rs = fact[n] * revfact[r];
  if (rs >= MOD) rs %= MOD;
  rs *= revfact[n - r];
  if (rs >= MOD) rs %= MOD;
  return rs;
}
void init() {
  fact[0] = 1, revfact[0] = 1;
  for (int i = 1; i < FACTSIZE; ++i) {
    fact[i] = fact[i - 1] * i % MOD;
    revfact[i] = binpow(fact[i], MOD - 2, MOD);
  }
}
signed ip[10000];
double w[(MAXSIZE / 4)];

void convolve2(double *a, signed n) {
  signed s = n + n - 1;  // 畳み込み結果のサイズ
  signed t = 1;          // FFT に使う配列のサイズ（2 の累乗）
  while (t < s) t *= 2;
  ip[0] = 0;
  cdft(2 * t, 1, a, ip, w);
  for (signed i = 0; i < t; i++) {
    signed addr = 2 * i;
    double tmpa = a[addr];
    a[addr] = a[addr] * a[addr] - a[addr + 1] * a[addr + 1];
    a[addr + 1] = 2 * tmpa * a[addr + 1];
  }
  cdft(2 * t, -1, a, ip, w);
  for (signed j = 0; j < t; j++) {
    a[j * 2] *= 1.0 / t;
    a[j] = a[j * 2];
  }
}

signed N;
signed main() {
  init();
  cin.tie(0), ios::sync_with_stdio(false);
  cout << fixed << setprecision(10);

  int N;
  cin >> N;

  signed maxV = 0;
  REP(i, N) {
    cin >> a[i] >> b[i];
    chmax(maxV, a[i]);
    chmax(maxV, b[i]);
  }
  maxV *= 2;
  maxV += 1;

  int maxAB = 0;
  REP(i, N) {
    int val = a[i] * maxV + b[i];
    ab[2 * val]++;
    chmax(maxAB, val);
  }
  maxAB++;
  convolve2(ab, maxAB);
  int ans = 0;
  int maxV2 = maxV * maxV;
  DUMP(maxAB, maxV2);
  REP(i, maxV2) {
    int i0 = i / maxV;
    int i1 = i % maxV;
    int val = round(ab[i]) % MOD;
    ans += val * ncr(i0 + i1, i0);
    ans %= MOD;
  }
  int cnt = 0;
  REP(i, N) cnt += ncr(2 * (a[i] + b[i]), 2 * a[i]), cnt %= MOD;
  cnt = MOD - cnt;
  ans += cnt;
  ans %= MOD;
  ans = ans * revfact[2];
  ans %= MOD;
  print(ans);
}

Compiler output:
a.cc: In function 'int main()':
a.cc:1511:28: error: invalid operands of types 'double' and 'const long long int' to binary 'operator%'
 1511 |     int val = round(ab[i]) % MOD;
      |               ~~~~~~~~~~~~ ^ ~~~
      |                    |         |
      |                    double    const long long int