HDU 5528 Count a - b 欧拉函数 素数

原题见HDU 5528

题意：$f(m)$表示$ab\mod m\neq 0$的$(a,b)$的个数，$a,b$为小于m的非负整数。
$G(n)=\sum\limits_{m|n}f(m)$
给出$n(1 \le n \le 10^9)$，求$G(n)\mod 2^{64}$.

分析

$\mod 2^{64}$

最后结果要 mod $2^{64}$，故结果用unsigned long long(64 bit)存，只要计算过程中只有加减乘，即是 mod $2^{64}$的结果。

求$f(m)$

由于$0\equiv m\mod m$，故
{$a,b$为小于m的非负整数得到的结果}$=${$a,b$为1到m的正整数得到的结果}
对于$a$，有$g=gcd(a,m)$，则第a行满足$ab\mod m = 0$的$b={m\over g},2{m\over g},...,g{m\over g}.共g个$。而对于1到m的a，满足$g=gcd(a,m)的a有\phi({m\over g})$个。

$$f(m)=m^2-\sum\limits_{g\mid m}g·\phi({m\over g})$$

求$G(n)$

$$G(n)=\sum\limits_{m|n}f(m)=\sum\limits_{m|n}(m^2-\sum\limits_{g\mid m}g·\phi({m\over g}))$$ $$=\sum\limits_{m|n}m^2-\sum\limits_{m|n}\sum\limits_{g\mid m}g·\phi({m\over g})$$ 令 $$X(n)=\sum\limits_{m|n}m^2,~~~~Y(n)=\sum\limits_{m|n}\sum\limits_{g\mid m}g·\phi({m\over g})$$
$X(n)$即所有n的约数的平方和。若对n作质因数分解 $$n=p_1^{a_1}p_2^{a_2}...p_s^{a_s}$$ 则 $$X(n)=(p_1^0+p_1^2+...+p_1^{2a_1})(p_2^0+p_2^2+...+p_2^{2a_2})...(p_s^0+p_s^2...+p_s^{2a_s})$$ $$=\large \Pi_{1\le i\le s}\small\sum\limits_{0\le j\le a_i}\large p_i^{2j}$$
注意求和部分不能用等比公式求出通项。因为这样会产生除法，如果分子在运算时已经超 $2^{64}$会造成无法整除等运算错误。

对于$Y(n)$，若给定n和m，则会产生数对$(g,{m\over g})$.且一个数对只对应唯一的m，可以得到所有的数对$(a,b)$满足$ab\mid n$，故

$$Y(n)=\sum\limits_{xy\mid n}x·\phi(y)=\sum\limits_{x\mid n}(x\sum\limits_{y\mid {n\over x}}\phi(y))$$ 已知 欧拉函数性质 $$\sum\limits_{d\mid n}\phi (d)=n$$ 则 $$Y(n)=\sum\limits_{x\mid n}(x·{n\over x})=n\sum\limits_{x\mid n}1=n\Pi_{1\le i\le s}(a_i+1)$$ 故 $$G(n)=\large \Pi_{1\le i\le s}\small\sum\limits_{0\le j\le a_i}\large p_i^{2j}-n\Pi_{1\le i\le s}(a_i+1)$$

附代码

注意：要估算好素数表的大小避免超时
sum[i][j]$=\sum\limits_{0\le k\le j} p_i^{2k}$
p[i]=素数在表中时保存序号，在表外时直接保存素数

/*--------------------------------------------
 * File Name: HDU 5528
 * Author: Danliwoo
 * Mail: Danliwoo@outlook.com
 * Created Time: 2016-07-10 14:05:24
--------------------------------------------*/

#include <bits/stdc++.h>
using namespace std;
#define LL unsigned long long
#define N 32000
int num[N], prim[N], pn;
#define M 32
LL sum[N][M], p[M];
int a[M];
void st(){
    pn = 0;
    memset(num, -1, sizeof(num));
    for(int i = 2;i < N;i++){
        if(num[i]) prim[pn++] = i;
        for(int j = 0;j < pn && 1LL*i*prim[j] < N;j++){
            num[i*prim[j]] = 0;
            if(i % prim[j] == 0) break;
        }
    }
    memset(sum, 0, sizeof(sum));
    for(int i = 0;i < pn;i++){
        sum[i][0] = 1;
        LL pm = 1LL*prim[i] * prim[i];
        for(int j = 1;j < M;j++)
            sum[i][j] = sum[i][j-1] * pm;
        for(int j = 1;j < M;j++)
            sum[i][j] += sum[i][j-1];
    }
}
int main()
{
    st();
    int T;
    scanf("%d", &T);
    while(T--){
        int n;
        scanf("%d", &n);
        memset(p, 0, sizeof(p));
        memset(a, 0, sizeof(a));
        int m = n, s = 0;
        for(int i = 0;i < pn;i++) if(m % prim[i] == 0){
            p[s++] = i;
            while(m % prim[i] == 0){
                a[s-1]++;
                m /= prim[i];
            }
            if(m == 1) break;   //大素数表的时候这个break很有用
        }
        if(m > 1){
            p[s++] = m;
            a[s-1] = 1;
        }
        LL x = 1, y = 1;
        for(int i = 0;i < s;i++){
            if(p[i] < pn)
                x *= sum[p[i]][a[i]];
            else 
                x *= (1+p[i]*p[i]);
        }
        for(int i = 0;i < s;i++)
            y *= 1LL*(a[i]+1);
        y *= 1LL*n;
        printf("%I64u\n", x-y);
    }
    return 0;
}