Mujin Programming Challenge 2018 G - 移動

問題ページ
数式を用いて問題を整理します．ベクトル $\vec{v_1} = (x_1, y_1), \vec{v_2} = (x_2, y_2), \vec{v_3} = (x_3, y_3)$ に対して， $a\vec{v_1} + b\vec{v_2} + c\vec{v_3} \quad (a+b+c \leq K)$ となる $(a,b,c)$ は何種類ですか？という問題になりました．

もし， $(a,b,c)$ が異なる場合は全て違う点に移動するならば， $a+b+c \leq K$ となるような $(a,b,c)$ が何種類か数える問題となり，これはcombinationを用いた計算で求めることができます．異なる $(a,b,c)$ であっても同じ点に移動することがあるのが難しいです．

移動する点が被るのはどのような状況か考えます． $p\vec{v_1} + q\vec{v_2} + r\vec{v_3} = \vec{0}$ となるような $p, q, r$ の組を求めます．この $p, q, r$ は定数倍を除いて一意に定まります．外積が0であるような2本のベクトルは一次独立であり， $\vec{v_1}, \vec{v_2}$ は基底となります．この2次元空間で $\vec{v_3}$ を表す方法は一意に定まることから， $p, q, r$ は定数倍を除いて一意に定まる，つまり比が常に一定であることがわかります．

線形方程式 $px_1 + qx_2 + rx_3 = 0, py_1 + qy_2 + ry_3 = 0$ の解を求めます． $r = 1$ としたとき， $p,q$ について連立方程式を解くと， $q = (x_1y_3-x_3y_1)/(x_2y_1-x_1y_2), p = (x_3y_2-x_2y_3)/(x_2y_1-x_1y_2)$ となります． $p, q, r$ は整数である場合を考えるので， $x_2y_1-x_1y_2$ を掛けることで分母をはらっておきます．すると， $p = x_3y_2-x_2y_3, q = x_1y_3-x_3y_1, r = x_2y_1-x_1y_2$ となります．また， $p, q, r$ のGCDでそれぞれ割っておくことで，もっとも簡単な形になるようにしておきます．

外積が0でない制約から $p,q,r \neq 0$ です．また定数倍しても問題がないことから， $p,q \gt 0, r \neq 0$ としても一般性を失いません．

$r \gt 0$
$a \geq p, b \geq q, c \geq r$ のとき $a\vec{v_1} + b\vec{v_2} + c\vec{v_3} = (a-p)\vec{v_1} + (b-q)\vec{v_2} + (c-r)\vec{v_3}$ です． $a+b+c \lt k$ and $(a \lt p$ or $b \lt q$ or $c \lt r)$ となる $(a, b, c)$ は何通り？という問題が解ければよいです．
- $p+q \lt -r$
  $c \geq -r$ のとき $a\vec{v_1} + b\vec{v_2} + c\vec{v_3} = (a+p)\vec{v_1} + (b+q)\vec{v_2} + (c+r)\vec{v_3}$ です． $a+b+c \lt k$ and $c \lt -r$ となる $(a, b, c)$ は何通り？という問題が解ければよいです．
- $p+q \geq -r$
  $a \geq p, b \geq q$ のとき $a\vec{v_1} + b\vec{v_2} + c\vec{v_3} = (a-p)\vec{v_1} + (b-q)\vec{v_2} + (c-r)\vec{v_3}$ です． $a+b+c \lt k$ and $(a \lt p$ or $b \lt q)$ となる $(a, b, c)$ は何通り？という問題が解ければよいです．

これらの数え上げは重複組み合わせを用いることで求められます．ユークリッドの互除法でgcdを求める部分がボトルネックで，各クエリについて $O(\log X)$ 程度で求められます．

#include <bits/stdc++.h>  
using namespace std;  
using ll = long long;  
using PII = pair<ll, ll>;  
#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)  
#define REP(i, n) FOR(i, 0, n)  
#define ALL(x) x.begin(), x.end()  
template<typename T> void chmin(T &a, const T &b) { a = min(a, b); }  
template<typename T> void chmax(T &a, const T &b) { a = max(a, b); }  
struct FastIO {FastIO() { cin.tie(0); ios::sync_with_stdio(0); }}fastiofastio;  
#ifdef DEBUG_   
#include "../program_contest_library/memo/dump.hpp"  
#else  
#define dump(...)  
#endif  
const ll INF = 1LL<<60;  
  
template<ll MOD>  
struct modint {  
    ll x;  
    modint(): x(0) {}  
    modint(ll y) : x(y>=0 ? y%MOD : y%MOD+MOD) {}  
    static constexpr ll mod() { return MOD; }  
    // e乗  
    modint pow(ll e) {  
        ll a = 1, p = x;  
        while(e > 0) {  
            if(e%2 == 0) {p = (p*p) % MOD; e /= 2;}  
            else {a = (a*p) % MOD; e--;}  
        }  
        return modint(a);  
    }  
    modint inv() const {  
        ll a=x, b=MOD, u=1, y=1, v=0, z=0;  
        while(a) {  
            ll q = b/a;  
            swap(z -= q*u, u);  
            swap(y -= q*v, v);  
            swap(b -= q*a, a);  
        }  
        return z;  
    }  
    // Comparators  
    bool operator <(modint b) { return x < b.x; }  
    bool operator >(modint b) { return x > b.x; }  
    bool operator<=(modint b) { return x <= b.x; }  
    bool operator>=(modint b) { return x >= b.x; }  
    bool operator!=(modint b) { return x != b.x; }  
    bool operator==(modint b) { return x == b.x; }  
    // Basic Operations  
    modint operator+(modint r) const { return modint(*this) += r; }  
    modint operator-(modint r) const { return modint(*this) -= r; }  
    modint operator*(modint r) const { return modint(*this) *= r; }  
    modint operator/(modint r) const { return modint(*this) /= r; }  
    modint &operator+=(modint r) {  
        if((x += r.x) >= MOD) x -= MOD;  
        return *this;  
    }  
    modint &operator-=(modint r) {  
        if((x -= r.x) < 0) x += MOD;  
        return *this;  
    }  
    modint &operator*=(modint r) {  
    #if !defined(_WIN32) || defined(_WIN64)  
        x = x * r.x % MOD; return *this;  
    #endif  
        unsigned long long y = x * r.x;  
        unsigned xh = (unsigned) (y >> 32), xl = (unsigned) y, d, m;  
        asm(  
            "divl %4; \n\t"  
            : "=a" (d), "=d" (m)  
            : "d" (xh), "a" (xl), "r" (MOD)  
        );  
        x = m;  
        return *this;  
    }  
    modint &operator/=(modint r) { return *this *= r.inv(); }  
    // increment, decrement  
    modint operator++() { x++; return *this; }  
    modint operator++(signed) { modint t = *this; x++; return t; }  
    modint operator--() { x--; return *this; }  
    modint operator--(signed) { modint t = *this; x--; return t; }  
    // 平方剰余のうち一つを返す なければ-1  
    friend modint sqrt(modint a) {  
        if(a == 0) return 0;  
        ll q = MOD-1, s = 0;  
        while((q&1)==0) q>>=1, s++;  
        modint z=2;  
        while(1) {  
            if(z.pow((MOD-1)/2) == MOD-1) break;  
            z++;  
        }  
        modint c = z.pow(q), r = a.pow((q+1)/2), t = a.pow(q);  
        ll m = s;  
        while(t.x>1) {  
            modint tp=t;  
            ll k=-1;  
            FOR(i, 1, m) {  
                tp *= tp;  
                if(tp == 1) { k=i; break; }  
            }  
            if(k==-1) return -1;  
            modint cp=c;  
            REP(i, m-k-1) cp *= cp;  
            c = cp*cp, t = c*t, r = cp*r, m = k;  
        }  
        return r.x;  
    }  
  
    template<class T>  
    friend modint operator*(T l, modint r) { return modint(l) *= r; }  
    template<class T>  
    friend modint operator+(T l, modint r) { return modint(l) += r; }  
    template<class T>  
    friend modint operator-(T l, modint r) { return modint(l) -= r; }  
    template<class T>  
    friend modint operator/(T l, modint r) { return modint(l) /= r; }  
    template<class T>  
    friend bool operator==(T l, modint r) { return modint(l) == r; }  
    template<class T>  
    friend bool operator!=(T l, modint r) { return modint(l) != r; }  
    // Input/Output  
    friend ostream &operator<<(ostream& os, modint a) { return os << a.x; }  
    friend istream &operator>>(istream& is, modint &a) {   
        is >> a.x;  
        a.x = ((a.x%MOD)+MOD)%MOD;  
        return is;  
    }  
    friend string to_frac(modint v) {  
        static map<ll, PII> mp;  
        if(mp.empty()) {  
            mp[0] = mp[MOD] = {0, 1};  
            FOR(i, 2, 1001) FOR(j, 1, i) if(__gcd(i, j) == 1) {  
                mp[(modint(i) / j).x] = {i, j};  
            }  
        }  
        auto itr = mp.lower_bound(v.x);  
        if(itr != mp.begin() && v.x - prev(itr)->first < itr->first - v.x) --itr;  
        string ret = to_string(itr->second.first + itr->second.second * ((int)v.x - itr->first));  
        if(itr->second.second > 1) {  
            ret += '/';  
            ret += to_string(itr->second.second);  
        }  
        return ret;  
    }  
};  
using mint = modint<998244353>;  
  
const mint inv6 = mint(6).inv();  
mint combi(ll N, ll K) {  
    if(N < K) return 0;  
    mint n(N);  
    return n*(n-1)*(n-2)*inv6;  
}  
  
void solve(ll x1, ll y1, ll x2, ll y2, ll x3, ll y3, ll k) {  
    // pv_1 + qv_2 + rv_3 = 0 となるような p,q,r  
    ll p = x3*y2-x2*y3, q = x1*y3-x3*y1, r = x2*y1-x1*y2;  
    ll g = __gcd(__gcd(p, q), r);  
    p /= g, q /= g, r /= g;  
    // p>0, q>0, r!=0 となるようにswapする  
    if((p<0) + (q<0) + (r<0) >= 2) {  
        p *= -1, q *= -1, r *= -1;  
    }  
    if(p < 0) {  
        swap(p, r);  
        swap(x1, x3);  
        swap(y1, y3);  
    } else if(q < 0) {  
        swap(q, r);  
        swap(x2, x3);  
        swap(y2, y3);  
    }  
  
    if(r > 0) {  
        // a>=p, b>=q, c>=r となっているときは a-p, b-q, c-r をそれぞれの係数にしたとき同じものがでてくる (pv_1 + qv_2 + rv_3 = 0 なので)  
        // a+b+c<k && (a<p || b<q || c<r) となるのは何通り？  
        mint tmp1 = combi(k+3, 3), tmp2 = combi(k-p-q-r+3, 3);  
        mint ret = tmp1 - tmp2;   
        dump(k+3, k-p-q-r+3, tmp1, tmp2, ret);  
        cout << ret << endl;   
    } else {  
        if(p+q < -r) {  
            mint ret = combi(k+3, 3) - combi(k+r+3, 3);  
            dump(k+3, k+r+3, ret);  
            cout << ret << endl;  
        } else {  
            mint tmp1 = combi(k+3, 3), tmp2 = combi(k-p-q+3, 3);  
            mint ret = tmp1 - tmp2;  
            dump(k+3, k-p-q+3, tmp1, tmp2, ret);  
            cout << ret << endl;  
        }  
    }  
}  
  
int main(void) {  
    ll q;  
    cin >> q;  
    while(q--) {  
        ll x1, y1, x2, y2, x3, y3, k;  
        cin >> x1 >> y1 >> x2 >> y2 >> x3 >> y3 >> k;  
        solve(x1, y1, x2, y2, x3, y3, k);  
    }  
  
    return 0;  
}

ferinの競プロ帳

競プロについてのメモ

Mujin Programming Challenge 2018 G - 移動