SRM 410 DIV1 Hard WifiPlanet - antaの競技プログラミング練習日記

問題

整数座標x,yで単純多角形が与えられる。整数denomが与えられる。
全ての整数a,bに対して座標(denom*a,denom*b)を考える。この座標のうち与えられる多角形に含まれる座標すべてにwifiを置きたい(wifi点と呼ぼう)。wifi点はいくつあるか求めよ。

3 ≦ |x| ≦ 50
1 ≦ x[i], y[i] ≦ 10^9
2 ≦ denom ≦ 10^9
i≠jなら(x[i],y[i])≠(x[j],y[j])
多角形は単純。つまり、自己交差も端点以外の自己接点も持たない。
wifi点は与えられる多角形の辺・頂点上にはない。

解答

Editorialを参考にした。
多角形をどのように対処するか。これは符号付き面積のように考えられる。「辺が右向きか左向きかによる±×ある多角形の辺の下にあるwifi点の数」の総和を取り、最後に絶対値を取れば求められる。ここで、"辺の下にある"の"下"は、全ての辺で共通してさえいればどこでもよいので、y=0の直線とする。
線分を(x1,y1),(x2,y2)とし、x1＜x2とする。重複させないために、x2の座標にある点は数えないとする(問題の制約によりこれでも正しい)。
まず単純に考える。x座標をceil(x1/denom)*denomからfloor( (x2-1)/d)*denomまでdenomごとに考えたとき、その下にある点の数はfloor( (y1+(y2-y1)/(x2-x1)*(x-x1))/denom)なので、それを足せば良い。これは単に直線の式を見ているだけである。
しかし、これでは時間がかかりすぎる。とりえあず、ループ内の式を展開してみる。

xl = ceil(x1/denom), xr = floor((x2-1)/d) とする。
  Σ_{x=xl}^xr (floor((y1 + (y2-y1)/(x2-x1)*(x*denom-x1)) / denom))
= Σ_{x=xl}^xr (floor((y1*(x2-x1) + (y2-y1)*(x*denom-x1)) / ((x2-x1)*denom)))
= Σ_{x=0}^{xr-xl} (floor((y1*(x2-x1) + (y2-y1)*((x+xl)*denom-x1)) / ((x2-x1)*denom)))
= Σ_{x=0}^{xr-xl+1-1} (floor(((y2-y1)*denom*x + (y2-y1)*denom*xl + y1*(x2-x1) - (y2-y1)*x1) / ((x2-x1)*denom)))
ここで、適切に変数にletすれば、
= ��_{x=0}^{n-1} (floor(a*x + b) / m)
と表せる。

さて、それなりに整った形にはなったが、この総和をどのように計算できるのか。これにとても賢いアルゴリズムが存在する。ユークリッドの互除法のような解析でO(log m)のアルゴリズムだ。

f(n,a,b,m) = ��_{x=0}^{n-1} (floor(a*x + b) / m)
とする。m＞0とする。
まず簡単に計算できる、mの倍数の部分を別に計算する。a' = a mod m, a'' = floor(a/m), b' = b mod m, b'' = floor(b/m)とすると、
= ��_{x=0}^{n-1} (floor(a'*x + b') / m) + ��_{x=0}^{n-1} (a''*x + b'')
足し算の右辺は等差数列なので、公式を使ってΘ(1)の演算回数で計算できる。
  ��_{x=0}^{n-1} (a''*x + b'')
= a'' * (n * (n-1) / 2) + b'' * n
左辺。b'＜mなのでa'が0の場合は0となることがわかるので、a'は0でないとする。

さて、ここで階段状の図をイメージする。
今までの総和はこれを"横"に分割して1行ごとに数えていたが、今度はこれを"縦"に分割して1列ごとに考える。つまり、floor((a'*x+b')/m)がある整数y以上になる回数をそれぞれ数える。ここで、a'＞0であって階段が右上向きであるためにこのようにできることに注意。
    y ≦ floor((a'*x + b') / m)
<-> y ≦ (a'*x + b') / m
<-> (m * y - b') / a' ≦ x (* 0＜a'であることに注意 *)
また、x≦n-1であるので、この区間の整数の数は、
  = n - ceil((m * y - b') / a')
となる。また、yの上限は
    ceil((m * y - b') / a') ≦ n - 1
<-> (m * y - b') / a' ≦ n - 1
<-> y ≦ ((n-1) * a' + b') / m
<-> y ≦ floor(((n-1) * a' + b') / m)
yの下限が1であることに注意してまとめる。
t = a' - b' - 1, u = floor(((n-1) * a' + b') / m) + 1 とすると、
  ��_{x=0}^{n-1} (floor(a'*x + b') / m)
= ��_{y=1}^{u-1} (n - floor((m * y + t) / a')))
= n * (u-1) - ��_{y=1}^{u-1} floor((m * y + t) / a'))
= n * (u-1) - ��_{y=0}^{u-1} floor((m * y + t) / a')) + floor(t / a')
この総和の部分は自分自身の関数の再帰で計算できる。

次にこのアルゴリズムの計算量を考える。

再帰以外の部分はΘ(1)の演算で計算できる。
引数の中でaとmに注目する。それ以外の引数は終了条件にも関わらず、aにもmにも影響しないので無視して良い。
T(f(a,m)) = Θ(1) +
            if a mod m == 0
            then Θ(1)
            else Θ(1) + T(f(m,a mod m))
これは実際にユークリッドの互除法とほとんど同じ形をしている。違うのは終了条件だけだが、ユークリッドの互除法で"a mod m == 0"のとき次のmが0になってΘ(1)で終了することを考えると同じ形になる。
ユークリッドの互除法はO(log b)なので、このアルゴリズムもO(log m)の計算量(整数演算回数)となる。

ところで、xが負かもしれない時にfloor(x/y)や(x mod y)を計算するには？
(x mod y)は"x - floor(x/y)*y"とすればよい。
floor(x/y)は

return x / y - (x % y < 0);

でよい。ただし、yが負かもしれない場合はもう少しだけコードを足す必要がある。

自分で式を立ててやってみたらEditorialに書かれているものと違う形になってしまった。Editorialのもののほうが少しきれいな気もする。
ユークリッドの互除法風のO(log n)のアルゴリズムは一般に頭いいと思う。なんというか、一般に、x軸とy軸を交互に交換してく感じになるのかな。
このアルゴリズムは他にも使えそう。様々な問題が作れそう。それでいて比較的あまり知られていないような？どうだろう？もしそうなら、チャンスだな。

コード

#include <vector>
#include <algorithm>
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
using namespace std;
typedef long long ll;   

template<typename T, typename U>
inline auto floordiv(T x, U y) -> decltype(x/y) {
	auto q = x / y, r = x % y;
	return q - ((r!=0) & ((r<0) ^ (y<0)));
}

//\Sigma_{x=0}^{n-1} floor((a * x + b) / m)
ll linear_floordiv_sum(ll n, ll a, ll b, ll m) {
	if(m < 0) return linear_floordiv_sum(n, -a, -b, -m);
	ll ad = floordiv(a, m), bd = floordiv(b, m);
	ll am = a - ad * m, bm = b - bd * m;
	ll r = ad * (n * (n-1) / 2) + bd * n;
	if(am == 0) return r;
	ll u = floordiv((n-1) * am + bm, m) + 1;
	ll t = am - bm - 1;
	r += n * (u-1);
	r -= linear_floordiv_sum(u, m, t, am);
	r += floordiv(t, am);
	return r;
}

struct WifiPlanet {
	//[x1,x2)
	ll count(int d, int x1, int x2, int y1, int y2) {
		if(x1 == x2) return 0;
		if(x1 > x2) return - count(d, x2, x1, y2, y1);
		ll r = 0;
		int xl = (x1+d-1)/d, xr = (x2-1)/d;
		int n = xr - xl + 1;
		ll a = (ll)(y2 - y1) * d;
		ll b = (ll)x2*y1 - (ll)x1*y2 + (ll)xl*d*y2 - (ll)xl*d*y1;
		ll m = (ll)(x2-x1) * d;
		r = linear_floordiv_sum(n, a, b, m);
//		r += floordiv((ll)y1 * (x2-x1) + (ll)(y2-y1) * (x-x1), (x2-x1) * d);	//(y1+(y2-y1)/(x2-x1)*(x-x1))/denom
		return r;
	}
	long long routersNeeded(vector <int> x, vector <int> y, int denom) {
		int n = x.size();
		ll r = 0;
		rep(i, n)
			r += count(denom, x[i], x[(i+1)%n], y[i], y[(i+1)%n]);
		return abs(r);
	}
};

問題

解答

コメント

コード