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

原题见HDU 5528

题意：f(m)表示abmodm≠0的(a,b)的个数，a,b为小于m的非负整数。

G(n)=∑m|nf(m)

给出n(1≤n≤109)，求G(n)mod264.

分析

mod264

最后结果要 mod 264，故结果用

unsigned long long

(64 bit)存，只要计算过程中只有加减乘，即是 mod 264的结果。

求f(m)

由于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 <bits/stdc++.h>
using namespace std;
#define LL unsigned long long
#define N 32000
int num
, prim
, pn;
#define M 32
LL sum
[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;
}