Submission #1550832

---

Source Code Expand

/**
 * code generated by JHelper
 * More info: https://github.com/AlexeyDmitriev/JHelper
 * @author
 */

#ifndef MAJK_LIB
#define MAJK_LIB

#include <vector>
#include <stack>
#include <iostream>
#include <unordered_map>
#include <map>
#include <iomanip>
#include <set>
#include <functional>
#include <fstream>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <string>
#include <sstream>
#include <queue>
#include <bitset>
using namespace std;

#define x first
#define y second
constexpr int MOD = 1000000007;

typedef std::pair<int,int> pii;
typedef long long ll;
typedef unsigned long long ull;
typedef unsigned int ui;

template <typename T, typename U> std::istream&operator>>(std::istream&i, pair<T,U>&p) {i >> p.x >> p.y; return i;}
template<typename T>std::istream&operator>>(std::istream&i,vector<T>&t) {for(auto&v:t){i>>v;}return i;}
template<typename T>std::ostream&operator<<(std::ostream&o,const vector<T>&t) {for(size_t i=0;i<t.size();++i){o<<t[i]<<" \n"[i == t.size()-1];}return o;}
template<typename T> using minheap = priority_queue<T, vector<T>, greater<T>>;
template<typename T> using maxheap = priority_queue<T, vector<T>, less<T>>;
auto fraclt = [](const pii&a,const pii&b) { return (ll)a.x * b.y < (ll)b.x * a.y; };
struct cmpfrac { bool operator()(const pii&a,const pii&b)const { return (ll)a.x * b.y < (ll)b.x * a.y; }};
template <typename T> bool in(T a, T b, T c) { return a <= b && b < c; }
int logceil(ll x) {int b=0;while(x){x>>=1;++b;}return b;}

namespace std {
    template<typename T,typename U>struct hash<pair<T,U>>{hash<T>t;hash<U>u;size_t operator()(const pair<T,U>&p)const{return t(p.x)^(u(p.y)<<7);}};
}
template<typename T,typename F>T bsh(T l,T h,const F&f){T r=-1,m;while(l<=h){m=(l+h)/2;if(f(m)){l=m+1;r=m;}else{h=m-1;}}return r;}
template<typename T,typename F>T bsl(T l,T h,const F&f){T r=-1,m;while(l<=h){m=(l+h)/2;if(f(m)){h=m-1;r=m;}else{l=m+1;}}return r;}
template<typename T> T gcd(T a,T b) { if (a<b) swap(a,b); return b?gcd(b,a%b):a; }

template<typename T>class vector2:public vector<vector<T>>{public:vector2(){} vector2(size_t a,size_t b,T t=T()):vector<vector<T>>(a,vector<T>(b,t)){}};
template<typename T>class vector3:public vector<vector<vector<T>>>{public:vector3(){} vector3(size_t a,size_t b,size_t c,T t=T()):vector<vector<vector<T>>>(a,vector<vector<T>>(b,vector<T>(c,t))){}};
template <typename T> struct bounded_priority_queue {
	inline bounded_priority_queue(ui X) : A(X), B(0) {}
	inline void push(ui L, T V) { B = max(B, L); A[L].push(V); }
	inline const T &top() const { return A[B].front(); }
	inline void pop() { A[B].pop(); while (B > 0 && A[B].empty()) --B; }
	inline bool empty() const { return A[B].empty(); }
	inline void clear() { B = 0; for (auto &a: A) a = queue<T>(); }
private:
	vector<queue<T>> A; ui B;
};


#endif
#ifndef MOD_H
#define MOD_H



template <unsigned int N> class Field {
    typedef unsigned int ui;
    typedef unsigned long long ull;
	inline ui pow(ui a, ui p){ui r=1,e=a;while(p){if(p&1){r=((ull)r*e)%N;}e=((ull)e*e)%N;p>>=1;}return r;}
	/*extended GCD(slow):ll t=0,nt=1,r=N,nr=a;while(nr){ll q=r/nr;t-=q*nt;swap(t,nt);r-=q*nr;swap(r,nr);}assert(r<=1);return(t<0)?t+N:t;*/
	inline ui inv(ui a){return pow(a,N-2);}
public:
    inline Field(int x = 0) : v(x) {}
	inline Field<N> pow(int p){return (*this)^p; }
	inline Field<N> operator^(int p){return {(int)pow(v,(ui)p)};}
    inline Field<N>&operator+=(const Field<N>&o) {if ((ll)v+o.v >= N) v = (ll)v+o.v-N; else v = v+o.v; return *this; }
    inline Field<N>&operator-=(const Field<N>&o) {if (v<o.v) v = N-o.v+v; else v=v-o.v; return *this; }
    inline Field<N>&operator*=(const Field<N>&o) {v=(ull)v*o.v % N; return *this; }
    inline Field<N>&operator/=(const Field<N>&o) { return *this*=inv(o.v); }
    inline Field<N> operator+(const Field<N>&o) const {Field<N>r{*this};return r+=o;}
    inline Field<N> operator-(const Field<N>&o) const {Field<N>r{*this};return r-=o;}
    inline Field<N> operator*(const Field<N>&o) const {Field<N>r{*this};return r*=o;}
    inline Field<N> operator/(const Field<N>&o) const {Field<N>r{*this};return r/=o;}
    inline Field<N> operator-() {if(v) return {(int)(N-v)}; else return {0};};
    inline Field<N>& operator++() { ++v; if (v==N) v=0; return *this; }
    inline Field<N> operator++(int) { Field<N>r{*this}; ++*this; return r; }
    inline Field<N>& operator--() { --v; if (v==-1) v=N-1; return *this; }
    inline Field<N> operator--(int) { Field<N>r{*this}; --*this; return r; }
    inline bool operator==(const Field<N>&o) const { return o.v==v; }
	inline bool operator!=(const Field<N>&o) const { return o.v!=v; }
	inline explicit operator ui() const { return v; }
	inline static vector<Field<N>>fact(int t){vector<Field<N>>F(t+1,1);for(int i=2;i<=t;++i){F[i]=F[i-1]*i;}return F;}
	inline static vector<Field<N>>invfact(int t){vector<Field<N>>F(t+1,1);Field<N> X{1};for(int i=2;i<=t;++i){X=X*i;}F[t]=1/X;for(int i=t-1;i>=2;--i){F[i]=F[i+1]*(i+1);}return F;}
private: ui v;
};
template<unsigned int N>istream &operator>>(std::istream&is,Field<N>&f){unsigned int v;is>>v;f=v;return is;}
template<unsigned int N>ostream &operator<<(std::ostream&os,const Field<N>&f){return os<<(unsigned int)f;}
template<unsigned int N>Field<N> operator+(int i,const Field<N>&f){return Field<N>(i)+f;}
template<unsigned int N>Field<N> operator-(int i,const Field<N>&f){return Field<N>(i)-f;}
template<unsigned int N>Field<N> operator*(int i,const Field<N>&f){return Field<N>(i)*f;}
template<unsigned int N>Field<N> operator/(int i,const Field<N>&f){return Field<N>(i)/f;}


typedef Field<1000000007> FieldMod;

struct Ring {
	static int div(int p, int q, int N) {
		ll t=0,nt=1,r=N,nr=q;
		while(nr){ ll q=r/nr;t-=q*nt;r-=q*nr;swap(t,nt);swap(r,nr); }
		t=(t<0)?t+N:t;
		r=(r<0)?r+N:r;
		if (gcd(p,N)%r) { return 0; }
		return (int)(((ll)t*(ll)p/r)%N);
	}
};
#endif



#include <array>

template <ui P> struct Group {};
template <> struct Group<4194304001u> { constexpr static ui bits = 25, root = 199, rootInverse = 758768563u; };
template <> struct Group<998244353u> { constexpr static ui bits = 23, root = 31, rootInverse = 128805723; };
template <> struct Group<104857601u> { constexpr static ui bits = 22, root = 21, rootInverse = 49932191; };
template <> struct Group<924844033u> { constexpr static ui bits = 21, root = 3597, rootInverse = 508059997; };

template <ui Prime> class FFT {
public:
	typedef Field<Prime> F;

	FFT() {
		OMEGA[Group<Prime>::bits] = Group<Prime>::root;
		OMEGA_INV[Group<Prime>::bits] = Group<Prime>::rootInverse;
		for (size_t i = Group<Prime>::bits; i > 0; --i) {
			OMEGA[i-1] = ui((ull(OMEGA[i])*ull(OMEGA[i]))%Prime);
			OMEGA_INV[i-1] = ui((ull(OMEGA_INV[i])*ull(OMEGA_INV[i]))%Prime);

			// cerr << OMEGA[i-1] << ' ' << OMEGA_INV[i-1] << endl;
		}
	}

	void fft(std::vector<F> &V, bool inverse = false) {
		int P = 1;
		while ((1u << P) < V.size()) ++P;
		V.resize(1u << P);
		fft(V, P, inverse);
	}

	void fft(std::vector<F> &V, const int P, bool inverse = false) {
		if (P == 0) return;
		if (P > Group<Prime>::bits) throw std::runtime_error("Number of bits in FFT too large");
		if (V.size() != (1 << P)) throw std::runtime_error("Vector has wrong size");
		for (int i = 1, j = 0; i < (1 << P); ++i) {
			int bit = (1 << (P - 1));
			for (; j >= bit; bit >>= 1) j -= bit;
			j += bit;
			if (i < j) std::swap(V[i], V[j]);
		}

		for (int B = 1; B <= P; ++B) {
			F omega = inverse ? OMEGA_INV[B] : OMEGA[B];
			for (int i = 0; i < (1 << P); i += (1 << B)) {
				F x = 1;
				for (int j = 0; j < (1 << (B - 1)); j++) {
					F u = V[i + j], v = V[i + j + (1 << (B - 1))] * x;
					V[i + j] = u + v;
					V[i + j + (1 << (B - 1))] = u - v;
					x *= omega;
				}
			}
		}

		if (inverse) {
			F q = 1 / F(1 << P);
			for (F &v:V) v *= q;
		}
	}

private:
	std::array<ui, Group<Prime>::bits+1> OMEGA, OMEGA_INV;
};


typedef Field<998244353> TF;
class F {
public:
    FFT<998244353u> Z;
    vector<TF> All, Conv1, Conv2, Ans;

    void perform(int pos, int size) {
        if (size <= 128) {
            for (int i = 1; i < size; ++i) {
                for (int j = 0; j < i; ++j) {
                    Ans[pos+i] += Ans[pos+j] * Conv1[i-j-1];
                }
            }
            return;
        }

        perform(pos, size/2);

        vector<TF> Cur(size, 0);
        vector<TF> Mul(size, 0);
        for (int i = 0; i < size/2; ++i) Cur[i] = Ans[i+pos];
        for (int i = 0; i < size; ++i) Mul[i] = Conv1[i];
        Z.fft(Cur);
        Z.fft(Mul);
        for (int i = 0; i < size; ++i) Cur[i] *= Mul[i];
        Z.fft(Cur, true);
        for (int i = 0; i < size/2; ++i) Ans[pos+size/2+i] += Cur[i+size/2-1];

        perform(pos+size/2, size/2);
    }


    void solve(istream& cin, ostream& cout) {
        int N, M; cin >> N >> M;
        int size = 1 << logceil(min(N, M));

        auto F = TF::fact(max(M+N,2*size));
        auto I = TF::invfact(max(M+N,2*size));

        All.clear(); Conv2.clear();
        All.resize(2*size); Conv1.resize(size); Conv2.resize(2*size);
        for (int i = 0; i < size; ++i) {
            All[i] = F[2*i] * I[i] * I[i];
            Conv1[i] = i>0 ? All[i] - F[2*i] * I[i+1] * I[i-1] : 1;
            Conv2[i] = Conv1[i] * (2*i+3);
            Conv1[i] *= 2;
        }

        Z.fft(All);
        Z.fft(Conv2);
        for (int i = 0; i < 2*size; ++i) Conv2[i] *= All[i];
        Z.fft(Conv2, true);
        Ans = vector<TF>(2*size, 0);
        for (int i = 1; i < size; ++i) {
            Ans[i] = Conv2[i - 1];
        }

        perform(0, size);
        if (N==M) {
            cout << Ans[N] * I[2*N] * F[N] * F[N] << endl;
        } else {
            if (M<N) swap(N,M);
            TF ret = Ans[N]+F[2*N]*I[N]*I[N]*(M-N);
            for (int i = 0; i < N; ++i) {
                TF ways = F[N-i + M-i-1] * I[N-i] * I[M-i-1] - F[N-i + M-i-1] * I[N-i-1] * I[M-i];
                ret += ways * Ans[i] + ways * F[2*i] * I[i] * I[i] * (M-i);
            }
            cout << ret * I[N+M] * F[N] * F[M] << endl;
        }
    }
};


int main() {
	ios_base::sync_with_stdio(false); cin.tie(nullptr); cout.tie(nullptr);
	F solver;
	std::istream& in(std::cin);
	std::ostream& out(std::cout);
	solver.solve(in, out);
    return 0;
}

Submission Info

Judge Result