51nod 1584 weighted divisor and Mobius inversion

meaning of the title

Last year, tangjz loved to do number theory problems, but a year later, tangjz was not so good at it.
In sorting out the previous test questions, he found such a question: "find σ σ (i), where 1 ≤ i ≤ N, σ (i) represents the sum of the divisors of i."
Now that he has grown up, the problem has become difficult, so please help him to solve a number problem.
He needs you to find the value of the following expression:
∑ni=1∑nj=1max(i,j)∗σ(i∗j)
Where max(i,j) is the maximum value in i and j, and σ (i ⋅ j) is the sum of the divisors of i ⋅ J.
For example, when N=2, the answer should be 1 ⋅ σ (1 ⋅ 1) + 2 ⋅ σ (1 ⋅ 2) + 2 ⋅ σ (2 ⋅ 1) + 2 ⋅ σ (2 ⋅ 1) + 2 ⋅ σ (2 ⋅ 2) = 27.
He found that the answer is a little big, so you just need to tell him the value of answer module 100000007.
1≤T≤50000,1≤N≤1000000

Analysis

Consider removing the annoying max first.

ans=2∗∑i=1ni∑j=1iσ(ij)−∑i=1ni∗σ(i2)

Let's continue pushing
sum=∑i=1ni∑j=1iσ(ij)

=∑i=1ni∑j=1i∑p|i∑q|j[(p,q)=1]pjq

We can get it by inversion
sum=∑d=1nμ(d)d2∑i=1⌊nd⌋i∑j=1i∑p|i∑q|jpq

=∑d=1nμ(d)d2∑i=1⌊nd⌋i∗σ(i)∑j=1nσ(j)

This is already a block form, but in this problem still can T.

Code

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

typedef long long LL;

const int N=1000005;
const int MOD=1000000007;

int sig[N],mu[N],s[N],sig2[N],tot,prime[N],low[N],tmp[N],ans[N];
bool not_prime[N];

int read()
{
    int x=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}

void get_prime(int n)
{
    sig[1]=mu[1]=sig2[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i]) prime[++tot]=i,sig[i]=i+1,mu[i]=-1,low[i]=i,sig2[i]=(1+i+(LL)i*i%MOD)%MOD;
        for (int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                low[i*prime[j]]=low[i]*prime[j];
                sig[i*prime[j]]=((LL)sig[i]*prime[j]+(LL)sig[i/low[i]])%MOD;
                sig2[i*prime[j]]=((LL)sig2[i]*prime[j]%MOD*prime[j]+(LL)sig2[i/low[i]]+(LL)sig2[i/low[i]]*prime[j])%MOD;
                break;
            }
            low[i*prime[j]]=prime[j];
            sig[i*prime[j]]=(LL)sig[i]*sig[prime[j]]%MOD;
            sig2[i*prime[j]]=(LL)sig2[i]*sig2[prime[j]]%MOD;
            mu[i*prime[j]]=-mu[i];
        }
    }
    for (int i=1;i<=n;i++) mu[i]+=mu[i]<0?MOD:0,s[i]=(s[i-1]+sig[i])%MOD,tmp[i]=((LL)tmp[i-1]+(LL)i*sig2[i])%MOD;
    for (int i=1;i<=n;i++)
        for (int j=i;j<=n;j+=i)
        {
            ans[j]+=(LL)mu[i]*i%MOD*j%MOD*sig[j/i]%MOD*s[j/i]%MOD;
            ans[j]-=ans[j]>=MOD?MOD:0;
        }
    for (int i=1;i<=n;i++) ans[i]+=ans[i-1],ans[i]-=ans[i]>=MOD?MOD:0;
}

int main()
{
    get_prime(1000000);
    int T=read();
    for (int i=1;i<=T;i++)
    {
        int x=read();
        printf("Case #%d: %d\n",i,(ans[x]*2%MOD+MOD-tmp[x])%MOD);
    }
    return 0;
}

Posted on Sun, 03 May 2020 23:04:19 -0700 by kilby