1. 程式人生 > >【HDOJ5528】Count a * b(積性函式)

【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 

 1 #include<cstdio>
 2
#include<cstring> 3 #include<iostream> 4 #include<algorithm> 5 #include<vector> 6 typedef long long ll; 7 typedef unsigned long long ull; 8 using namespace std; 9 #define N 40000 10 #define M 32 11 #define oo 10000000 12 #define MOD 105225319 13 14 ull s[N][M]; 15 int prime[N],isprime[N],tot;
16 17 ull calc1(int n) 18 { 19 int k=n; 20 ull ans=1; 21 for(int i=1;i<=tot;i++) 22 { 23 if(prime[i]*prime[i]>n) break; 24 if(k==1) break; 25 int t=0; 26 while(k%prime[i]==0) 27 { 28 t++; 29 k/=prime[i]; 30 } 31
ans*=s[i][t]; 32 } 33 if(k>1) ans*=((ull)k*k+1); 34 return ans; 35 } 36 37 ull calc2(int n) 38 { 39 int k=n; 40 ull ans=1; 41 for(int i=1;i<=tot;i++) 42 { 43 if(prime[i]*prime[i]>n) break; 44 if(k==1) break; 45 int t=0; 46 while(k%prime[i]==0) 47 { 48 t++; 49 k/=prime[i]; 50 } 51 ans*=(t+1); 52 } 53 if(k>1) ans*=2; 54 return ans; 55 } 56 57 int main() 58 { 59 tot=0; 60 for(int i=2;i<N;i++) 61 { 62 if(!isprime[i]) prime[++tot]=i; 63 for(int j=1;j<=tot;j++) 64 { 65 int t=prime[j]*i; 66 if(t>N) break; 67 isprime[t]=1; 68 if(i%prime[j]==0) break; 69 } 70 } 71 for(int i=1;i<=tot;i++) 72 { 73 ull t=1; s[i][0]=1; 74 for(int j=1;j<M;j++) 75 { 76 t*=prime[i]*prime[i]; 77 s[i][j]=s[i][j-1]+t; 78 } 79 } 80 int cas; 81 scanf("%d",&cas); 82 for(int v=1;v<=cas;v++) 83 { 84 int n; 85 scanf("%d",&n); 86 ull ans=calc1(n)-calc2(n)*n; 87 printf("%I64u\n",ans); 88 } 89 return 0; 90 } 91