リンク

http://arc012.contest.atcoder.jp/tasks/arc012_4

問題概要

N個の格子点(x_i,y_i)が与えられる。それぞれ時間1で上下左右に1ずつ動ける時、すべての点が時刻Tに一堂に会する動き方の組み合わせを求めよ。

制約

|x|,|y|<=10^6
T<=10^5
N<=10^5
1<=modulo<=1000000007

観察

  • 各点について独立に組み合わせを求めて全部かければよい。

部分点解法1

|x|,|y|<=10^2 のやつ
無し。
t=0,原点を1通りとしてスタートして、時間を増やしていってBFS的にDPする解法を想定としていたが、O(|x||y|T+N)かかり、明らかに間に合わない。どうしてこうなった\(^o^)/
なんか自分が送ったものが(制約含めて)ほとんどそのまま採用されていたみたいで完全に僕のミスですありがとうございました。

部分点解法2

modulo=1000000007のやつ
逆でもどうせ同じなので、(0,0)からスタートして(a,b)に時刻Tに着く方法の個数f(a,b)を求める。a>=0,b>=0としてよい。
T回のステップのうち、上下左右に進む回数をu,d,l,r回とおくと、
u+d+l+r=T, u-d=a, r-l=b
f(a,b)=\sum_{(u,d,l,r)} \frac{T!}{u!d!l!r!}. sumの(u,d,l,r)は上記の条件をすべて満たす。
d,lを消去して、
f(a,b)=\sum_{(u,r)} \frac{T!}{u!(u-a)!r!(r-b)!}.
u,rの満たす条件は、(2u-a)+(2r-b)=Tからu+r=(T+a+b)/2
これはrもuで表せることを意味する。r=(T+a+b)/2-uを代入して、
f(a,b)=\sum_u \frac{T!}{u!(u-a)!( (T+a+b)/2-u)!( (T+a+b)/2-u-b)!}.
1,3項目と2,4項目を足すとuが消えることに注目して、
f(a,b)=\left(\sum_u \binom{(T+a+b)/2}{u}\binom{(T-a-b)/2}{(T+a-b)/2-u}\right)\binom{T}{(T+a+b)/2}.
Vandermonde's identityから、
f(a,b)=\binom{(T+a+b)/2+(T-a-b)/2}{(T+a-b)/2}\binom{T}{(T+a+b)/2}
=\binom{T}{(T+a-b)/2}\binom{T}{(T+a+b)/2}.
T+a+bが2で割り切れない場合はf(a,b)=0となる。

あとはこれをN人分計算すれば良い。modulo=1000000007の場合、x!と(x!)^-1をあらかじめx<=10^5で求めておけば即座に計算できる。繰り返し数は10^5log 10^5+3N程度?

上記の複雑な計算をしなくても実験して推測できるんではないかと淡い期待を抱いていた。

満点解法

moduloの値が小さくなって、二項係数の計算がめんどくさくなる。割り算のところでmoduloの素因数が入っていると破綻するので、これを別口で計算する方針。
二項係数\binom{n}{r}~\text{mod}~moduloを亜光速で計算できるようにする。
moduloを素因数分解する。10^9+7以下なので、moduloを構成する相異なる素因数はたかだか9個である。これらをSとおく。Sの素因数以外で構成された数は、moduloに対して逆元が必ず存在するので、特に何も考えず混ぜてしまって良い。x!を列挙するときに、Sに属する素因数(p_i)の乗数(e_{x,i})と、x!からSに属する素因数を全て除いた積(x&#039;!)を保持しておく。後者に関しては逆元( (x&#039;!)^{-1})も、totient(modulo)-2を計算するか、拡張互除法を使った方法で求められるのでこれも保持しておく。
\binom{n}{r}を求めるときは、まず、n&#039;!(r&#039;!)^{-1}(n-r)&#039;^{-1}!を求めておく。次に、Sに属する乗数それぞれについて、\text{pow}(p_i, e_{n,i}-e_{r,i}-e_{n-r,i})~\text{mod}~moduloを求めて、かける。乗数が0未満になることはない。
pow(x,y)の計算にO(log y)かかるとしても、実行時間O( (Nlog T+T)log modulo)以下でいける。( (最初のk個の素数の積→k)のオーダーってどんなもんじゃろ)

所感

難しすぎましたね!
この問題は、自分がadvent calendarの二項係数の計算のところを見ていて、SegmentTreeの方法すげー無駄なことしているなーと思ったのがきっかけでできあがったものです。
前半の計算結果が割りと綺麗な式になるため、実は周知の事実ではないかと危惧してしまったのがいけなかったか・・。

コード

writer解。
+ ...
import java.io.ByteArrayInputStream;
import java.io.IOException;
import java.io.InputStream;
import java.io.PrintWriter;
import java.util.Arrays;
import java.util.InputMismatchException;
 
public class Main {
	static InputStream is;
	static PrintWriter out;
	static String INPUT = "";
 
	static void solve()
	{
		int n = ni(), t = ni(), mod = ni();
 
		int[] fs = enumFac(mod);
		int m = fs.length;
		int[] fact = new int[t+1];
		int[] ifact = new int[t+1];
		int[][] es = new int[t+1][m];
		fact[0] = ifact[0] = 1;
 
		// sieve
		for(int i = 1;i <= t;i++)fact[i] = i;
		for(int j = 0;j < m;j++){
			long p = fs[j];
			for(long q = p;q <= t;q *= p){
				for(long r = q;r <= t;r += q){
					fact[(int)r] /= p;
					es[(int)r][j]++;
				}
			}
		}
 
		for(int i = 1;i <= t;i++){
			fact[i] = (int)((long)fact[i-1] * fact[i] % mod);
			ifact[i] = (int)invl(fact[i], mod);
			for(int j = 0;j < m;j++)es[i][j] += es[i-1][j];
		}
 
		long ret = 1;
		int[] e = new int[m];
		for(int i = 0;i < n;i++){
			int a = Math.abs(ni()), b = Math.abs(ni());
			if(((t^a^b)&1)==1 || a+b > t){
				out.println(0);
				return;
			}
			// C(T, (T-a+b)/2), C(T, (T+a+b)/2)
			ret = ret
					*fact[t]%mod
					*ifact[(t-a+b)/2]%mod
					*ifact[(t+a-b)/2]%mod
					*fact[t]%mod
					*ifact[(t+a+b)/2]%mod
					*ifact[(t-a-b)/2]%mod;
			for(int j = 0;j < m;j++){
				e[j] += es[t][j]-es[(t-a+b)/2][j]-es[(t+a-b)/2][j];
				e[j] += es[t][j]-es[(t+a+b)/2][j]-es[(t-a-b)/2][j];
			}
		}
 
		for(int j = 0;j < m;j++){
			ret = ret * pow(fs[j], e[j], mod) % mod;
		}
 
		out.println(ret);
	}
 
	public static long pow(long a, long n, long mod)
	{
		long ret = 1;
		int x = 63-Long.numberOfLeadingZeros(n);
		for(;x >= 0;x--){
			ret = ret * ret % mod;
			if(n<<63-x<0)ret = ret * a % mod;
		}
		return ret;
	}
 
	public static long invl(long a, long mod)
	{
		long b = mod;
		long p = 1, q = 0;
		while(b > 0){
			long c = a / b;
			long d;
			d = a; a = b; b = d % b;
			d = p; p = q; q = d - c * q;
		}
		return p < 0 ? p + mod : p;
	}
 
	static int[] enumFac(int n)
	{
		int[] a = new int[10];
		int q = 0;
		for(int p = 2;p*p <= n;p++){
			if(n % p == 0)a[q++] = p;
			while(n%p == 0)n /= p;
		}
		if(n > 1)a[q++] = n;
		return Arrays.copyOf(a, q);
	}
 
	public static void main(String[] args) throws Exception
	{
		long S = System.currentTimeMillis();
		is = INPUT.isEmpty() ? System.in : new ByteArrayInputStream(INPUT.getBytes());
		out = new PrintWriter(System.out);
 
		solve();
		out.flush();
		long G = System.currentTimeMillis();
		tr(G-S+"ms");
	}
 
	private static byte[] inbuf = new byte[1024];
	static int lenbuf = 0, ptrbuf = 0;
 
	private static int readByte()
	{
		if(lenbuf == -1)throw new InputMismatchException();
		if(ptrbuf >= lenbuf){
			ptrbuf = 0;
			try { lenbuf = is.read(inbuf); } catch (IOException e) { throw new InputMismatchException(); }
			if(lenbuf <= 0)return -1;
		}
		return inbuf[ptrbuf++];
	}
 
	private static int ni()
	{
		int num = 0, b;
		boolean minus = false;
		while((b = readByte()) != -1 && !((b >= '0' && b <= '9') || b == '-'));
		if(b == '-'){
			minus = true;
			b = readByte();
		}
 
		while(true){
			if(b >= '0' && b <= '9'){
				num = num * 10 + (b - '0');
			}else{
				return minus ? -num : num;
			}
			b = readByte();
		}
	}
 
	private static void tr(Object... o) { if(INPUT.length() != 0)System.out.println(Arrays.deepToString(o)); }
}
 
最終更新:2013年02月11日 21:01