【jzoj3117】【WinterCamp 2013】【模积和】【分块】【扩展欧几里得】

题目大意

求∑∑((n mod i)*(m mod j))其中1<=i<=n,1<=j<=m,i≠j。答案mod 19940417。

解题思路

观察可知ans=∑∑((n mod i)(m mod j))1<=i<=n,1<=j<=m,i≠j。

=>ans=∑(n mod i)∑(m mod j)-∑((n mod i)(m mod i))

f(n)=∑(n mod i)

g(n,m)=∑(n mod i)(m mod i)

ans=f(n)f(m)-g(n,m)

$f(n)=\sum(n-\lfloor n/i \rfloor i)$

对于$\lfloor n/i \rfloor$我们可以使用分块，对于[i,j]满足n/i==n/j可以一起处理，对于n可以用乘法求，对于*i可以用等差数列求和。

$g(n,m)=\sum(nm+\lfloor n/i \rfloor\lfloor m/i \rfloor)i^2-\lfloor n/i \rfloor mi-\lfloor m/i \rfloor ni)$

对于第一项可以用乘法解决，对于第二项可以分块后用完全平方和求和公式求，第三四项用等差数列求和。

∑(1^2+2^2+……+n^2)==n(n+1)(2n+1)/6详细证明不再给出

观察可知模数不是质数，使用扩展欧几里得求出逆元。

aa’==1(%mo)

=>aa’+moy==1

由于a，mo互质，a’可以用扩展欧几里得求出。

code

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define LF double
#define LL long long
#define min(a,b) ((a<b)?a:b)
#define max(a,b) ((a>b)?a:b)
#define fo(i,j,k) for(int i=j;i<=k;i++)
#define fd(i,j,k) for(int i=j;i>=k;i--)
using namespace std;
int const mo=19940417;
int n,m,ni;
int calc(int n){
    int ans=1ll*n*n%mo;
    for(int i=1,j=1;i<=n;i=j+1){
        j=n/(n/i);
        ans=(ans-1ll*(i+j)*(j-i+1)/2%mo*(n/i)%mo)%mo;
    }
    return (ans+mo)%mo;
}
int calc3(int n){
    return 1ll*n*(n+1)%mo*(2*n+1)%mo*ni%mo;
}
int calc2(int n,int m){
    int mi=min(n,m),ans=1ll*n*m%mo*mi%mo;
    for(int i=1,j=1;i<=mi;i=j+1){
        j=min(n/(n/i),m/(m/i));
        ans=(ans+1ll*(calc3(j)-calc3(i-1))%mo*(n/i)%mo*(m/i)%mo
                      -1ll*(i+j)*(j-i+1)/2%mo*(n/i)%mo*m%mo
                      -1ll*(i+j)*(j-i+1)/2%mo*(m/i)%mo*n%mo)%mo;
    }
    return (ans+mo)%mo;
}
void exgcd(int a,int &x,int b,int &y){
    if(!b){x=1;y=0;return;}
    exgcd(b,y,a%b,x);
    y-=a/b*x;
}
int main(){
    freopen("d.in","r",stdin);
    freopen("d.out","w",stdout);
    scanf("%d%d",&n,&m);
    int y;exgcd(6,ni,mo,y);
    printf("%d",(1ll*calc(n)*calc(m)%mo-calc2(n,m)+mo)%mo);
    return 0;
}