首先观察$b*f(a,a+b)=(a+b)*f(a,b)$这个东西
可以化成$frac{f(a,a+b)}{a+b}=frac{f(a,b)}{b}$，发现这类似辗转相除求gcd
而我们两边同乘一个a就能得到$frac{f(a)}{a}$是个定值的这个结论
那么有$f(a,b)=frac{a*b}{gcd(a,b)^2}*f(gcd(a,b),gcd(a,b))$
为了方便现在设$gcd(i,j)=g$，现在把这个东西放进原来的式子里
$sumlimits_{i=1}^ksumlimits_{j=1}^kf(i,j)$
$=sumlimits_{i=1}^ksumlimits_{j=1}^kfrac{i*j}{g^2}*f(g,g)$
改为枚举$g$，把$f(g,g)$前提
$=sumlimits_{d=1}^kf(d,d)sum_{d|i}sum_{d|j}[gcd(i,j)==d]frac{i*j}{d^2}$
熟悉的，后面改为枚举$d$
$=sumlimits_{d=1}^kf(d,d)sumlimits_{i=1}^{leftlfloorfrac{k}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{k}{d} ight floor}[gcd(i,j)==1]i*j$
开始搞后面那个玩意
$sum_{i=1}^nsum_{j=1}^n[gcd(i,j)=1]i*j$
$=2sum_{i=1}^nsum_{j=1}^{i}[gcd(i,j)==1]i*j-sum_{i=1}^ni$
我们知道$n$以内$(n>=2)$和$n$互质的数的和是$frac{n*φ(n)}{2}$
，然后$i=1$的时候就是$1$
$=sum_{i=1}^ni^2φ(i)$
回到原来的式子
$=sumlimits_{d=1}^kf(d,d)sumlimits_{i=1}^{leftlfloorfrac{k}{d} ight floor}i^2φ(i)$
可做了，前面那个预处理，后面的分块
1 #include<cmath> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 using namespace std; 6 const int N=4000005,Sq=2005,mod=1e9+7; 7 int pri[N],npr[N],phi[N]; 8 int blo[N],stp[Sq],edp[Sq]; 9 int val[N],bsum[N],asum[Sq],func[N]; 10 int a,b,k,n,m,cnt,tot,sqr; long long x,ans; 11 int exGCD(int a,int b,int &x,int &y) 12 { 13 if(!b) {x=1,y=0; return a;} 14 int g=exGCD(b,a%b,y,x); y-=a/b*x; 15 return g; 16 } 17 int GCD(int a,int b) 18 { 19 return b?GCD(b,a%b):a; 20 } 21 int Inv(int nm,int md) 22 { 23 int xx,yy; 24 exGCD(nm,md,xx,yy); 25 return (xx%md+md)%md; 26 } 27 void Prework() 28 { 29 phi[1]=1,npr[1]=true; 30 sqr=sqrt(n)+5,stp[cnt=1]=1; 31 for(int i=2;i<=n;i++) 32 { 33 if(!npr[i]) 34 pri[++tot]=i,phi[i]=i-1; 35 for(int j=1;j<=tot&&1ll*i*pri[j]<=n;j++) 36 { 37 npr[i*pri[j]]=true; 38 phi[i*pri[j]]=phi[i]*pri[j]; 39 if(i%pri[j]) phi[i*pri[j]]-=phi[i]; 40 else break; 41 } 42 } 43 for(int i=1;i<=n;i++) 44 { 45 val[i]=1ll*i*i%mod,blo[i]=(i-1)/sqr+1; 46 if(i%sqr==0) edp[cnt++]=i,stp[cnt]=i+1; 47 func[i]=(func[i-1]+1ll*i*i%mod*phi[i]%mod)%mod; 48 } 49 (n%sqr)?edp[cnt]=n:cnt--; 50 for(int i=1;i<=cnt;i++) 51 { 52 bsum[stp[i]]=val[stp[i]]; 53 for(int j=stp[i]+1;j<=edp[i];j++) 54 bsum[j]=(bsum[j-1]+val[j])%mod; 55 asum[i]=(asum[i-1]+bsum[edp[i]])%mod; 56 } 57 } 58 void Change(int a,int b,int x) 59 { 60 int g=GCD(a,b),bl=blo[g],v=Inv(1ll*a*b%mod,mod); 61 val[g]=1ll*g*g%mod*x%mod*v%mod; 62 bsum[g]=(g==stp[bl])?val[g]:(bsum[g-1]+val[g])%mod; 63 for(int i=g+1;i<=edp[bl];i++) 64 bsum[i]=(bsum[i-1]+val[i])%mod; 65 for(int i=bl;i<=cnt;i++) 66 asum[i]=(asum[i-1]+bsum[edp[i]])%mod; 67 } 68 long long Query(int pos) 69 { 70 return pos?(asum[blo[pos]-1]+bsum[pos])%mod:0; 71 } 72 int main() 73 { 74 scanf("%d%d",&m,&n),Prework(); 75 while(m--) 76 { 77 scanf("%d%d%lld%d",&a,&b,&x,&k); 78 x%=mod,Change(a,b,x),ans=0; 79 for(int i=1,j;i<=k;i=j+1) 80 j=k/(k/i),ans+=1ll*(Query(j)-Query(i-1)+mod)%mod*func[k/i]%mod,ans%=mod; 81 printf("%lld ",ans); 82 } 83 return 0; 84 }