最小公倍数之和 V2 51Nod

给出2个数a, b，求LCM(a,b) + LCM(a+1,b) + .. + LCM(b,b)。
例如：a = 1, b = 6，1,2,3,4,5,6 同6的最小公倍数分别为6,6,6,12,30,6，加在一起 = 66。
由于结果可能很大，输出Mod 10^9 + 7的结果。（测试数据为随机数据，没有构造特别坑人的Test）
Input
第1行：一个数T，表示后面用作输入测试的数的数量。(1 <= T <= 50000)
第2 - T + 1行：每行2个数a, b，中间用空格分隔(1 <= a <= b <= 10^9)
Output
共T行，输出对应的最小公倍数之和Mod 10^9 + 7的结果。
Sample Input
3
1 6
10 15
41 90
Sample Output
66
675
139860
关键：f(n)=∑d|nμ(d)∗d是积性函数
解：

∑i=ablcm(i,b)=∑i=abi∗bgcd(i,b)=b∗∑d∑a<=i<=bid∗(gcd(i,b)=d)=b∗∑d|b∑⌈ad⌉<=i<=bdi∗(gcd(i,bd)=1)=b∗∑d|b∑⌈ad⌉<=i<=bdi∗∑d′|gcd(i,bd)μ(d′)=b∗∑d|b∑d′|bdμ(d′)∑⌈ad⌉<=i<=bd,d′|ii=b∗∑d|b∑d′|bdμ(d′)∗d′∗(⌈⌈ad⌉d′⌉+⌊⌊bd⌋d′⌋)∗(⌊⌊bd⌋d′⌋−⌈⌈ad⌉d′⌉+1)2=b∗∑d|b∑d′|bdμ(d′)∗d′∗(⌈add′⌉+⌊bdd′⌋)∗(⌊bdd′⌋−⌈add′⌉+1)2=b∗∑d|b(⌈ad⌉+⌊bd⌋)∗(⌊bd⌋−⌈ad⌉+1)2∑d′|dμ(d′)∗d′

bonus：不过在推导的过程中，出了一些叉子，但是也得出了一下额外结论hh，我假设a是从1开始的，那么设：

sum(n)=∑i=1nlcm(i,n)

所以：

sum(n)=n∑i=1nigcd(i,n)

我们从枚举gcd的角度看看，若gcd(i,n)=g,则i=ag,n=bg,gcd(a,b)=1,我们实际上要求的就是a的这个值，但是有什么要求呢与b互质，所有说求的就是1到b内与b互质的数的和，定义这个函数为h(b)
所以：

sum(n)=n∑d|nh(nd)=n∑d|nh(d)

现在关键是有没有这个函数h（n）的简单形式，的确有且h(n)=nϕ(n)2,证明比较简单反正法即可，所以：

sum(n)=n∑d|nh(nd)=n∑d|ndϕ(d)2,(ϕ为欧拉函数)

题目代码

#include<cstdio>#include<iostream>#include<algorithm>#include<cstring>#include<cmath>#include<bitset>#include<vector>#define N 50000#define INF 10000000000#define mod 1000000007using namespace std;typedef long long  ll;int a,b;bool tag[50000];vector<int> prime;void getPrime(){    for(int i=2;i<N;i++)    {        if(!tag[i])        {            prime.push_back(i);            for(int j=i+i;j<N;j+=i)                tag[j]=true;        }    }}int fac[100];int e[100];int cnt;void getFac(int x){    cnt=0;    for(int i=0;(ll)prime[i]*prime[i]<=x;i++)    {        if(x%prime[i]==0)        {            fac[cnt]=prime[i];            e[cnt]=0;            while(x%prime[i]==0)            {                x/=prime[i];                e[cnt]++;            }            cnt++;        }    }    if(x>1)    {        fac[cnt]=x;        e[cnt]=1;        cnt++;    }}int fac2[100];int cnt2=0;int cal(int x){    int now=1;    for(int i=0;i<cnt2;i++)        now*=1-fac2[i];    return now;}ll get(int x){    int l=a/x+(a%x?1:0);    int r=b/x;    return (ll)(l+r)*(r-l+1)/2%mod;}ll res;void dfs(int cur,int temp){    if(cur==cnt)    {        ll now=get(temp)*cal(temp)%mod;        res=(res+now+mod)%mod;        return;    }    dfs(cur+1,temp);    int c=temp;    for(int i=1;i<=e[cur];i++)    {        c*=fac[cur];        fac2[cnt2++]=fac[cur];        dfs(cur+1,c);        cnt2--;    }}int main(){    int t;    getPrime();    scanf("%d",&t);    while(t--)    {        res=0;        scanf("%d%d",&a,&b);        getFac(b);        dfs(0,1);        res=res*b%mod;        printf("%lld\n",res);    }    return 0;}