[CodeChef] Adi and the Matrix

发表于 2019-03-08 更新于 2020-04-23 分类于题解

题目链接

有$N*M\leq 550$

注意到数据范围，它的含义是说，N和M中较小的那个不会超过23

我们考虑置换，首先假设$N\leq M$

枚举N的划分方案，注意要从小到大枚举，这个数不会太大

这个划分方案中每一个集合对应着一个循环

如果元素个数是$n$，划分方案为$(a_i,b_i)$，即长度为$a_i$的循环有$b_i$个

那么它对应的置换总数为

$\frac{n!}{\prod a_i^{b_i}b_i!}$

这个式子的意思是说，将所有元素进行排列，每一个循环中都有$a_i$种等价方案，即旋转0次、旋转1次，……，旋转$a_i-1$次，长度相同的循环之间可以互换顺序

假如我们已经知道了M对应的循环，那么考虑如何计算不动点个数

对于一个格子，它在行内的循环假设长度为$a$，在列内的循环假设长度为$b$。那么我们知道，它在经过$a$次置换后会回到原来所在的行，在经过$b$次置换后会回到原来所在的列

所以，最少需要$lcm(a,b)$次置换才能回到原点

也就是说，每一种方案，都有$lcm(a,b)$种方案是与它等价的

对于这两个循环中，我们知道一共有$ab$个格子

所以，不动点的个数是$\frac{ab}{lcm(a,b)}=gcd(a,b)$

每个不动点都有两种染色方案，即0或1

我们考虑使用dp来计算贡献

dp[i][j]表示还剩下i列没有分配，还没有分配循环长度j时所有方案的不动点个数之和

我们枚举长度为j的循环个数

求出之前枚举出来的$N$的划分中每一个$a_i$与j的$gcd$之和

注意循环对应的置换个数是$\frac{M!}{\prod a_i^{b_i}b_i!}$，这里我们把$M$提出来，在dp的时候顺便转移分母

代码如下

#include <bits/stdc++.h>
using namespace std;
#define N 600
typedef long long LL;
int dp[N][N], a[N], top, n, m, frac[N], inv[N], g[N][N], Pow2[N];
const int mod = 1e9 + 7;
inline int Pow(int x, int y) {
    int res = 1;
    for (;y;y >>= 1, x = (LL)x * x % mod) if (y & 1) res = (LL)res * x % mod;
    return res;
}
inline void Inc(int &x, int y) {x += y, x -= x >= mod ? mod : 0;}
int gcd(int a, int b) {return b ? gcd(b, a % b) : a;}
inline int solve() {
    memset(dp, 0, sizeof(dp)), dp[m][m] = 1;
    for (int i = m;i >= 0;i--)
        for (int j = m;j >= 1;j--)
            if (dp[i][j]) {
                int val = 0, tmp = dp[i][j];
                for (int k = 1;k <= top;k++) 
                    val += g[j][a[k]];
                for (int k = 0;j * k <= i;k++) {
                    Inc(dp[i - j * k][j - 1], (LL)tmp * inv[k] % mod);
                    tmp = (LL)tmp * inv[j] % mod * frac[j - 1] % mod * Pow2[val] % mod;
                }
            }
    return (LL)dp[0][0] * frac[m] % mod;
}
inline int calc() {
    // for (int i = 1;i <= top;i++) cout << a[i] << ' ';
    // cout << endl;
    int res = frac[n];
    for (int i = 1;i <= top;i++) {
        int cnt = 1;
        while (i <= top && a[i] == a[i + 1]) i++, cnt++;
        res = (LL)res * Pow(Pow(a[i], cnt), mod - 2) % mod * inv[cnt] % mod;
    }
    return (LL)res * solve() % mod;
}
int ans = 0;
void dfs(int n, int cur, int dep) {
    if (n < 0) return;
    if (!n) {top = dep, Inc(ans, calc()); return;}
    if (n < cur) return;
    dfs(n, cur + 1, dep);
    a[++dep] = cur, dfs(n - cur, cur, dep);
}
int main() {
    frac[0] = inv[0] = Pow2[0] = 1;
    for (int i = 1;i <= N - 10;i++) frac[i] = (LL)frac[i - 1] * i % mod, inv[i] = Pow(frac[i], mod - 2), Pow2[i] = Pow2[i - 1] * 2 % mod;
    for (int i = 0;i <= N - 10;i++)
        for (int j = 0;j <= i;j++)
            g[i][j] = g[j][i] = gcd(i, j);
    scanf("%d%d", &n, &m);
    if (n > m) swap(n, m);
    dfs(n, 1, 0);
    printf("%d\n", (LL)ans * inv[n] % mod * inv[m] % mod);
    return 0;
}

其中，inv[i]表示i的阶乘的逆元，g[i][j]表示$gcd(i,j)$，Pow2[i]表示$2^i$，frac[i]表示$i!$

那个inv[j] * frac[j - 1]等价于除以j

注意最后需要将答案除以置换总数$n!m!$

A – Necklace of Beads

Burnside定理模板

假设有n个条件$S$，每个条件都是形如“$A$与$B$等价”的形式，要求集合$X$在这些条件下的方案数

将条件分解成循环的形式，那么方案数

$=\frac{1}{n}\sum_{g\in S}|X^g|$

其中，$X^g$是置换$g$作用于集合$X$之后的不动点个数，即不变的元素的个数

那么对于这道题，我们考虑两类置换

1.旋转

假设有一条长度为n的项链，旋转之后相同被视为相同方案，那么显然地，我们有n种对应的置换

即不旋转，旋转1次，旋转2次，…，旋转n-1次

考虑计算不动点个数

假设现在旋转了2次

如果要求旋转之后不变（不动点），那么1号点、3号点、5号点的颜色必须相同

因为1号点转一次可以转到3号点，而3号点转一次可以转到5号点

如果这三个点的颜色不相同，那么旋转之后就变了，不能再称之为不动点

同理，2号点、4号点、6号点的颜色必须相同

总结一下：

对于长度为n的项链，旋转i次之后得到的不动点个数为$gcd(i,n)$

不动点之间互不影响，假设有m种颜色，那么此时的染色方案为$m^{gcd(i,n)}$

2.翻转

分奇偶考虑

先考虑奇数

一条对称轴必定穿过一个顶点，也就是说，一个顶点对应一条对称轴

如果要满足变换之后不变，显然2、5号点必须相等，3、4号点必须相等，1号点随意

如果有n个点，m种颜色，且n为奇数，那么这种置换一共有n个，每种置换有$m^{\frac{n+1}{2}}$种方案

考虑偶数的情况

稍微麻烦一点

1.对称轴穿过一条边的中点，不穿过点

显然，此时有$m^{\frac{n}{2}}$种方案，共$\frac{n}{2}$种置换

2.对称轴穿过两个点

其它点两两配对，这两个点随便染色

共$m^{\frac{n-2}{2}+2}$种方案，共$\frac{n}{2}$种置换

最后除以置换总数2n即可

B – Let it Bead

和A题相同，将颜色数从3改为m即可

C – Color

N颗珠子，N种颜色，而且$N\leq 10^9$

显然枚举旋转方案i不行

考虑枚举gcd

假设gcd=x，那么实际上就是要计算有多少个i满足$i\leq n$且$gcd(i,n)=x$

首先，x必须是n的因数，而且是i的因数

假设$n=x*a,i=x*b$，那么a与b互质

所以，i的个数为$\varphi(\frac{n}{x})$

x只有$\sqrt n$种取值

筛一下2e5之内的质数即可

D – Magic Bracelet

n颗珠子，m种颜色，k种限制，每种限制形如“颜色a与颜色b不能放在一起”

先不管限制

考虑旋转i步之后的状态

假设虚线中是一个循环节，显然，1号点和5号点的颜色应当相同，而与2、3、4、5号点的颜色无关

如果要求满足限制，那么在上图中的含义等价为“放5颗珠子，第1颗珠子与第5颗珠子颜色相同并且满足限制的方案数”

循环节长度为i，那么珠子数就为i+1

将这个过程想象成一张图

如果a与b不能相邻，那么点a到点b之间没有边，否则有一条双向边

最开始的图是一张完全图，每个条件会删去一条双向边

答案就等于在这张图上走i+1步，且起点与终点相同的方案数

这张图的邻接矩阵是一个01矩阵，那么走i+1步就是取这个矩阵的i次方

答案为对角线上的数之和

由于又是$n\leq 10^9$，所以需要矩阵乘法

剩下的和上一道题一样，先筛质数然后求$\varphi$即可

E – Who’s Aunt Zhang

一个三阶魔方，给每个面、每个角、每条棱上色，共n中颜色，将魔方整体旋转之后相同的视为等价情况，问方案数

首先面数+棱数+角数=54+12+8=74

有4类置换

1.不动

方案数$n^{74}$

有1种置换

2.以某个面的中心为轴旋转

有3种方案：转90度、转180度、转270度

注意到第一种方案和第三种方案要求轴所在的那一面的四条棱都相等，而第二种方案仅要求两条棱相等

对于第一种和第三种方案，答案为$n^{20}$

正对着的两个面共有$3*2$)种方案

中心有1种，棱上的面有1种，角上的面有1种，这样的面有2个，所以有6种

剩下的四个面染色方案必须一样，但每个面中的颜色独立，共9种

正对着的两个面每个面上的棱只有1种方案，共2种

剩下的四条棱有1种

角分为两组，坐标一组，右边一组，每组颜色必须相同，共2种

所以共$6+9+3+2=20$种方案

置换个数为$2*3=6$种（以某两个面为轴，转90度或270度）

对于第二种方案，答案为$n^{38}$

正对着的两个面上，将面两两分组，每组颜色必须相同

共$(2+2+1)*2=10$种（角上的面两种，棱上的面两种，中心一种，共两个这样的面）

其余的四个面两两配对，每个面上颜色独立，共$9*2=18$种

将棱两两配对，每对的颜色必须相同，共6种

将角两两配对，每对的颜色必须相同，共4种

总方案数为$10+18+6+4=32$种

置换有3种（将面两两配对，每一对只有1种置换）

以相对的两个角为中心旋转

在这种情况下，旋转120度和旋转240度是等价的

都有$n^{26}$种方案

对着的两个角特殊考虑，其余的角、边、棱三个为一组，每组颜色互不影响，共$2+\frac{74-2}{3}=26$种

共$2*4=8$种置换

以相对的两条棱的中心为轴旋转

选中的两条棱特殊考虑，其余元素两两配对

共$2+\frac{74-2}{2}=38$种方案

有$\frac{12}{2}=6$种置换

所以总置换数为24，最后除以24即可

F – Toy

前置题目轮状病毒

首先不考虑等价的情况

有两种方法

矩阵树定理+打表(oeis)

矩阵是度数矩阵-邻接矩阵，再随便去掉某一行和某一列

然后打表

递推

先不考虑连成环的情况

设$f(n)$表示除了中心点，还有n个点的生成树个数

那么有

$f(n)=\sum_{i=1}^{n}f(n-i)\*i$

含义是，假设这次选取的连通块有i个节点，剩下的点共有$f(n-i)$种方案，这个连通块可以选择任意一个点连向中心点，共i种方案

化简这个式子

$f(n)=\sum_{i=0}^{n-1}f(i)\*(n-i)$ $f(n-1)=\sum_{i=0}^{n-2}f(i)\*(n-i-1)$ $f(n)=f(n-1)+\sum_{i=0}^{n-1}f(i)$ $f(n-1)=f(n-2)+\sum_{i=0}^{n-2}f(i)$ $f(n)-f(n-1)=f(n-1)-f(n-2)+f(n-1)$ $f(n)=3f(n-1)-f(n-2)$

注意这个式子仅当$n\geq 3$时成立

考虑连成环的情况

就是1号点和n号点连了起来

假设这个连通块有i个点

那么有i-1种连通块可以选择（包含1-n这条边的连通块数），随便选一个点连向中心点

设方案数为$g(n)$

那么

$g(n)=\sum_{i=2}^{n}i\*(i-1)\*f(n-i)$ $=\sum_{i=1}^{n}i\*(i-1)\*f(n-i)$

总方案数

$F(n)=f(n)+g(n)$ $=\sum_{i=1}^{n}f(n-i)\*i^2$

我们有$f(0)=1,f(1)=1,f(2)=3$

$\because \sum_{i=n-2}^nf(n-i)\*i^2=3\*(n-2)^2+(n-1)^2+n^2$ $=3\*(n-1)^2+2\*(n-2)^2+2$ $\therefore F(n)=\sum_{i=1}^{n-3}f(n-i)\*i^2+3\*(n-1)^2+2\*(n-2)^2+2$ $\because n-i\geq 3$ $\therefore F(n)=\sum_{i=1}^{n-3}[3\*f(n-i-1)-f(n-i-2)]\*i^2+3\*(n-1)^2\*f(0)+(n-2)^2\*[3\*f(1)-f(0)]+2$ $=\sum_{i=1}^{n-2}[3\*f(n-i-1)-f(n-i-2)]\*i^2+3\*(n-1)^2\*f(0)+2$ $=3\*\sum_{i=1}^{n-1}f(n-i-1)\*i^2-\sum_{i=1}^{n-2}f(n-i-2)+2$ $=3\*F(n-1)-F(n-2)+2$

写成矩阵的形式就是

$\begin{bmatrix}3&-1&2\\\1&0&0\\\0&0&1\end{bmatrix}^{n-2}\*\begin{bmatrix}F(2)\\\F(1)\\\1\end{bmatrix}$

然后再求欧拉函数，枚举gcd即可

G – Birthday Toy

考虑n个位置，m种颜色，要求相邻颜色不同，并且首位颜色一样的方案数（$n,m\leq 10^9$）

先考虑朴素的dp:dp[i][j][k]表示到了第i个位置，这一位颜色为j，最开始的颜色为k的方案数

转移很好转移

考虑优化：我们并不关心当前位置的颜色是什么，我们只关心它与首位颜色相不相同，所以可以简化状态

dp[i][0/1]表示到了第i位，与当前首位颜色相同/不相同的方案数

那么有

1 2	dp[i][0]=dp[i-1][1] dp[i][1]=dp[i-1][0](m-1)+dp[i-2][1](m-2)

矩阵乘法加速即可

I – Leonardo’s Notebook

结论：

两个长度为n的相同循环相乘，当n为奇数是结果是一个长为n的循环，否则是两个长度为$\frac{n}{2}$的循环的乘积

所以长度为奇数的循环一定可以被拆成两个相同循环的乘积，长度为偶数的循环需要两两配对，如果能配对上就可以

K – Find the Permutations

下标与序列构成了一个置换

将置换分解，注意到每个循环所需要的交换次数是循环长度-1

加起来就是n-循环节个数

所以dp，dp[i][j]表示到第i个位置，一共有j个循环的方案数

i这个数可以新开一个循环，也可以插入到前面任意一个数的后面

1	dp[i][j]=dp[i-1][j-1]+dp[i-1][j]*(i-1)

L – Necklace

每种颜色的珠子有限制个数

还是枚举gcd，但是要求$\frac{n}{gcd}$必须是每种珠子个数的因数

然后可重排列

翻转有点毒瘤，考虑两种情况：n为奇数或是偶数

奇数比较简单，偶数又要分两种情况

注意细节，还是比较模板的

M – Cubes

L + E

莫比乌斯反演与杜教筛专题

发表于 2019-03-08 更新于 2020-04-23 分类于题解

A – GCD

莫比乌斯反演模板，注意去重

当$i\leq min(n,m)$时，每个gcd都被计算了两遍，除以2即可

B – CA Loves GCD

见莫比乌斯反演 CA Loves GCD

C – Hillan and the girl

首先转化成总数-gcd是完全平方数的个数

$\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)\in sqrt\ number]$

按照套路，先枚举gcd，此处是gcd的平方根

$=\sum_{d=1}^{\sqrt n}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d^2]$

莫比乌斯反演

$=\sum_{d=1}^{\sqrt n}\sum_{l=1}^{\lfloor\frac{n}{d^2}\rfloor}\mu(l)\lfloor\frac{n}{d^2l}\rfloor\lfloor\frac{m}{d^2l}\rfloor$

然后换T

$T=d^2l$

改写式子

$=\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{l|T}\mu(l)$

最后那个可以预处理

对于每一个$\mu$，枚举一个完全平方数更新，最后算一下前缀和即可

D – Trick GCD

很经典的一道题

题目要求$gcd\geq 2$，转化成求$gcd=1$的方案数，再用总数去减

设$f(x)$为gcd恰好为x的方案数，$g(x)$为gcd为x的倍数的方案数

显然有

$g(x)=\sum_{x|d}f(d)$ $f(x)=\sum_{x|d}g(d)\mu(\frac{d}{x})$

此时f中的x是1

现在关键是要求出g

假设当前是$g(x)$

对于每一个$A_i$，有$\lfloor\frac{A_i}{x}\rfloor$的贡献，注意到一段区间内的$A_i$贡献相同，所以对A整体分块计算贡献，然后乘起来就好了

一部分代码

for (int i = 1, a;i <= n;i++) scanf("%d", &a), buk[a]++, mn = min(mn, a), mx = max(mx, a);
for (int i = 1;i <= mx * 2;i++) buk[i] += buk[i - 1]; 
for (int i = 2;i <= mn;i++) {
    int tmp = -mu[i];
    for (int j = i;j <= mx;j += i) 
        tmp = (LL)tmp * Pow(j / i, buk[j + i - 1] - buk[j - 1]) % mod;
        ans = (ans + tmp) % mod;
}

E,F,G – GCD Extreme(x)

直接用$\varphi$替换gcd即可

经典题目

H – Zap

莫比乌斯反演模板

上界同时除以d即可

I – 数表

求

$\sum_{i=1}^{n}\sum_{j=1}^{m}\sigma(gcd(i,j))[\sigma(gcd(i,j))\leq a]$

先莫比乌斯反演一波

$=\sum_{i=1}^{n}\sum_{j=1}^{m}\sigma(gcd(i,j))[\sigma(gcd(i,j))\leq a]$ $=\sum_{d=1}^{n}\sigma(d)[\sigma(d)\leq a]\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]$ $=\sum_{d=1}^{n}\sigma(d)[\sigma(d)\leq a]\sum_{l=1}^{\lfloor\frac{n}{d}\rfloor}\mu(l)\lfloor\frac{n}{dl}\rfloor\lfloor\frac{m}{dl}\rfloor$

设

$T=dl$

改写式子

$=\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{d|T}\sigma(d)\mu(\frac{T}{d})[\sigma(d)\leq a]$

约数和这个函数是可以直接求的

考虑将询问按照限制大小从小到大排序

将约数和也从小到大排序

当一个约数和达到限制以下后，枚举它的倍数更新

由于计算要用到前缀和，使用树状数组维护即可

J – Crash的数字表格

见莫比乌斯反演

K – DZY Loves Math

神仙题

先来反演一波

$\sum_{i=1}^{n}\sum_{j=1}^{m}f(gcd(i,j))$ $=\sum_{d=1}^{n}f(d)\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]$ $=\sum_{d=1}^{n}f(d)\sum_{l=1}^{\lfloor\frac{n}{d}\rfloor}\mu(l)\lfloor\frac{n}{dl}\rfloor\lfloor\frac{m}{dl}\rfloor$

设

$T=dl$

则有

$=\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{d|T}f(d)\mu(\frac{T}{d})$

关键就是要求最后的那个东西

假设有$T=\prod p_i^{k_i}$，要想最后的那个$\mu$不为0，d的每一个质数的次数最少不会少于T中所对应的次数-1

也就是说，假如$f(T)=x$，那么$x-1\leq f(d)\leq x$

进一步考虑，最后可能影响答案的只有那些次数最高的项。如果T中质数的次数不全相等，那么对于$f(d)=x$的每一种方案，其它的项可以乱选，方案数是2的这些项的个数次方（要么次数不减，要么次数减1）

注意到$\mu(\frac{T}{d})$的值在上述的这两种情况中取值互为相反数，仅当T中质数的次数全部相等时，剩余项的方案数是$2^0$，在这一步中不能被抵消，其余情况均可两两抵消，这个$\sum$的值永远是0

考虑所有次数相等的情况

假设T的质因子个数为k

每个数仍然可以选或不选，$\mu$的和为0

但是对于所有质数都不选的情况，$f(d)$的值并不是x，而是x-1，所以答案要减去$-1*(-1)^k$，那个$(-1)^l$是$\mu(\frac{T}{d})$

所以此时的答案为

$(-1)^{k+1}$

假设当前要求的是

$\sum_{d|T}f(d)\mu(\frac{T}{d})$

那么：

假设存在x，使得$\mu(x)\not=0$，且T是x的k次方，那么$\sum_{d|T}f(d)\mu(\frac{T}{d})=-\mu(x)]$

否则这个值为0

代码如下

#include <bits/stdc++.h>
using namespace std;
#define N 10000010
#define LL long long
int prime[N], ncnt, mu[N], sum[N], flag[N];
inline void sieve() {
	mu[1] = 1;
	for (int i = 2;i <= N - 10;i++) {
		if (!flag[i]) prime[++ncnt] = i, mu[i] = -1;
		for (int j = 1;j <= ncnt && i * prime[j] <= N - 10;j++) {
			flag[i * prime[j]] = 1;
			if (i % prime[j] == 0) break;
			mu[i * prime[j]] = -mu[i];
		} 
	}
	for (int i = 2;i <= N - 10;i++)
		if (mu[i] != 0) for (LL j = i;j <= N - 10;j *= i) 
			sum[j] = -mu[i];
	for (int i = 2;i <= N - 10;i++) sum[i] += sum[i - 1]; 
} 
inline LL solve(int n, int m) {
	if (n > m) swap(n, m);
	LL ans = 0;
	for (int l = 1, r;l <= n;l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		ans += (LL)(n / l) * (m / l) * (sum[r] - sum[l - 1]);
	}
	return ans;
}
int main() {
	int T; scanf("%d", &T), sieve();
	while (T--) {
		int a, b; scanf("%d%d", &a, &b);
		printf("%lld\n", solve(a, b));
	}
}

L – Sum

杜教筛模板

见杜教筛

M – 能量采集

讲一讲求

$\sum_{i=1}^{n}\sum_{j=1}^{m}gcd(i,j)$

的两种求法

以前已经写过了，直接贴上来

第二种方法

可以看出选择合适的解法的重要性

N – DZY Loves Math VI

还是套路，先设出gcd，然后反演，最后换T

精雕细琢即可

简单的数学题

发表于 2019-03-08 更新于 2020-04-23 分类于题解

求

$\sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j)$ $n\leq 10^{10}$

根据$\varphi$的性质

$\sum_{i|d}\varphi(i)=d$

可以改写上面那个式子

$=\sum_{i=1}^{n}\sum_{j=1}^{n}ij\sum_{d|gcd(i,j)}\varphi(d)$

把d提到前面来

$=\sum_{i=1}^{n}\sum_{j=1}^{n}ij\sum_{d|gcd(i,j)}\varphi(d)$

为了方便快捷，设

$F(n)=\sum_{i=1}^{n}i$

那么原式就改写为

$=\sum_{d=1}^{n}\varphi(d)\*d^2\*F(\lfloor\frac{n}{d}\rfloor)^2$

注意到后面的那个F括号内只有$\sqrt n$种取值，直接分块

现在关键就是要求前面那个的前缀和

但是有$n\leq 10^{10}$，直接筛显然不行，所以上杜教筛

设

$S(n)=\sum_{i=1}^{n}\varphi(i)\*i^2$

现在关键就是要找一个函数g来与$f(i)=\varphi(i)*i^2$卷积

观察到有$\sum_{i|d}\varphi(i)=d$，所以中间的那个$\varphi$可以不管

但是那个$i^2$很烦，要弄掉这个

所以考虑使用$g(x)=x^2$这个函数来卷积

根据杜教筛的公式

$g(1)S(n)=\sum_{i=1}^{n}(g\*f)(i)-\sum_{i=2}^{n}g(i)S(\lfloor\frac{n}{i}\rfloor)$

这个$(g*f)(i)$其实就是

$\sum_{j=1}^{i}\varphi(j)\*j^2\*\frac{n^2}{j^2}=n^3$

而$g(1)=1^2=1$

所以

$S(n)=\sum_{i=1}^{n}i^3-\sum_{i=2}^{n}i^2S(\lfloor\frac{n}{i}\rfloor)$

有神仙公式

$\sum_{i=1}^{n}i^3=(\sum_{i=1}^{n}i)^2$

以及

$\sum_{i=1}^{n}i^2=\frac{1}{6}n\*(n+1)\*(2n+1)$

记忆化一波就可以了

杜教筛

发表于 2019-03-08 更新于 2020-04-23 分类于算法

假设有一道题：求

$\sum_{i=1}^{n}\mu(i)$

然后$n\leq 10^9$

线筛大家都知道怎么做

考虑把这个东西转化

设

$S(n)=\sum_{i=1}^{n}\mu(n)$

我们知道

$\sum_{i|d}\mu(i)=[d=1]$

如何利用这个性质呢？考虑用另一个函数来与这个$\mu$卷积

我们使用1来卷，因为没有什么要抵消的

$\sum_{i|d}\mu(i)\*1 = [d=1]=(\mu\*1)(d)$ $1=\sum_{d=1}^{n}\sum_{i|d}\mu(\frac{d}{i})\*1$ $=\sum_{i=1}^{n}\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}\mu(d)=\sum_{i=1}^{n}S(\lfloor\frac{n}{i}\rfloor)$

我们已经将这个东西转化成了$S$，移一下项

$1=\sum_{i=1}^{n}S(\lfloor\frac{n}{i}\rfloor)$ $1-\sum_{i=2}^{n}S(\lfloor\frac{n}{i}\rfloor)=S(n)$

注意到$\lfloor\frac{n}{i}\rfloor$的取值不超过$\sqrt n$个，直接分块即可

然后DFS

记得要记忆化！！！

一般式

假设要求一个积性函数$f(x)$的前缀和$S(n)$

先找到另一个积性函数$g(x)$，与f做卷积

这个g函数的要求有几点：

1.前缀和很好求，要不然找g之后反而更慢

2.与f的卷积的前缀和很好求

开始

$(g*f)(d)=\sum_{i|d}g(i)f(\frac{d}{i})$ $\sum_{i=1}^{n}(g*f)(i)=\sum_{i=1}^{n}\sum_{d|i}g(d)f(\frac{i}{d})$

枚举d

$=\sum_{d=1}^{n}g(d)\sum_{d|i}f(\frac{i}{d})$ $=\sum_{d=1}^{n}g(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)$ $=\sum_{d=1}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)$

移项

$g(1)S(n)=\sum_{i=1}^{n}(g*f)(i)-\sum_{i=2}^{n}g(i)S(\lfloor\frac{n}{i}\rfloor)$

然后开始分块即可

代码其实很短

map<int, int> m;
// sum[N - 10] 是预处理的前缀和，一般大小为3e6
inline int calc(int n) { // 算 sum mu
    if (n <= N - 10) return sum[n];
    else if (m.count(n)) return m[n]; // 记忆化
    int res = 1;
    for (LL l = 2, r;l <= n;l = r + 1) {
        r = n / (n / l);
        res -= (r - l + 1) * calc1(n / l);
    }
    return m1[n] = res;
}

莫比乌斯反演

发表于 2019-03-07 更新于 2020-04-23 分类于算法

一些定义

$\mu(i)$ 莫比乌斯函数

设$i=p_1^{k_1} p_2^{k_2} p_3^{k_3} \cdots p_m^{k_m}$

当$max(k) > 1$时$\mu(i)=0$

否则$\mu(i)=(-1)^{\sum k}$

$\varphi(i)$ 欧拉函数，为$[1,i]$内与i互质的数的个数

设$n=p_1^{k_1} p_2^{k_2} p_3^{k_3} \cdots p_m^{k_m}$

则$\varphi(n)=n*\prod{(1-\frac{1}{p_i})}$

关于狄利克雷卷积

定义二元关系$*$，$f*g $的狄利克雷卷积为

$(f*g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})$

卷积满足交换律

$f\*g=g\*f$

结合律

$(f\*g)\*h=f\*(g\*h)$

单位元

$e(n)=[n=1]$ $\displaystyle f*e=f$

逆元

$f*f^{-1}=e$

一些性质

$\sum_{d|n}\mu(d)=[n=1]$ $\displaystyle \sum_{d|n}\varphi(d)=n$

$\mu$与$\varphi$都是积性函数，可以线性筛

#define N 10000010
int flag[N], prime[N], mu[N], ncnt;
inline void sieve() {
    mu[1] = 1;
    for (int i = 2;i <= N - 10;i++) {
        if (!flag[i]) mu[i] = -1, prime[++ncnt] = i;
        //i是质数，mu[i]为-1
        for (int j = 1;j <= ncnt && i * prime[j] <= N;j++) {
            flag[i * prime[j]] = 1;
            //筛质数用
            if (i % prime[j] == 0) break;
            //此时根据定义，mu为0
            mu[i * prime[j]] = -mu[i];
            //mu[i * prime[j]] = mu[i] * mu[prime[j]] = -mu[i]
        }
    }
}

积性函数都可以线性筛

莫比乌斯反演：

定理1

若

$F(n)=\sum_{d|n}f(d)$

则

$f(n)=\sum_{d|n}F(d)\mu(\frac{n}{d})$

定理2

若

$F(n)=\sum_{n|d}f(d)$

则

$f(n)=\sum_{n|d}F(d)\mu(\frac{n}{d})$

证明

$F=f*1$

有

$\mu(n)\*1=\sum_{d|n}\mu(d)\*1=[n=1]=e$

所以

$\mu^{-1}=1$ $F=f\*1=f\*\mu^{-1}$ $F*\mu=f$ $f(n)=(F*\mu)(n)=\sum_{d|n}F(d)\mu(\frac{n}{d})$

一些套路

以Crash的数字表格为例

求

$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)$

套路1.设出$gcd$，然后按照$gcd$分类

$\sum_{k=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{i*j}{k}[gcd(i,j)=k]$

方法1.凑出$[gcd(i,j)=1]$然后使用$\mu$函数的性质换元

将i,j同时除以k，注意$\frac{i*j}{k}$需要乘回去

$\sum_{k=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}i\*j\*k[gcd(i,j)=1]$

有

$\sum_{d|n}\mu(d)=[n=1]$

所以

$[gcd(i,j)=1]=\sum_{d|gcd(i,j)}\mu(d)$

方法2.使用莫比乌斯反演基本公式

设

$f(x)=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{i*j}{x}[gcd(i,j)=x]$ $\displaystyle F(x)=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{i*j}{x}[x|gcd(i,j)]$

改写一下$F$

由于$gcd$是x的倍数，所以n,m都是x的倍数

同时除以x

$F(x)=\sum_{x|d}f(x)=\sum_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{x}\rfloor}i\*j\*x$

套用反演公式

$f(x)=\sum_{x|d}F(d)\mu(\frac{d}{x})$

得到

$f(x)=\sum_{x|d}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}i\*j\*d\*\mu(\frac{d}{x})$

代入原式即可

个人推荐推式子的时候用第一种方法，使用莫比乌斯反演容斥的时候使用第二种方法

从方法1的最后一步继续

$\sum_{k=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}i\*j\*k[gcd(i,j)=1]$ $=\sum_{k=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}i\*j\*k\sum_{d|gcd(i,j)}\mu(d)$

套路2.枚举d，并放到前面去

此时i,j都应该是d的倍数

同样，注意要乘回去

$=\sum_{k=1}^{n}k\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{kd}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{kd}\rfloor}i\*j\*d^2$

整理一下

$=\sum_{k=1}^{n}k\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)*d^2\sum_{i=1}^{\lfloor\frac{n}{kd}\rfloor}i\sum_{j=1}^{\lfloor\frac{m}{kd}\rfloor}j$

为了方便快捷，设

$F(n)=\sum_{i=1}^{n}i$

将原式替换成

$=\sum_{k=1}^{n}k\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\*d^2\*F(\lfloor\frac{n}{kd}\rfloor)\*F(\lfloor\frac{m}{kd}\rfloor)$

套路3.换T大法

为了提升时间复杂度，我们需要设一个神奇的T

设

$T=kd$

所以原式等于

$\sum_{k=1}^{n}k\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\*d^2\*F(\lfloor\frac{n}{T}\rfloor)\*F(\lfloor\frac{m}{T}\rfloor)$

枚举T，放到最前面

$\sum_{T=1}^{n}F(\lfloor\frac{n}{T}\rfloor)\*F(\lfloor\frac{m}{T}\rfloor)\sum_{d|T}\mu(d)\*\frac{T}{d}\*d^2$ $\sum_{T=1}^{n}F(\lfloor\frac{n}{T}\rfloor)\*F(\lfloor\frac{m}{T}\rfloor)\*T\sum_{d|T}\mu(d)\*d$

考虑优化最后那一坨的求法

直接求显然是$n\ log\ n$的

套路4.设出积性函数

观察到$\mu$以及$f(x)=x$（它还有个名字叫id）都是积性函数

两个积性函数的狄利克雷卷积以及直接对应相乘的结果都是积性函数

也就是说，这个东西是可以线性筛的

设

$g(x)=\sum_{d|x}\mu(d)*d$

对x分类

x是质数

此时显然有$g(x)=1-x$

x被表示成了$i*prime[j]$，其中$prime[j]\bot i$

此时根据积性函数性质，$g(iprime[j])=g(i)g(prime[j])$，直接乘即可

x被表示成了$i*prime[j]$，其中$prime[j]|i$

考虑乘之后多出来的因数

由于多出来的因数prime[j]这个质数的次数至少是2，而根据$\mu$的定义此时$\mu$为0，对g没有贡献，直接忽略即可

所以此时$g(i*prime[j])=g(i)$

套路5.数论分块

预处理完了g，我们再来说说最后答案怎么求

答案等于

$\sum_{T=1}^{n}F(\lfloor\frac{n}{T}\rfloor)\*F(\lfloor\frac{m}{T}\rfloor)\*T\*g(T)$

考虑从$\lfloor\frac{n}{T}\rfloor$,$\lfloor\frac{m}{T}\rfloor$ 开始优化

可以发现对于某一段的T，上面的两个式子的值都是不变的

这种不变的值的段数一共有$\sqrt n$段

考虑一段一段地枚举，而不是一个一个的枚举

此时变化的就只有$g(T)$

对$g(T)$求一下前缀和，算的时候直接求区间的和就好了

代码如下

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
#define LL long long
const LL mod=20101009;
int flag[10000010],prime[10000010],cnt,f[10000010],sum[10000010];
// sum是前缀和，f就是g
void sieve() {
    f[1]=sum[1]=1;
    for (int i=2;i<=10000000;i++) {
        if (!flag[i]) prime[++cnt]=i,f[i]=(1-i+mod)%mod;
        for (int j=1;j<=cnt&&i*prime[j]<=10000000;j++) {
            flag[i*prime[j]]=1;
            if (i%prime[j]==0) {
                f[i*prime[j]]=f[i]%mod;
                break;
            }
            f[i*prime[j]]=(LL)f[i]*f[prime[j]]%mod;
        }
        sum[i]=(sum[i-1]+(LL)f[i]*i%mod)%mod;
    }
}
LL zjk(LL n,LL m) {
    LL ans=0;
    if (n>m) swap(n,m);
    for (LL l=1,r;l<=n;l=r+1) {
        r=min(n/(n/l),m/(m/l));
        // 分块：[l,r]
        ans=(ans+((n/l*(n/l+1)/2)%mod*((m/l*(m/l+1)/2)%mod))%mod*(LL)(sum[r]-sum[l-1]))%mod;
    }
    return (ans+mod)%mod;
}
int main() {
    sieve();
    int n,m;scanf("%d%d",&n,&m);
    printf("%lld\n",zjk(n,m));
}

%%% ZJK

下面讨论一下使用莫比乌斯反演来容斥的题

CA Loves GCD

题意：给出n个[1,1000]的数，问所有选法中所选择的数的gcd的和

$n \leq 1000$

首先，还是对gcd分类，$f(i)$表示gcd为i的方案数

发现f好像很不好求，考虑设辅助函数

设$g(i)$表示gcd是i的倍数的方案数

很显然，在一个gcd为i的集合中，每一个数都是i的倍数

那么我们假设在这n个数中，i的倍数一共有k个

每个数可以选或者不选，一共有$2^k$种方案

减去所有数都不选的情况，就这一种

之后得到的就是$g(i)$的值

将g转化成f

显然有

$g(i)=\sum_{i|d}f(i)$

反演得

$f(i)=\sum_{i|d}g(d)\mu(\frac{i}{d})$

对于每一个i，枚举每一个i的倍数转移即可

对于如何求i的倍数的个数，可以开一个桶。先预处理每个数的约数，当输入一个数a时，就把a的所有约数在桶中加一，最后桶中的第i个位置就是i的倍数的个数

代码如下

#include <bits/stdc++.h>
using namespace std;
#define N 1010
vector<int> Div[N];
int flag[N], prime[N], ncnt, mu[N], Pow[N];
#define LL long long
const int mod = 1e8 + 7;
inline void sieve() {
    Pow[0] = 1, Pow[1] = 2, mu[1] = 1;
    for (int i = 2;i <= N - 10;i++) {
        Pow[i] = Pow[i - 1] * 2 % mod;
        if (!flag[i]) mu[i] = -1, prime[++ncnt] = i;
        for (int j = 1;j <= ncnt && i * prime[j] <= N - 10;j++) {
            flag[i * prime[j]] = 1; if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i];
        }
        for (int j = i;j <= N - 10;j += i) Div[j].push_back(i);
    }
    for (int i = 1;i <= N - 10;i++) Div[i].push_back(1);
}
int buk[N], ans[N], tmp[N];
int main() {
    int T; scanf("%d", &T), sieve();
    while (T--) {
        int n; scanf("%d", &n), memset(buk, 0, sizeof(buk)), memset(tmp, 0, sizeof(tmp));
        for (int i = 1, a;i <= n;i++) scanf("%d", &a), tmp[a]++;
        for (int i = 1;i <= 1000;i++)
            if (tmp[i]) for (int j = 0;j < Div[i].size();j++) 
                buk[Div[i][j]] += tmp[i];
        memset(ans, 0, sizeof(ans));
        for (int i = 1;i <= 1000;i++) 
            for (int j = i;j <= 1000;j += i) 
                ans[i] = (ans[i] + (Pow[buk[j]] - 1) * mu[j / i] + mod) % mod;
        int res = 0;
        for (int i = 1;i <= 1000;i++) 
            res = (res + (LL)i * ans[i] % mod) % mod;
        printf("%d\n", res);
    }
    return 0;
}

后缀数组

发表于 2019-03-07 更新于 2020-04-23 分类于算法

一些约定

sa[i] 排名为i的后缀在原串中的编号
rk[i] 原串中第i个后缀所对应的排名
height[i] 排名为i的后缀与排名为i – 1的后缀的最长公共前缀

一些例子

字符串S = bacaab 中，sa，rk，height分别为

一些解释

在S中，所有的后缀分别为bacaab, acaab, caab, aab, ab, b

将这些串按照字典序排序，可以得到aab, ab, acaab, b, bacaab, caab

对应出来的rk值就是5, 3, 6, 1, 2, 4

根据上面的定义可以得出sa[rk[i]] = i, rk[sa[i]] = i

sa[1] = 4 对应后缀aab，sa[2] = 5 对应后缀ab

因此height[2] = lcp(sa[1], sa[2]) = lcp(aab, ab) = 1

注：height[1] 没有意义，约定其等于0

倍增

直接暴力求sa显然不行，由于涉及到后缀排序，所以考虑倍增

假设对于每个后缀的前k个字符，我们已经排好了序，考虑进一步的转移

如果直接比较每个后缀的第k + 1个字符，忽视了太多已经求出的信息

对于两个后缀s1 = S[i, …], s2 = S[j, …]，如果它们的前k个字母不全相同，那么我们在之前的比较中就已经确定好了这两个后缀的顺序

如果这两个后缀的前k个字母都相同，我们这次就比较s1[k + 1, 2k], s2[k + 1, 2k]的大小

显然，这两个串的长度都是k，这意味着我们已经比较过这两个串

图中括号中的数字表示当前层中此串的排名，注意相同串的排名相同。

第一次比较将串s1与s2拼在一起，将串s2与s3拼在一起…

第二次比较将串s1与s3拼在一起，将串s2与s4拼在一起…

第三次比较将串s1与s5拼在一起，将串s2与s6拼在一起…

直到不同的串的个数为原串长度

那么问题来了，如何高效地比较？

首先利用好已得出的信息，将这次比较前各个串的排名作为第一关键字

在这次比较中，每个字符串都将添加k个字符，对于串s[i, …]来说，将串s[i + k, …]的排名作为第二关键字

注意空串的字典序最小，假设原串长度为n，注意到在n – k + 1到n位置的第二关键字都是空串，我们约定此时n位置的第二关键字最小，其次是n – 1，然后以此类推

for (int i = n;i >= n - k + 1;i--) y[++num] = i; 
// 添加的串为空
for (int i = 1;i <= n;i++) if (sa[i] > k) y[++num] = sa[i] - k; 
// 此时的i就是串sa[i] - k的第二关键字，注意枚举的i代表排名
// sa[i] <= k时没有串会在后面添加串sa[i]，所以无视
// y[i] 表示排名为i的第二关键字对应的是哪个串

对此时的所有串重新排序，第一关键字小的排前面，如果第一关键字相同，则第二关键字小的排前面

不能直接$n \log n$排序，这样的总时间复杂度与直接用哈希构造无异，因此选择桶排序

我们开一个桶c，将所有的第一关键字放进去，然后做一遍前缀和，可以得出每个元素的最大排名

由于第一关键字相同时，第二关键字较小的排名较前，因此我们从大到小枚举第二关键字，更新每个桶此时的最大排名和当前的sa数组

for (int i = 1;i <= m;i++) c[i] = 0; 
// 先清空， m是当前第一关键字大小
for (int i = 1;i <= n;i++) c[x[i]]++; 
// 将每个第一关键字加入桶
for (int i = 1;i <= m;i++) c[i] += c[i - 1]; 
// 获得每个值的最大排名
for (int i = n;i >= 1;i--) sa[c[x[y[i]]]--] = y[i]; 
// 第一关键字相同时，第二关键字越大排名越靠后

最后，我们需要更新此时的第一关键字

用一个num表示当前不同元素个数，仅当sa[i], sa[i – 1]至少有一个关键字不同时，num才会累加1，更新x数组

用m更新num

memcpy(y, x, sizeof(x)), num = x[sa[1]] = 1;
// 这个memcpy只是临时存一下，后面那个是初值
for (int i = 2;i <= n;i++)
    x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k]) ? num : ++num; 
// 更新此时的第一关键字
if (num == n) break; else m = num; 
// 更新第一关键字的最大值
// 如果m等于n，则代表已经分清了所有后缀的排名，没有必要继续比较下去

最后总代码如下

inline void GetSA(int n, char s[]) {
    for (int i = 1;i <= n;i++) c[x[i] = s[i]]++;
    for (int i = 1;i <= m;i++) c[i] += c[i - 1];
    for (int i = n;i >= 1;i--) sa[c[x[i]]--] = i;
    for (int k = 1, num = 0;k <= n;k <<= 1, num = 0) {
        for (int i = n;i >= n - k + 1;i--) y[++num] = i;
        for (int i = 1;i <= n;i++) if (sa[i] > k) y[++num] = sa[i] - k;
        for (int i = 1;i <= m;i++) c[i] = 0;
        for (int i = 1;i <= n;i++) c[x[i]]++;
        for (int i = 1;i <= m;i++) c[i] += c[i - 1];
        for (int i = n;i >= 1;i--) sa[c[x[y[i]]]--] = y[i];
        memcpy(y, x, sizeof(x)), num = x[sa[1]] = 1;
        for (int i = 2;i <= n;i++)
            x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k]) ? num : ++num;
        if (num == n) break; else m = num;
    }
}

关于DC3

时间复杂度$O(n)$，但是实现较为繁琐，不推荐

~~其实是我不会233333~~

关于height

$LCP$ 的几条定理

定义$LCP(i, j)$为排名为i的后缀与排名为j的后缀的最长公共前缀长度

根据height数组的定义，$height[i] = LCP(i – 1, i)$

$LCP(i, j) = LCP(j, i)$ // 显然

$LCP(i, i) = LCP(i, i)$ // 同上

$LCP(i, j) = \min\{LCP(i, k), LCP(k, j)\}$ // 对于任意i ≤ k ≤ j，即$LCP$具有传递性

设$h[i] = height[rk[i]]$，即$height[i] = h[sa[i]]$，有$h[i] ≥ h[i – 1] – 1$

第4条定理的证明：

代码如下

inline void Getheight(int n, char s[]) {
    for (int i = 1;i <= n;i++) rk[sa[i]] = i;
// 根据定义，给rk赋值
    for (int i = 1, k = 0;i <= n;i++) {
// k表示h[i - 1]，注意i是编号
        if (rk[i] == 1) continue; if (k) k--;
// height[1]没有意义，且h[i] >= h[i - 1] - 1
        int at = sa[rk[i] - 1];
        while (at + k <= n && i + k <= n && s[at + k] == s[i + k]) k++;
// 暴力向后匹配
        height[rk[i]] = k;
    }
}

关于最长公共前缀

由$LCP$的第三条定理可以快速确定任意两个后缀的最长公共前缀，具体地

$LCP(rk[i], rk[j]) = \min\{height[rk[i] + 1], height[rk[i] + 2], \cdots ,height[rk[j]]\}$

使用ST表查找区间最小值即可

inline void GetST(int n) {
    memset(mn, 0x3f, sizeof(mn));
    for (int i = 2;i <= n;i++) lg[i] = lg[i >> 1] + 1;
    for (int i = 1;i <= n;i++) mn[0][i] = height[i]; 
    for (int i = 1;i <= 20;i++)
        for (int j = 1;j + (1 << i) - 1 <= n;j++) 
            mn[i][j] = min(mn[i - 1][j], mn[i - 1][j + (1 << i - 1)]); 
}
inline int GetLCP(int L, int R) { // rank 为 L, R 
    if (L > R) swap(L, R);
    L++; int len = R - L + 1;
    if (!len) return 1e9;
    return min(mn[lg[len]][L], mn[lg[len]][R - (1 << lg[len]) + 1]);
}