【BZOJ】2956: 模积和 反演

题目地址：https://www.lydsy.com/JudgeOnline/problem.php?id=2956

#include <iostream>
#include <stdio.h>
#include  <stdlib.h>
#include  <string.h>
#include  <math.h>
using namespace std;
#define mod 19940417
#define niyuan  3323403
long long sum(long long l,long long r)
{
    return (l+r)*(r-l+1)/2%mod;
}
long long sum2(long long n)
{
    return n*(n+1)%mod*(2*n+1)%mod*niyuan%mod;
}
long long he1(long long n)
{
    long long ans=0;
    for(long long l=1,r; l<=n; l=r+1)
        {
            r=n/(n/l);
            ans+=sum(l,r)*(n/l);
            ans=ans%mod;
        }
    return ans;
}
long long he11(long long n,long long m)
{
    long long ans=0;
    for(long long l=1,r; l<=n; l=r+1)
        {
            r=m/(m/l);
            ans+=sum(l,r)*(m/l);
            ans=ans%mod;
        }
    return ans;
}
long long he2(long long n)
{
    return n*n%mod;
}
long long he3(long long n,long long m)
{
    long long ans=0;
    for(long long l=1,r; l<=n; l=r+1)
        {
            r=min(n/(n/l),m/(m/l));
            ans+=(n/l)*(m/l)*((sum2(r)-sum2(l-1)+mod)%mod);///
            ans=ans%mod;
        }
    return ans;
}
long long he4(long long n,long long m)///为什么不能化简
{
    long long ans=0;
    long long l=1,r;
    while (l<=n)
        {
            r=min(n/(n/l),m / (m / l));
            long long  x=n / l,y=m / l;
            ans=(ans+sum(l,r)*(n*y%mod+m*x%mod))%mod;
                 l=r+1;
        }
        return ans;
}
int main()
{
    long long n,m,ans1,ans2,x1,x2,sum,sum1;
    scanf("%lld%lld",&n,&m);
    if(n>m)
        swap(n,m);
    ans1=he1(n);
    ans2=he1(m);
    x1=(he2(n)-ans1+mod)%mod;
    x2=(he2(m)-ans2+mod)%mod;
    sum=x1*x2%mod;

    sum1=he2(n)*m%mod+he3(n,m);
    sum1=sum1%mod;
    sum1=sum+he4(n,m)-sum1+mod;
    sum1=sum1%mod;
    printf("%lld\n",sum1);
    return 0;
}