bzoj 2956: 模积和 分块+数学

题意

求∑∑((n mod i)*(m mod j))其中1<=i<=n,1<=j<=m,i≠j。
$n,m<=10^9$

分析

一开始没看到i≠j的条件，然后就直接把原式拆成$(\sum (n mod i))*(\sum (m mod j))$，分别计算后直接加起来即可。
但是有了i≠j这个条件，我们可以选择容斥，就是用上述式子减去$\sum (n mod i)*(m mod i)$，又可以拆成$n*n-n*(m/i)*i-m*(n/i)*i+(n/i)*(m/i)*i^2$，分块计算即可。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
#define MOD 19940417
using namespace std;

int n,m,ny,ny6;

int solve(int n)
{
    int ans=(LL)n*n%MOD;
    for (int i=1,last;i<=n;i=last+1)
    {
        last=min(n/(n/i),n);
        ans=(ans-(LL)(n/i)*(i+last)%MOD*(last-i+1)%MOD*ny%MOD)%MOD;
    }
    return ans;
}

int get1(int x,int y)
{
    return (LL)(x+y)*(y-x+1)%MOD*ny%MOD;
}

int get2(int x,int y)
{
    x--;
    return ((LL)y*(y+1)%MOD*(y*2+1)%MOD*ny6%MOD-(LL)x*(x+1)%MOD*(x*2+1)%MOD*ny6%MOD+MOD)%MOD;
}

int work(int n,int m)
{
    if (n>m) swap(n,m);
    int ans=0;
    for (int i=1,last;i<=n;i=last+1)
    {
        last=min(n/(n/i),m/(m/i));
        int w1=n/i,w2=m/i,i1=get1(i,last),i2=get2(i,last);
        ans=(ans+(LL)n*m%MOD*(last-i+1)%MOD-(LL)n*w2%MOD*i1%MOD-(LL)m*w1%MOD*i1%MOD+(LL)w1*w2%MOD*i2%MOD)%MOD;
    }
    return ans;
}

int main()
{
    ny=9970209;ny6=3323403;
    scanf("%d%d",&n,&m);
    int ans=(LL)solve(n)*solve(m)%MOD-work(n,m);
    printf("%d",(ans+MOD)%MOD);
    return 0;
}