第二类Stirling数实际上是集合的一个拆分，表示将n个不同的元素拆分成m个集合的方案数
LL Stirling2(int n, int m) {
LL ans = 0;
for (int k=0; k<=m; ++k) {
LL del = C(m, k) * bin(m-k, n, MOD);
if (k % 2) {
ans = (ans - del) % MOD + MOD;
} else {
ans = (ans + del) % MOD;
}
}
return ans;
}- 整除
- 素因子分解
- 因数分解
namespace math {
int mu[N], prime[N], pcnt, phi[N];
int minp[N];
int smu[N];
LL sphi[N];
void init() {
mu[1] = 1;
for (int i = 2; i < N; ++i) {
if (!minp[i]) {
minp[i] = i;
prime[pcnt++] = i;
mu[i] = -1;
}
for (int j = 0; j < pcnt; ++j) {
LL nextp = 1LL * i * prime[j];
if (nextp >= N) break;
if (!minp[nextp]) minp[nextp] = prime[j];
if (i % prime[j] == 0) {
mu[nextp] = 0;
break;
} else {
mu[nextp] = -mu[i];
}
}
}
for (int i = 1; i < N; ++i) {
smu[i] = smu[i - 1] + mu[i];
sphi[i] = sphi[i - 1] + phi[i];
}
}
template <typename T>
vector<pair<T, int>> primeFac(T x) {
vector<pair<T, int>> res;
if (x < N)
while (x > 1) {
int cnt = 0, p = minp[x];
while (x % p == 0) {
++cnt;
x /= p;
}
res.emplace_back(p, cnt);
}
else {
for (int i = 0; i < pcnt; ++i) {
int p = prime[i];
if (1LL * p * p > x) break;
if (x % prime[i] == 0) {
int cnt = 0;
while (x % prime[i] == 0) {
++cnt;
x /= prime[i];
}
res.emplace_back(p, cnt);
}
}
if (x > 1) res.emplace_back(x, 1);
}
return res;
}
template <typename T>
vector<T> fac(T x) {
vector<T> res;
for (int i = 1; 1LL * i * i <= x; ++i) {
if (x % i == 0) {
res.push_back(i);
if (x / i != i) res.push_back(x / i);
}
}
return res;
}
int getmu(int x) {
if (x < N)
return mu[x];
int res = 1;
auto vp = primeFac(x);
for (const auto &p : vp) {
if (p.second > 1)
return 0;
else
res = -res;
}
return res;
}
} // namespace math- Weaver_zhu 版本，质因数个数
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<pc; ++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<pc; ++j) {
LL p = prime[j];
for (int i=0; i<m; ++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);
}
}求
构造一个积性函数
当然，要能够由此计算
-
$f*g$ 要能够快速求前缀和。 -
$g$ 要能够快速求分段和（前缀和）。 - 对于正常的积性函数
$g(1)=1$ ，所以不会有什么问题。
在预处理
- 前置： 快速乘、快速幂
- 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/
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<LL> 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;
}- 大整数分解
mt19937 mt(time(0));
LL pollard_rho(LL n, LL c) {
LL x = uniform_int_distribution<LL>(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);
}- 杜教版
namespace linear_seq {
const int N=10010;
ll res[N],base[N],_c[N],_md[N];
vector<int> 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<<pnt)<=n) pnt++;
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)<SZ(B)+m) C.pb(0);
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)<SZ(B)+m) C.pb(0);
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$ 可能是负数
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;
}-
$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$使得
void solve(LL p1,LL q1,LL p2,LL q2,LL &x,LL &y)
{
LL l=p1/q1+1;
if (l*q2<p2){x=l; y=1;return;}
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$ 不是素数，使用拓展欧几里得 -
前置模板：快速幂 / 扩展欧几里得
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 的逆元
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;
}- 如果数较小，模较大时使用逆元
- 前置模板：逆元-预处理阶乘及其逆元
- 卢卡斯定理
- 任意模数组合数（很慢慎用，基于CRT）
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：模数不固定下使用，无法单独使用。
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;
}LL Lucas(LL n, LL m) { // m >= n >= 0
return m ? C(n % MOD, m % MOD) * Lucas(n / MOD, m / MOD) % MOD : 1;
}- 前置模板：素数，CRT
bool nop[MAXN];
vector<int> prime;
void get_prime() {
// ...
}
vector<pair<int, int> > 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) {
// ...
}
LL bin(LL a, LL b, LL p) {
// ...
}
LL get_inv(LL a, LL M) {
// ...
}
LL mo[MAXN], r[MAXN];
LL calcT(LL n, pair<int, int> 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) {
// ...
}
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);
}- 绝对值是
$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)$
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;已经退出历史舞台，详见多项式大餐
$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的二进制子集
template<typename T>
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
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<n; d*=3) {
for (int x=0; x<n; x+=3*d) {
for (int i=x; i<x+d; ++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<n; ++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 子集卷积
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]
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)); }- O(1)
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;
}- 前置模板：素数筛，快速幂，分解质因数
- 要求 p 为质数
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<i} \varphi(i)\overset{t=\frac id}{=} \frac {i(i+1)}{2} - \sum_{t=2}^n S(\frac n t)$ - 利用
$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$。且
若
设
**但是：**佩尔方程千万不要去推（虽然推起来很有趣，但结果不一定好看，会是两个式子）。记住佩尔方程结果的形式通常是 $a_n=ka_{n−1}−a_{n−2}$（$a_{n−2}$ 前的系数通常是 $−1$）。暴力 / 凑出两个基础解之后加上一个
$|X/G|={\frac {1}{|G|}}\sum _{{g\in G}}|X^{g}|$
注：$X^g$ 是
$|Y^X/G| = \frac{1}{|G|}\sum_{g \in G} m^{c(g)}$
注：用
-
$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
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
- 前置模板：扩展欧几里得
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$ 。(心态崩了)
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(k^2)$ ，如果插值在$[0, k]$ 复杂度为$O(k)$ -
$k$ 次多项式需要$k+1$ 个插值 - 注意 除法逆元去模后为 0 的情况！！！（心态崩了
namespace interlopation {
LL go(LL y[], int n, LL x) {
if (x < n) return y[x];
x %= MOD;
LL tot = 1;
int hit = 0;
for (int i=0; i<n; ++i) {
if (x == i) {
hit = 1;
} else {
tot = tot * (x + MOD -i) % MOD;
}
}
LL ans = 0;
// printf("%lld %d\n", tot, hit);
for (int i=0; i<n; ++i) {
LL del = 1;
if (hit) {
if (i == x) del = tot;
else del = 0;
} else {
del = tot * get_inv(x+MOD-i, MOD) % MOD;
}
int sgn = (n-i-1) % 2 ? -1 : 1;
del = (del * sgn + MOD) * invf[n-i-1] % MOD * invf[i] % MOD;
ans = (ans + del * y[i]) % MOD;
}
return ans;
}
}- uoj 97分
simplex::init()和simplex::solve()
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<n; ++i) ans[i] = ans[i+1];
}
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];
}
}
- 针对多次询问，复杂度
$O(\sqrt{qn})$
namespace BSGS {
LL a, p, m, n, q;
unordered_map<LL, LL> 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;
}
}- 模数可以非素数
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<LL, LL> 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;
}for (LL l = 1, v, r; l <= N; l = r + 1) {
v = N / l; r = N / v;
}$O(sqrt(n))$
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
- 对局面进行离散化，然后严格按照 sg 定理（mex函数）暴力求出 sg 函数找规律。然后尝试周期，注意周期应该忽略开头一小段（开头的用dp数组，太大的掐掉开头的按周期找）
- MTT基于毛子黑科技的任意模数fft，链接https://blog.csdn.net/lvzelong2014/article/details/80156989
- 加减乘除，模数乘
- 求逆
- 求线性递推式$O(klogklogn)$
- bzoj4161
#include <cstdio>
#include <cmath>
#include <cctype>
#include <algorithm>
#include <vector>
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<int> 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("");
}注意取模的时候多项式为 0 的情况。
这种实现的方法多项式为空的表示是 vector 里面只有一个 0
inline Poly operator%(const Poly &lhs, const Poly &rhs) {
if (sz(lhs) < sz(rhs))
return lhs;
if (sz(rhs) == 1)
return Poly(1, 0);
Poly res = lhs;
while (sz(res) >= sz(rhs)) {
Poly r(sz(res) - sz(rhs), 0);
r.insert(r.end(), rhs.begin(), rhs.end());
int c = 1LL*res.back() * bin(r.back(), MOD-2, MOD) % MOD;
for (int &x : r) {
x = 1LL*x*c%MOD;
}
res = res - r;
normalize(res);
dbg(sz(r), sz(rhs), sz(res));
}
return res;
}
Poly pow(Poly a, int b, const Poly p) {
Poly res = {1};
for (a = a % p; b; b >>= 1, a = a * a % p) {
dbg(sz(a), sz(p));
myassert(sz(p) != sz(a));
if (b & 1) res = res * a % p;
}
return res;
}
Poly gcd(Poly a, Poly b) {
while (sz(a) > 1 && a.back() == 0) a.pop_back();
while (sz(b) > 1 && b.back() == 0) b.pop_back();
while (!(sz(b) == 1 && b.back() == 0)) {
Poly c = a % b;
while (sz(c) > 1 && c.back() == 0) c.pop_back();
a = b;
b = c;
}
return a;
}- 输出两分数之间最小分母的最简分数
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$ 的最小最简分数
Frac gen(LL a, LL b, int n) {
pair<LL, LL> 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
#include <bits/stdc++.h>
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<n; ++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<MAXM; ++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<<pcnt)-1) return -1;
LL res = 0;
for (int i=0; i<MAXM; ++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<n; ++i) {
scanf("%lld", &a[i]);
}
scanf("%d", &m);
for (int i=0; i<m; ++i) {
scanf("%lld", &q[i]);
}
lb.init(a, n);
for (int i=0; i<m; ++i) {
LL x = lb.findkth(q[i]);
printf("%lld\n", x);
}
}
return 0;
}线性基求交
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<<"--------"<<endl;
for (int i=31;i>=0;i--){
if (basis[i]){
cerr<<basis[i]<<endl;
}
}
}
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 无解
-
浮点数版本
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