zoukankan      html  css  js  c++  java
  • 洲阁筛详解

    你还真信了

    链接

    这筛对积性函数的要求不同于杜教筛,只消函数在自变量为质数或质数整数幂时是一个低阶多项式即可。以下n<=1e11。

    首先有一个性质:1~n的每个数,大于$sqrt{n}$的质因子只有一个。根据是否有大于$sqrt{n}$的质因子,再根据他是积性函数,得

    $sum_{i=1}^nf(i)=sum_{i=1}^{n}[i has no prime greater than sqrt{n}]f(i)(1+sum_{i=sqrt{n}}^{left lfloor frac{n}{i} ight floor}[i is prime]f(i))$

    好的!来观察他。首先$sum_{i=1}^{sqrt{n}}[i has no prime greater than sqrt{n}]$是可以线性筛的。

    于是我们就差两部分:Part1:$sum_{i=sqrt{n}}^{left lfloor frac{n}{i} ight floor}[i is prime]f(i)$,Part2:$sum_{i=sqrt{n}}^n[i has no prime greater than sqrt{n}]f(i)$,并希望用$leq sqrt{n}$的质数来处理他们。

    Part1:$g(i,j)$--$[1,j]$中与前$i$个质数互质的数的k次幂和。有递推式$g(i,j)=g(i-1,j)-p_i^kg(i-1,left lfloor frac{j}{p_i} ight floor)$。直接算的复杂度是$frac{n}{log(n)}$。

    注意到$p_{i+1}>j$时$g(i,j)=1$,紧接着$p_i>left lfloor frac{j}{p_i} ight floor$即$p_i^2>j$时$g(i-1,left lfloor frac{j}{p_i} ight floor)=1$,因此$g(i,j)=g(i-1,j)-p_i^k$。

    也就是说算到$p_i^2>j$时就不必再算了,后面要用某个$g(i,j)$时,若$p_{i+1}>j$则$g(i,j)=1$,否则用$g(i0_j,j)$减去一段$p_i^k$即可,其中$i0_j$表示最慢枚举到$j$的$i$。

    Part2:$h(i,j)$--$[1,j]$中用前$i$个质数凑起来的数的$f$值的和。有$h(i,j)=h(i-1,j)+sum_{c=1}^{p_i^c<=j}f(p_i^c)h(i-1,left lfloor frac{j}{p_i^c} ight floor)$

    然而这里没有上面那个用平方来卡$j$上限的性质!!咋整!!

    反过来。把前$i$个改成后$i$个,注意是$leq sqrt{n}$的后$i$个。递推式同理。但!

    $p_i>j$时$h(i,j)=1$,$p_i^2>j$时用后$i-1$个质数凑得$leq j$的部分只有1,因此$h(i,j)=h(i-1,j)+f(p_i)$。棒!用类似方法解决。

    注意$h$实际上回答了$sum_{i=1}^{n}[i has no prime greater than sqrt{n}]f(i)$,因此最后加上即可,括号里那个1就不必理了。然后$g$中都是包含1的,但是我不要,因此在算答案的时候记得-1。

    以下是那道万年不变的例题的代码。

     1 //#include<iostream>
     2 #include<cstring>
     3 #include<cstdlib>
     4 #include<cstdio>
     5 //#include<map>
     6 #include<math.h>
     7 //#include<time.h>
     8 //#include<complex>
     9 #include<algorithm>
    10 using namespace std;
    11 
    12 #define LL long long
    13 int T; LL n;
    14 int m=320000;
    15 #define maxn 640011
    16 LL prime[maxn],f[maxn],sf[maxn],s[maxn],i0[maxn],who[maxn],g[maxn],ff[maxn],mi[maxn],ci[maxn]; 
    17 int lp,lw; bool notprime[maxn];
    18 
    19 #define maxh 1000007
    20 struct Hash
    21 {
    22     int first[maxh],le; struct Edge{LL to,v; int next;}edge[maxn];
    23     void clear(LL n) {le=2; for (LL i=1;i<=n;i=n/(n/i)+1) first[n/i%maxh]=0;}
    24     void insert(LL x,LL y) {int z=x%maxh; Edge &e=edge[le]; e.to=x; e.v=y; e.next=first[z]; first[z]=le++;}
    25     LL find(LL x) {int z=x%maxh; for (int i=first[z];i;i=edge[i].next) if (edge[i].to==x) return edge[i].v; return -1;}
    26 }h;
    27 
    28 void pre(int &n)
    29 {
    30     lp=0; f[1]=sf[1]=1;
    31     for (int i=2;i<=n;i++)
    32     {
    33         if (!notprime[i]) {prime[++lp]=i; f[i]=4; mi[i]=i; ci[i]=1;}
    34         sf[i]=sf[i-1]+f[i]; s[i]=s[i-1]+!notprime[i];
    35         for (int j=1,tmp;j<=lp && 1ll*i*prime[j]<=n;j++)
    36         {
    37             notprime[tmp=i*prime[j]]=1;
    38             if (i%prime[j]) {f[tmp]=f[i]*4; mi[tmp]=prime[j]; ci[tmp]=1;}
    39             else {f[tmp]=f[i/mi[i]]*(3*ci[i]+4); mi[tmp]=mi[i]*prime[j]; ci[tmp]=ci[i]+1; break;}
    40         }
    41     }
    42     lp--; n=prime[lp+1]-1;
    43 }
    44 
    45 void mg()
    46 {
    47     lw=0; h.clear(n);
    48     for (LL i=1;i<=n;i=n/(n/i)+1) lw++,who[lw]=g[lw]=n/i,h.insert(who[lw],lw);
    49     memset(i0,0,sizeof(i0));
    50     for (int i=1;i<=lp;i++)
    51         for (int j=1;j<=lw && prime[i]*prime[i]<=who[j];j++)
    52         {
    53             int k=h.find(who[j]/prime[i]);
    54             g[j]-=g[k]-(i-1-i0[k]);
    55             i0[j]=i;
    56         }
    57 }
    58 
    59 void mff()
    60 {
    61     for (int j=1;j<=lw;j++) ff[j]=1;
    62     for (int i=lp;i;i--)
    63         for (int j=1;j<=lw && 1ll*prime[i]*prime[i]<=who[j];j++)
    64         {
    65             if (prime[i+1]>who[j]) ff[j]=1;
    66             else if (prime[i+1]*prime[i+1]>who[j]) ff[j]=(s[min((LL)m,who[j])]-s[prime[i+1]-1])*4+1;
    67             LL tmp=prime[i],k2;
    68             for (int c=1;tmp<=who[j];tmp*=prime[i],c++)
    69             {
    70                 LL k=who[j]/tmp;
    71                 if (prime[i+1]>k) k2=1;
    72                 else if (prime[i+1]*prime[i+1]>k) k2=(s[min((LL)m,k)]-s[prime[i+1]-1])*4+1;
    73                 else k2=ff[h.find(k)];
    74                 ff[j]+=k2*(3*c+1);
    75             }
    76         }
    77 }
    78 
    79 int main()
    80 {
    81     scanf("%d",&T);
    82     pre(m);
    83     while (T--)
    84     {
    85         scanf("%lld",&n);
    86         if (n<=m) {printf("%lld
    ",sf[n]); continue;}
    87         LL ans=0,last; mg(); mff();
    88         for (int i=1,id;i<=m;i=last+1)
    89         {
    90             LL tmp=n/i,kk;
    91             if (prime[lp+1]>tmp) kk=1;
    92             else id=h.find(tmp),kk=g[id]-(lp-i0[id]);
    93             ans+=(sf[last=min((LL)m,n/tmp)]-sf[i-1])*(kk-1)*4;
    94         }
    95         printf("%lld
    ",ans+ff[1]);
    96     }
    97     return 0;
    98 }
    View Code

    配俩图,表示$g$和$h$的递推顺序。

  • 相关阅读:
    怎么往mac中finder个人收藏里添加文件夹
    UIView 动画
    添加.pch文件
    声明属性的关键字
    创建app前的环境配置/AppIcon/启动图片
    修改动画的旋转支点
    实现自定义xib和storyboard的加载,
    Quartz2D绘图 及实例:下载进度
    帧动画
    在职研究生第一单元第二单元第三单元第四单元是什么?
  • 原文地址:https://www.cnblogs.com/Blue233333/p/8490272.html
Copyright © 2011-2022 走看看