3309: DZY Loves Math

Description

對于正整數n，定義f(n)為n所含質因子的最大冪指數。例如f(1960)=f(2^3 * 5^1 * 7^2)=3, f(10007)=1, f(1)=0。
給定正整數a,b，求sigma(sigma(f(gcd(i,j)))) (i=1..a, j=1..b)。

Input

第一行一個數T，表示詢問數。
接下來T行，每行兩個數a,b，表示一個詢問。

Output

對于每一個詢問，輸出一行一個非負整數作為回答。

Sample Input

4
7558588 9653114
6514903 4451211
7425644 1189442
6335198 4957

Sample Output

35793453939901
14225956593420
4332838845846
15400094813

HINT

【數據規模】

T<=10000

1<=a,b<=10^7

---

?　　這道題我足足寫了一白板。因為這實在是一道好題。

　　我們要求∑∑f(gcd(i,j))，可以枚舉d=gcd(i,j)

　　則有∑f(d)∑∑[gcd(i,j)==d]=∑f(d)∑∑e(gcd(i,j))=∑f(d)∑∑∑[k|i][k|j]mu(k)

　　即為∑f(d)∑mu(k)(n/dk)(m/dk)

　　令T=dk，則d|T,k=T/d，故有∑(n/T)(m/T)∑f(d)mu(T/d)。設g(T)=∑f(d)mu(T/d)。

　　如果我們能夠O（n）算出g(T)及其前綴和，那么按照套路便能直接用分塊。O（sqrt（n））單次詢問。

　　現在看如何預處理出g(T)=∑f(d)mu(T/d)：先唯一分解出T和d，再分析mu(T/d)與f(d)的性質。然后找最大項數目，分類討論，再進一步分類討論即得解。這個分類討論當時證的很復雜。有篇博客所言甚善：

　　如果存在ai≠aj(i≠j)，那答案只取決于最大的冪，對于其余部分通過組合數的性質，取得的數字個數為奇數和偶數的數量相等，因此和都是0，故如果存在ai≠aj(i≠j)，則g(T)=0

　　如果所有的a值都相等，我們假設對于任意選取方案，f值都不變，那么由于選取奇數個元素和偶數個元素的方案數相等，和仍然為0

　　但是有一種選取方案的f值=a-1 即d=p1*p2*p3*...pk時，因此要把那個減去，根據莫比烏斯函數性質，減去的是(-1)^k，因為是減去所以最終結果為(-1)^(k+1)，故如果不存在ai≠aj，則g(T)=(-1)^(k+1)

　　然后，可以使用歐拉篩O（n）預處理。但是怎樣處理呢？考慮到歐拉篩中i*prime[j]的最小質因子是prime[j]。于是我們可以直接使用i的最小質因子與j比較。所以對于每一個數，需要存g[i]和i的最小質因子的冪指數a[i]。

　　然后，當i%prime[j]==0，我們怎么做呢？常見的做法是把i中所有的prime[j]都用除法除去。但事實上，這樣效率太慢。考慮到歐拉篩中i*prime[j]的最小質因子是prime[j]（可用反證法證明），那么i的最小質因子也是prime[j]。如果我們存下minpa[i]表示i最小質因數的a[i]次方,那么直接除去即可。

　　眼見為實。代碼如下：

1 /************************************************************** 2 Problem: 3309 3 User: Doggu 4 Language: C++ 5 Result: Accepted 6 Time:11000 ms 7 Memory:166836 kb 8 ****************************************************************/ 9 10 #include <cstdio> 11 #include <algorithm> 12 const int N = 1e7+5; 13 int minpa[N], a[N], g[N], prime[N], ptot; 14 bool vis[N]; 15 void EULER(int n) { 16 for( int i = 2; i <= n; i++ ) { 17 if(!vis[i]) prime[++ptot]=i, minpa[i]=i, a[i]=1, g[i]=1; 18 for( int j = 1; j <= ptot; j++ ) { 19 if((long long)i*prime[j]>n) break; 20 vis[i*prime[j]]=1; 21 minpa[i*prime[j]]=prime[j]; 22 a[i*prime[j]]=1; 23 g[i*prime[j]]=(a[i]==1?-g[i]:0); 24 if(i%prime[j]==0) { 25 minpa[i*prime[j]]=minpa[i]*prime[j]; 26 a[i*prime[j]]=a[i]+1; 27 int temp=i/minpa[i]; 28 if(temp==1) g[i*prime[j]]=1; 29 else g[i*prime[j]]=(a[temp]==a[i*prime[j]]?-g[temp]:0); 30 break; 31 } 32 } 33 } 34 for( int i = 1; i <= n; i++ ) g[i]+=g[i-1]; 35 } 36 int main() { 37 int T, n, m; 38 EULER(1e7); 39 scanf("%d",&T); 40 while(T--) { 41 scanf("%d%d",&n,&m); 42 long long ans=0; 43 if(n>m) std::swap(n,m); 44 for( int a = 1, ed; a <= n; a = ed + 1 ) { 45 ed=std::min(n/(n/a),m/(m/a)); 46 ans+=(long long)(g[ed]-g[a-1])*(n/a)*(m/a); 47 } 48 printf("%lld\n",ans); 49 } 50 } 51 狄利克雷卷積+函數推演+歐拉篩+分塊

轉載于:https://www.cnblogs.com/Doggu/p/bzoj3309.html

總結

以上是生活随笔為你收集整理的BZOJ 3309 DZY Loves Math的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。