ZJT's Blog

Keep your determination.

[自选题 #125] Div

2018-1-1 20:292018-1-1 20:29
集训队作业自选题杜教筛

传送门

我们考虑$(a+bi)$整除$n$意味着什么。记$n/(a+bi)=(c+di)$,则$ac-bd=n,ad+bc=0$。容易得到当$b\ne0$时,$\frac ab=-\frac cd$。$b=0$的情况我们可以特别处理,剩下就只用考虑$b\ne0$的情况了。我们设$p,q$为满足$\frac ab=-\frac cd=\frac pq,p>0$且$p,q$互质的一对数,再设$a=px,b=qx,c=py,d=-qy$,这时,$p,q,x,y$就确定了唯一的$a,b,c,d$。题目要求所有满足$1\leq ac-bd\leq n$的$a$之和,我们把$p,q,x,y$代进去,得到$1\leq(p^2+q^2)xy\leq n$。显然$q$的符号在这里没有影响,所以我们只用考虑$q>0$的情况。

$$ \begin{aligned} ans&=\sum_{p,q,x,y>0}[(p^2+q^2)xy\leq n][\gcd(p,q)=1]px\\ &=\sum_k(\sum_{p^2+q^2=k}[\gcd(p,q)=1]p)(\sum_{xy\leq \lfloor\frac nk\rfloor}x)\\ &=\sum_k(\sum_{p^2+q^2=k}[\gcd(p,q)=1]p)(\sum_{y\leq \lfloor\frac nk\rfloor}\sigma(y))\\ &=\sum_kG(k)F(\lfloor\frac nk\rfloor) \end{aligned} $$

其中$\sigma(y)$表示$y$的因子和。这个式子显然可以按$\lfloor\frac nk\rfloor$的值分段求。我们现在的任务是对所有的$x=\lfloor\frac nk\rfloor$求$G(x)$和$F(x)$。推一推发现都可以杜教筛,然后就可以在$O(n^{\frac 23})$的时间内做完了。(然而我线性预处理的地方直接用了1个log的gcd,感觉复杂度不是很靠谱?)

#include <bits/stdc++.h>
#define LL long long
#define MAXN 5000000

const LL P=1004535809;

LL n;
bool flag[MAXN];
int prime[MAXN],nump,mind[MAXN];
LL miu[MAXN];
LL sd[MAXN];
LL f0[MAXN],g0[MAXN],h0[MAXN];
LL f1[MAXN],g1[MAXN],h1[MAXN];

LL getPow(LL x,LL y){
    LL res=1;
    while(y){
        if(y&1) res=res*x%P;
        x=x*x%P;
        y>>=1;
    }
    return res;
}

const LL inv2=getPow(2,P-2);

int gcd(int x,int y){
    while(y){
        int t=y;
        y=x%y;
        x=t;
    }
    return x;
}

void init(){
    miu[1]=1;
    sd[1]=1;
    for(int i=2;i<MAXN;i++){
        if(!flag[i]){
            prime[++nump]=i;
            mind[i]=i;
            miu[i]=-1;
            sd[i]=i+1;
        }
        for(int j=1;j<=nump && i*prime[j]<MAXN;j++){
            flag[i*prime[j]]=1;
            if(i%prime[j]){
                miu[i*prime[j]]=-miu[i];
                mind[i*prime[j]]=prime[j];
                sd[i*prime[j]]=sd[i]*(prime[j]+1);
            }else{
                miu[i*prime[j]]=0;
                mind[i*prime[j]]=mind[i]*prime[j];
                if(i==mind[i]) sd[i*prime[j]]=sd[i]*prime[j]+1;
                else sd[i*prime[j]]=sd[i]/sd[mind[i]]*sd[mind[i]*prime[j]];
                break;
            }
        }
    }
    for(LL i=1;i*i<MAXN;i++){
        LL i2=i*i;
        for(LL j=1;j*j+i2<MAXN;j++){
            h0[i2+j*j]+=i;
            if(gcd(i,j)==1) g0[i2+j*j]+=i;
        }
    }
    for(int i=1;i<MAXN;i++){
        f0[i]=(f0[i-1]+sd[i])%P;
        g0[i]=(g0[i-1]+g0[i])%P;
        h0[i]=(h0[i-1]+h0[i])%P;
    }
}

LL gets_f(LL x){
    if(x<MAXN) return f0[x];
    return f1[n/x];
}

LL gets_g(LL x){
    if(x<MAXN) return g0[x];
    return g1[n/x];
}

LL gets_h(LL x){
    if(x<MAXN) return h0[x];
    return h1[n/x];
}

void gao(LL x){
    LL sh=0,sg=0,sf=0;
    for(LL i=1;i*i<x;i++)
        sh=(sh+i*(LL)(sqrt(x-i*i)+1e-5))%P;
    h1[n/x]=sh;
    for(LL i=1;i*i<=x;i++)
        sg=(sg+miu[i]*i*gets_h(x/i/i))%P;
    g1[n/x]=sg;
    LL k=1;
    while(x/k){
        LL k2=x/(x/k);
        LL _k2=k2%P,_k=k%P;
        sf=(sf+(_k2*(_k2+1)-_k*(_k-1))%P*inv2%P*(x/k%P))%P;
        k=k2+1;
    }
    f1[n/x]=sf;
}

int main(){
#ifdef DEBUG
    freopen("125.in","r",stdin);
#endif
    scanf("%lld",&n);
    init();
    for(int i=100000;i;i--)
        if(n/i>=MAXN)
            gao(n/i);
    LL ans=0,lasts=0;
    LL k=1;
    while(n/k){
        LL k2=n/(n/k);
        LL nows=gets_f(k2);
        ans=(ans+(nows-lasts)*gets_g(n/k))%P;
        k=k2+1;
        lasts=nows;
    }
    ans=(ans*2+gets_f(n))%P;
    if(ans<0) ans+=P;
    printf("%lld\n",ans);
    return 0;
}

查看详细内容

[自选题 #134] Counting Divisors (square)

2017-12-29 21:102017-12-29 21:10
集训队作业自选题杜教筛

传送门

$S_2(n)=\sum_i\sum_d[d|i^2]$

我们考虑$d|i^2$的意义,发现对于每个质因子,$d$的次数一定小于等于$i$的次数的$2$倍,这等价于$i$的次数大于等于$d$的次数除以$2$(向上取整)。我们记$g(x)$表示$x$的每个质因子的次数都除以$2$向上取整后变成的数。设当前次数为$x$,则次数为$2x-1$和$2x$的除以$2$都会变成$x$,即每个质因子有两种选择。这时我们的式子可以变成:

$$ \begin{aligned} S_2(n)&=\sum_{d'}(\sum_d[g(d)=d'])(\sum_{i=1}^n[d|i])\\ &=\sum_d2^{p_d}\lfloor\frac nd\rfloor \end{aligned} $$

其中$p_d$表示$d$的质因子个数。这条式子还可以进一步转化:

$$ \begin{aligned} S_2(n)&=\sum_d(\sum_p\mu^2(p)[p|d])\lfloor\frac nd\rfloor\\ &=\sum_p\mu^2(p)\sum_{p|d}\lfloor\frac nd\rfloor\\ &=\sum_p\mu^2(p)\sum_{d=1}^{\lfloor\frac np\rfloor}\lfloor\frac {\lfloor\frac np\rfloor}d\rfloor\\ &=\sum_p\mu^2(p)\sum_{d=1}^{\lfloor\frac np\rfloor}\lfloor\frac {\lfloor\frac np\rfloor}d\rfloor\\ \end{aligned} $$

这时式子分成了两个部分,前面是求$f(n)=\sum_{k=1}^n\mu^2 (k)$,后面是求$h(n)=\sum_{d=1}^n\lfloor\frac nd\rfloor$。这两个都可以用类似杜教筛的方法在$O(n^\frac 23)$的时间内求出所有$\lfloor\frac nx\rfloor$的函数值(线性筛预处理出前$O(n^\frac 23)$项,然后后面跟杜教筛一样处理)。大概的处理方法是:

$$ \begin{aligned} f(n)=\sum_k\mu(k)\lfloor\frac n{k^2}\rfloor \end{aligned}\\ h(n)=\sum_{k=1}^n\sum_d[d|k] $$

(其实我做的时候并没有意识到到$h(n)$可以线性预处理,然后强行写的带个log的预处理,然后理所当然的T掉了,结果后来竟然卡常卡过去了2333

(所以我这份代码复杂度还要再多个log。。。

#include <bits/stdc++.h>
#define LL long long
#define MAXN_0 30000000

int MAXN;
LL n;
bool flag[MAXN_0];
int prime[MAXN_0],nump,numls;
LL miu[MAXN_0],numd[MAXN_0];
LL s1[MAXN_0],s2[MAXN_0];
LL ls[MAXN_0][3];
int temp[MAXN_0];
int debug=0;
LL q[10010];

void init(){
    miu[1]=1;
    for(int i=2;i<MAXN;i++){
        if(!flag[i]){
            prime[++nump]=i;
            miu[i]=-1;
            numd[i]=1;
        }
        for(int j=1;j<=nump && i*prime[j]<MAXN;j++){
            flag[i*prime[j]]=1;
            if(i%prime[j]){
                miu[i*prime[j]]=-miu[i];
                numd[i*prime[j]]=numd[i]+1;
            }else{
                numd[i*prime[j]]=numd[i];
                break;
            }
        }
    }
    for(int i=1;i<MAXN;i++){
        s2[i]=miu[i]?1:0;
        if(miu[i]){
            ls[++numls][0]=i;
            ls[numls][1]=(LL)i*i;
            ls[numls][2]=miu[i];
        }
        for(int j=i;j<MAXN;j+=i) temp[j]++;
    }
    for(int i=1;i<MAXN;i++)
        s1[i]=temp[i]+s1[i-1],s2[i]+=s2[i-1];
}

inline LL gets_div(LL x){
    if(x<MAXN) return s1[x];
    LL k=1;
    LL res=0;
    while(k*k<x){
        res+=x/k;
        k++;
    }
    while(x/k){
        LL k2=x/(x/k);
        res+=x/k*(k2-k+1);
        k=k2+1;
    }
    return res;
}

inline LL gets_miu2(LL x){
    if(x<MAXN) return s2[x];
    LL res=0;
    for(int i=1;i<=numls && ls[i][1]<=x;i++)
        res+=x/ls[i][1]*ls[i][2];
    return res;
}

LL solve(){
    LL k=1;
    LL lasts=0;
    LL ans=0;
    while(n/k){
        LL k2=n/(n/k);
        LL nows=gets_miu2(k2);
        LL divs=gets_div(n/k);
        ans+=(nows-lasts)*divs;
        lasts=nows;
        k=k2+1;
    }
    return ans;
}

int main(){
#ifdef DEBUG
    freopen("134.in","r",stdin);
#endif
    int T;
    scanf("%d",&T);
    LL maxq=0;
    for(int i=1;i<=T;i++){
        scanf("%lld",q+i);
        if(q[i]>maxq) maxq=q[i];
    }
    if(maxq<=1000000000LL) MAXN=10000000;
    else if(maxq<=10000000000LL) MAXN=20000000;
    else if(maxq<=100000000000LL) MAXN=20000000;
    else MAXN=30000000;
    init();
    for(int i=1;i<=T;i++){
        n=q[i];
        printf("%lld\n",solve());
    }
    return 0;
}

查看详细内容