测试地址：平均最小公倍数
做法：这一题需要用到杜教筛。这一题推式子的过程比较经典，是做杜教筛更难题的基础。
【数学-数论|【51Nod1227】平均最小公倍数-杜教筛】首先定义几个接下来要用到的符号和函数：
幂函数： idk(n)=nk （完全积性）
lcm(i,j) 和 gcd(i,j) ： i 和 j 的最小公倍数和最大公约数（这个不懂就……）
函数 f 和 g 的狄利克雷卷积： f?g （ (f?g)(n)=∑d|nf(d)g(n/d) ）
函数 f 和 g 的普通乘法： f?g （ (f?g)(n)=f(n)g(n) ）
[S] ：若式子 S 为真则 [S]=1 ，否则 [S]=0 。
我们用公式表达要求的函数： A(n)=1n∑ni=1lcm(n,i) ，求 F(n)=∑ni=1A(i) 。 A 不是积性函数，所以我们要想办法将式子化成积性函数前缀和的形式。
我们知道 lcm(i,j)×gcd(i,j)=ij ，所以 A(n)=1n∑ni=1nigcd(n,i)=∑ni=1igcd(n,i) 。但是 igcd(n,i) 仍然不是积性的，所以我们进行进一步化简：
A(n)=∑i=1nigcd(n,i)=∑i=1n∑d|n[gcd(n,i)=d]×id=∑d|n∑id=1nd[gcd(nd,id)=1]×id
注意到 ∑n/di/d=1[gcd(nd,id)=1]×id其实就是问小于等于 nd的数中和它互质的数的和，又因为如果 i和 x互质，那么 x?i也和 x互质（这个很显然，随便证证就可以了），推测出小于等于 x的与 x互质的数是成对出现的，而且每一对的贡献都是 x，所以小于等于 x的与 x互质的数的和为 x×φ(x)×12，特别地，小于等于 1的与 1互质的数的和为 1（若用上面公式算将算出 12），代入上面的式子中得：
A(n)=12(1+∑d|nd×φ(d))
我们发现这个式子已经有点积性函数的意思了，我们再用这个式子来算 F(n)：
F(n)=∑i=1nA(i)=12(n+∑i=1n∑d|id×φ(d))=12(n+∑i=1n∑d=1?ni?d×φ(d))
那么 G(n)=∑ni=1i×φ(i)就是一个积性函数前缀和了。
要求函数 id?φ的前缀和，考虑卷上一个好算前缀和的函数，从而得到另一个好算前缀和的函数。显然 (id?φ)?id=id2，因为 ∑d|nd×φ(d)×nd=n∑d|nφ(d)=n2，那么：
n(n+1)(2n+1)6=∑i=1ni2=∑i=1n∑d|id×φ(d)×id=∑id=1nid∑d=1?nid?d×φ(d)=∑i=1ni×G(?ni?)
看不懂上述式子的同学可以尝试从约数、倍数和贡献的方面来理解一下……
将 G(n)一项提出来，得：
G(n)=n(n+1)(2n+1)6?∑i=2ni×G(?ni?)
再结合这个式子：
F(n)=12(n+∑i=1nG(?ni?))
这样就可以使用杜教筛来计算了，使用哈希表判重，再预处理前 n23个 G值，就可以达到 O(n23)的时间复杂度了。注意因为答案要模 109+7，所以不能直接将结果乘上 12，而应该乘上 2关于模 109+7的逆元： 5×108+4。
以下是本人代码：
#include
#include
#include
#include
#include
#define hashsize 5000000
#define limit 1000000
#define mod 1000000007
#define ll long long
using namespace std;
ll a,b,phi[limit+5],sum[limit+5],h[hashsize+5]={0},f[hashsize+5];
bool prime[limit+5]={0};
int hash(ll x)
{
ll s=x%hashsize;
while(h[s]&&h[s]!=x) s=(s+1)%hashsize;
return s;
}void init()
{
for(int i=1;
i<=limit;
i++) phi[i]=i;
for(ll i=2;
i<=limit;
i++)
if (!prime[i])
{
for(ll j=1;
i*j<=limit;
j++)
{
prime[i*j]=1;
phi[i*j]=phi[i*j]/i*(i-1);
}
}
sum[0]=0;
for(int i=1;
i<=limit;
i++)
sum[i]=(sum[i-1]+i*phi[i])%mod;
}ll gcd(ll a,ll b)
{
return (b==0)?a:gcd(b,a%b);
}ll sum1(ll a,ll b)
{
return (b*(b+1)/2-a*(a-1)/2)%mod;
}ll sum2(ll x)
{
ll a=x,b=x+1,c=2*x+1,s=6,g;
g=gcd(a,s),a/=g,s/=g;
g=gcd(b,s),b/=g,s/=g;
g=gcd(c,s),c/=g,s/=g;
return (((a*b)%mod)*c)%mod;
}ll G(ll x)
{
if (x<=limit) return sum[x];
int pos=hash(x);
if (h[pos]==x) return f[pos];
ll i=2,next,s=0;
while(i<=x)
{
next=x/(x/i);
s=(s+sum1(i,next)*G(x/i))%mod;
i=next+1;
}
h[pos]=x,f[pos]=(sum2(x)-s+mod)%mod;
return f[pos];
}ll F(ll x)
{
ll i=1,next,s=0;
while(i<=x)
{
next=x/(x/i);
s=(s+(next-i+1)*G(x/i))%mod;
i=next+1;
}
return (500000004*(s+x))%mod;
}int main()
{
init();
scanf("%lld%lld",&a,&b);
printf("%lld",(F(b)-F(a-1)+mod)%mod);
return 0;
}
