zoukankan      html  css  js  c++  java
  • Min25筛

    Min25筛


    我是沙雕。。。

    yyb博客蒯的

    要求:(sum_{i=1}^nF(x))

    (F(x))是积性函数。

    (Min25)筛能用的前提:质数处的(f(p))值是关于(p)的多项式,质数次方处的(f(p^e))值 可以快速计算。

    预处理

    设完全积性函数(F'(x)),在质数处取值(F(p)=F'(p))

    预处理一个(g)函数。(g(n,j)=sum_{i=1}^nF'(i)[iinmathbb{P} ext{ or }min_{p|i}p>P_j])

    也就是(i)是质数或者(i)最小值因子(>P_j)的贡献(注意求的是(F')

    转移:

    (g(n,j)=egin{cases}g(n,j-1) (P_j^2>n)\g(n,j-1)-F'(P_j)[g(frac{n}{P_j},j-1)-g(P_{j-1},j-1)](P_jleq n)end{cases})

    考虑(g(n,j-1)-g(n,j))。如果(P_j^2>n)显然没有多余贡献。

    否则会多余的一定是(P_j)的倍数,那么因为(F')是完全积性函数,可以将多余的数全部除以(P_j),并乘上贡献(F'(P_j))。那么减掉的是(F'(P_j))乘【(frac{n}{P_j})以内最小质因数(geq P_j)的数的贡献和】。这个值就是(g(frac{n}{P_j},j-1)-g(P_{j-1},j-1))。前面多算了(<P_j)的质数,后面减掉了这一部分。

    可以看出只有(sqrt n)内的质数有用,线性筛到(sqrt n)就行了。

    后面的(g(P_{j-1},j-1))按照定义等于(sum_{i=1}^{P_{j-1}}F'(i))。线性筛时预处理一下就行了。

    这个(g)后面只用了质数处的(F'),所以是对的而且很妙

    求值

    (S(n,j)=sum_{i=1}^nF(i)[min_{p|i}p>P_j]),即最小质因数大于(P_j)的所有数的(F)和,那么要求的答案是(S(n,0))

    (S)的递推式为:(S(n,j)=g(n,|P|)-sum_{i=1}^{j-1}f(P_i)+sum_{p_k^eleq n,k>j}F(p_k^e)(S(frac{n}{p_k^e},k)+[e>1]))

    前面部分是质数的贡献,后面是合数的

    后面的(sum)意思就是枚举合数最小的质因子及次数,因为是积性函数可以直接乘起来。([e>1])意思是如果(e>1)那么有一个(p_k^e)没算进答案((S)不计算(1)的贡献,但(e=1)(p_k^e)就是质数已经在前面算过了),要算进答案。

    代码

    一些细节:

    • (n)一直整除一个数之后还是(n)整除一个数,也就是只需要用到所有(frac ni)。开两个数组(a1,a2),对于(sqrt n)的向(a1)存,下标为(frac ni),否则向(a2)存,下标为(i)。两个数组都是(sqrt n)的。
    • 由于(S)没有算(1),最后还要加上(1)的贡献。
    #include<bits/stdc++.h>
    #define il inline
    #define vd void
    #define mod 1000000007
    typedef long long ll;
    il ll gi(){
        ll x=0,f=1;
        char ch=getchar();
        while(!isdigit(ch))f^=ch=='-',ch=getchar();
        while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
        return f?x:-x;
    }
    ll n;
    int qt;
    int pr[100010],yes[100010],cnt,gp1[100010],gp2[100010];
    ll w[200010],sw;
    int g1[200010],g2[200010];
    int id1[100010],id2[100010];
    il int S(ll x,int y){
        if(pr[y]>=x)return 0;
        int p=x<=qt?id1[x]:id2[n/x];
        int ret=((0ll+g2[p]-g1[p]-(gp2[y]-gp1[y]))%mod+mod)%mod;
        for(int i=y+1;i<=cnt&&1ll*pr[i]*pr[i]<=x;++i){
            ll pe=pr[i];
            for(int e=1;pe<=x;++e,pe*=pr[i]){
                int o=pe%mod;
                ret=(ret+1ll*o*(o-1)%mod*(S(x/pe,i)+(e!=1)))%mod;
            }
        }
        return ret;
    }
    int main(){
    #ifdef XZZSB
        freopen("in.in","r",stdin);
        freopen("out.out","w",stdout);
    #endif
        n=gi();qt=sqrt(n);
        yes[1]=1;
        for(int i=2;i<=qt;++i){
            if(!yes[i])pr[++cnt]=i;
            for(int j=1;j<=cnt&&i*pr[j]<=qt;++j){
                yes[i*pr[j]]=1;
                if(i%pr[j]==0)break;
            }
        }
        for(int i=1;i<=cnt;++i)gp1[i]=(gp1[i-1]+pr[i])%mod,gp2[i]=(gp2[i-1]+1ll*pr[i]*pr[i])%mod;//递推时用的质数处F前缀和
        for(ll l=1,r;l<=n;l=r+1){//预处理(n/i) & 计算g的j=0边界,边界就是F'的前缀和,由于质数处F(p)是多项式所以可以快速算
            r=n/(n/l);w[++sw]=n/r;
            g1[sw]=w[sw]%mod;
            g2[sw]=(1ll*g1[sw]*(g1[sw]+1)%mod*(g1[sw]*2+1)%mod*166666668%mod-1)%mod;//1..w[sw]平方和
            g1[sw]=(1ll*g1[sw]*(g1[sw]+1)%mod*500000004-1)%mod;//1..w[sw]等差数列和
            if(n/r<=qt)id1[n/r]=sw;else id2[r]=sw;
        }
        //j从1到|P|递推g
        for(int i=1;i<=cnt;++i){
            ll sqr_pi=1ll*pr[i]*pr[i];
            for(int j=1;j<=sw&&sqr_pi<=w[j];++j){
                ll p=w[j]/pr[i];
                p=(p<=qt?id1[p]:id2[n/p]);//定位p的坐标
                g1[j]=(g1[j]-1ll*pr[i]*(g1[p]-gp1[i-1]+mod)%mod+mod)%mod;
                g2[j]=(g2[j]-1ll*pr[i]*pr[i]%mod*(g2[p]-gp2[i-1]+mod)%mod+mod)%mod;
            }
        }
        printf("%d
    ",(S(n,0)+1)%mod);
        return 0;
    }
    
  • 相关阅读:
    Spring基础知识
    Hibernate基础知识
    Struts2基础知识
    在eclipse里头用checkstyle检查项目出现 File contains tab characters (this is the first instance)原因
    java后台获取cookie里面值得方法
    ckplayer 中的style.swf 中的 style.xml 中的修改方法
    java hql case when 的用法
    Windows下Mongodb安装及配置
    Mongodb中经常出现的错误(汇总)child process failed, exited with error number
    Mac 安装mongodb
  • 原文地址:https://www.cnblogs.com/xzz_233/p/11085338.html
Copyright © 2011-2022 走看看