# 数学
### 矩阵运算
```cpp
struct Mat {
static const LL M = 2;
LL v[M][M];
Mat() { memset(v, 0, sizeof v); }
void eye() { FOR (i, 0, M) v[i][i] = 1; }
LL* operator [] (LL x) { return v[x]; }
const LL* operator [] (LL x) const { return v[x]; }
Mat operator * (const Mat& B) {
const Mat& A = *this;
Mat ret;
FOR (k, 0, M)
FOR (i, 0, M) if (A[i][k])
FOR (j, 0, M)
ret[i][j] = (ret[i][j] + A[i][k] * B[k][j]) % MOD;
return ret;
}
Mat pow(LL n) const {
Mat A = *this, ret; ret.eye();
for (; n; n >>= 1, A = A * A)
if (n & 1) ret = ret * A;
return ret;
}
Mat operator + (const Mat& B) {
const Mat& A = *this;
Mat ret;
FOR (i, 0, M)
FOR (j, 0, M)
ret[i][j] = (A[i][j] + B[i][j]) % MOD;
return ret;
}
void prt() const {
FOR (i, 0, M)
FOR (j, 0, M)
printf("%lld%c", (*this)[i][j], j == M - 1 ? '\n' : ' ');
}
};
```
### 分数类
```cpp
struct Frac {
LL a, b;
Frac(LL _a=0, LL _b=0) {
a = _a;
b = _b;
LL g = __gcd(a, b);
a /= g;
b /= g;
if (b < 0) {
a = -a;
b = -b;
}
}
Frac operator + (const Frac &other) {
LL newa = a*other.b + b*other.a,
newb = b*other.b;
return Frac(newa, newb);
}
Frac operator * (const Frac &other) {
return Frac(a*other.a, b*other.b);
}
};
```
## 筛
* 线性筛
```cpp
const LL p_max = 1E6 + 100;
LL pr[p_max], p_sz;
void get_prime() {
static bool vis[p_max];
FOR (i, 2, p_max) {
if (!vis[i]) pr[p_sz++] = i;
FOR (j, 0, p_sz) {
if (pr[j] * i >= p_max) break;
vis[pr[j] * i] = 1;
if (i % pr[j] == 0) break;
}
}
}
```
* 线性筛+欧拉函数
```cpp
const LL p_max = 1E5 + 100;
LL phi[p_max];
void get_phi() {
phi[1] = 1;
static bool vis[p_max];
static LL prime[p_max], p_sz, d;
FOR (i, 2, p_max) {
if (!vis[i]) {
prime[p_sz++] = i;
phi[i] = i - 1;
}
for (LL j = 0; j < p_sz && (d = i * prime[j]) < p_max; ++j) {
vis[d] = 1;
if (i % prime[j] == 0) {
phi[d] = phi[i] * prime[j];
break;
}
else phi[d] = phi[i] * (prime[j] - 1);
}
}
}
```
* 线性筛+莫比乌斯函数
```cpp
const LL p_max = 1E5 + 100;
LL mu[p_max];
void get_mu() {
mu[1] = 1;
static bool vis[p_max];
static LL prime[p_max], p_sz, d;
mu[1] = 1;
FOR (i, 2, p_max) {
if (!vis[i]) {
prime[p_sz++] = i;
mu[i] = -1;
}
for (LL j = 0; j < p_sz && (d = i * prime[j]) < p_max; ++j) {
vis[d] = 1;
if (i % prime[j] == 0) {
mu[d] = 0;
break;
}
else mu[d] = -mu[i];
}
}
}
```
## 亚线性筛
### min_25
+ zerol 版本, $f(p^c) = p \ xor \ c$
```cpp
namespace min25 {
const int M = 1E6 + 100;
LL B, N;
// g(x)
inline LL pg(LL x) { return 1; }
inline LL ph(LL x) { return x % MOD; }
// Sum[g(i),{x,2,x}]
inline LL psg(LL x) { return x % MOD - 1; }
inline LL psh(LL x) {
static LL inv2 = (MOD + 1) / 2;
x = x % MOD;
return x * (x + 1) % MOD * inv2 % MOD - 1;
}
// f(pp=p^k)
inline LL fpk(LL p, LL e, LL pp) { return p ^ e; }
// f(p) = fgh(g(p), h(p))
inline LL fgh(LL g, LL h) { return h - g; }
LL pc, sg[M], sh[M];
LL w[M];
LL id1[M], id2[M], h[M], g[M];
inline LL id(LL x) { return x <= B ? id1[x] : id2[N / x]; }
LL go(LL x, LL k) {
if (x <= 1 || (k >= 0 && prime[k] > x)) return 0;
LL t = id(x);
LL ans = fgh((g[t] - sg[k + 1]), (h[t] - sh[k + 1]));
FOR (i, k + 1, pc) {
LL p = prime[i];
if (p * p > x) break;
ans -= fgh(pg(p), ph(p));
for (LL pp = p, e = 1; pp <= x; ++e, pp = pp * p)
ans += fpk(p, e, pp) * (1 + go(x / pp, i)) % MOD;
}
return ans % MOD;
}
LL solve(LL _N) {
N = _N;
B = sqrt(N + 0.5);
pc = upper_bound(prime, prime+pcnt, B) - prime;
for (int i=1; i<=pc; ++i) {
sg[i] = (sg[i-1] + pg(prime[i-1])) % MOD;
sh[i] = (sh[i-1] + ph(prime[i-1])) % MOD;
}
int sz = 0;
for (LL l = 1, v, r; l <= N; l = r + 1) {
v = N / l; r = N / v;
w[sz] = v; g[sz] = psg(v); h[sz] = psh(v);
if (v <= B) id1[v] = sz; else id2[r] = sz;
sz++;
}
FOR (k, 0, pc) {
LL p = prime[k];
FOR (i, 0, sz) {
LL v = w[i]; if (p * p > v) break;
LL t = id(v / p);
g[i] = (g[i] - (g[t] - sg[k]) * pg(p)) % MOD;
h[i] = (h[i] - (h[t] - sh[k]) * ph(p)) % MOD;
}
}
return (go(N, -1) % MOD + MOD + 1) % MOD;
}
}
```
+ Weaver_zhu 版本，质因数个数
```cpp
namespace min25 {
const LL N = 1e6+5;
LL id1[N], id2[N], w[N];
LL n, B, m, pc;
inline LL id(LL x) { return x > B? id2[n/x]:id1[x];}
inline LL fp(LL x) { return 1;}
inline LL sp(LL x) { return x%MOD-1; }
inline LL fpk(LL p, LL k, LL pp) { return 0;}
LL g[N], sg[N];
LL S(LL n, int j) {
LL ans = (g[id(n)] - sg[j] + MOD) % MOD;
for (int k=j; k
LL p = prime[k]; if (p * p > n) break;
for (LL pp=p, e=1; pp*p<=n; ++e, pp*=p) {
LL del = S(n / pp, k+1);
ans = (ans + fpk(p, e, pp) * del + fpk(p, e+1, pp*p)) % MOD;
}
}
return ans;
}
LL solve(LL _n) {
n = _n;
B = sqrt(n + 0.5) + 1;
m = 0;
for (LL l=1, r, v; l<=n; l=r+1) {
v = n / l; r = n / v;
if (v <= B) id1[v] = m;
else id2[r] = m;
g[m] = sp(v);
w[m++] = v;
}
pc = upper_bound(prime, prime+pcnt, B) - prime;
for (int i=1; i<=pc; ++i) {
sg[i] = (sg[i-1] + fp(prime[i-1])) % MOD;
}
for (int j=0; j
LL p = prime[j];
for (int i=0; i
if (p * p > w[i]) break;
g[i] = (g[i] - fp(p) * (g[id(w[i]/p)] % MOD - sg[j])) % MOD + MOD;
}
}
return S(n, 0);
}
}
```
### min25 cheat sheets
#### 素数项等幂求和
+ 前置模板：等幂求和
+ 调用`init(n, k)`初始化
+ 调用`calc(n)`使用
```cpp
namespace min25 {
const int M = 1e6+100;
LL B, N;
vector c; LL k;
inline LL ph(LL x) { return bin(x, k, MOD);}
inline LL psh(LL x) {
x %= MOD;
LL res = 0;
int cn = c.size();
LL xx = 1;
for (int i=0; i
res = (res + c[i] * xx) % MOD;
xx = xx * x % MOD;
}
return res-1;
}
inline LL fpk(LL p, LL e, LL pp) {
return bin(pp, k, MOD);
}
inline LL fgh(LL g, LL h) { return h - g;}
LL pr[M], pc;
LL sh[M];
void get_prime(LL n) {
static bool vis[M];
memset(vis, 0, sizeof(bool) * (n+5));
pc = 0;
sh[0] = 0;
for (int i=2; i
if (!vis[i]) {
pr[pc++] = i;
sh[pc] = (sh[pc-1] + ph(i)) % MOD;
}
for (int j=0; j
if (1LL * pr[j] * i > n) break;
vis[pr[j] * i] = 1;
if (i % pr[j] == 0) break;
}
}
}
LL w[M];
LL id1[M], id2[M], h[M];
inline LL id(LL x) { return x <= B? id1[x] : id2[N / x];}
void init(LL n, LL _k) {
k = _k;
c = sumk::go(k);
B = sqrt(n);
N = n;
get_prime(B);
}
void init(int _k) {
k = _k;
c = sumk::go(k);
}
LL go(LL _N) {
N = _N;
int sz = 0;
for (LL l=1, v, r; l<=N; l=r+1) {
v = N / l;
r = N / v;
w[sz] = v;
h[sz] = psh(v);
if (v <= B) id1[v] = sz; else id2[r] = sz;
sz++;
}
for (int k=0; k
LL p = pr[k];
for (int i=0; i
LL v = w[i]; if (p * p > v) break;
LL t = id(v / p);
h[i] = (h[i] - (h[t] - sh[k]) * ph(p)) % MOD;
h[i] = (h[i] + MOD) % MOD;
}
}
return h[id(N)-1];
}
inline LL calc(LL x) { if (x == 0) return 0; return h[id(x)];}
}
```
### 杜教筛
求 $S(n)=\sum_{i=1}^n f(i)$，其中 $f$ 是一个积性函数。
构造一个积性函数 $g$，那么由 $(f*g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})$，得到 $f(n)=(f*g)(n)-\sum_{d|n,d
当然，要能够由此计算 $S(n)$，会对 $f,g$ 提出一些要求：
+ $f*g$ 要能够快速求前缀和。
+ $g$ 要能够快速求分段和(前缀和)。
+ 对于正常的积性函数 $g(1)=1$，所以不会有什么问题。
在预处理 $S(n)$ 前 $n^{\frac{2}{3}}$ 项的情况下复杂度是 $O(n^{\frac{2}{3}})$。
```cpp
namespace dujiao {
const int M = 5E6;
LL f[M] = {0, 1};
void init() {
static bool vis[M];
static LL pr[M], p_sz, d;
FOR (i, 2, M) {
if (!vis[i]) { pr[p_sz++] = i; f[i] = -1; }
FOR (j, 0, p_sz) {
if ((d = pr[j] * i) >= M) break;
vis[d] = 1;
if (i % pr[j] == 0) {
f[d] = 0;
break;
} else f[d] = -f[i];
}
}
FOR (i, 2, M) f[i] += f[i - 1];
}
inline LL s_fg(LL n) { return 1; }
inline LL s_g(LL n) { return n; }
LL N, rd[M];
bool vis[M];
LL go(LL n) {
if (n < M) return f[n];
LL id = N / n;
if (vis[id]) return rd[id];
vis[id] = true;
LL& ret = rd[id] = s_fg(n);
for (LL l = 2, v, r; l <= n; l = r + 1) {
v = n / l; r = n / v;
ret -= (s_g(r) - s_g(l - 1)) * go(v);
}
return ret;
}
LL solve(LL n) {
N = n;
memset(vis, 0, sizeof vis);
return go(n);
}
}
```
#### 欧拉函数前缀和
```cpp
namespace bigsum_phi {
bool noprime[MAXN];
int phi[MAXN];
int prime[MAXN], pcnt;
LL sumphi[MAXN];
unordered_map mp;
void getp() {
phi[1] = 1;
FOR (i, 2, MAXN) {
if (!noprime[i]) {
prime[pcnt++] = i;
phi[i] = i-1;
}
FOR (j, 0, pcnt) {
LL nextp = i * prime[j];
if (nextp >= MAXN) break;
noprime[nextp] = 1;
if (i % prime[j] == 0) {
phi[nextp] = phi[i] * prime[j];
} else {
phi[nextp] = phi[i] * phi[prime[j]];
}
}
}
FOR (i, 1, MAXN) {
sumphi[i] = (sumphi[i-1] + phi[i])%MOD;
}
}
LL bigsumphi(LL n) {
if (n < MAXN) return sumphi[n];
if (mp.count(n)) return mp[n];
LL ans = 0;
for (LL l=2, v, r; l<=n; l=r+1) {
v = n/l; r = n/v;
LL del = (r-l+1) * bigsumphi(v) % MOD;
ans = (ans + del) % MOD;
}
LL A = n*(n+1)/2 % MOD;
ans = (A+MOD-ans)%MOD;
mp[n] = ans;
return ans;
}
}
```
## 素数测试
+ 前置： 快速乘、快速幂
+ int 范围内只需检查 2, 7, 61
+ long long 范围 2, 325, 9375, 28178, 450775, 9780504, 1795265022
+ 3E15内 2, 2570940, 880937, 610386380, 4130785767
+ 4E13内 2, 2570940, 211991001, 3749873356
+ http://miller-rabin.appspot.com/
```cpp
bool checkQ(LL a, LL n) {
if (n == 2 || a >= n) return 1;
if (n == 1 || !(n & 1)) return 0;
LL d = n - 1;
while (!(d & 1)) d >>= 1;
LL t = bin(a, d, n); // 不一定需要快速乘
while (d != n - 1 && t != 1 && t != n - 1) {
t = mul(t, t, n);
d <<= 1;
}
return t == n - 1 || d & 1;
}
bool primeQ(LL n) {
static vector t = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
if (n <= 1) return false;
for (LL k: t) if (!checkQ(k, n)) return false;
return true;
}
```
## Pollard-Rho
+ 大整数分解
```cpp
mt19937 mt(time(0));
LL pollard_rho(LL n, LL c) {
LL x = uniform_int_distribution(1, n - 1)(mt), y = x;
auto f = [&](LL v) { LL t = mul(v, v, n) + c; return t < n ? t : t - n; };
while (1) {
x = f(x); y = f(f(y));
if (x == y) return n;
LL d = gcd(abs(x - y), n);
if (d != 1) return d;
}
}
LL fac[100], fcnt; // 结果
void get_fac(LL n, LL cc = 19260817) {
if (n == 4) { fac[fcnt++] = 2; fac[fcnt++] = 2; return; }
if (primeQ(n)) { fac[fcnt++] = n; return; }
LL p = n;
while (p == n) p = pollard_rho(n, --cc);
get_fac(p); get_fac(n / p);
}
```
## BM 线性递推
```cpp
namespace BerlekampMassey {
inline void up(LL& a, LL b) { (a += b) %= MOD; }
V mul(const V&a, const V& b, const V& m, int k) {
V r; r.resize(2 * k - 1);
FOR (i, 0, k) FOR (j, 0, k) up(r[i + j], a[i] * b[j]);
FORD (i, k - 2, -1) {
FOR (j, 0, k) up(r[i + j], r[i + k] * m[j]);
r.pop_back();
}
return r;
}
V pow(LL n, const V& m) {
int k = (int) m.size() - 1; assert (m[k] == -1 || m[k] == MOD - 1);
V r(k), x(k); r[0] = x[1] = 1;
for (; n; n >>= 1, x = mul(x, x, m, k))
if (n & 1) r = mul(x, r, m, k);
return r;
}
LL go(const V& a, const V& x, LL n) {
// a: (-1, a1, a2, ..., ak).reverse
// x: x1, x2, ..., xk
// x[n] = sum[a[i]*x[n-i],{i,1,k}]
int k = (int) a.size() - 1;
if (n <= k) return x[n - 1];
if (a.size() == 2) return x[0] * bin(a[0], n - 1, MOD) % MOD;
V r = pow(n - 1, a);
LL ans = 0;
FOR (i, 0, k) up(ans, r[i] * x[i]);
return (ans + MOD) % MOD;
}
V BM(const V& x) {
V a = {-1}, b = {233}, t;
FOR (i, 1, x.size()) {
b.push_back(0);
LL d = 0, la = a.size(), lb = b.size();
FOR (j, 0, la) up(d, a[j] * x[i - la + 1 + j]);
if (d == 0) continue;
t.clear(); for (auto& v: b) t.push_back(d * v % MOD);
FOR (_, 0, la - lb) t.push_back(0);
lb = max(la, lb);
FOR (j, 0, la) up(t[lb - 1 - j], a[la - 1 - j]);
if (lb > la) {
b.swap(a);
LL inv = -get_inv(d, MOD);
for (auto& v: b) v = v * inv % MOD;
}
a.swap(t);
}
for (auto& v: a) up(v, MOD);
return a;
}
}
```
+ 杜教版
```cpp
namespace linear_seq {
const int N=10010;
ll res[N],base[N],_c[N],_md[N];
vector Md;
void mul(ll *a,ll *b,int k) {
rep(i,0,k+k) _c[i]=0;
rep(i,0,k) if (a[i]) rep(j,0,k) _c[i+j]=(_c[i+j]+a[i]*b[j])%mod;
for (int i=k+k-1;i>=k;i--) if (_c[i])
rep(j,0,SZ(Md)) _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod;
rep(i,0,k) a[i]=_c[i];
}
int solve(ll n,VI a,VI b) { // a 系数 b 初值 b[n+1]=a[0]*b[n]+...
ll ans=0,pnt=0;
int k=SZ(a);
assert(SZ(a)==SZ(b));
rep(i,0,k) _md[k-1-i]=-a[i];_md[k]=1;
Md.clear();
rep(i,0,k) if (_md[i]!=0) Md.push_back(i);
rep(i,0,k) res[i]=base[i]=0;
res[0]=1;
while ((1ll<
for (int p=pnt;p>=0;p--) {
mul(res,res,k);
if ((n>>p)&1) {
for (int i=k-1;i>=0;i--) res[i+1]=res[i];res[0]=0;
rep(j,0,SZ(Md)) res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod;
}
}
rep(i,0,k) ans=(ans+res[i]*b[i])%mod;
if (ans<0) ans+=mod;
return ans;
}
VI BM(VI s) {
VI C(1,1),B(1,1);
int L=0,m=1,b=1;
rep(n,0,SZ(s)) {
ll d=0;
rep(i,0,L+1) d=(d+(ll)C[i]*s[n-i])%mod;
if (d==0) ++m;
else if (2*L<=n) {
VI T=C;
ll c=mod-d*powmod(b,mod-2)%mod;
while (SZ(C)
rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod;
L=n+1-L; B=T; b=d; m=1;
} else {
ll c=mod-d*powmod(b,mod-2)%mod;
while (SZ(C)
rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod;
++m;
}
}
return C;
}
int gao(VI a,ll n) {
VI c=BM(a);
c.erase(c.begin());
rep(i,0,SZ(c)) c[i]=(mod-c[i])%mod;
return solve(n,c,VI(a.begin(),a.begin()+SZ(c)));
}
};
```
## 扩展欧几里得
* 求 $ax+by=gcd(a,b)$ 的一组解
* 如果 $a$ 和 $b$ 互素，那么 $x$ 是 $a$ 在模 $b$ 下的逆元
* 注意 $x$ 和 $y$ 可能是负数
```cpp
LL ex_gcd(LL a, LL b, LL &x, LL &y) {
if (b == 0) { x = 1; y = 0; return a; }
LL ret = ex_gcd(b, a % b, y, x);
y -= a / b * x;
return ret;
}
```
+ 卡常欧几里得
```cpp
inline int ctz(LL x) { return __builtin_ctzll(x); }
LL gcd(LL a, LL b) {
if (!a) return b; if (!b) return a;
int t = ctz(a | b);
a >>= ctz(a);
do {
b >>= ctz(b);
if (a > b) swap(a, b);
b -= a;
} while (b);
return a << t;
}
```
## 类欧几里得
* $m = \lfloor \frac{an+b}{c} \rfloor$.
* $f(a,b,c,n)=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor$: 当 $a \ge c$ or $b \ge c$ 时，$f(a,b,c,n)=(\frac{a}{c})n(n+1)/2+(\frac{b}{c})(n+1)+f(a \bmod c,b \bmod c,c,n)$；否则 $f(a,b,c,n)=nm-f(c,c-b-1,a,m-1)$。
* $g(a,b,c,n)=\sum_{i=0}^n i \lfloor\frac{ai+b}{c}\rfloor$: 当 $a \ge c$ or $b \ge c$ 时，$g(a,b,c,n)=(\frac{a}{c})n(n+1)(2n+1)/6+(\frac{b}{c})n(n+1)/2+g(a \bmod c,b \bmod c,c,n)$；否则 $g(a,b,c,n)=\frac{1}{2} (n(n+1)m-f(c,c-b-1,a,m-1)-h(c,c-b-1,a,m-1))$。
* $h(a,b,c,n)=\sum_{i=0}^n\lfloor \frac{ai+b}{c} \rfloor^2$: 当 $a \ge c$ or $b \ge c$ 时，$h(a,b,c,n)=(\frac{a}{c})^2 n(n+1)(2n+1)/6 +(\frac{b}{c})^2 (n+1)+(\frac{a}{c})(\frac{b}{c})n(n+1)+h(a \bmod c, b \bmod c,c,n)+2(\frac{a}{c})g(a \bmod c,b \bmod c,c,n)+2(\frac{b}{c})f(a \bmod c,b \bmod c,c,n)$；否则 $h(a,b,c,n)=nm(m+1)-2g(c,c-b-1,a,m-1)-2f(c,c-b-1,a,m-1)-f(a,b,c,n)$。
求两直线之间最小的$x,y$使得 $\frac{p1}{q1} < \frac{x}{y} < \frac{p2}{q2}$
```cpp
void solve(LL p1,LL q1,LL p2,LL q2,LL &x,LL &y)
{
LL l=p1/q1+1;
if (l*q2
if (p1==0){x=1; y=q2/p2+1;return;}
if (p1<=q1 && p2<=q2){solve(q2,p2,q1,p1,y,x);return;}
LL t=p1/q1;
solve(p1-q1*t,q1,p2-q2*t,q2,x,y);
x+=y*t;
}
```
## 逆元
* 如果 $p$ 不是素数，使用拓展欧几里得
* 前置模板：快速幂 / 扩展欧几里得
```cpp
inline LL get_inv(LL x, LL p) { return bin(x, p - 2, p); }
LL get_inv(LL a, LL M) {
static LL x, y;
assert(ex_gcd(a, M, x, y) == 1);
return (x % M + M) % M;
}
```
* 预处理 1~n 的逆元
```cpp
LL inv[N];
void inv_init(LL n, LL p) {
inv[1] = 1;
FOR (i, 2, n)
inv[i] = (p - p / i) * inv[p % i] % p;
}
```
* 预处理阶乘及其逆元
```cpp
LL invf[M], fac[M] = {1};
void fac_inv_init(LL n, LL p) {
FOR (i, 1, n)
fac[i] = i * fac[i - 1] % p;
invf[n - 1] = bin(fac[n - 1], p - 2, p);
FORD (i, n - 2, -1)
invf[i] = invf[i + 1] * (i + 1) % p;
}
```
## 组合数
+ 如果数较小，模较大时使用逆元
+ 前置模板：逆元-预处理阶乘及其逆元
+ 卢卡斯定理
+ 任意模数组合数(很慢慎用，基于CRT)
```cpp
inline LL C(LL n, LL m) { // n >= m >= 0
return n < m || m < 0 ? 0 : fac[n] * invf[m] % MOD * invf[n - m] % MOD;
}
```
+ 如果模数较小，数字较大，使用 Lucas 定理
+ 前置模板可选1：求组合数 (如果使用阶乘逆元，需`fac_inv_init(MOD, MOD);`)
+ 前置模板可选2：模数不固定下使用，无法单独使用。
```cpp
LL C(LL n, LL m) { // m >= n >= 0
if (m - n < n) n = m - n;
if (n < 0) return 0;
LL ret = 1;
FOR (i, 1, n + 1)
ret = ret * (m - n + i) % MOD * bin(i, MOD - 2, MOD) % MOD;
return ret;
}
```
```cpp
LL Lucas(LL n, LL m) { // m >= n >= 0
return m ? C(n % MOD, m % MOD) * Lucas(n / MOD, m / MOD) % MOD : 1;
}
```
+ 前置模板：素数，CRT
```cpp
bool nop[MAXN];
vector prime;
void get_prime() {
FOR (i, 2, MAXN) {
if (!nop[i]) prime.push_back(i);
int siz = sz(prime);
FOR (j, 0, siz) {
LL nextp = 1LL*i*prime[j];
if (nextp >= MAXN) break;
nop[nextp] = 1;
if (i % prime[j] == 0) break;
}
}
}
vector > fac;
int cntfac(LL n, LL p) {
int ans = 0;
LL op = p;
for (; n>=p; p*=op) {
ans += n/p;
}
dbg("cntfac", n, op, ans);
return ans;
}
LL ex_gcd(LL a, LL b, LL &x, LL &y) {
if (b == 0) { x = 1; y = 0; return a; }
LL ret = ex_gcd(b, a % b, y, x);
y -= a / b * x;
return ret;
}
LL bin(LL a, LL b, LL p) {
LL res = 1;
for (; b; b>>=1, a=a*a%p)
if (b & 1) res = res*a%p;
return res;
}
LL get_inv(LL a, LL M) {
static LL x, y;
assert(ex_gcd(a, M, x, y) == 1);
LL res = (x % M + M) % M;
dbg("get_inv", a, M, res);
return res;
}
LL mo[MAXN], r[MAXN];
LL calcT(LL n, pair p, LL M) {
LL resT = 1;
for (int i=1; i<=M; ++i) {
if (i % p.first)
resT = resT * i % M;
}
LL cntT = n / M, res = bin(resT, cntT, M), remain = n % M;
dbg(cntT, resT, res, remain);
for (int i=1; i<=remain; ++i) {
if (i % p.first)
res = res * i % M;
}
if (n>=p.first)
res = res * calcT(n/p.first, p, M) % M;
dbg("calcT", n, p.first, p.second, M, res);
return res;
}
LL CRT(LL *m, LL *r, LL n) {
if (!n) return 0;
LL M = m[0], R = r[0], x, y, d;
FOR (i, 1, n) {
d = ex_gcd(M, m[i], x, y);
if ((r[i] - R) % d) return -1;
x = (r[i] - R) / d * x % (m[i] / d);
R += x * M;
M = M / d * m[i];
R %= M;
}
return R >= 0 ? R : R + M;
}
LL lucas(int n, int m, LL p) {
fac.clear();
LL tp = p;
int psiz = prime.size();
FOR (j, 0, psiz) {
if (1LL*prime[j]*prime[j] > tp) break;
int cnt = 0;
while (tp % prime[j] == 0) {
tp /= prime[j];
++cnt;
}
fac.push_back(make_pair(prime[j], cnt));
}
if (tp != 1) fac.push_back(make_pair(tp, 1));
int msiz = fac.size(); dbg(msiz);
FOR (i, 0, msiz) {
dbg(fac[i].first, fac[i].second);
LL M = bin(fac[i].first, fac[i].second, MOD); // dont let it module
mo[i] = M;
int cntp = cntfac(n, fac[i].first) - cntfac(m, fac[i].first) - cntfac(n-m, fac[i].first);
if (cntp >= fac[i].second) r[i] = 0;
else {
LL nM = bin(fac[i].first, fac[i].second - cntp, MOD); // dont let it module
r[i] = calcT(n, fac[i], M) * get_inv(calcT(m, fac[i], M), M) % M * get_inv(calcT(n-m, fac[i], M), M) % M;
dbg(r[i]);
r[i] = r[i] * bin(fac[i].first, cntp, MOD) % M; // dont let it module
dbg(nM, r[i]);
}
dbg(i, M, cntp, r[i]);
}
return CRT(mo, r, msiz);
}
```
* 组合数预处理
```cpp
LL C[M][M];
void init_C(int n) {
FOR (i, 0, n) {
C[i][0] = C[i][i] = 1;
FOR (j, 1, i)
C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % MOD;
}
}
```
## 斯特灵数
### 第一类斯特灵数
+ 绝对值是 $n$ 个元素划分为 $k$ 个环排列的方案数。
+ $s(n,k)=s(n-1,k-1)+(n-1)s(n-1,k)$
### 第二类斯特灵数
+ $n$ 个元素划分为 $k$ 个等价类的方案数
+ $S(n, k)=S(n-1,k-1)+kS(n-1, k)$
```cpp
S[0][0] = 1;
FOR (i, 1, N)
FOR (j, 1, i + 1) S[i][j] = (S[i - 1][j - 1] + j * S[i - 1][j]) % MOD;
```
## FFT & NTT & FWT
### NTT
```cpp
LL wn[N << 2], rev[N << 2];
int NTT_init(int n_) {
int step = 0; int n = 1;
for ( ; n < n_; n <<= 1) ++step;
FOR (i, 1, n)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (step - 1));
int g = bin(G, (MOD - 1) / n, MOD);
wn[0] = 1;
for (int i = 1; i <= n; ++i)
wn[i] = wn[i - 1] * g % MOD;
return n;
}
void NTT(LL a[], int n, int f) {
FOR (i, 0, n) if (i < rev[i])
std::swap(a[i], a[rev[i]]);
for (int k = 1; k < n; k <<= 1) {
for (int i = 0; i < n; i += (k << 1)) {
int t = n / (k << 1);
FOR (j, 0, k) {
LL w = f == 1 ? wn[t * j] : wn[n - t * j];
LL x = a[i + j];
LL y = a[i + j + k] * w % MOD;
a[i + j] = (x + y) % MOD;
a[i + j + k] = (x - y + MOD) % MOD;
}
}
}
if (f == -1) {
LL ninv = get_inv(n, MOD);
FOR (i, 0, n)
a[i] = a[i] * ninv % MOD;
}
}
```
### 任意模数NTT
+ 相当于NTT三次
+ 注意内存和时间
+ 前置模板bin, get_inv, ex_gcd
```cpp
struct NTT {
static const int N = 1e5+5;
const int G = 3;
int nn;
LL MOD;
LL wn[N << 2], rev[N << 2];
int init(int n_, LL M = 998244353, int G = 3) {
MOD = M;
int step = 0, n = 1;
for (; n
for (int i=1; i
rev[i] = (rev[i>>1] >> 1) | ((i & 1) << (step - 1));
}
int g = bin(G, (MOD-1) / n, MOD);
wn[0] = 1;
for (int i=1; i<=n; ++i)
wn[i] = wn[i-1] * g % MOD;
nn = n;
return n;
}
void go(LL a[], int n, int f) {
for (int i=0; i
if (i < rev[i]) swap(a[i], a[rev[i]]);
for (int k=1; k
for (int i=0; i
int t = n/(k<<1);
for (int j=0; j
LL w = f == 1? wn[t*j] : wn[n-t*j];
LL x = a[i + j];
LL y = a[i + j + k] * w % MOD;
a[i + j] = (x+y) % MOD;
a[i + j + k] = (x-y+MOD)%MOD;
}
}
}
if (f == -1) {
LL ninv = get_inv(n, MOD);
for (int i=0; i
a[i] = a[i] * ninv % MOD;
}
}
void Conv(LL a[], LL b[]) {
go(a, nn, 1);
go(b, nn, 1);
for (int i=0; i
a[i] = a[i] * b[i] % MOD;
go(a, nn, -1);
}
} ntt[3];
LL a[MAXN], b[MAXN];
LL ta[3][MAXN<<2], tb[3][MAXN<<2], ans[MAXN<<2];
int n, m, p;
const LL M[3] = {998244353, 469762049, 1004535809},
G[3] = {3, 3, 3};
int main() {
scanf("%d%d%d", &n, &m, &p); ++n; ++m;
for (int i=0; i
for (int i=0; i
for (int i=0; i<3; ++i) {
ntt[i].init(n+m, M[i], G[i]);
}
for (int i=0; i<3; ++i) {
for (int j=0; j
for (int j=0; j
ntt[i].Conv(ta[i], tb[i]);
}
n = n + m - 1;
LL A = M[0], B = M[1], C = M[2], AB = A*B, k[5], x[5];
for (int j=0; j
x[1] = ta[0][j];
x[2] = ta[1][j];
x[3] = ta[2][j];
LL g = x[2]-x[1];
LL res = ex_gcd(A, B, k[1], k[0]);
g *= A;
if (g<0)
g=g%AB+AB;
x[4] = fmul(k[1], g, AB);
x[4] = (x[4] + x[1]) % AB;
LL tmp = (x[3] + C - x[4]%C)%C;
tmp = tmp * get_inv(AB, C) % C;
k[4] = tmp;
k[4] %= p;
ans[j] = (k[4] * (AB%p)%p + x[4]%p)%p;
}
for (int j=0; j
printf("%lld ", ans[j]);
}
printf("%lld\n", ans[n-1]);
return 0;
}
```
### FFT
+ n 需补成 2 的幂 (n 必须超过 a 和 b 的最高指数之和)
```cpp
typedef double LD;
const LD PI = acos(-1);
struct C {
LD r, i;
C(LD r = 0, LD i = 0): r(r), i(i) {}
};
C operator + (const C& a, const C& b) {
return C(a.r + b.r, a.i + b.i);
}
C operator - (const C& a, const C& b) {
return C(a.r - b.r, a.i - b.i);
}
C operator * (const C& a, const C& b) {
return C(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
}
void FFT(C x[], int n, int p) {
for (int i = 0, t = 0; i < n; ++i) {
if (i > t) swap(x[i], x[t]);
for (int j = n >> 1; (t ^= j) < j; j >>= 1);
}
for (int h = 2; h <= n; h <<= 1) {
C wn(cos(p * 2 * PI / h), sin(p * 2 * PI / h));
for (int i = 0; i < n; i += h) {
C w(1, 0), u;
for (int j = i, k = h >> 1; j < i + k; ++j) {
u = x[j + k] * w;
x[j + k] = x[j] - u;
x[j] = x[j] + u;
w = w * wn;
}
}
}
if (p == -1)
FOR (i, 0, n)
x[i].r /= n;
}
void conv(C a[], C b[], int n) {
FFT(a, n, 1);
FFT(b, n, 1);
FOR (i, 0, n)
a[i] = a[i] * b[i];
FFT(a, n, -1);
}
```
### FWT
+ $C_k=\sum_{i \oplus j=k} A_i B_j$
+ FWT 完后需要先模一遍
+ 对于$xor$运算，$FWT[i] =\sum\limits_{|j\&i|=0}{a_j} - \sum\limits_{|j\&i|=1}{a_j}= \sum\limits_{j}(-1)^{|j\& i|}$，$|j\&i|$表示$i\&j$二进制1的个数的奇偶性
+ 对于$or$运算，$FWT[i] = \sum\limits_{j|i=i}a[j]$，j可以认为是i的二进制子集。
+ 对于$and$运算，$FWT[i] = \sum\limits_{j\&i=i}a[i]$, i可以认为是j的二进制子集
```cpp
template
void fwt(LL a[], int n, T f) {
for (int d = 1; d < n; d *= 2)
for (int i = 0, t = d * 2; i < n; i += t)
FOR (j, 0, d)
f(a[i + j], a[i + j + d]);
}
void AND(LL& a, LL& b) { a += b; }
void OR(LL& a, LL& b) { b += a; }
void XOR (LL& a, LL& b) {
LL x = a, y = b;
a = (x + y) % MOD;
b = (x - y + MOD) % MOD;
}
void rAND(LL& a, LL& b) { a -= b; }
void rOR(LL& a, LL& b) { b -= a; }
void rXOR(LL& a, LL& b) {
static LL INV2 = (MOD + 1) / 2;
LL x = a, y = b;
a = (x + y) * INV2 % MOD;
b = (x - y + MOD) * INV2 % MOD;
}
```
+ 三进制异或fwt
```cpp
struct RC {
int x, y;
RC(int x=0, int y=0):x(x), y(y) {}
inline RC operator + (const RC &other) {
return RC(add(x+other.x), add(y+other.y));
}
inline RC operator - (const RC &other) {
return RC(sub(x-other.x), sub(y-other.y));
}
/*
attention: unusual definition!
(a,b) = a+bw
w^2+w+1=0
therefore mul operator is different
*/
inline RC operator * (const RC &other) {
return RC(sub((1LL*x*other.x-1LL*y*other.y)%p), sub((1LL*x*other.y+1LL*y*other.x-1LL*y*other.y)%p));
}
static inline RC w31(const RC &a) {
return RC(sub(-a.y), sub(a.x-a.y));
}
static inline RC w32(const RC &a) {
return RC(sub(a.y-a.x), sub(-a.x));
}
bool operator < (const RC &other) const {
if (x != other.x) return x < other.x;
else return y < other.y;
}
} f[MAXN<<2], g[MAXN<<2];
void fwt(RC a[], bool rfwt = 0) { // n = 3 ^ m, wi 是单位负数根
for (int d=1; d
for (int x=0; x
for (int i=x; i
RC x = a[i], y = a[i+d], z = a[i+(d<<1)];
a[i] = (1LL*x+y+z)%p;
if (rfwt) {
a[i+d] = (1LL*x+1LL*W2*y+1LL*W1*z)%p;
a[i+(d<<1)] = (1LL*x+1LL*W1*y+1LL*W2*z)%p;
} else {
a[i+d] = (x+1LL*W1*y+1LL*W2*z);
a[i+(d<<1)] = (1LL*x+1LL*W2*y+1LL*W1*z)%p;
}
}
}
}
if (rfwt) {
for (int i=0; i
a[i] =1LL* a[i] * invn % p;
}
}
```
+ k进制异或fwt
+ 正变换矩阵 $A_{ij} = w_n^{ij}$ (0 based) $FWT[f] = Af$
+ 逆变换矩阵 $A_{ij} = w_n^{-ij}$ (0 based) $rFWT[f] = Af$
+ FWT 子集卷积
```text
a[popcount(x)][x] = A[x]
b[popcount(x)][x] = B[x]
fwt(a[i]) fwt(b[i])
c[i + j][x] += a[i][x] * b[j][x]
rfwt(c[i])
ans[x] = c[popcount(x)][x]
```
## simpson 自适应积分
```cpp
LD simpson(LD l, LD r) {
LD c = (l + r) / 2;
return (f(l) + 4 * f(c) + f(r)) * (r - l) / 6;
}
LD asr(LD l, LD r, LD eps, LD S) {
LD m = (l + r) / 2;
LD L = simpson(l, m), R = simpson(m, r);
if (fabs(L + R - S) < 15 * eps) return L + R + (L + R - S) / 15;
return asr(l, m, eps / 2, L) + asr(m, r, eps / 2, R);
}
LD asr(LD l, LD r, LD eps) { return asr(l, r, eps, simpson(l, r)); }
```
## 快速乘
```cpp
LL mul(LL a, LL b, LL m) {
LL ret = 0;
while (b) {
if (b & 1) {
ret += a;
if (ret >= m) ret -= m;
}
a += a;
if (a >= m) a -= m;
b >>= 1;
}
return ret;
}
```
+ O(1)
```cpp
LL mul(LL u, LL v, LL p) {
return (u * v - LL((long double) u * v / p) * p + p) % p;
}
LL mul(LL u, LL v, LL p) { // 卡常
LL t = u * v - LL((long double) u * v / p) * p;
return t < 0 ? t + p : t;
}
```
## 快速幂
* 如果模数是素数，则可在函数体内加上`n %= MOD - 1;`(费马小定理)。
```cpp
LL bin(LL x, LL n, LL MOD) {
LL ret = MOD != 1;
for (x %= MOD; n; n >>= 1, x = x * x % MOD)
if (n & 1) ret = ret * x % MOD;
return ret;
}
```
* 防爆 LL
* 前置模板：快速乘
```cpp
LL bin(LL x, LL n, LL MOD) {
LL ret = MOD != 1;
for (x %= MOD; n; n >>= 1, x = mul(x, x, MOD))
if (n & 1) ret = mul(ret, x, MOD);
return ret;
}
```
## 质因数分解
* 前置模板：素数筛
* 带指数
```cpp
LL factor[30], f_sz, factor_exp[30];
void get_factor(LL x) {
f_sz = 0;
LL t = sqrt(x + 0.5);
for (LL i = 0; pr[i] <= t; ++i)
if (x % pr[i] == 0) {
factor_exp[f_sz] = 0;
while (x % pr[i] == 0) {
x /= pr[i];
++factor_exp[f_sz];
}
factor[f_sz++] = pr[i];
}
if (x > 1) {
factor_exp[f_sz] = 1;
factor[f_sz++] = x;
}
}
```
* 不带指数
```cpp
LL factor[30], f_sz;
void get_factor(LL x) {
f_sz = 0;
LL t = sqrt(x + 0.5);
for (LL i = 0; pr[i] <= t; ++i)
if (x % pr[i] == 0) {
factor[f_sz++] = pr[i];
while (x % pr[i] == 0) x /= pr[i];
}
if (x > 1) factor[f_sz++] = x;
}
```
## 原根
* 前置模板：素数筛，快速幂，分解质因数
* 要求 p 为质数
```cpp
LL find_smallest_primitive_root(LL p) {
get_factor(p - 1);
FOR (i, 2, p) {
bool flag = true;
FOR (j, 0, f_sz)
if (bin(i, (p - 1) / factor[j], p) == 1) {
flag = false;
break;
}
if (flag) return i;
}
assert(0); return -1;
}
```
## 公式
### 一些数论公式
- 当 $x\geq\phi(p)$ 时有 $a^x\equiv a^{x \; mod \; \phi(p) + \phi(p)}\pmod p$
- $\mu^2(n)=\sum_{d^2|n} \mu(d)$
- $\sum_{d|n} \varphi(d)=n$
- $\sum_{d|n} 2^{\omega(d)}=\sigma_0(n^2)$，其中 $\omega$ 是不同素因子个数
- $\sum_{d|n} \mu^2(d)=2^{\omega(d)}$
### 一些狄利克雷卷积
+ 满足交换律和结合律，$f*g$也是积性函数
+ $id(n) = n， I(n)=1$
+ $\mu * 1 = [n==1]$ 枚举素因子种类
+ 除数函数$d(n)=\sum_{d|n}1=1*1$
+ $\sigma(n) = \sum_{d|n}d = id * 1$
+ $\phi * 1 = id$, $\phi = \mu * id$
+ $n=\prod_{i=1}^{t}p_i^{k_i},g(n)=f * 1 \Rightarrow g(n) = \prod_{i=1}^{t} \sum_{j=0}^{k_i}{f(p_i^j)}$ 相当于枚举素因数指数
### 一些数论函数求和的例子
+ $\sum_{i=1}^{n}\sum_{d|i}f(d) = \sum_{i=1}^{n}\sum_{i=1}^{\lfloor \frac{n}{i} \rfloor}f(d)$ 计算$f(d)$出现了多少次，杜教筛的核心思想
+ $\sum\limits_{i=1}^{n}{gcd(i, a)} = \sum\limits_{d|a}d\sum\limits_{i=1}^{n}{[gcd(i,a)==d]}=\sum\limits_{d|a}d\sum\limits_{i=1}^{n}[gcd(i/d, a/d)==1]=\sum\limits_{d|a}d\sum\limits_{i=1}^{n}u * 1(gcd(i/d,a/d))=\sum\limits_{d|a}\sum\limits_{i}\sum_{t|gcd(i/d, a/d)}\mu(t) 拿T=td替换,有d|T,=\sum\limits_d\sum\limits_{i}\sum_{T|gcd(a,d)}d\cdot \mu(T/d)=\sum\limits_{T|a}\frac{n}{T}\sum\limits_{d|T}d\mu{(\frac{T}{d})} = \sum\limits_{T|a}\frac{n}{T}\phi(T)$
+ $\sum_{i=1}^n i[gcd(i, n)=1] = \frac {n \varphi(n) + [n=1]}{2}$
+ $\sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=x]=\sum_d \mu(d) \lfloor \frac n {dx} \rfloor \lfloor \frac m {dx} \rfloor$
+ $\sum_{i=1}^n \sum_{j=1}^m gcd(i, j) = \sum_{i=1}^n \sum_{j=1}^m \sum_{d|gcd(i,j)} \varphi(d) = \sum_{d} \varphi(d) \lfloor \frac nd \rfloor \lfloor \frac md \rfloor$
+ $S(n)=\sum_{i=1}^n \mu(i)=1-\sum_{i=1}^n \sum_{d|i,d < i}\mu(d) \overset{t=\frac id}{=} 1-\sum_{t=2}^nS(\lfloor \frac nt \rfloor)$
+ 利用 $[n=1] = \sum_{d|n} \mu(d)$ 杜教筛
+ $S(n)=\sum_{i=1}^n \varphi(i)=\sum_{i=1}^n i-\sum_{i=1}^n \sum_{d|i,d
+ 利用 $n = \sum_{d|n} \varphi(d)$ 杜教筛
+ $S(n) = \sum\limits_{i=1}^{n}i^2\phi(i)=\sum\limits_{i=1}^{n}i^3 - \sum\limits_{i=2}^{n}i^2S(n/i)$ 杜教筛
+ $\sum_{i=1}^n \mu^2(i) = \sum_{i=1}^n \sum_{d^2|n} \mu(d)=\sum_{d=1}^{\lfloor \sqrt n \rfloor}\mu(d) \lfloor \frac n {d^2} \rfloor$
+ $\sum_{i=1}^n \sum_{j=1}^n gcd^2(i, j)= \sum_{d} d^2 \sum_{t} \mu(t) \lfloor \frac n{dt} \rfloor ^2 \\
\overset{x=dt}{=} \sum_{x} \lfloor \frac nx \rfloor ^ 2 \sum_{d|x} d^2 \mu(\frac tx)$
+ $\sum_{i=1}^n \varphi(i)=\frac 12 \sum_{i=1}^n \sum_{j=1}^n [i \perp j] - 1=\frac 12 \sum_{i=1}^n \mu(i) \cdot\lfloor \frac n i \rfloor ^2-1$
+ $\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}ij[gcd(i,j)==d] = \sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{n/d}ij[gcd(i,j)==1]=\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{n/d}i^2\phi(i)$ 然后后面的可以用杜教筛
### 斐波那契数列性质
- $F_{a+b}=F_{a-1} \cdot F_b+F_a \cdot F_{b+1}$
- $F_1+F_3+\dots +F_{2n-1} = F_{2n},F_2 + F_4 + \dots + F_{2n} = F_{2n + 1} - 1$
- $\sum_{i=1}^n F_i = F_{n+2} - 1$
- $\sum_{i=1}^n F_i^2 = F_n \cdot F_{n+1}$
- $F_n^2=(-1)^{n-1} + F_{n-1} \cdot F_{n+1}$
- $gcd(F_a, F_b)=F_{gcd(a, b)}$
- 模 $n$ 周期(皮萨诺周期)
- $\pi(p^k) = p^{k-1} \pi(p)$
- $\pi(nm) = lcm(\pi(n), \pi(m)), \forall n \perp m$
- $\pi(2)=3, \pi(5)=20$
- $\forall p \equiv \pm 1\pmod {10}, \pi(p)|p-1$
- $\forall p \equiv \pm 2\pmod {5}, \pi(p)|2p+2$
### 常见生成函数
+ $(1+ax)^n=\sum_{k=0}^n \binom {n}{k} a^kx^k$
+ $\dfrac{1-x^{r+1}}{1-x}=\sum_{k=0}^nx^k$
+ $\dfrac1{1-ax}=\sum_{k=0}^{\infty}a^kx^k$
+ $\dfrac 1{(1-x)^2}=\sum_{k=0}^{\infty}(k+1)x^k$
+ $\dfrac1{(1-x)^n}=\sum_{k=0}^{\infty} \binom{n+k-1}{k}x^k$
+ $e^x=\sum_{k=0}^{\infty}\dfrac{x^k}{k!}$
+ $\ln(1+x)=\sum_{k=0}^{\infty}\dfrac{(-1)^{k+1}}{k}x^k$
### 佩尔方程
若一个丢番图方程具有以下的形式：$x^2 - ny^2= 1$。且 $n$ 为正整数，则称此二元二次不定方程为**佩尔方程**。
若 $n$ 是完全平方数，则这个方程式只有平凡解 $(\pm 1,0)$(实际上对任意的 $n$，$(\pm 1,0)$ 都是解)。对于其余情况，拉格朗日证明了佩尔方程总有非平凡解。而这些解可由 $\sqrt{n}$ 的连分数求出。
$x = [a_0; a_1, a_2, a_3]=x = a_0 + \cfrac{1}{a_1 + \cfrac{1}{a_2 + \cfrac{1}{a_3 + \cfrac{1}{\ddots\,}}}}$
设 $\tfrac{p_i}{q_i}$ 是 $\sqrt{n}$ 的连分数表示：$[a_{0}; a_{1}, a_{2}, a_{3}, \,\ldots ]$ 的渐近分数列，由连分数理论知存在 $i$ 使得 $(p_i,q_i)$ 为佩尔方程的解。取其中最小的 $i$，将对应的 $(p_i,q_i)$ 称为佩尔方程的基本解，或最小解，记作 $(x_1,y_1)$，则所有的解 $(x_i,y_i)$ 可表示成如下形式：$x_{i}+y_{i}{\sqrt n}=(x_{1}+y_{1}{\sqrt n})^{i}$。或者由以下的递回关系式得到：
$\displaystyle x_{i+1} = x_1 x_i + n y_1 y_i$, $\displaystyle y_{{i+1}}=x_{1}y_{i}+y_{1}x_{i}$。
**但是：**佩尔方程千万不要去推(虽然推起来很有趣，但结果不一定好看，会是两个式子)。记住佩尔方程结果的形式通常是 $a_n=ka_{n−1}−a_{n−2}$($a_{n−2}$ 前的系数通常是 $−1$)。暴力 / 凑出两个基础解之后加上一个 $0$，容易解出 $k$ 并验证。
### Burnside & Polya
+ $|X/G|={\frac {1}{|G|}}\sum _{{g\in G}}|X^{g}|$
注：$X^g$ 是 $g$ 下的不动点数量，也就是说有多少种东西用 $g$ 作用之后可以保持不变。
+ $|Y^X/G| = \frac{1}{|G|}\sum_{g \in G} m^{c(g)}$
注：用 $m$ 种颜色染色，然后对于某一种置换 $g$，有 $c(g)$ 个置换环，为了保证置换后颜色仍然相同，每个置换环必须染成同色。
### 皮克定理
$2S = 2a+b-2$
+ $S$ 多边形面积
+ $a$ 多边形内部点数
+ $b$ 多边形边上点数
### 莫比乌斯反演
+ $g(n) = \sum_{d|n} f(d) \Leftrightarrow f(n) = \sum_{d|n} \mu (d) g( \frac{n}{d})$
+ $f(n)=\sum_{n|d}g(d) \Leftrightarrow g(n)=\sum_{n|d} \mu(\frac{d}{n}) f(d)$
### 二项式反演
+ s是非负整数
$a_n = \sum_{k=s}^{n}\binom{n}{k}b_k \Rightarrow b_n = \sum_{k=s}^{n}(-1)^{n-k}\binom{n}{k}a_k$
### 低阶等幂求和
+ $\sum_{i=1}^{n} i^{1} = \frac{n(n+1)}{2} = \frac{1}{2}n^2 +\frac{1}{2} n$
+ $\sum_{i=1}^{n} i^{2} = \frac{n(n+1)(2n+1)}{6} = \frac{1}{3}n^3 + \frac{1}{2}n^2 + \frac{1}{6}n$
+ $\sum_{i=1}^{n} i^{3} = \left[\frac{n(n+1)}{2}\right]^{2} = \frac{1}{4}n^4 + \frac{1}{2}n^3 + \frac{1}{4}n^2$
+ $\sum_{i=1}^{n} i^{4} = \frac{n(n+1)(2n+1)(3n^2+3n-1)}{30} = \frac{1}{5}n^5 + \frac{1}{2}n^4 + \frac{1}{3}n^3 - \frac{1}{30}n$
+ $\sum_{i=1}^{n} i^{5} = \frac{n^{2}(n+1)^{2}(2n^2+2n-1)}{12} = \frac{1}{6}n^6 + \frac{1}{2}n^5 + \frac{5}{12}n^4 - \frac{1}{12}n^2$
### 一些组合公式
+ 错排公式：$D_1=0,D_2=1,D_n=(n-1)(D_{n-1} + D_{n-2})=n!(\frac 1{2!}-\frac 1{3!}+\dots + (-1)^n\frac 1{n!})=\lfloor \frac{n!}e + 0.5 \rfloor$
+ 卡塔兰数($n$ 对括号合法方案数，$n$ 个结点二叉树个数，$n\times n$ 方格中对角线下方的单调路径数，凸 $n+2$ 边形的三角形划分数，$n$ 个元素的合法出栈序列数)：$C_n=\frac 1{n+1}\binom {2n}n=\frac{(2n)!}{(n+1)!n!}$
## 二次剩余
URAL 1132
```cpp
LL q1, q2, w;
struct P { // x + y * sqrt(w)
LL x, y;
};
P pmul(const P& a, const P& b, LL p) {
P res;
res.x = (a.x * b.x + a.y * b.y % p * w) % p;
res.y = (a.x * b.y + a.y * b.x) % p;
return res;
}
P bin(P x, LL n, LL MOD) {
P ret = {1, 0};
for (; n; n >>= 1, x = pmul(x, x, MOD))
if (n & 1) ret = pmul(ret, x, MOD);
return ret;
}
LL Legendre(LL a, LL p) { return bin(a, (p - 1) >> 1, p); }
LL equation_solve(LL b, LL p) {
if (p == 2) return 1;
if ((Legendre(b, p) + 1) % p == 0)
return -1;
LL a;
while (true) {
a = rand() % p;
w = ((a * a - b) % p + p) % p;
if ((Legendre(w, p) + 1) % p == 0)
break;
}
return bin({a, 1}, (p + 1) >> 1, p).x;
}
int main() {
int T; cin >> T;
while (T--) {
LL a, p; cin >> a >> p;
a = a % p;
LL x = equation_solve(a, p);
if (x == -1) {
puts("No root");
} else {
LL y = p - x;
if (x == y) cout << x << endl;
else cout << min(x, y) << " " << max(x, y) << endl;
}
}
}
```
## 中国剩余定理
+ 无解返回 -1
+ 前置模板：扩展欧几里得
```cpp
LL CRT(LL *m, LL *r, LL n) {
if (!n) return 0;
LL M = m[0], R = r[0], x, y, d;
FOR (i, 1, n) {
d = ex_gcd(M, m[i], x, y);
if ((r[i] - R) % d) return -1;
x = (r[i] - R) / d * x % (m[i] / d);
R += x * M;
M = M / d * m[i];
R %= M;
}
return R >= 0 ? R : R + M;
}
```
## 伯努利数和等幂求和
* 预处理逆元
* 预处理组合数
* $\sum_{i=0}^n i^k = \frac{1}{k+1} \sum_{i=0}^k \binom{k+1}{i} B_{k+1-i} (n+1)^i$.
* 也可以 $\sum_{i=0}^n i^k = \frac{1}{k+1} \sum_{i=0}^k \binom{k+1}{i} B^+_{k+1-i} n^i$。区别在于 $B^+_1 =1/2$。(心态崩了)
```cpp
namespace Bernoulli {
const int M = 100;
LL inv[M] = {-1, 1};
void inv_init(LL n, LL p) {
FOR (i, 2, n)
inv[i] = (p - p / i) * inv[p % i] % p;
}
LL C[M][M];
void init_C(int n) {
FOR (i, 0, n) {
C[i][0] = C[i][i] = 1;
FOR (j, 1, i)
C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % MOD;
}
}
LL B[M] = {1};
void init() {
inv_init(M, MOD);
init_C(M);
FOR (i, 1, M - 1) {
LL& s = B[i] = 0;
FOR (j, 0, i)
s += C[i + 1][j] * B[j] % MOD;
s = (s % MOD * -inv[i + 1] % MOD + MOD) % MOD;
}
}
LL p[M] = {1};
LL go(LL n, LL k) {
n %= MOD;
if (k == 0) return n;
FOR (i, 1, k + 2)
p[i] = p[i - 1] * (n + 1) % MOD;
LL ret = 0;
FOR (i, 1, k + 2)
ret += C[k + 1][i] * B[k + 1 - i] % MOD * p[i] % MOD;
ret = ret % MOD * inv[k + 1] % MOD;
return ret;
}
}
```
### 拉格朗日插值法
+ 输入点值表示和幂次
+ 输出多项式系数
+ 运算为`MOD`的剩余系
+ $O(n^2)$
```cpp
namespace interpolation {
/*
integer mod Lagrange interpolation formula
*/
const int MAXN = 1005;
typedef long long LL;
typedef vector Poly;
Poly operator * (const Poly &a, const Poly &b) {
/*
O(n^2) Polynomial multiplication method
*/
Poly c(a.size() + b.size() - 1, 0);
int an = a.size(), bn = b.size();
for (int i=0; i
for (int j=0; j
c[i+j] = (c[i+j] + a[i]*b[j])%MOD;
while (!c.empty() && c.back() == 0) c.pop_back();
return c;
}
Poly operator + (const Poly &a, const Poly &b) {
int an = a.size(), bn = b.size(), n = max(an, bn);
Poly c; c.resize(n);
for (int i=0; i
c[i] = 0;
if (i < an) c[i] += a[i];
if (i < bn) c[i] += b[i];
while (c[i] >= MOD) c[i] -= MOD;
}
while (!c.empty() && c.back() == 0) c.pop_back();
return c;
}
void operator *= (Poly &a, LL b) {
int n = a.size();
for (int i=0; i
a[i] = a[i] * b % MOD;
}
Poly g[MAXN];
LL c[MAXN];
Poly go(LL x[], LL y[], int n) {
/*
make sure x_i is unique
index start from 0
*/
for (int i=0; i
LL B = 1;
g[i] = Poly({1});
for (int j=0; j
if (i == j) continue;
B = B * (x[i] - x[j] + MOD) % MOD;
Poly p = Poly({ MOD-x[j], 1});
g[i] = g[i] * p;
}
g[i] *= (MOD + y[i]) * inv(B, MOD)%MOD;
}
Poly res; res.clear();
for (int i=0; i
res = res + g[i];
}
return res;
}
}
```
+ 点值 $x_i = i, 0 \le i < n$ 可$O(n)$求，否则带入 $x$ 可 $O(n^2)$ 单次求值
+ 预处理逆元
```cpp
namespace interpolation {
LL go(LL y[], int n, LL x) {
dbg(n, x);
if (x == -1) return 0;
if (x < n) return y[x];
LL res = 0;
for (int i=0; i
// dbg(i, f[x], bin(x-1, MOD-2, MOD), invf[x-n], invf[n-1-i], invf[i], y[i], x-n, n-1-i);
LL del = f[x] * invd[x-i] % MOD * invf[x-n] % MOD;
int sgn = (n-i-1)%2? -1:1;
del = del * sgn * invf[n-1-i] % MOD * invf[i] % MOD;
del = (del + MOD) % MOD * y[i] % MOD;
res = (res + MOD + del) % MOD;
// dbg(sgn, del, i);
}
return res;
}
}
```
### 拉格朗日插值法求等幂求和
```cpp
namespace sumk {
const int MAXK = 205;
LL a[MAXK], b[MAXK];
vector c;
int cn;
vector go(int k) {
int n = k+2;
LL res = 0;
for (int i=1; i<=n; ++i) {
a[i] = i;
res = (res + bin(i, k, MOD))%MOD;
b[i] = res;
}
c = interpolation::go(a, b, n);
cn = c.size();
return c;
}
LL calc(LL x) {
x %= MOD;
LL res = 0;
LL xx = 1;
for (int i=0; i
res = (res + xx * c[i]) % MOD;
xx = xx * x % MOD;
}
return res;
}
}
```
### 单纯形
+ uoj 97分
+ `simplex::init()`和`simplex::solve()`
```cpp
namespace simplex {
const int N = MAXN;
const int M = N << 1;
const double eps = 1e-10,
INF = 1e100;
double a[M][M], v, ans[N];
int b[N], id[N<<1];
int n, m;
void pivot(int x, int y) {
swap(id[n+x], id[y]);
double t = a[x][y]; a[x][y] = 1;
for (int i=0; i<=n; ++i) a[x][i] /= t;
for (int i=0; i<=m; ++i) {
if (i == x || sgn(a[i][y]) == 0) continue;
t = a[i][y]; a[i][y] = 0;
for (int j=0; j<=n; ++j) a[i][j] -= a[x][j] * t;
}
}
int simplex() {
while (1) {
int x = 0, y = 0;
for (int j=1; !y && j<=n; ++j) if (sgn(a[0][j]) == 1) y = j;
if (!y) return 1;
double minX = INF;
for (int i=1; i<=m; ++i)
if (sgn(a[i][y])==1 && minX > a[i][0] / a[i][y]) minX = a[i][0] / a[i][y], x = i;
if (!x) return 0; // unbound
pivot(x, y);
}
}
int solve() {
// 1 means sol, 0 for unbound(INF), -1 for infeasible
for (int i=1; i<=n; ++i) id[i] = i;
while (1) {
int x = 0, y = 0;
for (int i=1; i<=m; ++i) if (sgn(a[i][0]) == -1 && (!x || rand()&1)) x = i;
if (!x) {
int res = simplex();
if (res == 1) {
v = -a[0][0];
for (int i=1; i<=m; ++i) ans[id[n+i]] = a[i][0];
for (int i=0; i
}
return res;
}
for (int j=1; j<=n; ++j) if (sgn(a[x][j]) == -1 && (!y || rand()&1)) y = j;
if (!y) return -1; // infeasible
pivot(x, y);
}
}
void init(double _a[N][N], double b[N], double c[N], int _n, int _m) {
// _a[][], b[], c[] are 0-based
// m is the num of equations
// n is the num of variables, so m rows n cols
// _ax <= b / max(cx)
n = _n; m = _m;
for (int i=1; i<=n; ++i) a[0][i] = c[i-1];
for (int i=1; i<=m; ++i)
for (int j=1; j<=n; ++j)
a[i][j] = _a[i-1][j-1];
for (int i=1; i<=m; ++i)
a[i][0] = b[i-1];
}
}
```
## BSGS
+ 针对多次询问，复杂度 $O(\sqrt{qn})$
```cpp
namespace BSGS {
LL a, p, m, n, q;
unordered_map mp;
void init(LL _a, LL _p, LL _q = 1) {
a = _a; p = _p;
m = ceil(sqrt((double)p*_q+1.5));
// dbg(m);
n = ceil(1.0*p/m);
// swap(n, m);
// dbg(n, m);
mp.clear();
LL v = 1;
for (int i=1; i<=m; ++i) {
v = v * a % p;
mp[v] = i;
}
q = v;
dbg("init");
}
LL query(LL b) {
dbg(q, bin(a, m, p));
LL v = 1;
LL invb = inv(b, p);
dbg(invb);
for (int i=1; i<=n; ++i) {
v = v * q % p;
LL tar = v * invb % p;
if (mp.count(tar)) return i * m - mp[tar];
}
return -1;
}
}
```
## exBSGS
+ 模数可以非素数
```cpp
LL exBSGS(LL a, LL b, LL p) { // a^x = b (mod p)
a %= p; b %= p;
if (a == 0) return b > 1 ? -1 : b == 0 && p != 1;
LL c = 0, q = 1;
while (1) {
LL g = __gcd(a, p);
if (g == 1) break;
if (b == 1) return c;
if (b % g) return -1;
++c; b /= g; p /= g; q = a / g * q % p;
}
static map mp; mp.clear();
LL m = ceil(sqrt(p + 1.5))+10;
LL v = 1;
FOR (i, 1, m + 1) {
v = v * a % p;
mp[v * b % p] = i;
}
/* 要求正数解
FOR (i, 1, m + 1) {
v = v * a % p;
mp[v * b % p] = i;
}
*/
FOR (i, 1, m + 1) {
q = q * v % p;
auto it = mp.find(q);
if (it != mp.end()) return i * m - it->second + c;
}
return -1;
}
```
### 数论分块
$f(i) = \lfloor \frac{n}{i} \rfloor=v$ 时 $i$ 的取值范围是 $[l,r]$。
```cpp
for (LL l = 1, v, r; l <= N; l = r + 1) {
v = N / l; r = N / v;
}
```
### 约数个数和
+ $O(sqrt(n))$
```cpp
namespace dsum {
LL sum(LL i) {
LL x = i%MOD;
x = (1+x)*x/2%MOD;
return x;
}
LL solve(LL n) {
LL tar = sqrt(n);
LL ans = 0;
FOR (i, 1, tar+1) {
ans = (ans + i*((n/i)%MOD) %MOD) % MOD;
}
for (LL l=tar+1, r, v; l<=n; l=r+1) {
v = n/l;
r = n/v;
ans = (ans + v * (sum(r) - sum(l-1) + MOD)%MOD) %MOD;
}
return ans;
}
}
```
### 博弈
+ Nim 游戏：每轮从若干堆石子中的一堆取走若干颗。先手必胜条件为石子数量异或和非零。
+ 阶梯 Nim 游戏：可以选择阶梯上某一堆中的若干颗向下推动一级，直到全部推下去。先手必胜条件是奇数阶梯的异或和非零(对于偶数阶梯的操作可以模仿)。
+ Anti-SG：无法操作者胜。先手必胜的条件是：
+ SG 不为 0 且某个单一游戏的 SG 大于 1 。
+ SG 为 0 且没有单一游戏的 SG 大于 1。
+ Every-SG：对所有单一游戏都要操作。先手必胜的条件是单一游戏中的最大 step 为奇数。
+ 对于终止状态 step 为 0
+ 对于 SG 为 0 的状态，step 是最大后继 step +1
+ 对于 SG 非 0 的状态，step 是最小后继 step +1
+ 树上删边：叶子 SG 为 0，非叶子结点为所有子结点的 SG 值加 1 后的异或和。
尝试：
+ 打表找规律
+ 寻找一类必胜态(如对称局面)
+ 直接博弈 dp
### 多项式大餐 X 线性递推式
+ MTT基于毛子黑科技的任意模数fft，链接https://blog.csdn.net/lvzelong2014/article/details/80156989
+ 加减乘除，模数乘
+ 求逆
+ 求线性递推式$O(klogklogn)$
+ bzoj4161
```cpp
#include
#include
#include
#include
#include
inline int read() {
register int ret, cc, sign = 1;
while (!isdigit(cc = getchar())){ sign = cc == '-' ? -1 : sign; }ret = cc-48;
while ( isdigit(cc = getchar())) ret = cc-48+ret*10;
return ret * sign;
}
inline void write(int x, char ch = '\n') {
register int stk[20], tp;
stk[tp = !x] = 0;
while (x) stk[++tp] = x % 10, x /= 10;
while (tp) putchar(stk[tp--] + '0');
putchar(ch);
}
typedef std::vector Poly;
const double PI = acos(-1.0);
inline int getrev(int);
inline int Qpow(int, int, int);
inline Poly Inverse(const Poly&);
inline void Read(Poly&);
inline void Print(const Poly&);
inline Poly operator * (const Poly&, const Poly&);
inline Poly operator - (const Poly&, const Poly&);
inline Poly operator + (const Poly&, const Poly&);
inline Poly operator / (const Poly&, const Poly&);
inline Poly operator % (const Poly&, const Poly&);
const int MAXN = 300010;
const int MOD = 1e9 + 7;
int g[MAXN];
int f[MAXN];
int N, K;
inline void mul_mod(int* a, int* b) {
static int tmp[MAXN];
for (int i = 0; i < K << 1; ++i) tmp[i] = 0;
for (int i = 0; i < K; ++i)
for (int j = 0; j < K; ++j)
tmp[i+j] = (tmp[i+j] + 1ll*a[i]*b[j]%MOD) % MOD;
for (int i = (K<<1)-1; i >= K; --i)
for (int j = 1; j <= K; ++j)
tmp[i-j] = (tmp[i-j] + 1ll*tmp[i]*g[j]%MOD) % MOD;
for (int i = 0; i < K; ++i) a[i] = tmp[i];
}
Poly Tg, Rg;
inline Poly mul_mod(const Poly& lhs, const Poly& rhs) {
Poly A = lhs * rhs, B = A;
int len = K-1;
std::reverse(A.begin(), A.end());
A.resize(len);
Poly C = A * Rg;
C.resize(len);
std::reverse(C.begin(), C.end());
Poly D = B - Tg * C;
D.resize(K);
return D;
}
int a[MAXN];
int b[MAXN];
int main() {
#ifdef ARK
freopen("3.in", "r", stdin);
#endif
N = read(), K = read();
Tg.resize(K+1), Rg.resize(K+1);
for (int i = 1; i <= K; ++i) Tg[K-i] = Rg[K-i] = (MOD-read())%MOD;
for (int i = 0; i < K; ++i) f[i] = (read()+MOD)%MOD;
Tg[K] = Rg[K] = 1;
std::reverse(Rg.begin(), Rg.end());
Rg.resize(K-1);
Rg = Inverse(Rg);
Poly A(K), B(K);
A[0] = B[1] = 1;
while (N) {
if (N & 1) A = mul_mod(A, B);
B = mul_mod(B, B); N >>= 1;
//Print(B);
}
int ans = 0;
for (int i = 0; i < K; ++i) ans = (ans + 1ll*A[i]*f[i]%MOD)%MOD;
printf("%d\n", ans);
}
struct Complex {
double x, y;
Complex(double x = 0, double y = 0):x(x), y(y) { }
inline Complex operator + (const Complex& rhs) const { return Complex(x + rhs.x, y + rhs.y); }
inline Complex operator - (const Complex& rhs) const { return Complex(x - rhs.x, y - rhs.y); }
inline Complex operator * (const Complex& rhs) const { return Complex(x * rhs.x - y * rhs.y, x * rhs.y + y * rhs.x); }
inline Complex operator - () const { return Complex(-x, -y); }
inline Complex operator ! () const { return Complex( x, -y); }
void print() { printf("(%f, %f)\n", x, y); }
};
Complex W[2][MAXN];
int rev[MAXN];
int Inv[MAXN];
inline int getrev(int n) {
int bln = 1, bct = 0;
while (bln <= n) bln <<= 1, bct++;
for (int i = 0; i < bln; ++i)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bct - 1));
for (int i = 0; i < bln; ++i) {
W[0][i] = W[1][(bln-i)&(bln-1)] = Complex(cos(2*PI*i/bln), sin(2*PI*i/bln));
}
return bln;
}
inline void getinv(int n) {
Inv[1] = 1;
for (int i = 2; i < n; ++i)
Inv[i] = 1ll * (MOD - MOD / i) * Inv[MOD % i] % MOD;
}
inline void FFT(Complex* a, int n, int opt) {
for (int i = 0; i < n; ++i)
if (i < rev[i]) std::swap(a[i], a[rev[i]]);
for (int i = 1, t = n >> 1; i < n; i <<= 1, t >>= 1) {
for (int j = 0, p = (i << 1); j < n; j += p) {
Complex *w = W[opt];
for (int k = 0; k < i; ++k, w += t) {
Complex x = a[j + k];
Complex y = *w * a[j + k + i];
a[j + k] = x + y;
a[j + k + i] = x - y;
}
}
}
if (opt) for (int i = 0; i < n; ++i) a[i].x /= n, a[i].y /= n;
}
inline Poly MTT(const Poly& A, const Poly& B, int n) {
static Complex a[MAXN], b[MAXN], c[MAXN], d[MAXN];
for (size_t i = 0; i < A.size(); ++i) a[i] = Complex(A[i] & 32767, A[i] >> 15);
for (int i = A.size(); i < n; ++i) a[i] = Complex();
for (size_t i = 0; i < B.size(); ++i) b[i] = Complex(B[i] & 32767, B[i] >> 15);
for (int i = B.size(); i < n; ++i) b[i] = Complex();
FFT(a, n, 0), FFT(b, n, 0);
for (int i = 0; i < n; ++i) {
int j = (n - i) & (n - 1);
c[i] = Complex(a[i].x + a[j].x, a[i].y - a[j].y) * 0.5 * b[i];
d[i] = Complex(a[i].y + a[j].y, a[j].x - a[i].x) * 0.5 * b[i];
}
FFT(c, n, 1), FFT(d, n, 1);
Poly C(n);
for (int i = 0; i < n; ++i) {
long long u = llround(c[i].x) % MOD, v = llround(c[i].y) % MOD;
long long x = llround(d[i].x) % MOD, y = llround(d[i].y) % MOD;
C[i] = (u + ((v + x) << 15) % MOD + (y << 30)) % MOD;
}
return C;
}
inline int Qpow(int a, int p = MOD - 2, int mod = MOD) {
int ret = 1;
p += (p < 0) * (mod - 1);
while (p) {
if (p & 1) ret = 1ll * ret * a % mod;
a = 1ll * a * a % mod; p >>= 1;
}
return ret;
}
inline Poly operator * (const Poly& A, const Poly& B) {
int len = A.size() + B.size() - 1;
int bln = getrev(len);
Poly C = MTT(A, B, bln);
C.resize(len);
return C;
}
inline Poly operator + (const Poly& lhs, const Poly& rhs) {
Poly A(std::max(lhs.size(), rhs.size()));
for (size_t i = 0; i < A.size(); ++i)
A[i] = ((i < lhs.size() ? lhs[i] : 0) + (i < rhs.size() ? rhs[i] : 0)) % MOD;
return A;
}
inline Poly operator - (const Poly& lhs, const Poly& rhs) {
Poly A(std::max(lhs.size(), rhs.size()));
for (size_t i = 0; i < A.size(); ++i)
A[i] = ((i < lhs.size() ? lhs[i] : 0) - (i < rhs.size() ? rhs[i] : 0) + MOD) % MOD;
return A;
}
inline Poly Inverse(const Poly& A) {
Poly B(1, Qpow(A[0], MOD - 2));
int n = A.size() << 1;
for (int i = 2; i < n; i <<= 1) {
Poly C(A);
C.resize(i);
B = Poly(1, 2) * B - C * B * B;
B.resize(i);
}
B.resize(A.size());
return B;
}
inline Poly operator / (const Poly& lhs, const Poly& rhs) {
return lhs * Inverse(rhs);
}
inline Poly operator % (const Poly& lhs, const Poly& rhs) {
Poly A(lhs), B(rhs);
int len = A.size() - B.size() + 1;
std::reverse(A.begin(), A.end());
std::reverse(B.begin(), B.end());
A.resize(len), B.resize(len);
Poly C = A * Inverse(B);
C.resize(len);
std::reverse(C.begin(), C.end());
Poly D = lhs - rhs * C;
D.resize(rhs.size() - 1);
return D;
}
inline void Read(Poly& A) {
std::generate(A.begin(), A.end(), read);
}
inline void Print(const Poly& A) {
for (size_t i = 0; i < A.size(); ++i)
printf("%d ", A[i]);
puts("");
}
```
### 法雷序列
+ 输出两分数之间最小分母的最简分数
```cpp
Frac cal(const Frac &fa, const Frac &fb) {
LLL a = fa.a, b = fa.b, c = fb.a, d = fb.b;
LLL x = a/b+1;
if (x*d < c) return Frac(x, 1);
if (!a) return Frac(1, d/c+1);
if (a<=b && c<=d) {
Frac t = cal(Frac(d, c), Frac(b, a));
swap(t.a, t.b);
return t;
}
x = a/b;
Frac t = cal(Frac(a-b*x, b), Frac(c-d*x, d));
t.a += t.b*x;
return t;
}
```
+ 输出大于某分数，且分子分母小于 $n$ 的最小最简分数
```cpp
Frac gen(LL a, LL b, int n) {
pair l(0, 1), mid(1, 1), r(1, 0), res(-1, -1);
LL x=a, y=b;
while (x != y && max(mid.first, mid.second) <= n) {
if (a*mid.second < b*mid.first)
res = mid;
if (x < y) {
tie(l, mid, r) = make_tuple(l, make_pair(l.first+mid.first, l.second+mid.second), mid);
y -= x;
} else {
tie(l, mid, r) = make_tuple(mid, make_pair(mid.first+r.first, mid.second+r.second), r);
x -= y;
}
}
return Frac(res.first, res.second);
}
```
## 线性代数
### 线性基
+ 求第k大的异或子集
+ HDU 3949
```cpp
#include
using namespace std;
const int MAXN = 1e4+5;
typedef long long LL;
int n, m;
LL a[MAXN], q[MAXN];
const int MAXM = 63;
struct linear_basis {
int cnt;
bool zero;
LL p[MAXM], P[MAXN];
int pcnt;
static int highbit(LL x) {
for (int j=MAXM-1; j>=0; --j) {
if ((x>>j) & 1) return j;
}
return -1;
}
bool exist(LL x) {
for (int j=MAXM-1; j>=0; --j)
x = min(x, x^p[j]);
return x == 0;
}
void init(LL a[], int n) {
cnt = 0;
memset(p, -1, sizeof p);
zero = false;
for (int i=0; i
if (!a[i]) zero = true;
LL x = a[i];
for (int j=MAXM-1; j>=0; --j)
if ((x>>j) & 1) {
if (~p[j]) x ^= p[j];
else {
p[j] = x;
++cnt;
break;
}
}
if (x == 0) zero = true;
}
pcnt = 0;
for (int i=0; i
if (~p[i]) for (int j=i-1; ~j; --j) {
if (~p[j] && ((p[i]>>j) & 1))
p[i] ^= p[j];
}
if (~p[i]) P[pcnt++] = p[i];
}
// printf("%d\n", zero);
}
LL findkth(LL k) {
if (zero) --k;
if (k > (1LL<
LL res = 0;
for (int i=0; i
if ((k >> i) & 1) res ^= P[i];
return res;
}
} lb;
int main() {
int kk = 0;
int t;
scanf("%d", &t);
while (~scanf("%d", &n)) {
++kk;
printf("Case #%d:\n", kk);
for (int i=0; i
scanf("%lld", &a[i]);
}
scanf("%d", &m);
for (int i=0; i
scanf("%lld", &q[i]);
}
lb.init(a, n);
for (int i=0; i
LL x = lb.findkth(q[i]);
printf("%lld\n", x);
}
}
return 0;
}
```
线性基求交
```cpp
struct Linear_Basis{
//basis vector
ll basis[40];
bool droped = false;
void clear(){
memset(basis,0,sizeof basis);
}
void ins(ll x){
bool drop = true;
for (int i=31;i>=0;i--){
if (x & (1ll<< i)){
if (!basis[i]){basis[i] = x;drop = false;break;}
x ^= basis[i];
}
}
droped |= drop;
}
bool can(ll val){
for (int i=31;i>=0;i--){
if (val & (1ll << i)){
if (basis[i])val ^= basis[i];
else return false;
}
}
return val == 0;
}
void debug(){
cerr<
for (int i=31;i>=0;i--){
if (basis[i]){
cerr<
}
}
}
Linear_Basis (const Linear_Basis & other){
for (int i=0;i<=31;i++)basis[i] = other.basis[i];
}
Linear_Basis(bool full = false){
if (full){
for (int i=31;i>=0;i--){
basis[i] = 1ll << i;
}
}else clear();
}
Linear_Basis operator & (const Linear_Basis &other)const{
Linear_Basis ret = other;
Linear_Basis c,d;
for (int i=31;i>=0;i--){
d.basis[i] = 1ll << i;
}
for (int i=31;i>=0;i--){
if (basis[i]){
ll v = basis[i],k=0;
bool can = true;
for (int j=31;j>=0;j--){
if (v & (1ll << j)){
if (ret.basis[j]){
v ^= ret.basis[j];
k ^= d.basis[j];
}else{
can = false;
ret.basis[j] = v;
d.basis[j] = k;
break;
}
}
}
if (can){
ll v = 0;
for (int j=31;j>=0;j--){
if (k & (1ll << j)){
v ^= other.basis[j];
}
}
c.ins(v);
}
}
}
return c;
}
} dat[MAXN];
```
### 高斯消元
* n - 方程个数，m - 变量个数， a 是 n \* \(m + 1\) 的增广矩阵，free 是否为自由变量
* 返回自由变量个数，-1 无解
* 浮点数版本
```cpp
typedef double LD;
const LD eps = 1E-10;
const int maxn = 2000 + 10;
int n, m;
LD a[maxn][maxn], x[maxn];
bool free_x[maxn];
inline int sgn(LD x) { return (x > eps) - (x < -eps); }
int gauss(LD a[maxn][maxn], int n, int m) {
memset(free_x, 1, sizeof free_x); memset(x, 0, sizeof x);
int r = 0, c = 0;
while (r < n && c < m) {
int m_r = r;
FOR (i, r + 1, n)
if (fabs(a[i][c]) > fabs(a[m_r][c])) m_r = i;
if (m_r != r)
FOR (j, c, m + 1)
swap(a[r][j], a[m_r][j]);
if (!sgn(a[r][c])) {
a[r][c] = 0;
++c;
continue;
}
FOR (i, r + 1, n)
if (a[i][c]) {
LD t = a[i][c] / a[r][c];
FOR (j, c, m + 1) a[i][j] -= a[r][j] * t;
}
++r; ++c;
}
FOR (i, r, n)
if (sgn(a[i][m])) return -1;
if (r < m) {
FORD (i, r - 1, -1) {
int f_cnt = 0, k = -1;
FOR (j, 0, m)
if (sgn(a[i][j]) && free_x[j]) {
++f_cnt;
k = j;
}
if(f_cnt > 0) continue;
LD s = a[i][m];
FOR (j, 0, m)
if (j != k) s -= a[i][j] * x[j];
x[k] = s / a[i][k];
free_x[k] = 0;
}
return m - r;
}
FORD (i, m - 1, -1) {
LD s = a[i][m];
FOR (j, i + 1, m)
s -= a[i][j] * x[j];
x[i] = s / a[i][i];
}
return 0;
}
```
+ 数据
```
3 4
1 1 -2 2
2 -3 5 1
4 -1 1 5
5 0 -1 7
// many
3 4
1 1 -2 2
2 -3 5 1
4 -1 -1 5
5 0 -1 0 2
// no
3 4
1 1 -2 2
2 -3 5 1
4 -1 1 5
5 0 1 0 7
// one
```
一键复制
编辑
Web IDE
原始数据
按行查看
历史
1152
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
