[CQOI2017]小Q的表格——反演好题

// luogu-judger-enable-o2
#include<bits/stdc++.h>
#define reg register int
#define il inline
#define numb (ch^‘0‘)
using namespace std;
typedef long long ll;
il void rd(ll &x){
    char ch;x=0;bool fl=false;
    while(!isdigit(ch=getchar()))(ch==‘-‘)&&(fl=true);
    for(x=numb;isdigit(ch=getchar());x=x*10+numb);
    (fl==true)&&(x=-x);
}
namespace Miracle{
const int N=4e6+5;
const int mod=1e9+7;
const int blo=2e3;
bool vis[N];
ll n,m;
int pri[N],tot;
int phi[N];
int g[N];
int be[N];
int le[N],ri[N],cnt;//information of blo
int sum[N],pre[N];
int v[N];
int add(int x,int y){
    return (y>=0)?(x+y>=mod?x+y-mod:x+y):(x+y<0?x+y+mod:x+y);
}
int gcd(int a,int b){
    return b?gcd(b,a%b):a;
}
void sieve(){
    phi[1]=1;
    for(reg i=2;i<=n;++i){
        //if(i>3903333)cout<<" i "<<i<<endl;
        if(!vis[i]){
            pri[++tot]=i;
            phi[i]=i-1;
        }
        for(reg j=1;j<=tot;++j){
            if((ll)pri[j]*i>n) break;
            vis[pri[j]*i]=1;
            if(i%pri[j]==0){
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
    }
    for(reg i=1;i<=n;++i) g[i]=(ll)i*i%mod*phi[i]%mod;
    for(reg i=1;i<=n;++i) g[i]=add(g[i],g[i-1]);
}
int query(int l,int r){
    if(l!=le[be[l]]) return add(add(add(sum[be[r]-1],-sum[be[l]-1]),pre[r]),-pre[l-1]);
    else return add(add(sum[be[r]-1],-sum[be[l]-1]),pre[r]);
}
int qm(int x,int y){
    int ret=1;
    while(y){
        if(y&1) ret=(ll)ret*x%mod;
        x=(ll)x*x%mod;
        y>>=1;
    }
    return ret;
}
int main(){
    rd(m);rd(n);
    sieve();
    //cout<<" after sieve "<<endl;
    for(reg i=1;i<=n;++i){
        be[i]=(i-1)/blo+1;
        v[i]=(ll)i*i%mod;
        if(be[i]!=be[i-1])le[be[i]]=i;
        ri[be[i]]=max(ri[be[i]],i);
    }
    
    cnt=be[n];
    for(reg i=1;i<=cnt;++i){
        //cout<<i<<" "<<le[i]<<" "<<ri[i]<<endl;
        pre[le[i]]=(ll)le[i]*le[i]%mod;
        for(reg j=le[i]+1;j<=ri[i];++j){
            pre[j]=add(pre[j-1],(ll)j*j%mod);
        }
        sum[i]=add(sum[i-1],pre[ri[i]]);
    }
    //cout<<" after blo "<<endl;
    ll a,b,x,k;
    while(m--){
        rd(a);rd(b);rd(x);rd(k);
        int gc=gcd(a,b);
        x%=mod;
        x=(ll)gc*gc%mod*x%mod*qm((ll)a*b%mod,mod-2)%mod;
        v[gc]=x;
        if(gc==le[be[gc]]) pre[gc]=x;
        else pre[gc]=add(pre[gc-1],x);
        for(reg i=gc+1;i<=ri[be[gc]];++i){
            pre[i]=add(pre[i-1],v[i]);
        }
        for(reg i=be[gc];i<=cnt;++i){
            sum[i]=add(sum[i-1],pre[ri[i]]);
        }
        
        ll ans=0;
        for(reg i=1,x=0;i<=k;i=x+1){
            x=k/(k/i);
            ans=add(ans,(ll)query(i,x)*g[(k/i)]%mod);
        }
        printf("%lld\n",ans);
    }
    return 0;
}

}
signed main(){
    Miracle::main();
    return 0;
}

/*
   Author: *Miracle*
   Date: 2018/12/26 20:57:31
*/