/********************************************************** Yukicoder No.940 ワープ ε=ε=ε=ε=ε=│;p>д<│ https://yukicoder.me/problems/no/940 <解法> n=X+Y+Z とおく。 nより十分に大きいNをとり、 T=(1+x+x^2+..+x^(N-1))(1+y+y^2+..+y^(N-1))(1+z+z^2+..+z^(N-1))-1 =(1-x^N)/(1-x) (1-y^N)/(1-y) (1-z^N)/(1-z) -1 とおくと、 F=Σ(m=1..n)T^m の x^X y^Y z^Z の係数を求めればよい。 ※理由:x^i y^j z^k というのが、1回の移動でx方向にi、y方向にj、z方向にk 進むことを表していて、 m回移動したときに(X,Y,Z)に到達する場合の数が T^m の x^X y^Y z^Z の係数と等しいので。 F=Σ(m=1..n)T^m の x^X y^Y z^Z の係数を求めるときには FをxでX回、yでY回、zでZ回偏微分してx=y=z=0 を入れた値を X!Y!Z! で割ればよい。 この計算を行うときには、Nを十分大きくとっておけば、 T=(1-x^N)/(1-x) (1-y^N)/(1-y) (1-z^N)/(1-z) -1 のかわりに T=1/(1-x)(1-y)(1-z) -1 として計算してもよいことがわかる。 S=(1-x)^(-1) (1-y)^(-1) (1-z)^(-1) とおくと、 F=Σ(m=1..n)(S-1)^m =(S-1)(1-(S-1)^n)/(1-(S-1)) これは、(S-1)^nを二項定理で展開した後、多項式の掛け算、割り算をすれば、Σ(m=1..∞)Am S^m の形に書ける。 (係数Am を m<=n の範囲のみ計算することにすれば O(n)で求まる) あとは、各 S^m についてx^X y^Y z^Z の係数を求めるには S^m=(1-x)^(-m) (1-y)^(-m) (1-z)^(-m) を xでX回、yでY回、zでZ回偏微分してx=y=z=0 を入れればよいが、 これは m(m+1)..(m+X-1) と m(m+1)..(m+Y-1) と m(m+1)..(m+Z-1) の積になる。 **********************************************************/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #pragma warning(disable:4996) typedef long long ll; #define MIN(a, b) ((a)>(b)? (b): (a)) #define MAX(a, b) ((a)<(b)? (b): (a)) #define LINF 9223300000000000000 #define INF 2140000000 const long long MOD = 1000000007; //const long long MOD = 998244353; using namespace std; const int max_comb=2500000; vector fac(max_comb+1); //n! (mod M) vector ifac(max_comb+1); //k!^(-1) (mod M) ll mpow(ll x, ll n){ //x^n(mod M) ll ans = 1; while(n != 0){ if(n&1) ans = ans*x % MOD; x = x*x % MOD; n = n >> 1; } return ans; } ll minv(ll x){ return mpow( x, MOD-2 ); } ll comb(int a, int b){ // C(a,b) = a! * b!^(-1) * (a-b)^(-1) if(a == 0 && b == 0)return 1; if(a < b || a < 0)return 0; ll tmp = ifac[a-b]* ifac[b] % MOD; return tmp * fac[a] % MOD; } ll perm(int a, int b){ // P(a,b) = a! * (a-b)!^(-1) if(b == 0)return 1; if(a < b || a < 0)return 0; ll tmp = ifac[a-b] % MOD; return tmp * fac[a] % MOD; } void pre_comb() { fac[0] = 1; ifac[0] = 1; for(int i = 0; i S(n+1); // S[m]は、(1-x)^(-m) (1-y)^(-m) (1-z)^(-m) を xでX回、yでY回、zでZ回偏微分してx=y=z=0 を入れた値 int i; for(i=1; i<=n; i++) { S[i]=1; if(X>0) S[i]=S[i]*perm(i+X-1,X)%MOD; if(Y>0) S[i]=S[i]*perm(i+Y-1,Y)%MOD; if(Z>0) S[i]=S[i]*perm(i+Z-1,Z)%MOD; } vector A(n+1); // A[i]は、(1-(S-1)^n) の S^i の係数 for(i=0; i<=n; i++) { if(i==0) A[0]=1; ll tmp=comb(n,i); if((n-i)%2==0) tmp=-tmp; A[i]=(A[i]+tmp+MOD)%MOD; } vector B(n+1); // B[i]は、(S-1)(1-(S-1)^n) の S^i の係数 (nより大きい次数は無視) for(i=0; i<=n; i++) { if(i>0) B[i]=A[i-1]; B[i]=(B[i]-A[i]+MOD)%MOD; } vector C(n+1); // C[i]は、(S-1)(1-(S-1)^n)/(1-(S-1)) の S^i の係数 (nより大きい次数は無視 ll curr=0; for(i=0; i<=n; i++) { curr=(curr+B[i])%MOD; C[i]=curr*minv(2)%MOD; curr=C[i]; } ll ans=0; for(i=1; i<=n; i++) { ll tmp=S[i]; // S[i]の値を X!Y!Z! で割る if(X>0) tmp=tmp*minv(perm(X,X))%MOD; if(Y>0) tmp=tmp*minv(perm(Y,Y))%MOD; if(Z>0) tmp=tmp*minv(perm(Z,Z))%MOD; ans=(ans+tmp*C[i]%MOD)%MOD; } printf("%lld\n", ans); return 0; }