51nod 1965 奇怪的式子 —— min_25筛
题目：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+5; ll const mod=1e12+39,md=mod-1; 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<0)x+=mod; return x;} ll upt2(ll x){while(x>=md)x-=md; while(x<0)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=1; for(;b;b>>=1,a=mul(a,a))if(b&1)ret=mul(ret,a); return ret;} ll pw2(ll a,ll b){ll ret=1; for(;b;b>>=1,a=mul2(a,a))if(b&1)ret=mul2(ret,a); return ret;} void init(int mx) { for(int i=2;i<=mx;i++) { if(!vis[i])pri[++cnt]=i,ps[cnt]=upt2(ps[cnt-1]+i); for(int j=1;j<=cnt&&(ll)i*pri[j]<=mx;j++) { vis[i*pri[j]]=1; if(i%pri[j]==0)break; } } } int Id(ll x) { if(x>sqr)return n/x; return m-x+1; } ll F(ll x,int y) { if(pri[y]>x)return 0; ll ret=upt2(-h[Id(x)]+y-1); for(int i=y;i<=cnt&&(ll)pri[i]*pri[i]<=x;i++) ret=upt2(ret-F(n/pri[i],i+1)); return ret; } ll G(ll x,int y) { if(pri[y]>x)return 0; ll ret=upt2(-h[Id(x)]+y-1); for(int i=y;i<=cnt&&(ll)pri[i]*pri[i]<=x;i++) ret=upt2(ret-G(n/pri[i],i+1)-F(n/pri[i],i+1)); return ret; } int main() { int T; scanf("%d",&T); init(1e6); while(T--) { scanf("%lld",&n); sqr=sqrt(n); m=0;// int up; for(up=cnt;(ll)pri[up]*pri[up]>n;up--); for(ll i=1,j;i<=n;i=j+1) { w[++m]=n/i; j=n/w[m]; h[m]=w[m]-1; if(w[m]&1)sum[m]=mul2((w[m]+1)/2,w[m]); else sum[m]=mul2(w[m]/2,(w[m]+1)); s[m]=upt2(sum[m]-1); } s[m+1]=0; sum[m+1]=0;//!!! for(int j=1;j<=up;j++) for(int i=1;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-1)); s[i]=upt2(s[i]-mul2(s[k],pri[j])); s[i]=upt2(s[i]+mul2(ps[j-1],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=0; for(ll i=1,j;i<=n;i=j+1) { j=n/(n/i); res=upt2(res+mul2(sum[Id(n/i)],upt2(s[Id(j)]-s[Id(i-1)])));//s[m+1]! } for(int i=1;i<=up;i++) res=upt2(res-mul2(pri[i],sum[Id(n/pri[i])])); ll ans=pw(2,res); for(int i=1;i<=up;i++) for(ll k=1,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+1,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=1;i<=m;i++)f[i]=g[i]=upt2(-h[i]); 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]-j); g[i]=upt2(g[i]-g[k]-j-f[k]-j); } ans=mul(ans,pw(2,g[Id(n)])); printf("%lld\n",ans); } return 0; }
