題目描述

去年的Lucas非常喜歡數論題，但是一年以后的Lucas卻不那么喜歡了。

在整理以前的試題時，發現了這樣一道題目“求Sigma(f(i)),其中1<=i<=N”，其中 表示i的約數個數。他現在長大了，題目也變難了。 求如下表達式的值： 其中f(ij)表示ij的約數個數。 他發現答案有點大，只需要輸出模1000000007的值。

輸入

第一行一個整數n。

輸出

一行一個整數ans，表示答案模1000000007的值。

樣例輸入

2

樣例輸出

8

---

題解

莫比烏斯反演+杜教篩

首先有個神奇的結論：

證明：只考慮質數p，設n=a*p^x，m=b*p^y，那么等式左端p的貢獻顯然為x+y+1；而等式右邊p的貢獻為數對(i,j)：(p,1),(p^2,1),...,(p^x,1) , (1,p),(1,p^2),...,(1,p^y) , (1,1)，共x+y+1對，因此命題得證。

然后就有：

注意到n/d只有O(√n)種取值，因此我們可以枚舉這些n/d，找出對應的d的范圍，前邊的mu(d)求一下前綴和，后面的sigma，由于n/d/i只有O(√n/d)種取值，因此可以用O(√n/d)的復雜度求出。

由于這里面的n值過大，無法直接維護mu的前綴和，因此需要使用杜教篩的方法，參見 bzoj3944

因此總的時間復雜度貌似是O(n^3/4+n^2/3logn)。

#include <cstdio> #include <map> #define N 1000010 using namespace std; #define mod 1000000007 typedef long long ll; map<ll , ll> f; map<ll , ll>::iterator it; ll m = 1000000 , mu[N] , prime[N] , tot , sum[N]; bool np[N]; ll cal1(ll n) {if(n <= m) return sum[n];it = f.find(n);if(it != f.end()) return it->second;ll ans = 1 , i , last;for(i = 2 ; i <= n ; i = last + 1) last = n / (n / i) , ans = (ans - (last - i + 1) * cal1(n / i) % mod + mod) % mod;return f[n] = ans; } ll cal2(ll n) {ll ans = 0 , i , last;for(i = 1 ; i <= n ; i = last + 1) last = n / (n / i) , ans = (ans + n / i * (last - i + 1)) % mod;return ans * ans % mod; } int main() {ll n , i , j , last , ans = 0;mu[1] = sum[1] = 1;for(i = 2 ; i <= m ; i ++ ){if(!np[i]) mu[i] = -1 , prime[++tot] = i;for(j = 1 ; j <= tot && i * prime[j] <= m ; j ++ ){np[i * prime[j]] = 1;if(i % prime[j] == 0){mu[i * prime[j]] = 0;break;}else mu[i * prime[j]] = -mu[i];}sum[i] = (sum[i - 1] + mu[i] + mod) % mod;}scanf("%lld" , &n);for(i = 1 ; i <= n ; i = last + 1) last = n / (n / i) , ans = (ans + (cal1(last) - cal1(i - 1) + mod) % mod * cal2(n / i)) % mod;printf("%lld\n" , ans);return 0; }

?

轉載于:https://www.cnblogs.com/GXZlegend/p/6999146.html

《新程序員》：云原生和全面數字化實踐50位技術專家共同創作，文字、視頻、音頻交互閱讀

總結

以上是生活随笔為你收集整理的【bzoj4176】Lucas的数论 莫比乌斯反演+杜教筛的全部內容，希望文章能夠幫你解決所遇到的問題。

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