承接一下洛咕上的题解，这里基本就是谈谈优化，放个代码的

我们发现这里的常数主要来自于除法，那么我们优化除法次数，把所有的 \(n/1...n/s\) （\(s=\sqrt n\)）存下来，然后归并排（其实就是 merge 一下），最后 unique 去个重，然后就可以进行小常数的数论分块了

//by Judge
#pragma GCC optimize("Ofast")
#include<bits/stdc++.h>
#define Rg register
#define fp(i,a,b) for(Rg int i=(a),I=(b)+1;i<I;++i)
#define fd(i,a,b) for(Rg int i=(a),I=(b)-1;i>I;--i)
#define ll long long
using namespace std;
const int mod=1e9+7;
const int M=1e6+3;
typedef int arr[M*20];
#ifndef Judge
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
#endif
char buf[1<<21],*p1=buf,*p2=buf;
inline int read(){ Rg int x=0; char c=getchar(); for(;!isdigit(c);c=getchar());
    for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x;
} char sr[1<<21],z[21];int CCF=-1,Z;
inline void Ot(){fwrite(sr,1,CCF+1,stdout),CCF=-1;}
inline void print(int x,char chr='\n'){
    if(CCF>1<<20)Ot(); while(z[++Z]=x%10+48,x/=10);
    while(sr[++CCF]=z[Z],--Z); sr[++CCF]=chr;
} int cnt; ll phi[M]; arr v,p,A,B,C;
#define get(x,y) (x/y?x/(x/y):c)
#define swap(a,b) (a^=b^=a^=b)
inline int Min(Rg int x,Rg int y){return x<y?x:y;}
int main(){
    int T=read(),n=read(),m=read(); if(n>m) swap(n,m); v[1]=phi[1]=1;
    for(int i=2;i<=n;++i){ if(!v[i]) p[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*p[j]<=n;++j){ v[i*p[j]]=1;
            if(!(i%p[j])){ phi[i*p[j]]=phi[i]*p[j]; break;
            } phi[i*p[j]]=phi[i]*(p[j]-1);
        } phi[i]+=phi[i-1];
    } ++T;
    while(--T){ Rg int a=read()-1,b=read()-1,c=read(),d=read();
        if(c>d) swap(c,d),swap(a,b); Rg ll ans=0;
        *A=0; fd(i,sqrt(a),1) A[++*A]=a/i;
        *B=0; fd(i,sqrt(b),1) B[++*B]=b/i;
        merge(A+1,A+*A+1,B+1,B+*B+1,C+1),*C=*A+*B;
        *B=0; fd(i,sqrt(c),1) B[++*B]=c/i;
        merge(C+1,C+*C+1,B+1,B+*B+1,A+1),*A=*C+*B;
        *B=0; fd(i,sqrt(d),1) B[++*B]=d/i;
        merge(A+1,A+*A+1,B+1,B+*B+1,C+1),*C=*A+*B;
        *B=0; fp(i,1,sqrt(d)) B[++*B]=i;
        merge(C+1,C+*C+1,B+1,B+*B+1,A+1),*A=*C+*B;
        *A=unique(A+1,A+*A+1)-1-A,*C=*A,*A=0;
        Rg int l,r; fp(i,1,*C) l=A[i-1]+1,r=A[i],
            ans+=(phi[r]-phi[l-1])*(c/l-a/l)*(d/l-b/l); print(ans%mod);
    } return Ot(),0;
}

话说 min25 也太神仙了，交了一发，结果快成那个样子，惊