51nod 1227 average minimum common multiple Mobius inversion + Dujiao screen

Keywords: REST

meaning of the title

Lcm(a,b) represents the minimum common multiple of a and B, and A(n) represents the average number of Lcm(n,i) (1 < = I < = n),
For example: A(4) = (Lcm(1,4) + Lcm(2,4) + Lcm(3,4) + Lcm(4,4)) / 4 = (4 + 4 + 12 + 4) / 4 = 6.
F(a, b) = A(a) + A(a + 1) + ...... A(b). (F(a,b) = ∑A(k), a <= k <= b)
For example: F(2, 4) = A(2) + A(3) + A(4) = 2 + 4 + 6 = 12.
Give a,b and calculate f (a,b). Because the result may be very large, output the result of F (a,b)% 1000000007.
1 <= a <= b <= 10^9

Analysis

What we need is

∑i=1n∑j=1ijgcd(i,j)

=∑d=1n∑i=1⌊nd⌋∑j=1ij∗[gcd(i,j)=1]

So far, there are two methods. One is to use φ (d) to do it. It will be much better. But I didn't expect that, so I used the inversion method.
It's easy to get a wave by routine inversion
ans=12∗∑k=1nμ(k)∗k∗∑d=1⌊nk⌋∑i=1⌊nkd⌋(i2+i)

Let s(d) = μ (d) * d,f(d)=d * (d+1) * 2 (2d+1)6+n * (n+1)2
that
ans=12∗∑k=1ns(k)∗∑d=1⌊nk⌋f(⌊nkd⌋)

Let T=kd to get
ans=12∗∑T=1nf(⌊nT⌋)∗∑d|Ts(d)

Now the problem is how to find the prefix sum of g(T) = ∑ d|T μ (d) * D.
It's not hard to find a picture
∑d=1ng(d)=∑d=1nS(⌊nd⌋)

Where S(d) = ∑ ni=1s(i).
It is noted that both s(d) and g(d) are product functions, and a part of them can be pretreated by linear sieve, and the rest can be pretreated by Dujiao sieve.

Code

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

typedef long long LL;

const int MOD=1000000007;
const int ny2=500000004;
const int ny6=166666668;
const int N=10000005;

int l,r,prime[N],tot,s[N],g[N];
bool not_prime[N];
map<int,int> w,wg;

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

int get_s(int n)
{
    if (n<=10000000) return s[n];
    if (w[n]) return w[n];
    int ans=1;
    for (int i=2,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        (ans-=(LL)((LL)last*(last+1)/2%MOD-(LL)i*(i-1)/2%MOD)*get_s(n/i)%MOD)%=MOD;
    }
    ans+=ans<0?MOD:0;
    return w[n]=ans;
}

int get_g(int n)
{
    if (n<=10000000) return g[n];
    if (wg[n]) return wg[n];
    int ans=0;
    for (int i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        ans+=(LL)(last-i+1)*get_s(n/i)%MOD;
        ans-=ans>=MOD?MOD:0;
    }
    return wg[n]=ans;
}

int get_f(int n)
{
    int ans=(LL)n*(n+1)%MOD*(n*2+1)%MOD*ny6%MOD+(LL)n*(n+1)/2%MOD;
    ans-=ans>=MOD?MOD:0;
    return ans;
}

int solve(int n)
{
    int ans=0;
    for (int i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        ans+=(LL)(get_g(last)+MOD-get_g(i-1))*get_f(n/i)%MOD;
        ans-=ans>=MOD?MOD:0;
    }
    ans+=ans<0?MOD:0;
    return (LL)ans*ny2%MOD;
}

int main()
{
    get_prime(10000000);
    scanf("%d%d",&l,&r);
    printf("%d",(solve(r)-solve(l-1)+MOD)%MOD);
    return 0;
}

Posted by rocksolidhq on Fri, 01 May 2020 20:23:57 -0700