Description
求\(\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j)\)，答案模1e9+9輸出，多組詢問

Input
一個正整數T表示數據組數
接下來T行 每行兩個正整數 表示N、M

Output
T行 每行一個整數 表示第i組數據的結果

Sample Input
1
4 5

Sample Output
122

HINT
T <= 10000
N, M<=10000000

---

我們令n<m，然后將柿子化簡
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j)\]
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m \dfrac{i\times j}{\gcd(i,j)}\]
\[\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^m\dfrac{i\times j}ze8trgl8bvbq[\gcd(i,j)=d]\]
\[\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor\frac{n}ze8trgl8bvbq\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}ze8trgl8bvbq\rfloor}\dfrac{d^2\times i\times j}ze8trgl8bvbq[\gcd(i,j)=1]\]
\[\sum\limits_{d=1}^n d\sum\limits_{i=1}^{\lfloor\frac{n}ze8trgl8bvbq\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}ze8trgl8bvbq\rfloor}i\times j\sum\limits_{x|i,x|j}\mu(x)\]
\[\sum\limits_{d=1}^n d\sum\limits_{x=1}^{\lfloor\frac{n}ze8trgl8bvbq\rfloor}\mu(x)\times x^2\sum\limits_{i=1}^{\lfloor\frac{n}{dx}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{dx}\rfloor}i\times j\]
我們發現最后那個是等差數列，繼續化簡
\[\sum\limits_{d=1}^n d\sum\limits_{x=1}^{\lfloor\frac{n}ze8trgl8bvbq\rfloor}\mu(x)\times x^2\dfrac{\lfloor\frac{n}{dx}\rfloor(\lfloor\frac{n}{dx}\rfloor+1)\lfloor\frac{m}{dx}\rfloor(\lfloor\frac{m}{dx}\rfloor+1)}{4}\]
然后我們令\(T=dx\)，那么得到
\[\sum\limits_{T=1}^n\dfrac{\lfloor\frac{n}{T}\rfloor(\lfloor\frac{n}{T}\rfloor+1)\lfloor\frac{m}{T}\rfloor(\lfloor\frac{m}{T}\rfloor+1)}{4}T\sum\limits_{x|T}\mu(x)x\]
我們設\(f(T)=\sum\limits_{x|T}\mu(x)x\)，預處理出f，就可以分塊了

設\(g(x)=\mu(x)x\)，當a,b互質，\(g(a)\times g(b)=ab\mu(a)\mu(b)=ab\mu(ab)=g(ab)\)，所以g是積性函數，根據莫比烏斯反演的性質，f也是積性函數

令\(T=\prod\limits_{i=1}^k P_i^{x_i}\)，\(f(P_i^{x_i})=(1-P_i)\)，那么
\[f(T)=\prod\limits_{i=1}^k(1-P_i)\]
這樣子我們可以在\(O(n)\)時間內線篩出來，然后維護一下\(T\times f(T)\)的前綴和，然后就可以在\(O(\sqrt N)\)的時間內完成每次詢問

/*program from Wolfycz*/ #include<cmath> #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #define inf 0x7f7f7f7f using namespace std; typedef long long ll; typedef unsigned int ui; typedef unsigned long long ull; inline int read(){int x=0,f=1;char ch=getchar();for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;for (;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+ch-'0';return x*f; } inline void print(int x){if (x>=10) print(x/10);putchar(x%10+'0'); } const int N=1e7,p=1e8+9; int prime[N+10],f[N+10],tot; bool inprime[N+10]; void prepare(){f[1]=1;for (int i=2;i<=N;i++){if (!inprime[i]) prime[++tot]=i,f[i]=1-i+p;for (int j=1;j<=tot&&i*prime[j]<=N;j++){inprime[i*prime[j]]=1;if (i%prime[j]==0){f[i*prime[j]]=f[i];break;}f[i*prime[j]]=1ll*f[i]*f[prime[j]]%p;}}for (int i=1;i<=N;i++) f[i]=(f[i-1]+1ll*i*f[i]%p)%p; } int get(int x){return (1ll*x*(x+1)>>1)%p;} int main(){prepare();for (int Data=read();Data;Data--){int n=read(),m=read(),pos,Ans=0;if (n>m) swap(n,m);for (int T=1;T<=n;T=pos+1){pos=min(n/(n/T),m/(m/T));Ans=(Ans+1ll*get(n/T)*get(m/T)%p*(f[pos]-f[T-1]+p)%p)%p;}printf("%d\n",Ans);}return 0; }

轉載于:https://www.cnblogs.com/Wolfycz/p/9489894.html

總結

以上是生活随笔為你收集整理的[BZOJ2693]jzptab的全部內容，希望文章能夠幫你解決所遇到的問題。

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