杂题选做：数论（一）

标签：le gcd 数论 sum sqrt int 杂题 sim

早期 shitpost 重修。

P1835 素数密度

求区间 \([L,R]\) 内的素数个数。

\(1\le L\le R<2^{31},R-L\le 10^6.\)

如果 \(x\) 是合数，那么 \([2,\sqrt x]\) 范围内一定存在 \(x\) 的因子。

所以用 \([2,\sqrt R]\) 内的素数筛一下就好了。

时间复杂度 \(O((R-L)\log \log R).\) （忽略了比较小的 \(\sqrt R\) 相关项）

为了卡常，代码中对素数 2 做了特殊处理。可能不太便于阅读。仅作展示。

#include<stdio.h>
#include<math.h>
#define int unsigned

const int N=1e6+5,M=4.7e5+5;
int l,r,ans; bool v[M],V[N];

signed main() {
    scanf("%u%u",&l,&r);
    const int n=sqrt(r);
    for (int i=2; i<=n; ++i) if (!v[i])
        for (int j=i+i; j<=n; j+=i) v[j]=1;
    for (int i=3; i<=n; i+=2) if (!v[i]) {
        const int mn=(l>i+i?(l+i-1)/i*i:i+i);
        for (int j=mn; j<=r; j+=i) V[j-l]=1;
    }
    const int len=r-l;
    for (int i=!(l&1); i<=len; i+=2)
        if (!V[i]) ++ans;
    printf("%u\n",ans+(l==2));
    return 0;
}

SP26045 GCDMAT2 - GCD OF MATRIX (hard)

给定 \(T,n,m\) ，你要回答 \(T\) 组询问。

每组询问给定 \(i1,j1,i2,j2\) ，求 \(\displaystyle \sum_{i=i1}^{i2}\sum_{j=j1}^{j2}\gcd(i,j)\) 。

答案对 \(P=10^9+7\) 取模。

\(1\le T\le 50000,1\le i1\le i2\le n\le 10^6,1\le j1\le j2\le m\le 10^6.\)

先推式子。这种直接对 \(\gcd\) 求和的情况用欧拉反演比较方便：

\[\sum_{i=i1}^{i2}\sum_{j=j1}^{j2}\sum_{d|i,d|j}\varphi(d) \]

交换求和号，化简：

\[\def\div#1#2{\left\lfloor\dfrac{#1}{#2}\right\rfloor} \sum_d \varphi(d)\left(\div{i2}d-\div{i1-1}d\right)\left(\div{j2}d-\div{j1-1}d\right)\]

预处理 \(\varphi\) ，整除分块回答询问。时间复杂度 \(O(n+T\sqrt n)\) 。

但是这个题严重卡常。

如果你的整除分块是这样写的：

    for (int l=1,r; l<=n; l=r+1) {
        r=n/(n/l);
        /* ... */
    }

那就直接暴毙了。

最大的问题在于除法的效率太低。甚至不如浮点乘快。
所以可以预处理 \(1\sim n\) 的倒数，用 double 存起来，然后化除为乘。

优化之后大概就可以通过了，但还是会跑得很慢。

为了继续优化，我们需要一个更好的整除分块写法：

    const int mid=sqrt(n);
    for (int i=1; i<=mid; ++i)
        /* 暴力计算 */
    int t=n/mid;
    for (int l=mid+1,r; l<=n; l=r+1) {
        r=n/(--t);
        /* ... */
    }

原理：

\(1\sim \sqrt n\) 用整除分块处理的话，每个点就是一个区间，显然太蠢了。不如直接暴力。
\(\sqrt n\sim n\) 范围内，\(n/l\) 的变化是连续的，也就是说 \((n/l)_{now}=(n/l)_{last}-1\) 。
因此可以直接维护 \(t=n/l\) 而避开除法。

以上讨论的是只有一个 \(n\) 的简单情况。
这题要复杂一些，整除分块中有 \(4\) 个不同的 \(n\) 。但上述优化仍然可以使用。

最后总结一下本题中有用的卡常技巧：

除法优化：预处理倒数，化除为乘；
整除分块优化：\(1\sim \sqrt n\) 暴力计算；
（其实本题中可以把这个范围适当扩大几倍，会跑得更快）
整除分块优化：不使用除法而直接维护 \(n/l\) ；
IO 优化：使用 getchar 快读；
（实测 fread/putchar/fwrite 都是负优化，很魔幻。也许是因为 SPOJ 的某些奇怪特性）
取模优化：注意到答案用 long long 存得下，所以计算时可以不取模，最后输出时再模。

放一下我目前最快的代码（2.00s）。

#include<stdio.h>
#include<math.h>
#include<ctype.h>
typedef long long ll;

const int N=1e6+5,M=8e4+5,P=1e9+7;
const double eps=1e-14;
int T,n,m,t,k,x,y,p[M],phi[N];
double I[N]; ll S[N]; bool v[N];

inline void read(int &x) {
    x=0; char ch=getchar();
    while (!isdigit(ch)) ch=getchar();
    while (isdigit(ch)) x=x*10+(ch^48),ch=getchar();
}

int main() {
    read(T),read(n),read(m);
    if (n>m) t=n,n=m,m=t;
    phi[1]=1;
    for (int i=2; i<=n; ++i) { //线性筛
        if (!v[i]) p[++k]=i,phi[i]=i-1;
        for (int j=1; j<=k; ++j) {
            const int t=i*p[j];
            if (t>n) break;
            v[t]=1;
            if (i%p[j]==0) { phi[t]=phi[i]*p[j]; break; }
            phi[t]=phi[i]*(p[j]-1);
        }
    }
    for (int i=1; i<=n; ++i) S[i]=S[i-1]+phi[i];
    for (int i=1; i<=m; ++i) I[i]=1./i+eps; //预处理倒数
    const auto dv=[](int x,int y) { return y?int(x*I[y]):N; };
    const auto min=[](int x,int y) { return x<y?x:y; };
    while (T--) {
        read(x),read(y),--x,--y,read(n),read(m);
        if (n>m) t=x,x=y,y=t,t=n,n=m,m=t;
        const int mid=min(n,sqrt(m*4)); ll s=0;
        for (int i=1; i<=mid; ++i) //暴力计算
            s+=1ll*phi[i]*(dv(n,i)-dv(x,i))*(dv(m,i)-dv(y,i));
        int tx=dv(x,mid),ty=dv(y,mid),tn=dv(n,mid),tm=dv(m,mid);
        int t1=dv(x,tx),t2=dv(y,ty),t3=dv(n,tn),t4=dv(m,tm);
        for (int l=mid+1,r; l<=n; l=r+1) { //整除分块
            if (l>t1) t1=dv(x,--tx);
            if (l>t2) t2=dv(y,--ty);
            if (l>t3) t3=dv(n,--tn);
            if (l>t4) t4=dv(m,--tm);
            r=min(min(t1,t2),min(t3,t4));
            s+=(S[r]-S[l-1])*(tn-tx)*(tm-ty);
        }
        printf("%lld\n",s%P);
    }
    return 0;
}

P3172 [CQOI2015]选数

给定 \(n,k,L,H\) ，求满足条件的长为 \(n\) 的序列 \(a\) 的个数：

\(L\le a_i\le H~~(i=1,2,\cdots,n)\)

\(\gcd(a_1,a_2,\cdots,a_n)=K\)

答案对 \(P=10^9+7\) 取模。

\(1\le n,k\le 10^9,1\le L\le H\le 10^9,H-L\le 10^5.\)

设 \(l=\lceil\frac Lk\rceil,r=\lfloor\frac Rk\rfloor\) ，可将问题简化为：
在 \([l,r]\) 范围内选取 \(a_1\sim a_n\) ，且满足它们的 \(\gcd\) 为 \(1\) ，求方案数。

莫反 + 杜教筛可做，但这样太野蛮了。下面是容斥做法。

注意到如果 \(a_1\sim a_n\) 不全相同，那么它们的 \(\gcd\) 一定不会超过 \(len=r-l\) 。

设 \(f_d\) 表示：在 \([l,r]\) 范围内选取 \(a_1\sim a_n\) ，满足它们不全相同，且 \(\gcd\) 为 \(d\) 的方案数。
那么 \(\forall d>len\) 都有 \(f_d=0\) ，所以只需考虑 \(f_1\sim f_{len}\) 。

要满足 \(\gcd(a_1,\cdots,a_n)=d\) ，首先 \(a_1\sim a_n\) 都必须是 \(d\) 的倍数。
\([l,r]\) 内 \(d\) 的倍数共有 \(cnt=\lfloor\frac rd\rfloor-\lfloor\frac{l-1}d\rfloor\) 个。
那么方案数为 \(cnt^n-cnt\) 。注意排除 \(a_1\sim a_n\) 全相等的情况。

但这样只能保证 \(d\) 是 \(a_1\sim a_n\) 的公约数。
而最大公约数可能是 \(d,2d,3d,\cdots\)
把多余的方案数都减掉就好了。

\[f_d=cnt^n-cnt-f_{2d}-f_{3d}-\cdots \]

倒推即可。

最后记得把 \(a_1\sim a_n\) 全相同的方案加回来。
如果 \(l=1\) ，则 \(a_1=a_2=\cdots=a_n=1\) 也是合法方案。
所以答案是 \(f_1+[l=1]\) 。

时间复杂度 \(O((r-l)\log (r-l))\) 。

#include<stdio.h>

const int P=1e9+7,N=1e5+5;
int n,k,l,r,f[N];

inline int power(int x,int y) {
    int s=1;
    while (y) {
        if (y&1) s=1ll*s*x%P;
        x=1ll*x*x%P; y>>=1;
    }
    return s;
}
inline int inc(int x,int y) { return (x+=y)<P?x:x-P; }
inline int dec(int x,int y) { return (x-=y)<0?x+P:x; }

int main() {
    scanf("%d%d%d%d",&n,&k,&l,&r);
    l=(l+k-1)/k,r/=k;
    const int len=r-l;
    for (int i=len; i; --i) {
        const int cnt=r/i-(l-1)/i;
        f[i]=dec(power(cnt,n),cnt);
        for (int j=i+i; j<=len; j+=i)
            f[i]=dec(f[i],f[j]);
    }
    printf("%d\n",inc(f[1],l==1));
    return 0;
}

SP5971 LCMSUM - LCM Sum / P1891 疯狂 LCM

（双倍经验）

给定 \(n\) ，求 \(\displaystyle \sum_{i=1}^n{\rm lcm}(i,n)\) 。\(T\) 组数据。

\(1\le T\le 3\times10^5,1\le n\le 10^6.\)

首先 \(\rm lcm\) 肯定要化成 \(\gcd\) 才能做。

\[n\sum_{i=1}^n\dfrac{i}{\gcd(i,n)} \]

枚举 \(d=\gcd(i,n)\) 。

\[n\sum_{d|n}\sum_{i=1}^n[\gcd(i,n)=d]\dfrac id \]

\[n\sum_{d|n}\sum_{i=1}^{n/d}[\gcd(i,n/d)=1]i \]

枚举因数时，\(d\) 和 \(n/d\) 是对称的，可以直接替换。

\[n\sum_{d|n}\sum_{i=1}^d[\gcd(i,d)=1]i \]

注意到 \(\gcd(i,d)=1\Longleftrightarrow \gcd(d-i,d)=1\) 。
所以被求和的 \(i\) 是成对出现的。

它们的平均数是 \(\dfrac d2\) ，项数则为欧拉函数 \(\varphi(d)\) 。

因此求和结果应当是 \(\dfrac{\varphi(d)d}2\) 。

注意 \(d=1\) 时求和结果为 \(1\) ，上面的分析不适用，所以需要特判。

为了方便，暂时先忽略这个特判，那么答案为

\[\dfrac12n\sum_{d|n}\varphi(d)d \]

至此显然已经有了 \(O(n\log n+T)\) 的做法，足以通过。

可用 Dirichlet 前缀和的技巧优化至 \(O(n\log \log n+T)\) 。

也可以用更加硬核的线性筛做到 \(O(n+T)\) 。（见代码）

#include<stdio.h>
typedef long long ll;

const int N=1e6+5,M=8e4+5,K=3e5+5;
int T,mx,k,p[M],n[K]; bool v[N]; ll f[N];

int main() {
    scanf("%d",&T);
    for (int i=1; i<=T; ++i)
        scanf("%d",n+i),mx=(mx>n[i]?mx:n[i]);
    f[1]=1;
    for (int i=2; i<=mx; ++i) {
        if (!v[i]) p[++k]=i,f[i]=1ll*i*i-i+1;
        for (int j=1,t; j<=k&&(t=i*p[j])<=mx; ++j) {
            v[t]=1;
            if (i%p[j]==0) {
                f[t]=f[i]+(f[i]-f[i/p[j]])*p[j]*p[j];
                break;
            }
            f[t]=f[i]*f[p[j]];
        }
    }
    for (int i=1; i<=T; ++i)
        printf("%lld\n",(f[n[i]]+1)*n[i]>>1);
    return 0;
}

P3312 [SDOI2014]数表

\(Q\) 组数据，给定 \(n,m,a\) ，求

\[\sum_{i=1}^n\sum_{j=1}^m[\sigma(\gcd(i,j))\le a]\sigma(\gcd(i,j)) \]
其中 \(\sigma(n)\) 表示 \(n\) 的正因子之和。
答案对 \(2^{31}\) 取模。

\(1\le Q\le 2\times10^4,1\le n,m\le 10^5,|a|\le 10^9.\)

首先是常规推式子。过程就省略了。

\[\sum_{d=1}^n[\sigma(d)\le a]\sigma(d)\sum_{t=1}^{n/d}\mu(t)\left\lfloor\dfrac n{dt}\right\rfloor\left\lfloor\dfrac m{dt}\right\rfloor \]

替换 \(T=dt\) ：

\[\sum_{T=1}^n\left\lfloor\dfrac nT\right\rfloor\left\lfloor\dfrac mT\right\rfloor\sum_{d|T}[\sigma(d)\le a]\sigma(d)\mu(T/d) \]

如果忽略 \(\sigma(d)\le a\) 的限制，那么可以预处理 \(F(T)=\sum\limits_{d|T}\sigma(d)\mu(T/d)\) ，然后整除分块求解。

加上限制之后，每次询问中只有一部分 \(d\) 能对 \(F\) 产生贡献。
考虑 \(a\) 增大的过程，此时满足 \(\sigma(d)\le a\) 的 \(d\) 只增不减。

将所有询问离线，并按照 \(a\) 从小到大排序。
每次询问时，相比上一次询问，只是多了一部分满足限制的 \(d\) 。
将所有 \(d\) 按照 \(\sigma(d)\) 从小到大排序，那么很容易找到这部分 \(d\) 。然后累加贡献，更新 \(F\) 即可。

但是到这里还没有结束，因为整除分块的过程中，还要用到 \(F\) 的区间和。
处理单点修改和区间求和，当然是用树状数组了。

共 \(O(n\log n)\) 次单点修改和 \(O(Q\sqrt n)\) 次区间求和。
因此时间复杂度是 \(O(n\log ^2n+Q\sqrt n\log n)\) 。

离线询问的时候有个 \(O(Q\log Q)\) 的排序，但这项好像比较小，所以就丢了。

#include<stdio.h>
#include<algorithm>
#include<math.h>
#define rep(i,l,r) for (int i=l; i<=r; ++i)

const int N=1e5+5,Q=2e4+5,inf=1e9;
const double eps=1e-14;
int T,mx,Mx,t,k,tp=1,n[Q],m[Q],a[Q],id[Q],p[Q];
bool v[N]; int mu[N],s[N],g[N],d[N];
double I[N]; int f[N],F[N],ans[Q];

inline void upd(int i,int x) {
    f[i]+=x;
    while (i<=mx) F[i]+=x,i+=(i&-i);
}
inline int qry(int i) {
    int s=0;
    while (i) s+=F[i],i&=i-1;
    return s;
}

int main() {
    scanf("%d",&T);
    const auto max=[](int x,int y) { return x>y?x:y; };
    rep(i,1,T) {
        scanf("%d%d%d",n+i,m+i,a+i);
        if (n[i]>m[i]) t=n[i],n[i]=m[i],m[i]=t;
        mx=max(mx,n[i]),Mx=max(Mx,m[i]);
    }
    rep(i,1,T) id[i]=i;
    std::sort(id+1,id+T+1,[](int i,int j) { return a[i]<a[j]; }); //离线 排序
    mu[1]=s[1]=g[1]=1,s[0]=inf;
    rep(i,2,mx) { //筛mu,sigma
        if (!v[i]) p[++k]=i,mu[i]=-1,s[i]=g[i]=i+1;
        for (int j=1; j<=k&&(t=i*p[j])<=mx; ++j) {
            v[t]=1;
            if (i%p[j]==0) {
                s[t]=s[i]/g[i]*(g[t]=g[i]*p[j]+1);
                break;
            }
            mu[t]=-mu[i],s[t]=s[i]*(g[t]=p[j]+1);
        }
    }
    rep(i,1,mx) d[i]=i;
    std::sort(d+1,d+mx+1,[](int i,int j) { return s[i]<s[j]; });
    rep(i,1,Mx) I[i]=1./i+eps;
    const auto dv=[](int x,int y) { return int(x*I[y]); };
    const auto min=[](int x,int y) { return x<y?x:y; };
    rep(i,1,T) {
        const int tn=n[id[i]],tm=m[id[i]],ta=a[id[i]];
        while (s[t=d[tp]]<=ta) { //枚举d
            for (int j=1,k=t; k<=mx; ++j,k+=t)
                if (mu[j]) upd(k,mu[j]*s[t]); //更新F
            ++tp;
        }
        const int mid=min(tn,sqrt(tm<<1));
        int S=0,ts1=qry(mid),ts2;
        rep(i,1,mid) S+=f[i]*dv(tn,i)*dv(tm,i);
        int t1=dv(tn,mid),t2=dv(tm,mid);
        int nt=dv(tn,t1),mt=dv(tm,t2);
        for (int l=mid+1,r; l<=tn; l=r+1) { //整除分块
            if (l>nt) nt=dv(tn,--t1);
            if (l>mt) mt=dv(tm,--t2);
            ts2=qry(r=min(nt,mt));
            S+=(ts2-ts1)*t1*t2,ts1=ts2;
        }
        ans[id[i]]=(S<0?S+2147483648u:S);
    }
    rep(i,1,T) printf("%d\n",ans[i]);
    return 0;
}

标签：le,gcd,数论,sum,sqrt,int,杂题,sim
From： https://www.cnblogs.com/icyM3tra/p/17087216.html