zoukankan      html  css  js  c++  java
  • 【51NOD1965】奇怪的式子 min_25筛

    题目描述

      给你(n),求

    [prod_{i=1}^n{sigma_0(i)}^{i+mu(i)} ]

      对({10}^{12}+39)取模。

      (sigma_0(i))表示约数个数。

    题解

      把式子拆成两部分:

    [prod_{i=1}^n{sigma_0(i)}^{i+mu(i)}=prod_{i=1}^n{sigma_0(i)}^{i}prod_{i=1}^n{sigma_0(i)}^{mu(i)} ]

      先看前面这部分

    [egin{align} S(n)&=sum_{i=1}^ni=frac{n(n+1)}{2}\ prod_{i=1}^nsigma_0(i)^i&=prod_{p}prod_{k=1}{(k+1)}^{p^kS(frac{n}{p^k})-p^{k+1}S(frac{n}{p^{k+1}})}\ &=prod_{pleq sqrt n}prod_{k=1}{(k+1)}^{p^kS(frac{n}{p^k})-p^{k+1}S(frac{n}{p^{k+1}})} imes prod_{p>sqrt n}2^{pS(frac{n}{p})} end{align} ]

      前面那个式子可以直接计算,后面那个式子可以用min_25筛筛出质数和,然后套上一个数论分块解决。

      再看后面那部分:

      因为(i)的每个质数只出现了一次,所以(sigma_0(i)=2^{i ext{的质因子个数}})

    [egin{align} prod_{i=1}^nsigma_0(i)^{mu(i)}&=prod_{i=1}^n2^{i ext{的质因子个数} imes mu(i)} end{align} ]

      然后设

    [egin{align} F1_{n,j}&=sum_{i=2}^n[i ext{为质数或}i ext{的最小因子都不小于}p_j]mu(i) imes i ext{的质因子个数}\ F2_{n,j}&=sum_{i=2}^n[i ext{为质数或}i ext{的最小因子都不小于}p_j]mu(i)\ F1_{n,j}&=-([p_j,n] ext{之间的质数个数})-sum_{igeq j,p_i^2leq n}(F1_{frac{n}{p_i},i+1}+F2_{frac{n}{p_i},i+1})\ F2_{n,j}&=-([p_j,n] ext{之间的质数个数})-sum_{igeq j,p_i^2leq n}F2_{frac{n}{p_i},i+1}\ end{align} ]

      直接上min_25筛就好了。

      乘法取模可以用黑科技,或者__int128

      时间复杂度:(O(frac{n^frac{3}{4}}{log n}))

    代码

    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<cstdlib>
    #include<ctime>
    #include<utility>
    #include<cmath>
    #include<functional>
    using namespace std;
    typedef long long ll;
    typedef unsigned long long ull;
    typedef pair<int,int> pii;
    typedef pair<ll,ll> pll;
    void sort(int &a,int &b)
    {
    	if(a>b)
    		swap(a,b);
    }
    void open(const char *s)
    {
    #ifndef ONLINE_JUDGE
    	char str[100];
    	sprintf(str,"%s.in",s);
    	freopen(str,"r",stdin);
    //	sprintf(str,"%s.out",s);
    //	freopen(str,"w",stdout);
    #endif
    }
    int rd()
    {
    	int s=0,c,b=0;
    	while(((c=getchar())<'0'||c>'9')&&c!='-');
    	if(c=='-')
    	{
    		c=getchar();
    		b=1;
    	}
    	do
    	{
    		s=s*10+c-'0';
    	}
    	while((c=getchar())>='0'&&c<='9');
    	return b?-s:s;
    }
    void put(int x)
    {
    	if(!x)
    	{
    		putchar('0');
    		return;
    	}
    	static int c[20];
    	int t=0;
    	while(x)
    	{
    		c[++t]=x%10;
    		x/=10;
    	}
    	while(t)
    		putchar(c[t--]+'0');
    }
    int upmin(int &a,int b)
    {
    	if(b<a)
    	{
    		a=b;
    		return 1;
    	}
    	return 0;
    }
    int upmax(int &a,int b)
    {
    	if(b>a)
    	{
    		a=b;
    		return 1;
    	}
    	return 0;
    }
    const ll p=1000000000039;
    const ll p1=1000000000038;
    ll mul(ll a,ll b)
    {
    	return (__int128)a*b%p;
    //	a%=p;
    //	b%=p;
    //	return (a*b-(ll)((long double)a/p*b+1e-8)*p)%p;
    }
    ll mul1(ll a,ll b)
    {
    	return (__int128)a*b%p1;
    //	a%=p1;
    //	b%=p1;
    //	return (a*b-(ll)((long double)a/p1*b+1e-8)*p1)%p1;
    }
    ll fp(ll a,ll b){ll s=1;for(;b;b>>=1,a=mul(a,a))if(b&1)s=mul(s,a);return s;}
    ll fp1(ll a,ll b){ll s=1;for(;b;b>>=1,a=mul1(a,a))if(b&1)s=mul1(s,a);return s;}
    const int M=350010;
    const int m=350000;
    ll n;
    ll sq;
    int b[M],pri[M],cnt;
    ll f1[M],f2[M],g1[M],g2[M];
    void init()
    {
    	for(int i=2;i<=m;i++)
    	{
    		if(!b[i])
    			pri[++cnt]=i;
    		for(int j=1;j<=cnt&&i*pri[j]<=m;j++)
    		{
    			b[i*pri[j]]=1;
    			if(i%pri[j]==0)
    				break;
    		}
    	}
    	pri[cnt+1]=m+1;
    }
    ll s[50];
    void gao()
    {
    	for(int i=2;i<=m;i++)
    	{
    		f1[i]=(ll)(i+2)*(i-1)/2%p1;
    		g1[i]=i-1;
    	}
    	for(int i=1;n/i>m;i++)
    	{
    		f2[i]=((n/i)&1?mul1((n/i-1)/2,(n/i+2)):mul1((n/i+2)/2,(n/i-1)));
    		g2[i]=n/i-1;
    	}
    	for(int i=1;i<=cnt;i++)
    	{
    		int j;
    		ll x1=f1[pri[i]-1];
    		ll x2=g1[pri[i]-1];
    		for(j=1;n/j/pri[i]>m&&n/j>=(ll)pri[i]*pri[i];j++)
    		{
    			f2[j]=(f2[j]-pri[i]*(f2[j*pri[i]]-x1))%p1;
    			g2[j]-=g2[j*pri[i]]-x2;
    		}
    		for(;n/j>m&&n/j>=(ll)pri[i]*pri[i];j++)
    		{
    			f2[j]=(f2[j]-pri[i]*(f1[n/j/pri[i]]-x1))%p1;
    			g2[j]-=g1[n/j/pri[i]]-x2;
    		}
    		for(j=m;j>=(ll)pri[i]*pri[i];j--)
    		{
    			f1[j]=(f1[j]-pri[i]*(f1[j/pri[i]]-x1))%p1;
    			g1[j]-=g1[j/pri[i]]-x2;
    		}
    	}
    }
    pll getf(ll x,int y)
    {
    	if(x<=1||x<pri[y])
    		return pll();
    	ll s1=-((x<=m?g1[x]:g2[n/x])-g1[pri[y]-1]);
    	ll s2=s1;
    	for(int i=y;i<=cnt&&(ll)pri[i]*pri[i]<=x;i++)
    	{
    		pll v=getf(x/pri[i],i+1);
    		s1=(s1-v.first-v.second)%p1;
    		s2-=v.second;
    	}
    	return pll(s1,s2);
    }
    ll get(ll x)
    {
    	return x<=m?f1[x]:f2[n/x];
    }
    ll sum(ll x)
    {
    	return x&1?mul1((x+1)/2,x):mul1(x/2,x+1);
    }
    void solve()
    {
    	scanf("%lld",&n);
    	sq=0;
    	while((sq+1)*(sq+1)<=n)
    		sq++;
    	memset(s,0,sizeof s);
    	gao();
    	pll res=getf(n,1);
    	s[2]=(s[2]+res.first)%p1;
    	for(ll i=1,j;i<=n;i=j+1)
    	{
    		j=n/(n/i);
    		if(j>sq)
    			s[2]=(s[2]+mul1(get(j)-get(max(i-1,sq)),sum(n/i)))%p1;
    	}
    	for(int i=1;i<=cnt&&(ll)pri[i]*pri[i]<=n;i++)
    	{
    		ll s1=pri[i],s2=(ll)pri[i]*pri[i];
    		for(int j=1;s1<=n;s1=s2,s2*=pri[i],j++)
    			s[j+1]=(s[j+1]+mul1(s1,sum(n/s1))-mul1(s2,sum(n/s2)))%p1;
    	}
    	ll ans=1;
    	for(int i=2;i<=50;i++)
    	{
    		s[i]=(s[i]+p1)%p1;
    		ans=mul(ans,fp(i,s[i]));
    	}
    	ans=(ans+p)%p;
    	printf("%lld
    ",ans);
    }
    int main()
    {
    	open("51nod1965");
    	init();
    	int t;
    	scanf("%d",&t);
    	while(t--)
    		solve();
    	return 0;
    }
    
  • 相关阅读:
    POJ3693 Maximum repetition substring —— 后缀数组 重复次数最多的连续重复子串
    SPOJ
    POJ2774 Long Long Message —— 后缀数组 两字符串的最长公共子串
    POJ3261 Milk Patterns —— 后缀数组 出现k次且可重叠的最长子串
    POJ1743 Musical Theme —— 后缀数组 重复出现且不重叠的最长子串
    SPOJ
    AC自动机小结
    HDU3247 Resource Archiver —— AC自动机 + BFS最短路 + 状压DP
    POJ1625 Censored! —— AC自动机 + DP + 大数
    Herding
  • 原文地址:https://www.cnblogs.com/ywwyww/p/9113284.html
Copyright © 2011-2022 走看看