#1025750 (C++23) No.2959 Dolls' Tea Party

提出ソース
結果

問題	No.2959 Dolls' Tea Party
ユーザー	nouka28
提出日時	2024-11-09 07:17:09
言語	C++23 (gcc 13.3.0 + boost 1.87.0)
結果	TLE
実行時間	-
コード長	17,495 bytes
コンパイル時間	10,503 ms
コンパイル使用メモリ	339,948 KB
実行使用メモリ	129,996 KB
最終ジャッジ日時	2024-11-09 07:17:32
合計ジャッジ時間	21,867 ms
ジャッジサーバーID （参考情報）	judge3 / judge1
このコードへのチャレンジ
（要ログイン）
ファイルパターン	結果
sample	AC * 4
other	AC * 5 TLE * 1 -- * 27
権限があれば一括ダウンロードができます
ソースコード

raw source code
#include<bits/stdc++.h>
using namespace std;

#include<atcoder/all>
using namespace atcoder;
using mint=atcoder::modint998244353;

#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")

#define int long long

using ld=long double;

template<class mint>
struct FormalPowerSeries:vector<mint>{
	using vector<mint>::vector;
	using vector<mint>::operator=;
	using fps=FormalPowerSeries;
	using sfps=vector<pair<int,mint>>;

	fps& operator+=(const fps& g){
		if(g.size()>this->size())this->resize(g.size());
		for(int i=0;i<(int)g.size();i++)(*this)[i]+=g[i];
		return *this;
	}

	fps& operator+=(const mint& v){
		if(this->empty())this->resize(1);
		(*this)[0]+=v;
		return *this;
	}

	fps operator+(const fps& g)const{return fps(*this)+=g;}
	fps operator+(const mint& v)const{return fps(*this)+=v;}
	friend fps operator+(const mint& v,const fps& f){return f+v;}
	fps& operator+=(const int& v){*this+=mint(v);return *this;}
	fps operator+(const int& v){return fps(*this)+=v;;}
	friend fps operator+(const int& v,const fps& f){return f+v;}

	fps& operator-=(const fps& g){
		if(g.size()>this->size())this->resize(g.size());
		for(int i=0;i<(int)g.size();i++)(*this)[i]-=g[i];
		return *this;
	}

	fps& operator-=(const mint& v){
		if(this->empty())this->resize(1);
		(*this)[0]-=v;
		return *this;
	}
	fps operator-(const fps& g)const{return fps(*this)-=g;}
	fps operator-(const mint& v)const{return fps(*this)-=v;}
	friend fps operator-(const mint& v,const fps& f){return -(f-v);}
	fps& operator-=(const int& v){*this-=v;return *this;}
	fps operator-(const int& v){return fps(*this)-=v;}
	friend fps operator-(const int& v,const fps& f){return -(f-v);}

	fps operator-()const{return fps(*this)*=-1;}

	fps& operator*=(const mint& v){for(auto&e:*this)e*=v;return *this;}
	fps operator*(const mint& v)const{return fps(*this)*=v;}
	friend fps operator*(const mint& v,const fps& f){return f*v;}
	fps& operator*=(const int& v){*this*=mint(v);return *this;}
	fps operator*(const int& v)const{return fps(*this)*=v;}
	friend fps operator*(const int&v,const fps& f){return f*v;}

	fps& operator<<=(const int d){
		this->insert(this->begin(),d,0);
		return *this;
	}

	fps operator<<(const int d)const{return fps(*this)<<=d;}
	
	fps& operator>>=(const int d){
		this->erase(this->begin(),this->begin()+min((int)this->size(),d));
		return *this;
	}

	fps operator>>(const int d)const{return fps(*this)>>=d;}

	//fast
	fps& operator*=(const fps& g){
		*this=atcoder::convolution(*this,g);
		return *this;
	}
	
	//naive
	// fps& operator*=(const fps& g){
	// 	this->resize(this->size()+g.size()-1);
	// 	for(int i=(int)this->size()-1;i>=0;i--){
	// 		for(int j=1;j<=(int)g.size();j++){
	// 			if(i+j>=(int)this->size())break;
	// 			(*this)[i+j]+=(*this)[i]*g[j];
	// 		}
	// 		(*this)[i]*=g[0];
	// 	}
	// 	return *this;
	// }

	fps operator*(const fps& g)const{return fps(*this)*=g;}

	fps inv(int d)const{
		fps g={(*this)[0].inv()};
		for(int k=1;k<d;k*=2){
			g=(2-*this*g)*g;
			g.resize(2*k);
		}
		g.resize(d+1);
		return g;
	}

	fps& operator/=(const fps& g){return *this=fps(*this*=g.inv(this->size())).pre(this->size());}
	fps& operator/=(const mint& v){for(auto&e:*this)e/=v;return *this;}
	fps operator/(const fps& g)const{return fps(*this)/=g;}
	fps operator/(const mint& v)const{return fps(*this)/=v;}
	
	fps quotient(const fps& g)const{
		if(this->size()<g.size())return fps();
		return (fps(this->rev()/g.rev()).pre(this->size()-g.size()+1)).rev();
	}
	fps reminder(const fps& g)const{return fps(*this-this->quotient(g)*g).pre(g.size()-1);}
	
	pair<fps,fps> quotient_reminder(const fps& g)const{
		pair<fps,fps> res;
		res.first=this->quotient(g);
		res.second=fps(*this-res.first*g).pre(g.size()-1);
		return res;
	}
	
	void shrink(){
		while(this->size()&&this->back()==mint(0))this->pop_back();
	}
	fps rev()const{fps g(*this);reverse(g.begin(),g.end());return g;}
	fps dot(fps g)const{
		fps res(min(this->size(),g.size()));
		for(int i=0;i<(int)res.size();i++)res[i]=(*this)[i]*g[i];
		return res;
	}
	fps pre(int d)const{
		fps res(begin(*this),begin(*this)+min((int)this->size(),d));
		if((int)res.size()<d)res.resize(d);
		return res;
	}

	fps& operator*=(const sfps& g){
		auto it0=g.begin();
		mint g0=0;
		if(it0->first==0){
			g0=it0->second;
			it0++;
		}
		for(int i=this->size()-1;i>=0;i--){
			for(auto it=it0;it!=g.end();it++){
				auto[j,gj]=*it;
				if(i+j>=this->size())break;
				(*this)[i+j]+=(*this)[i]*gj;
			}
			(*this)[i]*=g0;
		}
		return *this;
	}
	fps operator*(const sfps& g)const{return fps(*this)*=g;}

	fps& operator/=(const sfps& g){
		auto it0=g.begin();
		assert(it0->first==0&&it0->second!=0);
		mint g0_inv=it0->second.inv();
		it0++;
		for(int i=0;i<(int)this->size();i++){
			(*this)[i]*=g0_inv;
			for(auto it=it0;it!=g.begin();it++){
				auto[j,gj]=*it;
				if(i+j>=this->size())break;
				(*this)[i+j]-=(*this)[i]*gj;
			}
		}
		return *this;
	}
	fps operator/(const sfps& g)const{return fps(*this)/=g;}

	fps pow(long long d,const fps& g)const{
		fps res={1},pow2=*this;
		while(d>0){
			if(d&1)res=(res*pow2).reminder(g);
			pow2=(pow2*pow2).reminder(g);
			d>>=1;
		}
		return res;
	}

	fps derivative()const{
		fps res;
		for(int i=1;i<(int)this->size();i++)res.push_back((*this)[i]*i);
		return res;
	}

	fps integral()const{
		fps res={0};
		for(int i=0;i<(int)this->size();i++)res.push_back((*this)[i]/(i+1));
		return res;
	}

	fps log(int d)const{
		return fps(this->derivative()*this->inv(d)).integral().pre(d);
	}

	fps exp(int d)const{
		fps g={1};
		for(int k=1;k<d;k*=2){
			g=g*(*this+1-g.log(2*k));
			g.resize(2*k);
		}
		return g.pre(d);
	}

	fps pow(long long k,int d)const{
		if(k==0){
			fps res(d,mint(0));
			if(d)res[0]=1;
			return res;
		}
		int i0=0;
		while(i0<(int)this->size()&&(*this)[i0]==mint(0))i0++;
		if(i0==(int)this->size())return fps(d,mint(0));
		mint c0=(*this)[i0];
		fps fs=(*this>>i0)/c0;
		if(i0>=(d+k-1)/k)return fps(d,mint(0));
		int ds=(int)(d-k*i0);
		fps gs=fps(mint(k)*fs.log(ds)).exp(ds);
		fps g=fps(gs*c0.pow(k))<<(int)(k*i0);
		return g;
	}

	friend istream& operator>>(istream& is,fps&f){
		for(auto&e:f)cin>>e;
		return is;
	}
	friend ostream& operator<<(ostream& os,const fps& f){
		if((int)f.size()==0)os<<0;
		else{
			for(int i=0;i<(int)f.size();i++){
				os<<f[i].val();
				if(i<(int)f.size()-1)os<<" ";
			}
			return os;
		}
		return os;
	}
};

template<long long mod,long long MAX_N>
struct factional_prime{
	long long inv_[MAX_N+1];
    long long fac_[MAX_N+1];
    long long fac_inv_[MAX_N+1];

    factional_prime(){
        inv_[0]=0;inv_[1]=fac_[0]=fac_[1]=fac_inv_[0]=fac_inv_[1]=1;
        for(long long i=2;i<=MAX_N;i++){
            inv_[i]=((mod-mod/i)*inv_[mod%i])%mod;
            fac_[i]=(fac_[i-1]*i)%mod;
            fac_inv_[i]=(fac_inv_[i-1]*inv_[i])%mod;
        }
    }
    long long inv(long long n){
        if(n<0)return 0;
        return inv_[n];
    }
    long long fac(long long n){
        if(n<0)return 0;
        return fac_[n];
    }
    long long finv(long long n){
        if(n<0)return 0;
        return fac_inv_[n];
    }
    long long nCr(long long n,long long r){
        if(n<r||n<0||r<0)return 0;
        return ((fac_[n]*fac_inv_[n-r])%mod*fac_inv_[r])%mod;
    }
    long long nPr(long long n,long long r){
        if(n<r||n<0||r<0)return 0;
        return (fac_[n]*fac_inv_[n-r])%mod;
    }
};

factional_prime<998244353,5000000> fp;

using fps=FormalPowerSeries<mint>;
using sfps=vector<pair<int,mint>>;

//math
namespace nouka28{
	random_device rnd; 
	mt19937 mt(rnd());
	const long long MT_MAX=(1LL<<62)-1;
	uniform_int_distribution<long long> rd(0,MT_MAX);
	double randd(){
		return 1.0*rd(mt)/MT_MAX;
	}
	long long randint(long long a,long long b){
		// [a,b]の乱数を生成
		return a+rd(mt)%(b-a+1);
	}

    template<class T=long long>
    vector<T> Quotients(T n){
        vector<T> retl,retr;
        for(int i=1;i*i<=n;i++){
            retl.push_back(i);
            if(i<n/i)retr.push_back(n/i);
        }
        reverse(retr.begin(),retr.end());
        retl.insert(retl.end(),retr.begin(),retr.end());
        return retl;
    }
    template<class T=long long>T ceil_sqrt(T n){
        T l=-1,r=n;
        while(r-l>1){
            T m=(l+r)>>1;
            if(m*m>=n)r=m;
            else l=m;
        }
        return r;
    }

    //ceil(a/b)
    template<class T=long long>T ceil(T a,T b){
        if(a>=0){
            return (a+b-1)/b;
        }else{
            return (a)/b;
        }
    };

    //floor(a/b)
    template<class T=long long>T floor(T a,T b){
        if(a>=0){
        return a/b;
        }else{
            return -(-a+b-1)/b;
        }
    };

    //x^y mod m
    template<class T=long long>T modpow(T x,T y,T m){
        T res=1%m;x%=m;
        while(y){
            if(y%2)res=(res*x)%m;
            x=(x*x)%m;
            y>>=1;
        }
        return res;
    }

    //a^0+a^1+..+a^(n-1) (mod m)
    template<class T>
    T geometric_progression(T a,T n,T m){
        if(n==0)return 0;
        if(n%2==1){
            return (geometric_progression(a,n-1,m)*a+1)%m;
        }else{
            return (geometric_progression(a*a%m,n/2,m)*(1+a))%m;
        }
    };

    //素数判定(高速)
    bool is_prime(long long n){
        if(n<=1)return 0;
        if(n==2)return 1;
        if(n%2==0)return 0;
        long long s=0,d=n-1;
        while(d%2==0)d/=2,s++;
        if(n<4759123141LL){
            for(long long e:{2,7,61}){
                if(n<=e)break;
                long long t,x=modpow<__int128_t>(e,d,n);
                if(x!=1){
                    for(t=0;t<s;t++){
                        if(x==n-1)break;
                        x=__int128_t(x)*x%n;
                    }
                    if(t==s)return 0;
                }
            }
            return 1;
        }else{
            for(long long e:{2,325,9375,28178,450775,9780504,1795265022}){
                if(n<=e)break;
                long long t,x=modpow<__int128_t>(e,d,n);
                if(x!=1){
                    for(t=0;t<s;t++){
                        if(x==n-1)break;
                        x=__int128_t(x)*x%n;
                    }
                    if(t==s)return 0;
                }
            }
            return 1;   
        }
    }

    //Xor Shift
    unsigned xor_shift_rng(){
        static unsigned tx=123456789,ty=362436069,tz=521288629,tw=88675123;
        unsigned tt=(tx^(tx<<11));
        tx=ty,ty=tz,tz=tw;
        return (tw=(tw^(tw>>19))^(tt^(tt>>8)));
    }

    //ロー法 Nを割り切る素数を見つける
    long long pollard(long long n){
        if(n%2==0)return 2;
        if(is_prime(n))return n;
        long long step=0;
        while(true){
            long long r=(long long)xor_shift_rng();
            auto f=[&](long long x)->long long {return ((__int128_t(x)*x)%n+r)%n;};
            long long x=++step,y=f(x);
            while(true){
                long long p=__gcd(abs(y-x),n);
                if(p==0||p==n)break;
                if(p!=1)return p;
                x=f(x);
                y=f(f(y));
            }
        }
    }

    //internal fast factrize vector
    void _internal_factrize_vector(long long n,vector<long long>&v){
        if(n==1)return;
        long long p=pollard(n);
        if(p==n){v.push_back(p);return;}
        _internal_factrize_vector(p,v);
        _internal_factrize_vector(n/p,v);
    }

    //fast factrize vector
    vector<long long> factrize_vector(long long n){
        vector<long long> res;
        _internal_factrize_vector(n,res);
        sort(res.begin(),res.end());
        return res;
    }

    //internal fast factrize map
    void _internal_factrize_map(long long n,map<long long,long long>&v){
        if(n==1)return;
        long long p=pollard(n);
        if(p==n){v[p]++;return;}
        _internal_factrize_map(p,v);
        _internal_factrize_map(n/p,v);
    }

    //fast factrize map
    map<long long,long long> factrize_map(long long n){
        map<long long,long long> res;
        _internal_factrize_map(n,res);
        return res;
    }

    //fast factor
    vector<long long> factor(long long n){
        map<long long,long long> fm;_internal_factrize_map(n,fm);
        vector<long long> res={1};
        for(auto[i,j]:fm){
            vector<long long> tmp;
            int p=1;
            for(long long k=0;k<=j;k++){
                for(auto e:res){
                    tmp.push_back(e*p);
                }
                p*=i;
            }
            swap(res,tmp);
        }
        return res;
    }

    //euler phi function
    long long euler_phi(long long n){
        vector<long long> ps=factrize_vector(n);
        ps.erase(unique(ps.begin(),ps.end()),ps.end());
        for(long long p:ps){
            n/=p;n*=(p-1);
        }
        return n;
    }

    //ax+by=__gcd(a,b)
    template<class T=long long>
    tuple<T,T,T> extgcd(T a,T b){
        T x1=1,y1=0,d1=a,x2=0,y2=1,d2=b;
        while(d2!=0){
            T q=d1/d2,u=d1-d2*q,v=x1-x2*q,w=y1-y2*q;
            d1=d2;d2=u;x1=x2;x2=v;y1=y2;y2=w;
        }
        if(d1<0){
            d1=-d1;x1=-x1;y1=-y1;
        }
        return {d1,x1,y1};
    }

    //x inverse (mod m)
    long long modinv(long long a,long long m){
        long long b=m,u=1,v=0;
        while(b){
            long long t=a/b;
            a-=t*b;swap(a,b);
            u-=t*v;swap(u,v);
        }
        u%=m;
        if(u<0)u+=m;
        return u;
    }

	//find primitive root
	long long primitive_root(long long p){
		vector<long long> f=factrize_vector(p-1);
		f.erase(unique(f.begin(),f.end()),f.end());
		while(1){
			long long x=randint(1,p-1);
			bool flg=1;
			for(auto e:f)if(modpow<__int128_t>(x,(p-1)/e,p)==1){flg=0;break;}
			if(flg)return x;
		}
	}

    //x^k=y (mod m) __gcd(x,m)=1 k>=0
    long long discrete_logarithm_coprime_mod(long long x,long long y,long long m){
        x%=m;y%=m;
        if(y==1||m==1){
            return 0;
        }
        if(x==0){
            if(y==0)return 1;
            else return -1;
        }
        long long M=ceil_sqrt(m),a=modpow(modinv(x,m),M,m);
        unordered_map<long long,long long> mp;
        long long pow_x=1;
        for(long long i=0;i<M;i++){
            if(!mp.count(pow_x))mp[pow_x]=i;
            pow_x=pow_x*x%m;
        }
        long long ya=y;
        for(long long i=0;i<M;i++){
            if(mp.count(ya))return M*i+mp[ya];
            ya=ya*a%m;
        }
        return -1;
    }

    //x^k=y (mod m) __gcd(x,m)=1 k>=1
    long long discrete_Nlogarithm_coprime_mod(long long x,long long y,long long m){
        if(m==1){
            if(x==1)return 1;
            else return -1;
        }
        if(x==0){
            if(y==0)return 1;
            else return -1;
        }
        long long M=ceil_sqrt(m),a=modpow(modinv(x,m),M,m);
        unordered_map<long long,long long> mp;
        long long pow_x=1;
        for(long long i=0;i<=M;i++){
            if(!mp.count(pow_x))mp[pow_x]=i;
            pow_x=pow_x*x%m;
        }
        long long ya=y;
        for(long long i=0;i<M;i++){
            if(ya==1&&i>0)return i*M;
            else if(mp.count(ya))return M*i+mp[ya];
            ya=ya*a%m;
        }
        return -1;
    }

    //x^k=y (mod m) k>=0
    long long discrete_logarithm_arbitrary_mod(long long x,long long y,long long m){
        if(m==1){
            return 0;
        }
        x%=m;y%=m;
        long long d,pow_x=1;
        for(d=0;;d++){
            if(!(m>>d))break;
            if(pow_x==y){
                return d;
            }
            pow_x=pow_x*x%m;
        }
        long long g=__gcd(pow_x,m);
        if(y%g!=0){
            return -1;
        }
        m/=g;
        long long z=y*modinv(pow_x,m),t=discrete_logarithm_coprime_mod(x,z,m);
        if(t==-1)return -1;
        else return d+t;
    }

    //x^k=y (mod m) k>=1
    long long discrete_Nlogarithm_arbitrary_mod(long long x,long long y,long long m){
        if(m==1){
            if(x==1)return 1;
            else return -1;
        }
        x%=m;y%=m;
        long long d,pow_x=1;
        for(d=0;;d++){
            if(!(m>>d))break;
            if(pow_x==y&&d){
                return d;
            }
            pow_x=pow_x*x%m;
        }
        long long g=__gcd(pow_x,m);
        if(y%g!=0){
            return -1;
        }
        m/=g;
        long long z=y*modinv(pow_x,m),t;
        if(d)t=discrete_logarithm_coprime_mod(x,z,m);
        else t=discrete_Nlogarithm_coprime_mod(x,y,m);
        if(t==-1)return -1;
        else return d+t;
    }
}

signed main(){
	int N,K;cin>>N>>K;
	vector<int> A(N);for(auto&&e:A)cin>>e;

	vector<int> fa=nouka28::factor(K);

	int sz=fa.size();

	vector<mint> dp(sz);

	for(int S=0;S<sz;S++){

		int k=fa[S];

		int us=K/k;

		vector<int> cnt(us);

		int x=0;
		
		for(auto e:A){
			e/=k;
			if(e>=us){
				x++;
			}else{
				cnt[e]++;
			}
		}

		fps f={1};

		fps g={1};

		for(int i=1;i<us;i++){
			g.push_back(fp.finv(i));
			if(cnt[i]==0)continue;

			fps t=g.pow(cnt[i],us+1);

			f*=t;
			f.resize(us+1);
		}

		f.resize(us+1);

		mint ans=0;

		for(int i=0;i<=us;i++){
			mint tmp=f[i];

			tmp*=fp.finv(us-i);
			tmp*=mint(x).pow(us-i);

			ans+=tmp;
		}

		ans*=fp.fac(us);

		dp[S]=ans;
	}

	// for(auto&&e:fa)cout<<e<<" ";
	// cout<<endl;

	// for(auto&&e:dp)cout<<e.val()<<" ";
	// cout<<endl;

	mint ans=0;

	for(int i=sz-1;i>=0;i--){
		for(int j=i-1;j>=0;j--){
			if(fa[i]%fa[j]==0){
				dp[j]-=dp[i];
			}
		}

		ans+=dp[i]/(K/fa[i]);
	}

	cout<<ans.val()<<endl;
}
yukicoder

結果

ソースコード
