最小公倍数之和 V2 51Nod - 1190（莫比乌斯反演）

给出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)=\sum_{d|n}\mu(d)*d$是积性函数
解：

$$\sum_{i=a}^blcm(i,b) =\sum_{i=a}^b{{i*b}\over{gcd(i,b)}} =b*\sum_{d}\sum_{a<=i<=b}{{i}\over{d}}*(gcd(i,b)=d)\\ =b*\sum_{d|b}\sum_{\lceil{a\over{d}}\rceil<=i<={b\over{d}}}{i}*(gcd(i,{b\over{d}})=1)\\ =b*\sum_{d|b}\sum_{\lceil{a\over{d}}\rceil<=i<={b\over{d}}}i*\sum_{d'|gcd(i,{b\over{d}})}\mu(d')\\ =b*\sum_{d|b}\sum_{d'|{b\over{d}}}\mu(d')\sum_{\lceil{a\over{d}}\rceil<=i<={b\over{d}},d'|i}i\\ =b*\sum_{d|b}\sum_{d'|{b\over{d}}}\mu(d')*d'*{{({\lceil{{\lceil{a\over{d}}\rceil}\over{d'}}\rceil}+{\lfloor{{\lfloor{b\over{d}}\rfloor}\over{d'}}\rfloor})*({\lfloor{{\lfloor{b\over{d}}\rfloor}\over{d'}}\rfloor}-{\lceil{{\lceil{a\over{d}}\rceil}\over{d'}}\rceil}+1)}\over{2}}\\ =b*\sum_{d|b}\sum_{d'|{b\over{d}}}\mu(d')*d'*{{({{\lceil{a\over{dd'}}\rceil}}+{{\lfloor{b\over{dd'}}\rfloor}})*({{\lfloor{b\over{dd'}}}\rfloor}-{{\lceil{a\over{dd'}}\rceil}}+1)}\over{2}}\\ =b*\sum_{d|b}{{({{\lceil{a\over{d}}\rceil}}+{{\lfloor{b\over{d}}\rfloor}})*({{\lfloor{b\over{d}}}\rfloor}-{{\lceil{a\over{d}}\rceil}}+1)}\over{2}}\sum_{d'|d}\mu(d')*d'\\$$

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

$$sum(n)=\sum_{i=1}^nlcm(i,n)$$
所以： $$sum(n)=n\sum_{i=1}^n\frac{i}{gcd(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\sum_{d|n}h(\frac{n}{d})=n\sum_{d|n}h(d)$$
现在关键是有没有这个函数h（n）的简单形式，的确有且 $h(n)=\frac{n\phi(n)}{2}$,证明比较简单反正法即可，所以： $$sum(n)=n\sum_{d|n}h(\frac{n}{d})=n\sum_{d|n}\frac{d\phi(d)}{2},(\phi为欧拉函数)$$
题目代码

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<bitset>
#include<vector>
#define N 50000
#define INF 10000000000
#define mod 1000000007
using 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;
}