51nod 1965 奇怪的式子 —— min_25筛

时间:2023-11-11 12:47:32

题目:http://www.51nod.com/Challenge/Problem.html#!#problemId=1965

推式子就同这里:https://www.cnblogs.com/yoyoball/p/9196092.html

一开始想设 \( g(n,j) = \sum\limits_{i=1}^{n} [ min(i) >= p_{j} ] f(i) \),其中 \( f(i) = d(i)  \mu(i) \) 或 \( f(i) = mu(i) \),\( d(i) \) 是质因子个数;

结果发现枚举最小质因子还好(但是太慢了会TLE啊!),改成循环来筛的话很不好赋初值啊!而且筛的过程也感觉怪怪的,弄了半天还是不对...

于是就也把 \( g(n,j) \) 设成 \( g(n,j) = \sum\limits_{i=1}^{n} [ i \in prime || min(i) >= p_{j} ] f(i) \) 了,果然就对了...

这里因为 \( d(i) \mu(i) \) 不是积性函数,所以要考虑 \( d(i) \) 少加的部分,再补回来即可(就是加一个 \( \sum \mu(i) \) );

long double 快速乘;

注意多组数据!虽然也不用清空什么数组,但因为我在过程中用到了 \( s[m+1] \) 和 \( sum[m+1] \),所以要把这两个清零。

代码如下:

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
typedef long double ld;
int const xn=1e6+;
ll const mod=1e12+,md=mod-;
int m,pri[xn],cnt,sqr;
ll n,w[xn],h[xn],s[xn],ps[xn],sum[xn],f[xn],g[xn];
bool vis[xn];
ll upt(ll x){while(x>=mod)x-=mod; while(x<)x+=mod; return x;}
ll upt2(ll x){while(x>=md)x-=md; while(x<)x+=md; return x;}
/*
ll mul(ll a,ll b){ll ret=0; for(;b;b>>=1,a=upt(a+a))if(b&1)ret=upt(ret+a); return ret;}
ll mul2(ll a,ll b){ll ret=0; for(;b;b>>=1,a=upt2(a+a))if(b&1)ret=upt2(ret+a); return ret;}
*/
ll mul(ll a,ll b){ll r=(ld)a*b/mod; return a*b-r*mod;}
ll mul2(ll a,ll b){ll r=(ld)a*b/md; return a*b-r*md;}
ll pw(ll a,ll b){ll ret=; for(;b;b>>=,a=mul(a,a))if(b&)ret=mul(ret,a); return ret;}
ll pw2(ll a,ll b){ll ret=; for(;b;b>>=,a=mul2(a,a))if(b&)ret=mul2(ret,a); return ret;}
void init(int mx)
{
for(int i=;i<=mx;i++)
{
if(!vis[i])pri[++cnt]=i,ps[cnt]=upt2(ps[cnt-]+i);
for(int j=;j<=cnt&&(ll)i*pri[j]<=mx;j++)
{
vis[i*pri[j]]=;
if(i%pri[j]==)break;
}
}
}
int Id(ll x)
{
if(x>sqr)return n/x;
return m-x+;
}
ll F(ll x,int y)
{
if(pri[y]>x)return ;
ll ret=upt2(-h[Id(x)]+y-);
for(int i=y;i<=cnt&&(ll)pri[i]*pri[i]<=x;i++)
ret=upt2(ret-F(n/pri[i],i+));
return ret;
}
ll G(ll x,int y)
{
if(pri[y]>x)return ;
ll ret=upt2(-h[Id(x)]+y-);
for(int i=y;i<=cnt&&(ll)pri[i]*pri[i]<=x;i++)
ret=upt2(ret-G(n/pri[i],i+)-F(n/pri[i],i+));
return ret;
}
int main()
{
int T; scanf("%d",&T); init(1e6);
while(T--)
{
scanf("%lld",&n); sqr=sqrt(n); m=;//
int up; for(up=cnt;(ll)pri[up]*pri[up]>n;up--);
for(ll i=,j;i<=n;i=j+)
{
w[++m]=n/i; j=n/w[m];
h[m]=w[m]-;
if(w[m]&)sum[m]=mul2((w[m]+)/,w[m]);
else sum[m]=mul2(w[m]/,(w[m]+));
s[m]=upt2(sum[m]-);
}
s[m+]=; sum[m+]=;//!!!
for(int j=;j<=up;j++)
for(int i=;i<=m&&(ll)pri[j]*pri[j]<=w[i];i++)
{
int k=Id(w[i]/pri[j]);
h[i]=upt2(h[i]-h[k]+(j-));
s[i]=upt2(s[i]-mul2(s[k],pri[j]));
s[i]=upt2(s[i]+mul2(ps[j-],pri[j]));//pri[j]
}
//for(int i=1;i<=m;i++)printf("sum[%d]=%lld ",i,sum[i]); puts("");
//for(int i=1;i<=m;i++)printf("s[%d]=%lld ",i,s[i]); puts("");
ll res=;
for(ll i=,j;i<=n;i=j+)
{
j=n/(n/i);
res=upt2(res+mul2(sum[Id(n/i)],upt2(s[Id(j)]-s[Id(i-)])));//s[m+1]!
}
for(int i=;i<=up;i++)
res=upt2(res-mul2(pri[i],sum[Id(n/pri[i])]));
ll ans=pw(,res);
for(int i=;i<=up;i++)
for(ll k=,p0=pri[i];p0<=n;p0*=pri[i],k++)
{
ll r=upt2(mul2(p0,sum[Id(n/p0)])-mul2(p0*pri[i],sum[Id(n/p0/pri[i])]));//sum[m+1]!
ans=mul(ans,pw(k+,r));
}
//ans=mul(ans,pw(2,G(n,1)));
/*
for(int i=1;i<=m;i++)
{
//if(w[i]>sqr)f[i]=g[i]=upt2(-(h[Id(n)]-h[i]+1));
//else f[i]=g[i]=upt2(-(h[Id(n)]-up));
if(h[i]<=up)f[i]=g[i]=upt2(-(h[Id(n)]-up));
else f[i]=g[i]=upt2(-(h[Id(n)]-h[i]+1));
}
for(int j=up;j;j--)
for(int i=1;i<=m&&(ll)pri[j]*pri[j]<=w[i];i++)
{
int k=Id(w[i]/pri[j]);
f[i]=upt2(f[i]-f[k]-1);//mu[pj]
g[i]=upt2(g[i]-g[k]-1-f[k]-1);
}
*/
for(int i=;i<=m;i++)f[i]=g[i]=upt2(-h[i]);
for(int j=up;j;j--)
for(int i=;i<=m&&(ll)pri[j]*pri[j]<=w[i];i++)
{
int k=Id(w[i]/pri[j]);
f[i]=upt2(f[i]-f[k]-j);
g[i]=upt2(g[i]-g[k]-j-f[k]-j);
}
ans=mul(ans,pw(,g[Id(n)]));
printf("%lld\n",ans);
}
return ;
}