【HDOJ5528】Count a * b（积性函数）

题意：设f(i)为0<=x,y<=i-1且xy%i=0的(x,y)对数，g(i)为sigma f(j) [i%j==0]

给定n，求g(n)，答案对2^64取模

T<=2e4，n<=1e9

思路：这题坚定了我要找一个专业数学手的决心……

x,y从[0,i-1]等价于从[1,i]

From Gold_7

最右边那个符号为约数个数

ANS=n所有约数的平方和-n*约数个数

设s[i][j]表示p[i]^0+p[i]^2+...+p[i]^2*j，欧拉筛之后预处理出来

中间有关于答案的变量全部用unsigned long long 

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<vector>
typedef long long ll;
typedef unsigned long long ull;
using namespace std;
#define N   40000 
#define M   32
#define oo  10000000
#define MOD 105225319

ull s[N][M];
int prime[N],isprime[N],tot; 

ull calc1(int n)
{
    int k=n;
    ull ans=1;
    for(int i=1;i<=tot;i++)
    {
        if(prime[i]*prime[i]>n) break;
        if(k==1) break;
        int t=0;
        while(k%prime[i]==0)
        {
            t++;
            k/=prime[i];
        }
        ans*=s[i][t];
    }
    if(k>1) ans*=((ull)k*k+1);
    return ans;
}

ull calc2(int n)
{
    int k=n;
    ull ans=1;
    for(int i=1;i<=tot;i++)
    {
        if(prime[i]*prime[i]>n) break;
        if(k==1) break;
        int t=0;
        while(k%prime[i]==0)
        {
            t++;
            k/=prime[i];
        }
        ans*=(t+1);
    }
    if(k>1) ans*=2;
    return ans;
}

int main()
{ 
     tot=0;
     for(int i=2;i<N;i++)
     {
         if(!isprime[i]) prime[++tot]=i;
         for(int j=1;j<=tot;j++)
         {
             int t=prime[j]*i;
             if(t>N) break;
             isprime[t]=1;
             if(i%prime[j]==0) break;
         }
    }
    for(int i=1;i<=tot;i++)
    {
        ull t=1; s[i][0]=1;
        for(int j=1;j<M;j++)
        {
            t*=prime[i]*prime[i];
            s[i][j]=s[i][j-1]+t;
        }
    }
    int cas; 
    scanf("%d",&cas);
    for(int v=1;v<=cas;v++)
    {
        int n;
        scanf("%d",&n);
        ull ans=calc1(n)-calc2(n)*n;
        printf("%I64u
",ans); 
    }
    return 0;
}