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;
    }
    
  • 相关阅读:
    查看tls指纹
    并行流
    方法引入2
    方法引入
    Optional.ofNullable
    stream.filter
    stream.skip limit
    反射
    Optional orElseGet
    nginx 预压缩(gzip)
  • 原文地址:https://www.cnblogs.com/ywwyww/p/9113284.html
Copyright © 2011-2022 走看看