首页 > 代码库 > bzoj 3309 DZY Loves Math - 莫比乌斯反演 - 线性筛

bzoj 3309 DZY Loves Math - 莫比乌斯反演 - 线性筛

对于正整数n,定义f(n)为n所含质因子的最大幂指数。例如f(1960)=f(2^3 * 5^1 * 7^2)=3, f(10007)=1, f(1)=0。 给定正整数a,b,求sigma(sigma(f(gcd(i,j)))) (i=1..a, j=1..b)。

 

Input

第一行一个数T,表示询问数。 接下来T行,每行两个数a,b,表示一个询问。

 

Output

对于每一个询问,输出一行一个非负整数作为回答。

 

Sample Input

47558588 96531146514903 44512117425644 11894426335198 4957

Sample Output

3579345393990114225956593420433283884584615400094813

Hint

【数据规模】
T<=10000
1<=a,b<=10^7


  题目大意 定义函数技术分享为数x所含质因子的最大幂指数。求技术分享

  继续用之前的套路统计每个f被计算的次数,然后将gcd化为莫比乌斯函数之和。

  技术分享

  如果这样,似乎并没有什么好方法求值,所以试试把那几个向下取整的分数搞到外层来,这样就可以找机会分段处理。于是得到了下面这个式子:

技术分享

  好了问题来了,如何把后面的值快速预处理出来,前面的是可以通过分段降低时间复杂度。

  为了方便描述,所以我们设技术分享

  为了更好地理解证明,你需要清楚一个东西,就是对于一个非空有限集合S,满足技术分享。至于它的证明,考虑S的k元素子集的个数,然后二项式定理搞成(1-1)的某次幂。

  由于我比较笨,不会大佬们的分类讨论做法。就用集合论乱搞骗结论。

  首先对n进行质因数分解,得到技术分享。然后设技术分享技术分享(技术分享是P‘的补集)。

  根据莫比乌斯函数的性质,对于任何一个数有一个质因子的幂指数大于1,则它的莫比乌斯函数值为0,即为了使式子对答案有贡献,就应当满足,当d分解后,pi的幂指数要么是ai要么是ai - 1然后有

技术分享

  式子中A的意义大概就是它里面的质因子pi的幂指数就是ai - 1.

  然后根据这个又可以得到一个事实,中间的max一坨的取值要么是a要么是(a - 1)。所以我们可以把它拆成两部分求和。

  一部分是包含整个P‘,所以此时那一坨的取值为a-1,另一部分是只包含P‘的一部分,所以此时那一坨的取值为a。

  技术分享

  然后前一部分事实上可以把A拆成技术分享的子集和P‘的并集,后一部分也可以把A拆成P‘的真子集和技术分享的子集。

技术分享

  然后发现左右两边的求和符号可以合并。于是得到了

技术分享

  接下来考虑技术分享

  1) 如果它不为空的话,那么意味着求和符号下的值为0。换言之,就是n中的所有幂指数不相等,那么就有技术分享

  2) 如果它为空的话,意味着求和符号的值为1,P‘的大小就是k,那么就有技术分享

  然后有什么用呢?

  考虑预处理,由于1e7的凶残数据范围所以不能够O(nlogn)预处理。就考虑线性筛。

  这里要用到线性筛一个很重要的性质,每个合数保证被其最小的质因子筛掉。所以,我们只需要记录每个数的最小质因子的次数pos。

  这样在prime[j]不整除i的时候我们就可以轻松处理掉(判断最小质因子的次数是否为1,g(i)是否为0)。

  现在考虑prime[j]整除i的时候,通过常用套路,我们考虑i把prime[j]除尽后的数x和 i / x * prime[j]。此时显然可以判断一下pos[x]和pos[i] + 1是否相等,更新g值就好了。

  因此,我们还需要记录这个数的最小质因子去除这个数,直到不能再除后得到的数。这些记录的东西都可以在线性筛中得到(假设读者都会,就直接贴代码了)。

  我表示这是有史以来我写过最长的题解和最短的模板(bits/std++.h直接抵掉我15行左右的include)。

Code

 1 /** 2  * bzoj 3  * Problem#3309 4  * Accepted 5  * Time:11804ms 6  * Memory:170040k 7  */ 8 #include <bits/stdc++.h> 9 #ifndef WIN3210 #define Auto "%lld"11 #else12 #define Auto "%I64d"13 #endif14 using namespace std;15 typedef bool boolean;16 17 #define LL long long18 19 const int limit = 1e7;20 21 int num = 0;22 int prime[700100];23 boolean vis[limit + 1];24 int last[limit + 1];25 int posm[limit + 1];26 LL g[limit + 1];27 inline void Euler() {28     memset(vis, false, sizeof(vis));29     g[0] = g[1] = -1;30     for(int i = 2; i <= limit; i++) {31         if(!vis[i])    prime[num++] = i, last[i] = 1, posm[i] = 1, g[i] = 1;32         for(int j = 0; j < num && i * 1LL * prime[j] <= limit; j++) {33             int c = i * prime[j];34             vis[c] = true;35             if(i % prime[j]) {36                 last[c] = i, posm[c] = 1;37                 g[c] = (posm[i] == 1) ? (-g[i]) : (0);38             } else {39                 last[c] = last[i];40                 posm[c] = posm[i] + 1;41                 g[c] = (posm[last[c]] == posm[c] || last[c] == 1) ? (-g[last[c]]) : (0);42                 break;43             }44         }45     }46     for(int i = 2; i <= limit; i++)47         g[i] += g[i - 1];48 }49 50 int a, b;51 inline void init() {52     scanf("%d%d", &a, &b);53 }54 55 inline void solve() {56     if(a > b)    swap(a, b);57     long long res = 0;58     for(int i = 1, j; i <= a; i = j + 1) {59         j = min(a / (a / i), b / (b / i));60         res += (a / i) * 1LL * (b / i) * (g[j] - g[i - 1]);61     }62     printf(Auto"\n", res);63 }64 65 int T;66 int main() {67     Euler();68     scanf("%d", &T);69     while(T--) {70         init();71         solve();72     }73     return 0;74 }

 

bzoj 3309 DZY Loves Math - 莫比乌斯反演 - 线性筛