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;
}