数学的な奴

abc130 c

```
int main() {
	ll w, h, x, y;
	cin >> w >> h >> x >> y;
	double ans = (double)(w)* h / 2;
	if (x*2==w && y*2==h) {
		cout << fixed << setprecision(10) << ans << " " << 1 << endl;
	}
	else {
		cout << fixed << setprecision(10) << ans << " " << 0 << endl;
	}
	return 0;
}
```
これのwについてるdoubleをちゃんと書かないとバグる。反省に書いておく


重複組み合わせ
abc 021-D 多重ループ

nHk=n+k-1Ck-1

```
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<map>
#include<math.h>
#include<queue>
#include<deque>
#include<stack>
#include<cstdio>
#include<utility>
#include<set>
#include<list>
#include<cmath>
#include<stdio.h>
#include<string.h>
#include<iomanip>
using namespace std;
#define FOR(i, a, b) for (ll i = (a); i <= (b); i++)
#define REP(i, n) FOR(i, 0, n - 1)
#define NREP(i, n) FOR(i, 1, n)
using ll = long long;
using pii = pair<int, int>;
using piii = pair<pii, pii>;
const int dx[4] = { 0, -1, 1, 0 };
const int dy[4] = { -1, 0, 0, 1 };
const int INF = 1e9 + 7;
int gcd(int x, int y) {
	if (x < y)swap(x, y);
	if (y == 0)return x;
	return gcd(y, x%y);
}

template<class T1, class T2> void chmin(T1 &a, T2 b) { if (a>b)a = b; }

template<class T1, class T2> void chmax(T1 &a, T2 b) { if (a<b)a = b; }
template<class T>
void Add(T &a, const T &b, const T &mod = 1000000007) {
	int val = ((a % mod) + (b % mod)) % mod;
	if (val < 0) { val += mod; }
	a = val;
}
////////////////////////////////////////
ll fac[510000], finv[510000], inv[510000];
void COMinit() {
	fac[0] = fac[1] = 1;
	finv[0] = finv[1] = 1;
	inv[1] = 1;
	for (int i = 2; i < 510000; ++i) {
		fac[i] = fac[i - 1] * i%INF;//i!のmod INF
		inv[i] = INF - inv[INF%i] * (INF / i) % INF;
		finv[i] = finv[i - 1] * inv[i] % INF;
	}
}

ll COM(int n, int k) {
	if (n < k)return 0;
	if (n < 0 || k < 0)return 0;
	return fac[n] * (finv[k] * finv[n - k] % INF) % INF;
}
int main() {
	int n, k;
	cin >> n >> k;
	COMinit();
	cout << COM(n + k - 1, k) << endl;
	return 0;
}
```

D
整数N,Mが与えられる。
a1*a2*.....aN=Mを満たす整数の組(a1,a2,,,,aN)が何通りあるか,1000000007で割った余りを求めよ。
まずは素数べきで考えるというのが１つの定石ではある。pを素数として
a1*a2*,,,,aN=pk
を満たす(a1,a2,,,,aN)が何通りあるかを考える。これは実はpの指数に着目してai=p*xiとしてx1+x2+,,,+xN=kを満たす0以上の整数が何通りあるかを数え上げる問題と考えることができる
これを各素因子ごとに独立に行って得られた結果を掛け算していけばよい

```
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<map>
#include<math.h>
#include<queue>
#include<deque>
#include<stack>
#include<cstdio>
#include<utility>
#include<set>
#include<list>
#include<cmath>
#include<stdio.h>
#include<cstdio>
using namespace std;
using ll = long long;
using pii = pair<int, int>;
const int INF = 1e9 + 7;
vector<pair<ll, ll>>prime_factorize(ll n) {
    vector<pair<ll, ll>>res;
    for (ll p = 2; p*p <= n; ++p) {
        if (n%p != 0)continue;
        int num = 0;
        while (n%p == 0) { ++num; n /= p; }
        res.push_back(make_pair(p, num));
    }
    if (n != 1)res.push_back(make_pair(n, 1));
    return res;
}
const int MAX = 210000;
ll fac[MAX], finv[MAX], inv[MAX];
void COMinit() {
    fac[0] = fac[1] = 1;
    finv[0] = finv[1] = 1;
    inv[1] = 1;
    for (int i = 2; i < MAX; ++i) {
        fac[i] = fac[i - 1] * i%INF;
        inv[i] = INF - inv[INF%i] * (INF / i) % INF;
        finv[i] = finv[i - 1] * inv[i] % INF;
    }
}
ll COM(int n, int k) {
    if (n < k)return 0;
    if (n < 0 || k < 0)return 0;
    return fac[n] * (finv[k] * finv[n - k] % INF) % INF;
}
int main() {
    int N, M;
    COMinit();
    cin >> N >> M;
    auto vec = prime_factorize(M);
    ll res = 1;
    for (auto pa : vec) {
        int num = pa.second;
        ll tmp = COM(num + N - 1, N - 1);
        res = (res*tmp) % INF;
    }
    cout << res << endl;
    return 0;
}
```
二項係数について

```
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<map>
#include<math.h>
#include<queue>
#include<deque>
#include<stack>
#include<cstdio>
#include<utility>
#include<set>
#include<list>
#include<cmath>
#include<stdio.h>
#include<cstdio>
using namespace std;
using ll = long long;
using pii = pair<int, int>;
const int INF = 1e9 + 7;
const int MAX = 510000;
ll fac[MAX], finv[MAX], inv[MAX];
//テーブルをつくる前処理
void COMinit() {
    fac[0] = fac[1] = 1;
    finv[0] = finv[1] = 1;
    inv[1] = 1;//mod pにおける1,2,,,nの逆元
    for (int i = 2; i < MAX; ++i) {
        fac[i] = fac[i - 1] * i%INF;
        inv[i] = INF - inv[INF%i] * (INF / i) % INF;
        finv[i] = finv[i - 1] * inv[i] % INF;
    }
}

ll COM(int n, int k) {
    if (n < k)return 0;
    if (n < 0 || k < 0)return 0;
    return fac[n] * (finv[k] * finv[n - k] % INF) % INF;
}

int main() {
    COMinit();//前処理
    cout << COM(100000, 50000) << endl;//計算例
}
```
使用可能場面1<=k<=n<=10*7
使用原理
nCk=(n!)(k!)^(-1)*((n-k)!)^(-1)

COMinit()でa!(fac[a])と(a!)^(-1)(finv[a])のテーブルを予め作っている。

Euclidの互除法による逆元計算
a^(-1)mod pを計算するとはすなわち
a*x+p*y=1
を満たすxを求めたいということになる。Euclidの互除法を利用する。すなわちpをaで割ってみる。
p=q*a+r
これを代入すると
a*x+(q*a+r)y=1<=>r*y+a(x+q*y)=1
になる。これによって(a,p)に関する問題がそれより数値の小さな(r,a)に関する問題に帰着できた。これを再帰的に解くのが拡張Euclidの互除法です。具体的には(r,a)に関する小問題を解いてr*s+a*t=1と解(s,t)が再帰的に得られたとすると、y=s,x+q*y=t<=>x=t-q*s,y=s
という風に元の問題の解を構成できます

最小公倍数
abc 070 c
a*b=gcd(a,b)*lcm(a,b)を利用する

```
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<map>
#include<math.h>
#include<queue>
#include<deque>
#include<stack>
#include<cstdio>
#include<utility>
#include<set>
#include<list>
#include<cmath>
#include<stdio.h>
#include<string.h>
#include<iomanip>
#include<cstdio>
#include<cstdlib>
#include<cstring>
using namespace std;
#define FOR(i, a, b) for (ll i = (a); i <= (b); i++)
#define REP(i, n) FOR(i, 0, n - 1)
#define NREP(i, n) FOR(i, 1, n)
using ll = long long;
using pii = pair<int, int>;
using piii = pair<pii, pii>;
const ll dx[4] = { 0, -1, 1, 0 };
const ll dy[4] = { -1, 0, 0, 1 };
const int INF = 1e9 + 7;
int gcd(int x, int y) {
	if (x < y)swap(x, y);
	if (y == 0)return x;
	return gcd(y, x%y);
}
void mul(ll a, ll b) {
	a = a * b % INF;
}
///////////////////////////////////////
ll Gcd(ll x, ll y) {
	if (x < y)swap(x, y);
	if (y == 0)return x;
	return Gcd(y, x%y);
}
ll lcm(ll a, ll b) {
	ll g = Gcd(a, b);
	return a / g * b;
}
int N;
ll T[110];
int main() {
    cin >> N;
	ll ans = 1;
	for (int i = 0; i < N; ++i) {
		ll T;
		cin >> T;
		ans = lcm(ans, T);
	}
	cout << ans << endl;
	return 0;
}
```
gcdもlcmも二つの数ごとに考える

3つ以上の整数の最小公倍数の求め方を考える
120=2*2*2*3*5,63=3*3*7,250=2*5*5*5ならそれぞれ素因数分解したときの指数の最大をとっていくことで求められるL=2*2*2*3*3*5*5*5*7=63000となる
agc 018 A

A1,,,ANすべての最大公約数をGとする。A1,,,,ANの最大値をMとする。このときkの書かれたボールを箱に入れることができる必要十分条件はKがgで割り切れるM以下の整数であることである。
M=max{Ai}として、k>Mのときは成り立たないことは自明
k<=Mの時を考える。
G=gcd{Ai}とする。MはAiのうちどれかなのでGの倍数。KがGの倍数ならばMからGを繰り返し引けばつくれる。逆にKがGの倍数でないならつくれない。Gの倍数同士の差は必ずGの倍数となるので。

```
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<map>
#include<math.h>
#include<queue>
#include<deque>
#include<stack>
#include<cstdio>
#include<utility>
#include<set>
#include<list>
#include<cmath>
#include<stdio.h>
#include<string.h>
#include<iomanip>
#include<cstdio>
#include<cstdlib>
#include<cstring>
using namespace std;
#define FOR(i, a, b) for (ll i = (a); i <= (b); i++)
#define REP(i, n) FOR(i, 0, n - 1)
#define NREP(i, n) FOR(i, 1, n)
using ll = long long;
using pii = pair<int, int>;
using piii = pair<pii, pii>;
const ll dx[4] = { 0, -1, 1, 0 };
const ll dy[4] = { -1, 0, 0, 1 };
const int INF = 1e9 + 7;
int gcd(int x, int y) {
	if (x < y)swap(x, y);
	if (y == 0)return x;
	return gcd(y, x%y);
}
void mul(ll a, ll b) {
	a = a * b % INF;
}
///////////////////////////////////////
ll Gcd(ll x, ll y) {
	if (x < y)swap(x, y);
	if (y == 0)return x;
	return Gcd(y, x%y);
}
ll lcm(ll a, ll b) {
	ll g = Gcd(a, b);
	return a / g * b;
}
int N, K;
int A[100100];
int main() {
	cin >> N >> K;
	int M = 0;
	for (int i = 0; i < N; ++i) {
		cin >> A[i];
		M = max(M, A[i]);
	}
	int ans = A[0];
	for (int i = 1; i < N; ++i) {
		ans = gcd(ans, A[i]);
	}
	if (K%ans == 0 && K <= M) {
		cout << "POSSIBLE" << endl;
	}
	else {
		cout << "IMPOSSIBLE" << endl;
	}
	return 0;
}
```
agc 001 B　mysterious lights

```
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<map>
#include<math.h>
#include<queue>
#include<deque>
#include<stack>
#include<cstdio>
#include<utility>
#include<set>
#include<list>
#include<cmath>
#include<stdio.h>
#include<string.h>
#include<iomanip>
#include<cstdio>
#include<cstdlib>
#include<cstring>
using namespace std;
#define FOR(i, a, b) for (ll i = (a); i <= (b); i++)
#define REP(i, n) FOR(i, 0, n - 1)
#define NREP(i, n) FOR(i, 1, n)
using ll = long long;
using pii = pair<int, int>;
using piii = pair<pii, pii>;
const ll dx[4] = { 0, -1, 1, 0 };
const ll dy[4] = { -1, 0, 0, 1 };
const int INF = 1e9 + 7;
int gcd(int x, int y) {
	if (x < y)swap(x, y);
	if (y == 0)return x;
	return gcd(y, x%y);
}
void mul(ll a, ll b) {
	a = a * b % INF;
}
///////////////////////////////////////

int main() {
	ll N, X; cin >> N >> X;
	ll res = N;
	ll a = max(N - X, X), b = min(N - X, X);
	while (b) {
		ll q = a / b;
		ll r = a % b;
		res += (b * 2)*q;
		if (r == 0) {
			res -= b;//最後だけ対角線のみになる
		}
		a = b, b = r;//次の平行四辺形は(b,r)
	}
	cout << res << endl;
	return 0;
}
```
いや頭いいな、、

JOI　2006本戦　最軽量のモビール

拡張ユークリッドの互除法
双六があります。整数が書かれておりこれは前後に無限に続いている。０がスタートで1がゴール。サイコロの目にはa,b,-a,-bしかない。4つの目をそれぞれ何回出せばゴールにたどり着けますか？複数解が存在する場合はどれか一つを出力せよ解がなければ-1を出力せよ。

この問題を式で表すとax+by=1となる整数x,yを求めよとなる。この方程式の解を求める関数をint extgcd(int a,int b,int& x,int& y)とし、返り値はgcd(a,b)とする。再帰的に定義すると
b*x'+(a%b)y'=gcd(a,b)の整数解x',y'が定まっているとして、このとき
a%b=a-(a/b)*bなので、これを代入すると
a*y'+b(x'-(a/b)*y')=gcd(a,b)となる。b=0のときは
a*1+b*0=a=gcd(a,b)です。


```
int extgcd(int a, int b, int& x, int& y) {
	int d = a;
	if (b != 0) {
		d = extgcd(b, a%b, y, x);
		y -= (a / b)*x;
	}
	else {
		x = 1; y = 0;
	}
	return d;
}
int main() {
	int x, y;
	cin >> x >> y;
	extgcd(111, 30, x, y);
	cout << x << "," << y << endl;//3,-11
	return 0;
}
```
素数判定
ARC017 A

```
bool is_prime(int n) {
	for (int i = 2; i*i <= n; ++i) {
		if (n%i == 0)return false;
	}
	return n != 1;
}
int main() {
	int N; cin >> N;
	if (is_prime(N)) {
		cout << "YES" << endl;
	}
	else {
		cout << "NO" << endl;
	}
	return 0;
}
```

```
//素因数分解
map<int, int>prime_factor(int n) {
	map<int, int>res;
	for (int i = 2; i*i <= n; ++i) {
		while (n%i == 0) {
			++res[i];
			n /= i;
		}
	}
	if (n != 1)res[n] = 1;
	return res;
}
```

```
//約数列挙
vector<int>divisor(int n) {
	vector<int>res;
	for (int i = 1; i*i <= n; ++i) {
		if (n%i == 0) {
			res.push_back(i);
			if (i != n / i) {
				res.push_back(n / i);
			}
		}
	}
	return res;
}
```

ABC 052 factors of　factorial

```
vector<int>key(1010, 0);
void prime(int n) {
	for (int i = 2; i*i <= n; ++i) {
		while (n%i == 0) {
			++key[i];
			n /= i;
		}
	}
	if (n != 1)key[n]++;
	return;
}
int main() {
	int n; cin >> n;
	ll ans = 1;
	for (int i = 2; i <= n; ++i) {
		prime(i);
	}
	for (int i = 2; i <= n; ++i) {
		ans *= (key[i] + 1);
		ans %= INF;
	}
	cout << ans << endl;
	return 0;
}
```
エラトステネスの篩
与えられた整数nに対してn以下の素数はいくつ存在しますか
n以下の素数を列挙するにはエラトステネスの篩と呼ばれるアルゴリズムがある。まず2以上n以下の整数をすべて書き出す。その中の最小の数である2は素数。２の倍数をすべて取り除く。残った最小の数である3はそれ以下の数で割り切れなかったので素数。表から3の倍数をすべて取り除く。、、、このような操作を繰り返す。

```
int prime[1001000];//i番目の素数
bool is_prime[1001010];//is_prime[i]がtrueならiは素数
//n以下の素数の数を返す
int sieve(int n) {
	int p = 0;
	for (int i = 0; i <= n; ++i) {
		is_prime[i] = true;
	}
	is_prime[0] = false;
	is_prime[1] = false;
	for (int i = 2; i <= n; ++i) {
		if (is_prime[i]) {
			prime[p] = i;
			p++;
			for (int j = 2 * i; j <= n; j += i) {
				is_prime[j] = false;
			}
		}
	}
	return p;
}
```
区間内の素数の個数
与えられた整数a,bに対して区間[a,b)に素数はいくつ存在しますか
b未満の素数でない整数の最小の素因数はたかだか√bなので√以下の素数の表があればエラトステネスの篩を[a,b)に用いることができる。つまり[2,√b)の表と[a,b)の表を別々に作っておいて[2,b)の表から素数が得られるたびに、その素数の倍数を[a,b)にある素数が列挙できる。

```
ll prime[1001024];//[a,b)の素数のうちi番目の素数
bool is_prime[1001024];//整数iが素数であるかどうか
bool is_prime_ab[1001024];//整数i+aが素数であるかどうか

ll segment_sieve(ll a, ll b) {
	fill(is_prime, is_prime + 1001024, true);
	fill(is_prime_ab, is_prime_ab + 1001024, true);
	for (ll i = 2; i*i < b; ++i) {
		if (!is_prime[i])continue;
		for (ll j = 2 * i; j*j < b; j += i)is_prime[j] = false;
		for (ll j = a - a % i; j < b; j += i) {
			if (j < a)continue;
			if (is_prime_ab[j - a])is_prime_ab[j - a] = false;
		}
	}
	ll res = 0;
	for (ll i = a; i < b; ++i) {
		if (is_prime_ab[i - a]) {
			prime[res++] = i;
		}
	}
	return res;
}
int main() {
	ll a, b; cin >> a >> b;
	ll res = segment_sieve(a,b);
	cout << res << endl;
	return 0;
}
```
ネットで調べた再帰の繰り返し2乗法

```
ll PowMod_RepeatSquaring(ll N, ll P, ll M) {
    if (P == 0)return 1;
    if (P % 2 == 0) {
        ll t = PowMod_RepeatSquaring(N, P / 2, M);
        return t * t%M;
    }
    return N * PowMod_RepeatSquaring(N, P - 1, M);
}
int main() {
    cout << PowMod_RepeatSquaring(2, 1000000000, INF) << endl;
    return 0;
}
```
カーマイケル数と繰り返し二乗法
任意の1<x<nに対して,x^n≡x(mod n)が成り立つような素数でない整数nをカーマイケル数といいます。与えられたnに対してそれがカーマイケル数かどうか判定しなさい

n=2^kのときはx^n=((x^2)^2)...となるので2乗をk回することで簡単に求まる。これを念頭に置くとnを2のべき乗の和で表すとよさそう。
n=2^k1+2^k2+2^k3...とすると
x^n=x^(2^k1)*x^(2^k2)*x^(2^k3)

```
bool is_prime(ll n) {
	for (ll i = 2; i*i <= n; ++i) {
		if (n%i == 0)return false;
	}
	return n != 1;
}

ll mod_pow(ll x, ll n, ll mod) {
	if (n == 0)return 1;
	ll res = mod_pow(x*x%mod, n / 2, mod);
	if (n & 1)res = res * x%mod;
	return res;
}

int main() {
    ll n; cin >> n;
	bool ok = true;
	if (is_prime(n)) {
		ok = false;
	}
	for (ll i = 2; i < n; ++i) {
		if (i%n != mod_pow(i, n, n)) {
			ok = false;
			break;
		}
	}
	if (ok)cout << "Yes" << endl;
	else cout << "No" << endl;
	return 0;
}
```
ゲーム系

caddi 2018D harlequin

```
int main() {
	int N; cin >> N;
	bool ok = true;
	REP(i, N) {
		int a; cin >> a;
		if (a % 2 != 0) {
			ok = false;
		}
	}
	if (ok)cout << "second" << endl;
	else cout << "first" << endl;
	
	//どの色のリンゴも偶数個という状態をE、いずれかの色のリンゴが奇数個という状態をOと呼ぶことにする
	//最後に自分が食べて全部０個にしてEの状態にすれば勝ち
	//現在の状態がEの場合、すでに全色０個でなければどのようにリンゴを食べても
	//食べた色のリンゴが奇数個になり必ずOに移る。一方で現在の状態がOの場合、
	//奇数個の色のりんごをだけを一個ずつ食べることでEに必ず移れる
	//よって初期状態がEの場合後手が先手の真似をすると先手はEに閉じ込められ
	//いずれ全色０個に至り負ける。
	//初期状態がOの場合は先手は最初に奇数個の色のリンゴだけ食べて後手にEを
	//押しつけそれ以降は後手の真似をすることで勝てる

	//状態で考える！！！！
	return 0;
}
```
3分探索
ムーアの法則
pが与えられた時、f(x)=x+p/2**(x/1.5)(0<=x<=100)を求めよ
これは三分探索で求められる
三分探索とは？→凸関数の極値を求めるのに用いる
探索領域(x0,x3)を三等分するx1,x2を選ぶ(x0,x1,x2,x3)でf(x1)とf(x2)を比べ、f(x)を比べ、f(x)が下に凸な関数なら値が大きいほう、上に凸な関数なら値が小さいほうの外側を捨てる。これを続ける

```cpp
double p;
const double EPS = 1e-10;
double f(double x, double p) {
	return x + (p / pow(2.0, x / 1.5));
}
int main() {
	cin >> p;
	double ans = 1e18;
	double left = 0.0;
	double right = 1000.0;
	while (abs(right - left) > EPS) {
		double x1 = left + (right - left) / 3.0;
		double x2 = left + 2.0*(right - left) / 3.0;
		double tmp1 = f(x1, p);
		double tmp2 = f(x2, p);
		ans = min(ans, tmp1);
		ans = min(ans, tmp2);
		if (tmp1 <= tmp2)right = x2;
		else left = x1;
	}
	cout << fixed << setprecision(10) << ans << endl;
	return 0;
}
```
約数包除
avc 162 E

```c++
//最大公約数がぴったりxであるための必要十分条件はxの倍数でありかつ2x,3x,,,,ではない

//x^y
ll mod_pow(ll x, ll n, ll mod) {
    if (n == 0)return 1;
    ll res = mod_pow(x*x%mod, n / 2, mod);
    if (n & 1)res = res * x%mod;
    return res;
}
int main() {
	int n, k; cin >> n >> k;
	ll res[100010];//key[x]で,x以上で
	for (int i = 0; i < 100010; ++i)res[i] = 0;
	ll ans = 0, cnt;//cntは最大公約数がiの数
	for (ll i = k; i >=1; --i) {
		cnt = mod_pow(k / i, n,(ll)INF);
		for (int j = i; j <= k; j+=i) {
			cnt -= res[j];
			if (cnt < 0)cnt += INF;
		}
		res[i] = cnt;
		ans += cnt * i;
		ans %= INF;
	}
	cout << ans << endl;
	return 0;
}
```
逆元をかけるのは
例えばAの逆元をINFで割った余りは
*=beki(A[i],INF-2)で表される

```c++
 for(int i = 0; i < N; i++) {
        ll now = LCM;
        now *= beki(A[i], mod - 2);
        now %= mod;
        ans += now;
        ans %= mod;
    }
```

```c++
ll powmod(ll a, ll k){
    ll ap=a, ans=1;
    while(k){
        if(k&1){
            ans*=ap;
            ans%=MOD;
        }
        ap=ap*ap;
        ap%=MOD;
        k>>=1;
    }
    return ans;
}
ll inv(ll a){
    return powmod(a, MOD-2);
}
```