Astral-23 / competitive-programing

競プロ用のライブラリ。 Utility/template.hpp を include しないと動かない物が多数。
Creative Commons Zero v1.0 Universal
0 stars 0 forks source link

math 1d function #1

Open Astral-23 opened 2 months ago

Astral-23 commented 2 months ago

TT T safe_floor(T a, T b) {
    if (b < 0) a = -a, b = -b;
    return a >= 0 ? a / b : (a + 1) / b - 1;
}
TT T safe_ceil(T a, T b) {
    if (b < 0) a = -a, b = -b;
    return a > 0 ? (a - 1) / b + 1 : a / b;
}

vec<pair<ll, ll>> factorize(ll x) {
    vec<pair<ll, ll>> res;
    for (ll i = 2; i * i <= x; i++) {
        if (i != 2 && !(i & 1)) continue;
        int c = 0;
        while (x % i == 0) {
            x /= i;
            c++;
        }

        if (c) res.emplace_back(i, c);
    }
    if (x != 1) res.emplace_back(x, 1);
    return res;
}

ll extgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }

    ll d = extgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

ll modinv(ll a, ll MOD) {
    ll x, y;
    extgcd(a, MOD, x, y);
    return (x % MOD + MOD) % MOD;
}

TT struct f1d {
    T a, b;

    f1d(T A = 0, T B = 0) : a(A), b(B) {}

    f1d composition(f1d g) {
        f1d res;
        res.a = a * g.a;
        res.b = a * g.b + b;
        return res;
    }

    T operator()(T x) {return a * x + b;}

    friend bool operator==(f1d f, f1d g) { return f.a == g.a && f.b == g.b; }
    friend bool operator!=(f1d f, f1d g) { return !(f == g); }
    friend bool operator<(f1d f, f1d g) {
        if (f.a != g.a) {
            return f.a < g.a;
        } else {
            return f.b < g.b;
        }
    }
    friend bool operator>(f1d f, f1d g) {
        if (f.a != g.a) {
            return f.a > g.a;
        } else {
            return f.b > g.b;
        }
    }
    friend bool operator<=(f1d f, f1d g) { return f < g || f == g; }
    friend bool operator>=(f1d f, f1d g) { return f > g || f == g; }
    friend ostream &operator<<(ostream &os, f1d f) {
        return os << f.a << "x + " << f.b;
    }
    friend f1d operator+(f1d f, f1d g) {
        f.a += g.a;
        f.b += g.b;
        return f;
    }
    friend f1d operator-(f1d f, f1d g) {
        f.a -= g.a;
        f.b -= g.b;
        return f;
    }
    friend f1d operator+=(f1d &f, f1d g) { return f = f + g; }
    friend f1d operator-=(f1d &f, f1d g) { return f = f - g; }
    template<typename S> friend f1d operator*(f1d f, S s) {
        f.a *= s;
        f.b *= s;
        return f;
    }
    template<typename S> friend f1d operator/(f1d f, S s) {
        f.a /= s;
        f.b /= s;
        return f;
    }
    template<typename S> friend f1d operator*=(f1d &f, S s) { return f = f * s; }
    template<typename S> friend f1d operator/=(f1d &f, S s) { return f = f / s; }
};

TT bool general_solution(ll a, f1d<T> &x, T b, f1d<T> &y, ll c) {
    // ax + by = c の x, yそれぞれの一般解f, g。存在する...true しない...false
    // LLONG_MIN <= a, b, c <= LLONG_MAX
    if (a == 0 && b == 0 && c == 0) {
        x = f1d<T>(1, 0);
        y = f1d<T>(1, 0);
        return true;
    }

    ll gc = abs(gcd(a, gcd(b, c)));
    a /= gc;
    b /= gc;
    c /= gc;

    ll s, t;
    ll d = extgcd(a, b, s, t);
    if(d == -1) {
        d = -d; s = -s; t = -t;
    }
    if (d != 1) {
        return false;
    } else {
        x = f1d<T>(-b, T(s) * c);
        y = f1d<T>(a, T(t) * c);
        return true;
    }
}

TT T count_x(T l, f1d<T> fx, T r) {
    // l <= f(x) <= rとなるxの個数。xは任意の時、-1。
    if (fx.a == 0) {
        if (l <= fx.b && fx.b <= r)
            return -1;
        else
            return 0;
    }

    l -= fx.b;
    r -= fx.b;
    if (fx.a < 0) {
        l = -l;
        r = -r;
        swap(l, r);
        fx.a = -fx.a;
    }
    T li = safe_ceil<T>(l, fx.a);
    T ri = safe_floor<T>(r, fx.a);

    return max<T>(T(0), ri - li + 1);
}
/*
l <= f(x) <= r : count_x(l, fx, r);
l <= f(x) < r  : count_x(l, fx, r - 1);
l < f(x) <= r  : count_x(l + 1, fx, r);
l < f(x) < r   : count_x(l + 1, fx, r - 1);
*/

/*
   C(n) = n * (1 - 1/P_1) * (1 - 1/P_2) * ...
   n, m互いに素 : C(nm) = C(n)C(m)

*/

https://atcoder.jp/contests/abc186/submissions/55989132 verified : https://jag-icpc.org/?2022%2FPractice%2F%E6%A8%A1%E6%93%AC%E5%9B%BD%E5%86%85%E4%BA%88%E9%81%B8%2F%E5%95%8F%E9%A1%8C%E6%96%87%E3%81%A8%E3%83%87%E3%83%BC%E3%82%BF%E3%82%BB%E3%83%83%E3%83%88

Astral-23 commented 2 months ago

使用例

int main() {

    //---- 1次関数 f1d-----------
    f1d<ll> f(3, 2);    // f(x) = (ll)(3) x + (ll)(2) を宣言。
    cout << f << endl;  // 3x + 2
    cout << f.a << " " << f.b << endl; // 3  2      直接メンバ変数へのアクセス。

    f1d<ll> g(10, 3);  // g(x) = 10x + 3;

    if (f == g) cout << "ERROR" << endl;
    if (f != g) cout << "f != g" << endl;
    if (f < g) cout << "辞書順比較" << endl;
    if (f > g) cout << "ERROR" << endl;
    if (f <= g) cout << "辞書順比較" << endl;
    if (f >= g) cout << "ERROR" << endl;

    f1d<ll> h = f.composition(g); //f(g(x))
    cout << h << endl; //30x + 11

    // ---- T general_solution<T>(ll a,ll b,f1d<T> xf,  f1d<T> yf, ll c)--------
    f1d<ll> xf, yf;
    bool res = general_solution<ll>(-3, 2, xf, yf, -5); 
    /* -3x + 2y = -5の、x, yの一般解について
    戻り値 :
      x, yの一般解が存在するならtrue
      しなければfalse
    xf, yf : 
     それぞれx, yの一般解が(存在するなら)代入される。

    注意点:
    もしa, b, cがlong long相当の大きさなら、f1dの型は__int128_tにする。
     また、テンプレート引数(=型)は省略しないで、一番大きいやつに合わせる。
        理由:
        ax + by = cとする。
        xf.a, yf.a それぞれについて、大きさのオーダーは O(max(|a||c|, |b||c|))
        よって、そこでオーバーフロウが発生し得る。
        O(log(max(a, b))) ぐらい? 
    */

   cout << xf << endl; // x = -2k + (-5) として存在した様だ。
   cout << yf << endl; // y = -3k + (-5) として存在した様だ。
   assert(res == true);

   // ---- T count_x(T l, f1d<T> fx, T r) --- 
   cout << xf << endl; // xf(k) = -2k - 5

   cout << count_x<ll>(-6, xf, 5) << endl; // 6
   /*
      -6 <= -2k - 5 <= 5 
      を満たすkの個数を返す。

      特にオーバーフロウは起きない。
      ただし、渡すf1bの型とl, rの型が異なる場合、より大きい方の型に合わせなければならない。
      O(1)
   */
}
Astral-23 commented 2 months ago
TT struct f1d {
    T a, b;

    f1d(T A = 0, T B = 0) : a(A), b(B) {}

    f1d composition(f1d g) {
        f1d res;
        res.a = a * g.a;
        res.b = a * g.b + b;
        return res;
    }

    friend bool operator==(f1d f, f1d g) { return f.a == g.a && f.b == g.b; }
    friend bool operator!=(f1d f, f1d g) { return !(f == g); }
    friend bool operator<(f1d f, f1d g) {
        if (f.a != g.a) {
            return f.a < g.a;
        } else {
            return f.b < g.b;
        }
    }
    friend bool operator>(f1d f, f1d g) {
        if (f.a != g.a) {
            return f.a > g.a;
        } else {
            return f.b > g.b;
        }
    }
    friend bool operator<=(f1d f, f1d g) { return f < g || f == g; }
    friend bool operator>=(f1d f, f1d g) { return f > g || f == g; }
    friend ostream &operator<<(ostream &os, f1d f) {
        return os << f.a << "x + " << f.b;
    }
    friend f1d operator+(f1d f, f1d g) {
        f.a += g.a;
        f.b += g.b;
        return f;
    }
    friend f1d operator-(f1d f, f1d g) {
        f.a -= g.a;
        f.b -= g.b;
        return f;
    }
    friend f1d operator+=(f1d &f, f1d g) { return f = f + g; }
    friend f1d operator-=(f1d &f, f1d g) { return f = f - g; }
    template<typename S> friend f1d operator*(f1d f, S s) {
        f.a *= s;
        f.b *= s;
        return f;
    }
    template<typename S> friend f1d operator/(f1d f, S s) {
        f.a /= s;
        f.b /= s;
        return f;
    }
    template<typename S> friend f1d operator*=(f1d &f, S s) { return f = f * s; }
    template<typename S> friend f1d operator/=(f1d &f, S s) { return f = f / s; }
};

・初期化の際、何も入れなかったら(0, 0)にするようにした(以前は未定義だった) ・スカラーに対する演算・f1d同士の演算を追加した https://atcoder.jp/contests/abc189/submissions/55445203