SRM613 Div1 Medium - RandomGCD

http://compro.tsutajiro.com/archive/181015_incexc.pdf
(包除原理 解ける数え上げの範囲を広げよう - tsutaj)
今の俺が求めていたものはこれだった
包除原理を極め鯛

問題

L以上r以下の整数を重複を許してn個選んで数列を作る
全要素のgcdがkになるような数列はいくつあるか

1 <= n, k, L <= 10^9
L <= r <= 10^9
0 <= r-L <= 10^5

解法

数列の要素として使えるのはL以上r以下のkの倍数のみ
そこで L = ceil(L/k), r = floor(r/k) とし、
「L以上r以下の整数のみからなる長さnの数列のうち、全要素のgcdが1になるようなものの個数」
を求めることにする

f(g): L以上r以下の整数のみからなる長さnの数列のうち、全要素のgcdがg以上になるようなものの個数
とすると、答えは包除原理の要領で
Σ f(g) * (-1)^(gの素因数の個数) (g: 同じ素因数で2度以上割り切れないような整数)
となる

素因数で2度以上割り切れないような整数 っていう条件はどういうことかというと
例えばg=4のケースは、g=1で足されてg=2で引かれてちょうど相殺されているから、足したり引いたりする必要がない
(g=12なんかも同様に、g=1で+、g=2で-、g=3で-、g=6で+ と既に消えている)

するとgの上限はどこまで考えればいいかが問題になる
1からRまで全部試せればそれでいいけど制約的に無理
そこで↑のΣで足し合わせるとき、L以上r以下の範囲に条件を満たす数が2個以上存在するようなもののみ考えることにする
(∵ 1個しか存在しない場合は当然同じ数字xがn個続く数列になり、gcdはxになるため、(x=1でない限りは)どうせ答えに影響しないから)
こうするとgは1からr-Lまで考えれば十分になる
(∵ L以上r以下の範囲にr-L+1以上の整数の倍数が2度以上現れることはないため)
L=1のときのみ答えに+1してやればok

実装

pf(n)はnの素因数一覧を重複なしで返すようになってる
だからこれの積がnと一致しなかったら素因数の重複があるってことになるから、
↑で言ったように足しも引きもしない

accumulateがあるなら積バージョンもあるんじゃね？ってずっと思ってたんだけど、
ググってみたらコードの通り
accumulate(a.begin(), a.end(), s, multiplies＜int＞())
とか書くとs * (指定範囲の各要素の積)を計算してくれるらしい
滅多に使わないけど覚えておこう

足し引きする値は、L以上r以下のうちiで割り切れる整数の個数を
m = floor(r/k) - floor((L-1)/k)
として、
(m個の整数から1個選ぶのをn回繰り返す) - (同じ整数をn回繰り返す数列の個数)
= n^m - m
となる

class RandomGCD {
public:
	vector<int> pf(int n) {
		vector<int> res;
		for (int i = 2; i*i <= n; i++) if (n%i==0) {
			res.push_back(i);
			while (n%i==0) n /= i;
		}
		if (n!=1) res.push_back(n);
		return res;
	}

	int countTuples(int n, int k, int l, int r) {
		l = (l+k-1)/k, r = r/k;
		if (r<l) return 0;
		ll res = (l==1);
		reps(i,1,r-l+1) {
			vector<int> p = pf(i);
			if (i!=1 && accumulate(all(p),1,multiplies<int>())!=i) continue;
			int m = r/i - (l-1)/i;
			(res += (mop(m,n)-m) * (i==1 || p.size()%2==0 ? 1 : -1)) %= mod;
		}
		(res += mod) %= mod;
		return (int)res;
	}
};

https://community.topcoder.com/stat?c=problem_statement&pm=13016&rd=15846

Code Festival 2018 予選B D - Sushi Restaurant

なんか本番中ずっと不適合度って言葉が気になってた
社会不適合度

本戦はtokumini君と寿司食ってきます✨

問題

略

解法

まずn人の空腹度が a[1],a[2],...,a[n] (昇順) で各皿に乗った寿司の個数が b[1],b[2],...,b[n] (昇順) である場合の不適合度は
Σ |a[i]-b[i]|
となる
このことからは「空腹度がi番目に少ない人」は「i番目に寿司が少ない皿」としか関係しないことがわかる

「空腹度がi番目に低い人」の空腹度がx[j]である確率をs(i,j)とする
またi人目が取る皿の寿司の個数をcとすると、この人の不満度の期待値は
Σ s(i,j) * |x[j]-c| (1 <= j <= m)
となる
これはよくあるやつで、cがx[K] (K: Σs(i,j)が0.5に最も近くなるi) のときに最小になる(中央値のやつ、証明略)
これを各iについて計算して足し合わせればO(nm)で答えがわかる
あとはs(i,j)を事前計算として求められればok

s(i,j)は単純に考えると、L <= i かつ i < r として
ΣΣ [(空腹度がx[j]未満の人がL-1人いる確率) * (空腹度がx[j]の人がr-L人いる確率)
* (空腹度がx[j]より大きい人がn-r+1人いる確率)] * n! / (L-1)! / (r-L)! / (n-r+1)! (並べ方の個数)
となる
Lとrを全通り試すとs(i,j)1個につきO(n^2)かかり、iがn通りでjがm通りあるからトータルでO(n^3 m)
ただ L <= i < r を満たす全てのiでs(i,j)は同じ値になるから、Lとrを全通り試しながら配列のs[L][j]からs[r-1][j]までいもす法(累積和)で足していくようにすると、iを全通り試す必要がなくなり、トータルO(n^2 m)になる
けどこれでもまだ間に合わない

空腹度がi番目に低い人の空腹度がx[j]以上である確率をt(i,j)とする
するとt(i,j)は、L <= iとして
Σ [(空腹度がx[j]未満の人がL-1人いる確率) * (空腹度がx[j]以上の人がn-L+1人いる確率)]
* n! / (L-1)! / (n-L+1)!
となる
これはさっきと同様にLを全通り試しつつL <= iの範囲に足していけばO(nm)で求められる
また定義より
s(i,j) = t(i,j) - t(i,j+1)
だからこれで解けた🤠

実装

コード中のsが↑で言うt

2000の階乗とかをdoubleでまともにやるとオーバーフローするから、途中はlogで計算して最後にexpする
logだから掛け算が足し算に、割り算が引き算に、べき乗が掛け算になってることに注意
あと0除算同様log(0)に気をつける必要がある
↓の書き方だとj=1のとき面倒なことになるから特別扱いしてやった
これ初めてやったけど便利だな

int n, m, q, x[2020], p[2020], ps[2020];
double f[2020], s[2020][2020], ans;

signed main() {
	reps(i,2,2020) f[i] = f[i-1] + log(i);
	cin >> n >> m >> q;
	reps(i,1,m+1) {
		cin >> x[i] >> p[i];
		ps[i] = p[i] + ps[i-1];
	}

	reps(i,1,n+1) {
		s[i][1] = (i==1 ? 1 : 0);
		reps(j,2,m+1) {
			s[i][j] += log(1.*ps[j-1]/q) * (i-1);
			s[i][j] += log(1.*(ps[m]-ps[j-1])/q) * (n-i+1);
			s[i][j] += f[n] - f[i-1] - f[n-i+1];
			s[i][j] = exp(s[i][j]);
		}
	}
	reps(i,1,n+1) reps(j,1,m+1) s[i][j] += s[i-1][j];

	reps(i,1,n+1) {
		double sum = 0;
		int id = -1;
		reps(j,1,m+1) {
			s[i][j] -= s[i][j+1];
			if (sum<0.5) id = j;
			sum += s[i][j];
		}
		reps(j,1,m+1) ans += s[i][j] * abs(x[j]-x[id]);
	}
	printf("%.14lf\n",ans);
}

模範解答だとloglistを事前に作ることで高速化とlog0避けを同時にしててスマート
条件を復元できる範囲でゆるくして変数を1個消すってテクニックは覚えておきたい

https://beta.atcoder.jp/contests/code-festival-2018-qualb/tasks/code_festival_2018_qualb_d

AGC028 B - Removing Blocks

問題

n個のブロックが一列に並んでおり、左からi番目のブロックの重さはa[i]である
これらのブロックに対して「まだ取り除かれていないブロックを1つ取り除き、そのブロックと連結なブロックの重さの総和をスコアに加える」という操作をn回行う
操作の列n!通り全てについて、スコアの合計を求め、その総和を求めよ

解法

p(i,j): ブロックiを取り除くときにブロックiとjが連結である確率
とする
ブロックiを取り除くときブロックiとjが連結であるということは、ブロックiとjの間にあるブロックの中で一番最初に取り除かれるブロックがブロックiであるということ
よって
p(i,j) = 1 / (abs(i-j)+1)
となる
あとは累積和をとっておいて
Σ a[i] * (Σs(i,j)) * n!
を求めればok

n!通り全ての和みたいな問題は確率を考えるってのが典型なのかな
何かの企業コンテストで同じ発想の問題を見た記憶がある

#define int long long

int n, a, f[100100], s[100100], res;

signed main() {
	f[0] = 1;
	reps(i,1,100100) f[i] = f[i-1]*i %mod;
	reps(i,1,100100) s[i] = (s[i-1] + mop(i,mod-2)) %mod;
	cin >> n;
	rep(i,n) {
		cin >> a;
		(res += (s[i+1] + s[n-i] - 1) * a) %= mod;
	}
	(res *= f[n]) %= mod;
	cout << res << endl;
}

鮮やかだなぁ・・

https://beta.atcoder.jp/contests/agc028/tasks/agc028_b

CF198 Div1 C - Iahub and Permutations

包除原理って良いなぁ😇

問題

Iahubは長さnの順列aを持っていたが、席を外している間に会社の同僚によってaのいくつかの要素を-1に書き換えられてしまった
書き換えられる前のaは、a[i]=iとなるような要素が1つもなかったことだけは覚えている
元のaとしてあり得るものの個数を求めよ

2 <= n <= 2000
aは-1になっている要素が少なくとも2つあり、条件を満たす順列が少なくとも1つは存在することが保証される

解法

-1の個数をk、消された数字xのうちa[x]が-1になっているようなものの個数をmとし、
f(i): a[x]=xとなるようなxがi個以上存在するような順列aの個数
とすると、包除原理より答えは
Σ f(i) * (-1)^i (0 <= i <= m)
となり、
f(i) = (m個からi個選ぶ選び方) * (残りk-i個を好きに当てはめる方法の個数)
= mCi * (k-i)!

#define int long long

int n, a[2020], b[2020], f[2020], m, k, res;

const int MAXS = 2020;
ll fact[MAXS+1], factr[MAXS+1], inv[MAXS+1];
ll comb(ll n, ll r) {
	if (fact[0]==0) {
		fact[0] = factr[0] = inv[1] = 1;
		for (int i=2; i<=MAXS; i++) inv[i] = inv[mod%i] * (mod-mod/i) %mod;
		for (int i=1; i<=MAXS; i++) fact[i] = fact[i-1]*i %mod, factr[i] = factr[i-1]*inv[i] %mod;
	}
	if (r<0 || n<r) return 0;
	return fact[n]*factr[r] %mod *factr[n-r] %mod;
}

signed main() {
	f[0] = 1;
	reps(i,1,2020) f[i] = f[i-1]*i %mod;
	cin >> n;
	rep(i,n) {
		cin >> a[i];
		if (a[i]!=-1) b[a[i]-1] = 1;
		else k++;
	}
	rep(i,n) if (!b[i] && a[i]==-1) m++;
	rep(i,m+1) (res += comb(m,i) * f[k-i] * (i%2 ? -1 : 1)) %= mod;
	(res += mod) %= mod;
	cout << res << endl;
}

これくらいの問題なら瞬殺できるようになってきた

http://codeforces.com/problemset/problem/340/E

yukicoder 391 - CODING WAR

包除原理は便利だなぁ！☺️

問題

n人の競技プログラマとm個の競技プログラミング問題がある
・各競技プログラマはどれかの問題を1問だけ解く
・各問題は必ず1人以上のプログラマによって解かれないといけない
・同じ問題を担当するプログラマが複数いてもいい
各プログラマに担当させる問題の組み合わせであって、これらの条件を全て満たすようなものの個数を求めよ

1 <= n <= 10^12
1 <= m <= 10^5

解法

f(k): 担当者がいない問題がk問以上あるような担当者の決め方の個数
とすると答えは包除原理より
Σ f(k) * (-1)^k (0 <= k <= m)
となる
f(k)はm問のうちどのk問を担当者不在にするかでmCk通りあり、その上でn人が残りm-k問のどれかを解くわけだから
f(k) = mCk * (m-k)^n

実装

mop = modpow

#define int long long

int n, m, res;

const int MAXS = 200200;
ll fact[MAXS+1], factr[MAXS+1], inv[MAXS+1];
ll comb(ll n, ll r) {
	if (fact[0]==0) {
		fact[0] = factr[0] = inv[1] = 1;
		for (int i=2; i<=MAXS; i++) inv[i] = inv[mod%i] * (mod-mod/i) %mod;
		for (int i=1; i<=MAXS; i++) fact[i] = fact[i-1]*i %mod, factr[i] = factr[i-1]*inv[i] %mod;
	}
	if (r<0 || n<r) return 0;
	return fact[n]*factr[r] %mod *factr[n-r] %mod;
}

signed main() {
	cin >> n >> m;
	rep(i,m+1) {
		int p = comb(m,i) * mop(m-i,n) %mod;
		(res += p * (i%2 ? -1 : 1)) %= mod;
	}
	(res += mod) %= mod;
	cout << res << endl;
}

https://yukicoder.me/problems/no/391

AOJ 2446 - Enumeration

スーパー包除原理マン

問題

n個の整数aとn個の整数pと整数mが与えられる
i番目の整数a[i]をp[i]%の確率で選ぶ、という操作を各iについて行い、0個以上n個以下の整数を選ぶ
1以上m以下の整数の中で、選ばれた整数の少なくとも1つで割り切れるものの個数の期待値を求めよ

1 <= n <= 20
1 <= m <= 10^18
1 <= a[i] <= 10^18
1 <= p[i] <= 99

解法

n個の整数aのある部分集合Sについて、
f(S): 1以上m以下の整数のうち、Sの全ての要素で割り切れるようなものの個数の期待値
とすると、包除原理より答えは
Σ f(S) * (-1)^(|S|+1)
となる (ここまでテンプレ)

Sの全ての要素で割り切れる数ってのは、Sの各要素のLCMで割り切れる数ってこと
Sを2^n通り全部試して、包除原理の符号に気を付けつつ
個数) 各要素のLCMでmを割った値
と
確率) 各要素を選ぶ確率の積
の積(=f(S))を足し合わせていけばいい

実装

LCMを求めるメソッドをちょっと改良して、オーバーフローする場合は2^60を返すようにした
今回はm <= 10^18 < 2^60だからこうすればオーバーフローしたときは答えが変動しないようになる

ll gcd(ll a,ll b) {return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b) {
	ll g = gcd(a,b);
	ll res = a/g*b;
	if (res/b != a/g) return 1ll<<60;
	return res;
}

ll n, m, a[22], x[22];
double res;

signed main() {
	cin >> n >> m;
	rep(i,n) cin >> a[i];
	rep(i,n) cin >> x[i];
	reps(b,1,1<<n) {
		double p = 1;
		ll lc = 1;
		rep(i,n) if (b>>i&1) {
			p *= 1.*x[i]/100;
			lc = lcm(lc,a[i]);
		}
		p *= m/lc;
		if (__builtin_parity(b)) res += p;
		else res -= p;
	}
	printf("%.14lf\n",res);
}

みんなのプロコン2018 決勝 A - Uncommon

包除原理マスターになりたい

問題

n個の相異なる整数aと整数mが与えられる
1以上m以下のそれぞれの整数iについて、aのうちiと互いに素であるものの個数を求めよ

1 <= n, m <= 10^5
1 <= a[i] <= 10^5

解法

xと互いに素な数 ってのは xの素因数のどれでも割り切れないような数 と言い換えられる
「xの相異なる素因数の集合」のある部分集合をSとし、
f(S): Sの全要素で割り切れるような数の個数
とすると、包除原理より答えは
Σ f(S) * (-1)^|S|
となる
f(S)はO(nlogn)で事前計算可能
10^5以下の自然数の相異なる素因数はたかだか6個しかないから、Sは全通り試しても十分高速に計算できる

int n, m, a[100100], d[100100];
set<int> s;

vector<int> pf(int n) {
	vector<int> res;
	for (int i = 2; i*i <= n; i++) if (n%i==0) {
		res.push_back(i);
		while (n%i==0) n /= i;
	}
	if (n!=1) res.push_back(n);
	return res;
}

signed main() {
	cin >> n >> m;
	rep(i,n) {
		cin >> a[i];
		s.insert(a[i]);
	}
	reps(i,1,m+1) {
		for(int j=i; j<100100; j+=i) if (s.count(j)) d[i]++;
		vector<int> p = pf(i);
		int res = 0;
		rep(j,1<<p.size()) {
			int num = 1;
			rep(k,p.size()) if (j>>k&1) num *= p[k];
			if (__builtin_parity(j)) res -= d[num];
			else res += d[num];
		}
		cout << res << endl;
	}
}