【数论分块】[BZOJ2956、LuoguP2260] 模积和
十年OI一场空，忘记取模见祖宗
题目：
求$$\sum_{i=1}^{n}\sum_{j=1}^{m} (n \bmod i)(m \bmod i)$$
(其中i,j不相等)
暴力拆式子：
$$ANS=\sum_{i=1}^{n}\sum_{j=1}^{m} (n- \left \lfloor \frac{n}{i} \right \rfloor*i)(m- \left \lfloor \frac{m}{i} \right \rfloor*i)-\sum_{i=1}^{min(n,m)} (n- \left \lfloor \frac{n}{i} \right \rfloor *i)(m- \left \lfloor \frac{m}{i} \right \rfloor *i)$$
令$f(n)=\sum_{i=1}^{n} (n- \left \lfloor \frac{n}{i} \right \rfloor *i)$
令$g(n)=\sum_{i=1}^{n}(n- \left \lfloor \frac{n}{i} \right \rfloor *i)(m- \left \lfloor \frac{m}{i} \right \rfloor *i)$
不妨设n<=m
则有
$$ANS=f(n)*f(m)-g(n)$$
其中$$g(n)=\sum_{i=1}^{n} n*m-n*\sum_{i=1}^{n} \left \lfloor \frac{m}{i} \right \rfloor *i-m*\sum_{i=1}^{n} \left \lfloor \frac{n}{i} \right \rfloor *i+\sum_{i=1}^{n} \left \lfloor \frac{n}{i} \right \rfloor* \left \lfloor \frac{m}{i} \right \rfloor *i^2$$
且易有$$\sum_{i=1}^{n} i^2=\frac{n*(n+1)*(2*n+1)}{6}$$
预处理6在模19940417意义下的逆元（我用了exgcd）
然后用数论分块把上面一堆东西算一下即可
1 #include<bits/stdc++.h> 2 #define int long long 3 #define writeln(x) write(x),puts("") 4 #define writep(x) write(x),putchar(' ') 5 using namespace std; 6 inline int read(){ 7 int ans=0,f=1;char chr=getchar(); 8 while(!isdigit(chr)){if(chr=='-') f=-1;chr=getchar();} 9 while(isdigit(chr)){ans=(ans<<3)+(ans<<1)+chr-48;chr=getchar();} 10 return ans*f; 11 }void write(int x){ 12 if(x<0) putchar('-'),x=-x; 13 if(x>9) write(x/10); 14 putchar(x%10+'0'); 15 }const int mod = 19940417; 16 int n,m,k; 17 inline void Add(int &x,int y){x+=y;x%=mod;} 18 void exgcd(int a,int b,int &x,int &y){ 19 if(b==0)return x=1,y=0,void(); 20 exgcd(b,a%b,x,y); 21 int t=x;x=y,y=t-a/b*y; 22 }int inv(int x){ 23 int xx,y; 24 exgcd(6,mod,xx,y); 25 xx=(xx%mod+mod)%mod; 26 return xx; 27 }const int inv6=inv(6); 28 int sum(int x){return (x)*(x+1)%mod*(2*x%mod+1)%mod*inv6%mod;} 29 int query1(int l,int r){return ((sum(r)-sum(l-1))%mod+mod)%mod;} 30 int query2(int l,int r){int ans=(r-l+1)*(l+r)/2;return ans%mod;} 31 int calc1(int n){ 32 int ans=0; 33 for(int i=1,j,t;i<=n;i=j+1){ 34 j=n/(n/i); 35 t=n/i*(i+j)*(j-i+1)/2; 36 t%=mod; 37 Add(ans,t); 38 }ans=n*n%mod-ans; 39 ans=(ans%mod+mod)%mod; 40 return ans; 41 }int calc2(int k){ 42 int ans=0; 43 for(int i=1,j,t;i<=n;i=j+1){ 44 j=min(n/(n/i),m/(m/i)); 45 int s1=n*(m/i)%mod*query2(i,j)%mod; 46 int s2=m*(n/i)%mod*query2(i,j)%mod; 47 int s3=(n/i)*(m/i)%mod*query1(i,j)%mod; 48 Add(s1,s2); 49 Add(ans,s1); 50 ans-=s3; 51 ans=((ans)%mod+mod)%mod; 52 }return ans; 53 } 54 signed main(){ 55 n=read(),m=read(); 56 if(n>m)swap(n,m); 57 int ans=calc1(n)*calc1(m)%mod; 58 ans-=n*m%mod*n%mod; 59 ans=(ans%mod+mod)%mod; 60 ans+=calc2(n); 61 ans=(ans%mod+mod)%mod; 62 cout<<ans<<endl; 63 return 0; 64 }
