作者: KSkun

[APIO2014]回文串 题解

[APIO2014]回文串 题解

题目地址:洛谷:【P3649】[APIO2014]回文串 – 洛谷、BZOJ: 

[洛谷4142]洞穴遇险 题解

[洛谷4142]洞穴遇险 题解

题目地址:洛谷:【P4142】洞穴遇险 – 洛谷 题目描述 前方有一片沼泽地. 

[ZJOI2009]对称的正方形 题解

[ZJOI2009]对称的正方形 题解

题目地址:洛谷:【P2601】[ZJOI2009]对称的正方形 – 洛谷、BZOJ:Problem 1414. — [ZJOI2009]对称的正方形

题目描述

Orez很喜欢搜集一些神秘的数据,并经常把它们排成一个矩阵进行研究。最近,Orez又得到了一些数据,并已经把它们排成了一个n行m列的矩阵。通过观察,Orez发现这些数据蕴涵了一个奇特的数,就是矩阵中上下对称且左右对称的正方形子矩阵的个数。 Orez自然很想知道这个数是多少,可是矩阵太大,无法去数。只能请你编个程序来计算出这个数。

题意简述

求上下左右对称的正方形子矩阵个数。

输入输出格式

输入格式:
文件的第一行为两个整数n和m。接下来n行每行包含m个正整数,表示Orez得到的矩阵。

输出格式:
文件中仅包含一个整数answer,表示矩阵中有answer个上下左右对称的正方形子矩阵。

输入输出样例

输入样例#1:

5 5
4 2 4 4 4 
3 1 4 4 3 
3 5 3 3 3 
3 1 5 3 3 
4 2 1 2 4 

输出样例#1:

27

说明

数据范围
对于30%的数据 n,m≤100
对于100%的数据 n,m≤1000 ,矩阵中的数的大小≤10^9

题解

参考资料:题解 P2601 【[ZJOI2009]对称的正方形】 – Miaomiao’s Home – 洛谷博客
本题可以先对矩阵元素间插入数字0以便处理奇偶回文的情况,对每行每列跑一遍Manacher,统计出回文半径。
对于每个格子,上下左右找一遍以该格子为顶点的区间回文半径不大于区间长度的最大区间,这些区间长度的最小值就是这个格子为中心的最大对称正方形,而消除插入0的影响后的值就是这个格子对答案的贡献。
思路很简单也很暴力,算法$O(n^2 \log n)$,但是写出来常数挺大的,而且其实不太好写。

代码

// Code by KSkun, 2018/6
#include <cstdio>
#include <cctype>
#include <cstring>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2)
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1; register char c = fgc();
    for(; !isdigit(c); c = fgc()) if(c == '-') neg = -1;
    for(; isdigit(c); c = fgc()) res = (res << 1) + (res << 3) + c - '0';
    return res * neg;
}

inline int min(int a, int b) {
    return a < b ? a : b;
}

const int MAXN = 2005;

int n, m, a[MAXN][MAXN], lh[MAXN][MAXN], lv[MAXN][MAXN], t[MAXN];

inline void manacher(int n, int *arr, int *len) {
    int id = 0, mx = 0;
    for(register int i = 1; i <= n; i++) {
        if(i < mx) len[i] = min(len[(id << 1) - i], mx - i);
        else len[i] = 1;
        while(arr[i - len[i]] == arr[i + len[i]] 
            && i - len[i] >= 1 && i + len[i] <= n) len[i]++;
        if(mx < i + len[i]) {
            id = i; mx = i + len[i];
        }
    }
}

int Log[MAXN], mn[MAXN][12], ans[MAXN][MAXN];

inline void initst(int n, int x, int arr[][MAXN]) {
    for(register int i = 1; i <= n; i++) {
        mn[i][0] = arr[i][x] - 1;
    }
    for(register int j = 1; j < 12; j++) {
        for(register int i = 1; i <= n; i++) {
            if(i + (1 << j) - 1 > n) break;
            mn[i][j] = min(mn[i][j - 1], mn[i + (1 << (j - 1))][j - 1]);
        }
    }
}

inline int query(int x, int y) {
    int lg = Log[y - x + 1];
    return min(mn[x][lg], mn[y - (1 << lg) + 1][lg]);
}

int main() {
    for(register int i = 2; i < MAXN; i++) {
        Log[i] = Log[i >> 1] + 1;
    }
    n = readint(); m = readint();
    for(register int i = 1; i <= n; i++) {
        for(register int j = 1; j <= m; j++) {
            a[(i << 1) - 1][(j << 1) - 1] = readint();
        }
    }
    n = (n << 1) - 1; m = (m << 1) - 1;
    for(register int i = 1; i <= n; i++) {
        for(register int j = 1; j <= m; j++) {
            t[j] = a[i][j];
        }
        manacher(m, t, lh[i]);
    }
    for(register int i = 1; i <= m; i++) {
        for(register int j = 1; j <= n; j++) {
            t[j] = a[j][i];
        }
        manacher(n, t, lv[i]);
    }
    memset(ans, 0x3f, sizeof(ans));
    for(register int i = 1; i <= n; i++) {
        initst(m, i, lv);
        int t = 1;
        for(register int j = 1; j <= m; j++) {
            while(t < j && query(t, j) < j - t) t++;
            ans[i][j] = min(ans[i][j], j - t);
        }
        t = m;
        for(register int j = m; j >= 1; j--) {
            while(t > j && query(j, t) < t - j) t--;
            ans[i][j] = min(ans[i][j], t - j);
        }
    }
    for(register int i = 1; i <= m; i++) {
        initst(n, i, lh);
        int t = 1;
        for(register int j = 1; j <= n; j++) {
            while(t < j && query(t, j) < j - t) t++;
            ans[j][i] = min(ans[j][i], j - t);
        }
        t = n;
        for(register int j = n; j >= 1; j--) {
            while(t > j && query(j, t) < t - j) t--;
            ans[j][i] = min(ans[j][i], t - j);
        }
    }
    int res = 0;
    for(register int i = 1; i <= n; i++) {
        for(register int j = 1; j <= m; j++) {
            if((i & 1) && (j & 1)) res += (ans[i][j] >> 1) + 1;
            else if(!(i & 1) && !(j & 1)) res += (ans[i][j] + 1) >> 1;
        }
    }
    printf("%d", res);
    return 0;
}
[CF916E]Jamie and Tree 题解

[CF916E]Jamie and Tree 题解

题目地址:Codeforces:Problem – 916E –  

[WC2010]重建计划 题解

[WC2010]重建计划 题解

题目地址:洛谷:【P4292】[WC2010]重建计划 – 洛谷、BZOJ:P 

后缀自动机原理及实现

后缀自动机原理及实现

概述

后缀自动机(Suffix Automaton,简称SAM)是一种用于字符串处理的有限状态自动机(DFA),它可以接受其母串的任何一个后缀,且构造算法为线性算法。下面对其原理及实现进行介绍。
本文中所有SAM示意图皆为SAM Builder生成的图片,构造法处略图来自Menci博客。本文部分内容有摘录参考资料或摘录并改编的部分,参考资料均已附链接。

结构

abc串的SAM示意图
如图为串abc的SAM示意图,其中,实线黑色箭头表示节点的转移边,虚线蓝色箭头表示后缀链接。

有向无环单词图(Directed Acyclic Word Graph,简称DAWG)

DAWG是SAM的一部分,对于一个DAWG,每个节点(除起始节点外)都可以代表一个或多个子串(的结束点),每条转移边上都有一个字符,从起始节点沿转移边走,每条路径都对应一个子串,即路径上所有边的字符相连构成的。
SAM是一个DFA,它可以接受一些字符串,可以接受的字符串是指从起始节点开始,沿着每一个字符对应的转移边在DAWG上转移,最后到达某个节点的字符串。显然,可以接受的字符串一定是母串的一个子串。而如果这个字符串是母串的一个后缀,则称它的终止节点为可接受的节点。

节点性质

  • 每个节点所表示的所有字符串,一定是母串的某个前缀的若干个长度相差1的后缀。
    例如,对于字符串abcabcabc,它的一个前缀是abcab,某个节点所表示的子串可能是cabbcababcab
  • 节点$v$所表示的所有字符串中,定义最短的子串为$\min(v)$,最长的子串为$\max(v)$。
    例如,对于字符串abcabcabc,某个节点$v$所表示的子串可能是cabbcababcab,因此,$\min(v) = \mathtt{cab}, \max(v) = \mathtt{abcab}$。
  • 节点$v$所表示的最长字符串在母串中的所有出现位置的结尾下标构成该节点的$\mathrm{endpos}(v)$集合。
    例如,对于字符串abcabcabc,某个节点$v$所表示的子串可能是cabbcababcab,因此,$\mathrm{endpos}(v) = \{ 5, 8 \}$。
  • 任意两个节点的$\mathrm{endpos}(v)$集合不同。
    如果存在两个节点的$\mathrm{endpos}(v)$集合相同,则说明这两个节点的出边是等价的,这两个节点可以合并使得DAWG更简化。
  • 只有表示整个母串的节点沿后缀链接向上到起始节点的路径上的点的$\mathrm{endpos}(v)$集合包含母串的末尾,所有只有这些节点才是接受节点,分别代表母串的一个或一些后缀。

转移边

节点$u$到$v$有转移边,表示$u$代表的所有子串加上转移边对应的字符后,得到的集合是$v$代表的子串集合的子集,但是,并不保证两个集合相等。

后缀链接与前缀树(又称parent树)

  • $u$的后缀链接指向$v$,记作$\mathrm{next}(u)=v$,当且仅当$|\min(u)|=|\max(v)|+1$且$v$代表的子串均为$u$子串的后缀。
    例如,对于字符串abcabcabc,一个节点$u$可能代表子串abcabcaca,一个节点$v$可能代表子串a,则存在后缀链接$\mathrm{next}(u)=v$。
  • 任意节点沿后缀链接的方向走,最终都会走到DAWG的起始节点,以后缀链接为边,所有的节点构成了一棵向根的有向树,称为前缀树(parent树),DAWG的起始节点是前缀树的根。
  • 前缀树中,子节点的$\mathrm{endpos}(v)$集合一定是其父节点的真子集,即$\mathrm{endpos}(v) \subsetneq \mathrm{endpos}(\mathrm{next}(v))$。
    因为$\max(\mathrm{next}(v))$一定是$\max(v)$的后缀,所以$\max(v)$出现的地方也会出现$\max(\mathrm{next}(v))$,所以$\mathrm{endpos}(v) \subseteq \mathrm{endpos}(\mathrm{next}(v))$。如果$\mathrm{endpos}(v) = \mathrm{endpos}(\mathrm{next}(v))$,则两个点等价可以被合并,因此$\mathrm{endpos}(v) \subsetneq \mathrm{endpos}(\mathrm{next}(v))$。
  • 一个节点的$\mathrm{endpos}(v)$集合包含其所有子节点$\mathrm{endpos}(v)$集合的并,如果这个节点表示了母串的一个前缀,则还包含这个前缀的位置。

构造

流程

在字符集为常数且较小时,SAM的构造法为线性复杂度的。这是一个增量算法,因此,实际上应该解释为,每在母串后插入一个字符,更新SAM的复杂度是$O(1)$的。
这里,我们将主要使用略图结合实际例子的形式描述增量过程。
SAM构造1
上图描述了一个串的SAM略图,我们要在这个串末尾加上字符$c$,图中,灰色的边表示字符$c$的转移边,黑色的边表示后缀链接,$\text{start}$为起始节点$\text{last}$为代表整个母串的节点。

  1. 从$\text{last}$点开始向上跳前缀树,直到遇到一个有字符$c$出边的节点。
    图中,由于第一个有字符$c$出边的节点为$v_3$,说明对应的后缀加上字符$c$后出现在原母串中了,因此,$v_4, v_5, v_6$这些代表$v_3$代表的字符串的后缀的节点对应的后缀,加上字符$c$后也会在原母串中出现。
    对于末尾那些没有字符$c$出边的节点,他们对应的较长后缀加上字符$c$后会形成新的后缀,这些后缀需要一个点来表示,因此,这些点均需要引出一条字符$c$的转移边,指向新建的后缀节点。对于新节点$u$,有$\max(u)=\max(v_1)+c, \min(u)=\min(v_2)+c$。
    而对于前面那些有字符$c$出边的节点,对应的后缀已经被已存在的节点$d$所代表,因此它们不需要新建节点表示新后缀。
    SAM构造2
  2. 如果不存在有字符$c$出边的节点,即第1步跳到了null节点,则新建节点$u$的后缀链接指向起始节点$\text{start}$,即$\mathrm{next}(u) = \text{start}$。
    因为新出现的这些后缀都用节点$u$表示了,即$|\min(u)|=1$,因此需要找到一个$\max(\mathrm{next}(u))=0$这样的$\mathrm{next}(u)$,自然就是起始节点$\text{start}$了。
  3. 如果$\max(d) = \max(v_3) + c$,则有$\mathrm{next}(u) = d$。
    因为$|\max(d)|=|\max(v_3)|+1, |\min(u)|=|\min(v_2)|+1, |\max(v_3)|=|\min(v_2)|+1$,所以可以推出$|\max(d)|=|\min(v_2)|+2=|\min(u)|+1$,得出结论$u$的后缀链接指向$d$。
    SAM构造3
  4. 如果2、3的情况都不符合,此时一定有$|\max(d)|>|\max(v_3)|+1$。因为字符串$\max(v_3)+c$一定存在于$d$中,并且存在另一个异于$v_i$的节点$t$,满足$\mathrm{next}(t)=v_3$且$t$有连向$d$的转移边。
    SAM构造4
    此时$d$不能直接作为$u$的后缀链接,因为$d$中的字符不全是$u$的后缀,这是由$t$带来的影响。因此,我们需要考虑把原来的$d$拆成两个点$o$和$n$,$n$表示长度不大于$|\max(v_3)|+1$的子串,另一个点代表$t$带来的更长的子串。显然,$n$可以继承$d$的出边,且$v_3$、$v_4$都可以转移到$n$,因此$n$表示的就是这两个节点的字符串加上字符$c$转移而来的,且$\mathrm{next}(o)=n$,$\mathrm{next}(n)=\mathrm{next}(d)$。这里借用Menci的例子,图示如下:
    SAM构造5
    SAM构造6
    此时拆出来的的$n$点就可以作为$u$的后缀链接了。

例子

到这里,算法的流程就结束了,下面我们用若干个例子来加深理解。

abc中加入a

abc的SAM结构abca的SAM结构
abc的SAM结构abca的SAM结构

对于这个串,从代表母串的节点4上跳到第一个有字符a出边的节点即是起始节点1,因此节点4需要在下面新建节点5代表新后缀abca,并将节点5的后缀链接指向节点2,即节点1的a出边对应的节点,这是由于$len_2 = len_1 + 1$。在图中,$len$指的就是$|\max|$。

abca中加入c

abca的SAM结构abca的SAM结构
abca的SAM结构abcac的SAM结构

首先从结尾节点5上跳,节点5、2没有字符c的出边,因此新建节点6来代表新后缀abcacac。节点1有字符c的出边,对应节点4,但是$len_4 \neq len_1 + 1$,因此需要对4拆点,新点为7,继承了4的出边,且从4上跳前缀树路径上所有指向4的边都应改为指向7。此时,6和4的后缀链接指向7,7继承了原来4的后缀链接指向1,完成构造。

实现

实现可以参考后文例题中的代码。
实现中,只需记下每个节点的$|\max|$即可,因为$|\min(u)| = |\max(\mathrm{next}(u))| + 1$。

例题:【P3804】【模板】后缀自动机 – 洛谷

题目描述

给定一个只包含小写字母的字符串 S ,请你求出 S 的所有出现次数不为 1 的子串的出现次数乘上该子串长度的最大值。

题解思路

出现次数不为1的子串,可以求出每个节点的$\mathrm{endpos}$集合大小,即上文中提到的,从前缀树的子节点加和,这个过程可以使用拓扑序来做。由于我们知道前缀树的节点越深$|\max|$越大,因此可以对这个值排序后直接转移统计。对于每个$\mathrm{endpos}$集合大小大于1的节点计算其$|\mathrm{endpos}||\max|$值更新答案即可。

代码

// Code by KSkun, 2018/6
#include <cstdio>
#include <cctype>
#include <cstring>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2)
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1; register char c = fgc();
    for (; !isdigit(c); c = fgc()) if (c == '-') neg = -1;
    for (; isdigit(c); c = fgc()) res = (res << 1) + (res << 3) + c - '0';
    return res * neg;
}

const int MAXN = 1000005;

struct SuffixAutomaton {
    int last, tot, ch[MAXN << 1][26], nxt[MAXN << 1], len[MAXN << 1], siz[MAXN << 1];
    SuffixAutomaton() {
        last = tot = 1;
    }
    inline void extend(int c) {
        int u = ++tot, v = last; 
        len[u] = len[v] + 1;
        for(; v && !ch[v][c]; v = nxt[v]) ch[v][c] = u;
        if(!v) {
            nxt[u] = 1;
        } else if(len[ch[v][c]] == len[v] + 1) {
            nxt[u] = ch[v][c];
        } else {
            int n = ++tot, o = ch[v][c]; len[n] = len[v] + 1;
            memcpy(ch[n], ch[o], sizeof(ch[o]));
            nxt[n] = nxt[o]; nxt[o] = nxt[u] = n;
            for(; v && ch[v][c] == o; v = nxt[v]) ch[v][c] = n;
        }
        last = u; siz[u] = 1;
    } 
} sam;

char str[MAXN];
int a[MAXN << 1], c[MAXN << 1];

int main() {
    scanf("%s", str + 1);
    for(int i = 1; str[i]; i++) {
        sam.extend(str[i] - 'a');
    }
    LL ans = 0;
    for(int i = 1; i <= sam.tot; i++) c[sam.len[i]]++;
    for(int i = 1; i <= sam.tot; i++) c[i] += c[i - 1];
    for(int i = 1; i <= sam.tot; i++) a[c[sam.len[i]]--] = i;
    for(int i = sam.tot; i; i--) {
        int p = a[i];
        sam.siz[sam.nxt[p]] += sam.siz[p];
        if(sam.siz[p] > 1) ans = std::max(ans, 1ll * sam.siz[p] * sam.len[p]);
    }
    printf("%lld", ans);
    return 0;
}

其他

大字符集、不确定字符集的处理方法

可以用一个map<int, int>来代替二维数组ch,这样,构造算法的复杂度就会变为$O(n \log n)$。

拓扑序

在例题中提到了,按照$|\max|$从大到小排序可以得到前缀树的拓扑序。这个排序可以使用计数排序,实现线性复杂度。

常见应用

本质不同子串数量

每个节点$u$表示的字符串长度在$[|\min(u)|, |\max(u)|]$范围内,而且每个节点代表的字符串又各不相同,因此对$\max(u)-\min(u)+1$求和即可。

字符串的最小表示法

对于字符串$S$,它的最小表示定义为,将$S$的一个前缀切去并连接在末尾,这样得到的字典序最小的字符串。
求最小表示可以先对$S+S$建立SAM,并从起始节点开始,每次沿着存在的最小字符转移边走,走$|S|$步后,路径对应的字符串即为最小表示。

求两个串的LCS

对其中一个字符串建SAM,记录一个当前匹配的长度$L$和当前走到的节点$v$,枚举另一个字符串的每个字符$c$,如果当前节点有$c$的出边,则沿出边走一步,对$L$加1;否则,跳至$\mathrm{next}(v)$,并将$L$重置为$|\max(\mathrm{next}(v))$,再次检查是否有$c$的出边,做出选择即可。

参考资料

[SPOJ-LCS2]Longest Common Substring II 题解

[SPOJ-LCS2]Longest Common Substring II 题解

题目地址:SPOJ:SPOJ.com – Problem LCS2、洛谷:【S 

[BZOJ3091]城市旅行 题解

[BZOJ3091]城市旅行 题解

题目地址:BZOJ:Problem 3091. — 城市旅行 题目描述 输入输 

[BZOJ3561]DZY Loves Math VI 题解

[BZOJ3561]DZY Loves Math VI 题解

题目地址:BZOJ:Problem 3561. — DZY Loves Math VI

题目描述

给定正整数n,m,求 \displaystyle \sum_{i=1}^n \sum_{i=1}^m \mathrm{lcm}(i, j)^{\mathrm{gcd}(i, j)}

输入输出格式

输入格式:
一行两个整数n,m。

输出格式:
一个整数,为答案模1000000007后的值。

输入输出样例

输入样例#1:

5 4

输出样例#1:

424

说明

数据规模:
1<=n,m<=500000,共有3组数据。

题解

参考资料:BZOJ3561 DZY Loves Math VI 数论 快速幂 莫比乌斯反演 – -zhouzhendong- – 博客园
我们知道有\mathrm{lcm}(i, j)\mathrm{gcd}(i, j) = ij,因此题目给的式子可以搞成
\sum_{i=1}^n \sum_{i=1}^m \left( \frac{ij}{\mathrm{gcd}(i, j)} \right)^{\mathrm{gcd}(i, j)}
类似反演经典题目那样,把gcd=k的问题转化成gcd=1的问题,这样就可以直接判定是否互质了
\sum_{d=1}^n \sum_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor} \sum_{j=1}^{\left\lfloor \frac{m}{d} \right\rfloor} \left( \frac{ijd^2}{d} \right)^d [\mathrm{gcd}(i, j) == 1]
然后你发现这个gcd判定可以用\mu(n)来换掉,有利于之后改变求和顺序
\sum_{d=1}^n \sum_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor} \sum_{j=1}^{\left\lfloor \frac{m}{d} \right\rfloor} \left( \frac{ijd^2}{d} \right)^d \sum_{p|i, p|j} \mu(p)
接着我们枚举gcd,改变求和顺序,在这个过程中还可以把不变的常数项提到前面来
\sum_{d=1}^n d^d \sum_{p=1}^{\left\lfloor \frac{n}{d} \right\rfloor} \mu(p) \sum_{i=1}^{\left\lfloor \frac{n}{pd} \right\rfloor} \sum_{j=1}^{\left\lfloor \frac{m}{pd} \right\rfloor} (ijp^2)^d
最后的整理
\sum_{d=1}^n d^d \sum_{p=1}^{\left\lfloor \frac{n}{d} \right\rfloor} \mu(p)p^{2d} \sum_{i=1}^{\left\lfloor \frac{n}{pd} \right\rfloor} i^d \sum_{j=1}^{\left\lfloor \frac{m}{pd} \right\rfloor} j^d
乘方求前缀和可以预处理,\mu(n)可以筛出来,d^d快速幂搞一下就好。总复杂度O(n \log n)

代码

// Code by KSkun, 2018/6
#include <cstdio>
#include <cctype>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2)
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1; register char c = fgc();
    for (; !isdigit(c); c = fgc()) if (c == '-') neg = -1;
    for (; isdigit(c); c = fgc()) res = (res << 1) + (res << 3) + c - '0';
    return res * neg;
}

const int MAXN = 500005, MO = 1e9 + 7;

bool notprime[MAXN];
int mu[MAXN], prime[MAXN], ptot;

inline void sieve() {
    notprime[1] = true; mu[1] = 1;
    for(int i = 2; i < MAXN; i++) {
        if(!notprime[i]) {
            prime[++ptot] = i; mu[i] = -1;
        }
        for(int j = 1; j <= ptot && prime[j] * i < MAXN; j++) {
            int k = prime[j] * i;
            notprime[k] = true;
            if(i % prime[j]) {
                mu[k] = -mu[i]; 
            } else {
                mu[k] = 0; break;
            }
        }
    }
}

inline LL fpow(LL x, LL k) {
    LL res = 1;
    for(; k; k >>= 1) {
        if(k & 1) res = res * x % MO;
        x = x * x % MO;
    }
    return res;
}

int n, m;
LL powd[MAXN], psum[MAXN];

int main() {
    sieve();
    n = readint(); m = readint();
    if(n > m) std::swap(n, m);
    for(int i = 1; i <= m; i++) {
        powd[i] = 1;
    }
    LL ans = 0;
    for(int d = 1; d <= n; d++) {
        LL res = 0;
        for(int i = 1; i <= m / d; i++) {
            powd[i] = powd[i] * i % MO; psum[i] = (psum[i - 1] + powd[i]) % MO;
        }
        for(int p = 1; p <= n / d; p++) {
            res = (res + mu[p] * powd[p] % MO * powd[p] % MO 
                * psum[n / p / d] % MO * psum[m / p / d] % MO) % MO;
        }
        res = (res % MO + MO) % MO;
        ans = (ans + res * fpow(d, d) % MO) % MO;
    }
    printf("%lld", ans);
    return 0;
}
[BZOJ3560]DZY Loves Math V 题解

[BZOJ3560]DZY Loves Math V 题解

题目地址:BZOJ:Problem 3560. — DZY Loves Mat