【51nod 1190】最小公倍数之和 V2

Description

给出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）

Solution

$$ans=\sum_{i=a}^{b}{i*b\over gcd(i,b)}=b*\sum_{d|b}\sum_{i=\lceil a/d\rceil}^{\lfloor b/d \rfloor}i*\sum_{d'|gcd(i,b/d)}\mu(d')=b*\sum_{d|b} \sum_{d'|\lfloor b/d \rfloor}\mu(d')\sum_{i=\lceil a/d\rceil,d'|i}^{\lfloor b/d \rfloor}i$$ $$ans=b*\sum_{d|b}\sum_{d'|\lfloor b/d \rfloor}\mu(d'){(\lceil{a\over dd'}\rceil+\lfloor{b\over dd'}\rfloor)(\lfloor{b\over dd'}\rfloor-\lceil{a\over dd'}\rceil+1)\over 2}*d'$$
令T=d*d’
$$ans=\sum_{T|b}{(\lceil{a\over T}\rceil+\lfloor{b\over T}\rfloor)(\lfloor{b\over T}\rfloor-\lceil{a\over T}\rceil+1)\over 2}\sum_{d'|T}\mu(d')*d'$$ 设 $f(T)=\sum_{d'|T}\mu(d')*d'$。
由于 $\mu(d)$和d都是积性函数，所以当gcd（p,q）=1时，f(p*q)=f(p)*f(q)。
并且由 $\mu(d)$的性质可知，T里的每个质数指数最多为1，否则 $\mu(d)=0$。
我们可以对b分解质因数，枚举b的约数T，同时就能算出 $\sum_{d'|T}\mu(d')*d'$，就如有 $x*p^k=T$p这个质因子时，你会有选与不选这两个决策，不选的话贡献为1，选的话贡献为-p,然后就 $f(T)=f(x)*(1-p)$。 $\sum_{d'|T}\mu(d')*d'$的计算类似于(a+b)*(d+c)=ac+ad+bc+bd。时间复杂度O( $2^8*test$).

Code

#include<iostream>
#include<math.h>
#include<string.h>
#include<stdio.h>
#include<algorithm>
#define ll long long
using namespace std;
const ll maxn=1e5,mo=1e9+7,mo2=5e8+4;
int d[maxn],bz[maxn+5],b1[maxn],c[maxn];
ll n,i,t,j,k,l,a,b,z,ans,p;
void dg(int x,int y,ll z){
    ll t=1,k;
    if (y>b1[0]){
        t=(a%x)?a/x+1:a/x;
        k=b/x;
        ans+=(t+k)*(k-t+1)*z;
        return;
    }
    dg(x,y+1,z);
    for (k=1;k<=c[y];k++)
        t*=b1[y],dg(x*t,y+1,z*(1-b1[y]));
}
int main(){
    scanf("%lld",&n);
    for (i=2;i<=maxn;i++){
        if (!bz[i]) d[++d[0]]=i;
        for (j=1;j<=d[0];j++){
            if (i*d[j]>maxn) break;
            bz[i*d[j]]=1;
            if (i%d[j]==0) break;
        }
    }
    for (i=1;i<=n;i++){
        scanf("%lld%lld",&a,&b);
        z=b;b1[0]=0;
        for (j=1;d[j]*d[j]<=z;j++){
            if (z%d[j]) continue;
            b1[++b1[0]]=d[j];c[b1[0]]=0;
            while (z%d[j]==0) z/=d[j],c[b1[0]]++;
        }
        if (z>1)b1[++b1[0]]=z,c[b1[0]]=1;
        ans=0;
        dg(1,1,1);ans=(ans%mo+mo)%mo*mo2%mo;
        ans=b%mo*ans%mo;
        printf("%lld\n",ans);
    }
}