• luogu P5325 Min_25筛


    LINK:Min_25筛

    新版感觉有点鬼畜 而且旧版的也够用了至少.

    这个并不算很简单也不算很困难的知识点 学起来还是很麻烦的。

    (误入了很多dalao的blog 说的云里雾里的 甚是懵逼 这里推荐几个blog一起看 能看出很多门道

    网上资源辣么多 我自然也不会去写一个非常正常的学习笔记辣.. 只会写几个容易疑惑的地方。

    注意 学会 和会写代码是两码事 因为代码中有一些细节需要细细揣摩。

    关于g数组的求出 其转移静下心来理解还是可以看懂的这里不再赘述。

    注意 为了方便(f(1))最后考虑.

    (maxx=sqrt n)

    (g_{n,j}=sum_{i=1}^{n}f(i)[iin P|min_p(i)>p_j])

    关于g数组 第二维的范围显然只有(maxx)的大小。

    因为大于maxx的转移都是承接上一个的状态 所以不必要求 或者可以理解欧拉筛的思想 筛到maxx大小的质数之后 所有数字都被筛完了。

    考虑第一维 这是根据 我们后面的S组数来定的 但是注意到后面的S每次只会用到(frac{n}{p_i})

    虽然我们不知道有多大但是数量级还是 是固定的 因为(frac{n}{i})最多只有(maxx)数量级 注意 真正的大小应该是(2cdot maxx)

    从状态转移方程 来看 第二维可以滚动。看似状态数量很多 实际上关于 (p_jcdot p_j>n)的那些状态都是没有必要的 所以状态数量在一个可控范围。

    一个笔者还未搞懂的问题 以后可能会完成回答:实际上在求g数组的时候利用到了完全积性函数的性质 当f不是完全积性函数的时候

    不过很多blog上说可以这样做 且 只有(g_{n,|P|})的位置上的值有效 这个地方还不太清楚为什么。

    下面关于求答案的部分。

    (S(n,j)=sum_{i=1}^{n}f(i)[min_p(i)>=p_j])

    容易得到(S(n,1)+f(1)) 且有(S(i,j)=g_{i,|P|}-sum_{j-1}+sum_{k>=j}sum_{e}f(p^e)(S(frac{i}{p^e},k+1)+[e eq 1]))

    然后就没了 复杂度的大概就是1s 1e10,3s 1e11的样子。

    code
    //#include<bits/stdc++.h>
    #include<iostream>
    #include<iomanip>
    #include<cstdio>
    #include<ctime>
    #include<cstdlib>
    #include<cctype>
    #include<cstring>
    #include<cmath>
    #include<string>
    #include<utility>
    #include<queue>
    #include<vector>
    #include<algorithm>
    #include<deque>
    #include<stack>
    #include<list>
    #include<bitset>
    #include<set>
    #include<map>
    #define INF 1000000000000000000ll
    #define rep(p,n,i) for(int i=p;i<=n;++i)
    #define fep(n,p,i) for(int i=n;i>=p;--i)
    #define vep(p,n,i) for(int i=p;i<n;++i)
    #define db double
    #define get(x) x=read()
    #define put(x) printf("%d
    ",x)
    #define pb push_back
    #define ll long long
    #define db double
    #define mod 1000000007
    #define putl(x) printf("%lld
    ",x)
    using namespace std;
    char *fs,*ft,buf[1<<15];
    inline char getc()
    {
    	return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;
    }
    inline ll read()
    {
    	ll x=0,f=1;char ch=getc();
    	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getc();}
    	while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getc();}
    	return x*f;
    }
    const int MAXN=200010,N=5,inv3=333333336;
    ll n;
    int maxx,cnt,top;
    ll g1[MAXN],g2[MAXN],w[MAXN];
    int sum1[MAXN],sum2[MAXN],id1[MAXN],id2[MAXN];
    int v[MAXN],p[MAXN];
    inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
    inline void prepare()
    {
    	rep(2,maxx,i)
    	{
    		if(!v[i])
    		{
    			v[i]=p[++top]=i;
    			sum1[top]=add(sum1[top-1],i);
    			sum2[top]=add(sum2[top-1],(ll)i*i%mod);
    		}
    		rep(1,top,j)
    		{
    			if(p[j]>maxx/i)break;
    			v[p[j]*i]=p[j];
    			if(v[i]==p[j])break;
    		}
    	}
    }
    inline void calc()
    {
    	for(ll i=1,j;i<=n;i=j+1)
    	{
    		w[++cnt]=n/i;
    		j=n/w[cnt];
    		g1[cnt]=w[cnt]%mod;
    		g2[cnt]=g1[cnt]*(g1[cnt]+1)/2%mod*(2*g1[cnt]+1)%mod*inv3%mod-1;//平方和.
    		g1[cnt]=(1+g1[cnt])*g1[cnt]/2%mod-1;
    		if(w[cnt]<=maxx)id1[w[cnt]]=cnt;
    		else id2[j]=cnt;
    	}
    	rep(1,top,i)
    	{
    		rep(1,cnt,j)
    		{
    			if((ll)p[i]*p[i]>w[j])break;
    			ll cc=w[j]/p[i],k;
    			if(cc<=maxx)k=id1[cc];else k=id2[n/cc];
    			g1[j]=(g1[j]-p[i]*(g1[k]-sum1[i-1]))%mod;
    			g2[j]=(g2[j]-p[i]*(g2[k]-sum2[i-1])%mod*p[i])%mod;
    		}
    	}
    }
    inline int S(ll x,int y)
    {
    	if(x<=1||p[y]>x)return 0;
    	int k=x<=maxx?id1[x]:id2[n/x];
    	ll ans=(g2[k]-g1[k]-sum2[y-1]+sum1[y-1])%mod;
    	for(int i=y;i<=top;++i)
    	{
    		if((ll)p[i]*p[i]>x)break;
    		ll pe=p[i],pp=pe*pe;
    		for(int e=1;pp<=x;++e,pe=pp,pp*=p[i])
    		{
    			ans=(ans+pe%mod*((pe-1)%mod)%mod*S(x/pe,i+1)+pp%mod*((pp-1)%mod)%mod)%mod;
    		}
    	}
    	return ans;
    }
    signed main()
    {
    	freopen("1.in","r",stdin);
    	get(n);maxx=(int)sqrt(n*1.0)+1;
    	prepare();calc();put((S(n,1)+1+mod)%mod);
    	return 0;
    }
    
  • 相关阅读:
    windows7通过Dns.GetHostAddresses(Dns.GetHostName())获得ipv6地址转换到ipv4
    题解 P3829 【[SHOI2012]信用卡凸包】
    点积与叉积
    点分治
    珂朵莉树
    NOIP2020模拟赛(二十五)7.26 结题报告
    树连剖分
    NOIP2020模拟赛(拾)解题报告
    题解 P2538 【[SCOI2008]城堡】
    模拟退火
  • 原文地址:https://www.cnblogs.com/chdy/p/13204742.html
Copyright © 2020-2023  润新知