原题见HDU 5528 题意： f(m) 表示 abmodm≠0 的 (a,b) 的个数， a,b 为小于m的非负整数。
G(n)=∑m|nf(m)
给出 n(1≤n≤109) ，求 G(n)mod264 .
分析
mod264
最后结果要 mod264 ，故结果用unsigned long long(64 bit)存，只要计算过程中只有加减乘，即是 mod264 的结果。
求 f(m)
【ACM-解题报告|HDU 5528 Count a - b 欧拉函数 素数】由于 0≡mmodm ，故
{ a,b 为小于m的非负整数得到的结果} = { a,b 为1到m的正整数得到的结果}
对于 a ，有 g=gcd(a,m) ，则第a行满足 abmodm=0 的 b=mg,2mg,...,gmg.共g个 。而对于1到m的a，满足 g=gcd(a,m)的a有?(mg) 个。
f(m)=m2?∑g∣mg??(mg)
求 G(n)
G(n)=∑m|nf(m)=∑m|n(m2?∑g∣mg??(mg))=∑m|nm2?∑m|n∑g∣mg??(mg)
令 X(n)=∑m|nm2,Y(n)=∑m|n∑g∣mg??(mg)
X(n)即所有n的约数的平方和。若对n作质因数分解 n=pa11pa22...pass
则 X(n)=(p01+p21+...+p2a11)(p02+p22+...+p2a22)...(p0s+p2s...+p2ass)=Π1≤i≤s∑0≤j≤aip2ji
注意求和部分不能用等比公式求出通项。因为这样会产生除法，如果分子在运算时已经超 264会造成无法整除等运算错误。
对于 Y(n) ，若给定n和m，则会产生数对 (g,mg) .且一个数对只对应唯一的m，可以得到所有的数对 (a,b) 满足 ab∣n ，故
Y(n)=∑xy∣nx??(y)=∑x∣n(x∑y∣nx?(y))
已知 欧拉函数性质 ∑d∣n?(d)=n
则 Y(n)=∑x∣n(x?nx)=n∑x∣n1=nΠ1≤i≤s(ai+1)
故 G(n)=Π1≤i≤s∑0≤j≤aip2ji?nΠ1≤i≤s(ai+1)
附代码 注意：要估算好素数表的大小避免超时
sum[i][j] =∑0≤k≤jp2ki
p[i]=素数在表中时保存序号，在表外时直接保存素数
/*--------------------------------------------
* File Name: HDU 5528
* Author: Danliwoo
* Mail: Danliwoo@outlook.com
* Created Time: 2016-07-10 14:05:24
--------------------------------------------*/#include
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;
} 
