#line 2 "library/template/template.hpp" #include using namespace std; #line 2 "library/template/util.hpp" using uint = unsigned int; using ll = long long int; using ull = unsigned long long; using i128 = __int128_t; using u128 = __uint128_t; template S SUM(const vector &a) { return accumulate(ALL(a), S(0)); } template inline bool chmin(T &a, T b) { if (a > b) { a = b; return true; } return false; } template inline bool chmax(T &a, T b) { if (a < b) { a = b; return true; } return false; } template int popcnt(T x) { return __builtin_popcountll(x); } template int topbit(T x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); } template int lowbit(T x) { return (x == 0 ? -1 : __builtin_ctzll(x)); } #line 6 "library/template/template.hpp" #line 2 "library/template/macro.hpp" #define rep(i, a, b) for (int i = (a); i < (int)(b); i++) #define rrep(i, a, b) for (int i = (int)(b) - 1; i >= (a); i--) #define ALL(v) (v).begin(), (v).end() #define UNIQUE(v) sort(ALL(v)), (v).erase(unique(ALL(v)), (v).end()) #define SZ(v) (int)v.size() #define MIN(v) *min_element(ALL(v)) #define MAX(v) *max_element(ALL(v)) #define LB(v, x) int(lower_bound(ALL(v), (x)) - (v).begin()) #define UB(v, x) int(upper_bound(ALL(v), (x)) - (v).begin()) #define YN(b) cout << ((b) ? "YES" : "NO") << "\n"; #define Yn(b) cout << ((b) ? "Yes" : "No") << "\n"; #define yn(b) cout << ((b) ? "yes" : "no") << "\n"; #line 8 "library/template/template.hpp" #line 2 "library/template/inout.hpp" struct Fast { Fast() { cin.tie(nullptr); ios_base::sync_with_stdio(false); cout << fixed << setprecision(15); } } fast; template istream &operator>>(istream &is, pair &p) { return is >> p.first >> p.second; } template ostream &operator<<(ostream &os, const pair &p) { return os << p.first << " " << p.second; } template istream &operator>>(istream &is, vector &a) { for (auto &v : a) is >> v; return is; } template ostream &operator<<(ostream &os, const vector &a) { for (auto it = a.begin(); it != a.end();) { os << *it; if (++it != a.end()) os << " "; } return os; } template ostream &operator<<(ostream &os, const set &st) { os << "{"; for (auto it = st.begin(); it != st.end();) { os << *it; if (++it != st.end()) os << ","; } os << "}"; return os; } template ostream &operator<<(ostream &os, const map &mp) { os << "{"; for (auto it = mp.begin(); it != mp.end();) { os << it->first << ":" << it->second; if (++it != mp.end()) os << ","; } os << "}"; return os; } void in() {} template void in(T &t, U &...u) { cin >> t; in(u...); } void out() { cout << "\n"; } template void out(const T &t, const U &...u) { cout << t; if (sizeof...(u)) cout << sep; out(u...); } #line 10 "library/template/template.hpp" #line 2 "library/template/debug.hpp" #ifdef LOCAL #define debug 1 #define show(...) _show(0, #__VA_ARGS__, __VA_ARGS__) #else #define debug 0 #define show(...) true #endif template void _show(int i, T name) { cerr << '\n'; } template void _show(int i, const T1 &a, const T2 &b, const T3 &...c) { for (; a[i] != ',' && a[i] != '\0'; i++) cerr << a[i]; cerr << ":" << b << " "; _show(i + 1, a, c...); } #line 2 "library/math/util.hpp" namespace Math { template T safe_mod(T a, T b) { assert(b != 0); if (b < 0) a = -a, b = -b; a %= b; return a >= 0 ? a : a + b; } template T floor(T a, T b) { assert(b != 0); if (b < 0) a = -a, b = -b; return a >= 0 ? a / b : (a + 1) / b - 1; } template T ceil(T a, T b) { assert(b != 0); if (b < 0) a = -a, b = -b; return a > 0 ? (a - 1) / b + 1 : a / b; } long long isqrt(long long n) { if (n <= 0) return 0; long long x = sqrt(n); while ((x + 1) * (x + 1) <= n) x++; while (x * x > n) x--; return x; } // return g=gcd(a,b) // a*x+b*y=g // - b!=0 -> 0<=x<|b|/g // - b=0 -> ax=g template T ext_gcd(T a, T b, T& x, T& y) { T a0 = a, b0 = b; bool sgn_a = a < 0, sgn_b = b < 0; if (sgn_a) a = -a; if (sgn_b) b = -b; if (b == 0) { x = sgn_a ? -1 : 1; y = 0; return a; } T x00 = 1, x01 = 0, x10 = 0, x11 = 1; while (b != 0) { T q = a / b, r = a - b * q; x00 -= q * x01; x10 -= q * x11; swap(x00, x01); swap(x10, x11); a = b, b = r; } x = x00, y = x10; if (sgn_a) x = -x; if (sgn_b) y = -y; if (b0 != 0) { a0 /= a, b0 /= a; if (b0 < 0) a0 = -a0, b0 = -b0; T q = x >= 0 ? x / b0 : (x + 1) / b0 - 1; x -= b0 * q; y += a0 * q; } return a; } template T inv_mod(T x, T m) { x %= m; if (x < 0) x += m; T a = m, b = x; T y0 = 0, y1 = 1; while (b > 0) { T q = a / b; swap(a -= q * b, b); swap(y0 -= q * y1, y1); } if (y0 < 0) y0 += m / a; return y0; } template T pow_mod(T x, T n, T m) { x = (x % m + m) % m; T y = 1; while (n) { if (n & 1) y = y * x % m; x = x * x % m; n >>= 1; } return y; } constexpr long long pow_mod_constexpr(long long x, long long n, int m) { if (m == 1) return 0; unsigned int _m = (unsigned int)(m); unsigned long long r = 1; unsigned long long y = x % m; if (y >= m) y += m; while (n) { if (n & 1) r = (r * y) % _m; y = (y * y) % _m; n >>= 1; } return r; } }; // namespace Math #line 2 "library/math/floor-monoid-product.hpp" template T FloorMonoidProduct(U n, U m, U a, U b, T x, T y, function pow = nullptr) { if (!pow) { pow = [](T t, U n) { T p = e(); while (n) { if (n & 1) p = op(p, t); t = op(t, t); n >>= 1; } return p; }; } assert(m != 0); T pl = e(), pr = e(); while (true) { if (a >= m) { U q = a / m; x = op(x, pow(y, q)); a -= m * q; } if (b >= m) { U q = b / m; pl = op(pl, pow(y, q)); b -= m * q; } U c = a * n + b; if (c < m) { pl = op(pl, pow(x, n)); break; } pr = op(op(y, pow(x, c % m / a)), pr); n = c / m - 1; b = m + a - b - 1; swap(a, m); swap(x, y); } return op(pl, pr); } /** * @brief モノイド版 Floor Sum * @docs docs/math/floor-monoid-product.md */ #line 4 "library/math/polynomial-floor-sum.hpp" template struct PolynomialFloorSumMonoid { static constexpr int D = max(D1, D2); using M = PolynomialFloorSumMonoid; using P = array, D1 + 1>; T x, y; P f; static M e() { return {0, 0, P{}}; } static M op(M a, M b) { static T binom[D + 1][D + 1]; if (binom[0][0] != T(1)) { binom[0][0] = T(1); for (int i = 0; i < D; i++) for (int j = 0; j <= i; j++) { binom[i + 1][j] += binom[i][j]; binom[i + 1][j + 1] += binom[i][j]; } } T pow_x[D1 + 1], pow_y[D2 + 1]; pow_x[0] = 1, pow_y[0] = 1; for (int i = 0; i < D1; i++) pow_x[i + 1] = pow_x[i] * a.x; for (int i = 0; i < D2; i++) pow_y[i + 1] = pow_y[i] * a.y; for (int i = 0; i <= D1; i++) for (int j = 0; j <= D2; j++) { T v = b.f[i][j]; for (int k = i + 1; k <= D1; k++) b.f[k][j] += v * pow_x[k - i] * binom[k][i]; } for (int i = 0; i <= D1; i++) for (int j = 0; j <= D2; j++) { T v = b.f[i][j]; for (int k = j; k <= D2; k++) a.f[i][k] += v * pow_y[k - j] * binom[k][j]; } return {a.x + b.x, a.y + b.y, a.f}; } static M elm_x() { M t = e(); t.x = 1, t.f[0][0] = 1; return t; } static M elm_y() { M t = e(); t.y = 1; return t; } }; // enumerate sum{k=l}^{r-1}k^i*floor((a*k+b)/m))^j template array, D1 + 1> PolynomialFloorSum(I l, I r, I m, I a, I b) { assert(l <= r && m != 0); using U = conditional_t, __uint128_t, uint64_t>; using M = PolynomialFloorSumMonoid; using P = array, D1 + 1>; if (m < 0) m = -m, a = -a, b = -b; if (a < 0) { P c = PolynomialFloorSum(-r + 1, -l + 1, m, -a, b); for (int i = 1; i <= D1; i += 2) for (int j = 0; j <= D2; j++) c[i][j] = -c[i][j]; return c; } b += a * l; I q = Math::floor(b, m); b -= q * m; M t = M::e(); t.x = l, t.y = q; M z = FloorMonoidProduct(r - l, m, a, b, M::elm_x(), M::elm_y()); return M::op(t, z).f; } // P(x,y): 2 variables polynomial, deg_x(P)<=D1, deg_y(Q)<=D2 // find sum{k=0}^{n-1}P(k,floor((a*k+b)/m)) template T PolynomialFloorSum(array, D1 + 1> poly, I l, I r, I m, I a, I b) { using M = PolynomialFloorSumMonoid; auto c = PolynomialFloorSum(l, r, m, a, b); T res = 0; for (int i = 0; i <= D1; i++) for (int j = 0; j <= D2; j++) res += poly[i][j] * c[i][j]; return res; } /** * @brief 多項式版 floor sum * @docs docs/math/polynomial-floor-sum.md */ #line 3 "main.cpp" void solve() { int n, m, x, y; in(n, m, x, y); auto c1 = PolynomialFloorSum(0, n, m, x, y); auto c2 = PolynomialFloorSum(0, n, m, x, 0); auto ans = 2 * c1[1][1] - (n - 1) * c1[0][1] + c2[1][1] - n * c2[0][1]; out(ans); } int main() { int t = 1; cin >> t; while (t--) solve(); }