SPOJ22549 DIVFACT4 - Divisors of factorial (extreme) 题解

代码\(:\)

#include <bits/stdc++.h>
#define LL long long 
#define db long double
using namespace std;
const int M = 1000050;
int p[78498+50],cntp,pre[M];
inline void sieve(int n){
	static bool vis[M]; static int i,j; memset(vis,0,sizeof(vis));
	for (i = 2; i <= n; ++i){
		pre[i] = pre[i-1]; if (!vis[i]) p[++cntp] = i,++pre[i];
		for (j = 1; i * p[j] <= n; ++j){ vis[i*p[j]] = 1; if (i % p[j] == 0) break; }
	}
}
LL P;
inline LL mul(LL x,LL y){
	static LL ans; x %= P,y %= P,ans = x*y-P*(LL)((db)(x)*y/P);
	if (ans < 0) ans += P; else if (ans >= P) ans -= P; return ans;
}
inline LL power(LL x,LL y){ static LL r; r = 1; while (y){ if (y&1) r = mul(r,x); x = mul(x,x),y>>=1; } return r; }
inline LL calc(LL n,int p){ static LL r; r = 0; while (n >= p) n /= p,r += n; return r; }
int id1[M],id2[M],cntid; LL n,m,num[M<<1],f[M<<1];
inline int Id(LL x){ return x <= m ? id1[x] : id2[n/x]; }
inline void Min25(){
	LL l = 1,v; cntid = 0; m = min(n,(LL)(sqrt(n)) + 1);
	if (n <= 1000000) return;
	while (l <= n){
		v = n/l; if (v <= m) id1[v] = ++cntid; else id2[n/v] = ++cntid;
		num[cntid] = v,f[cntid] = v-1,l = n/v+1;
	}
	register int i,j;
	for (i = 1; i <= cntp && (LL)p[i] * p[i] <= n; ++i)
	for (v = (LL)p[i] * p[i],j = 1; j <= cntid && num[j] >= v; ++j)
		f[j] -= f[Id(num[j]/p[i])]-(i-1);
}
inline LL ask(LL x){ return (x <= 1000000) ? pre[x] : f[Id(x)]; } 
inline void solve(){
	int i; LL ans = 1%P,l = 1,r,v,lst = 0,now; Min25();
	for (i = 1; i <= cntp && (LL)p[i] * p[i] <= n; ++i) ans = mul(ans,1+calc(n,p[i])),++lst,l = p[i]+1;
	while (l <= n) v = n/l,r = n/v,now = ask(r),ans = mul(ans,power(v+1,now-lst)),l = r+1,lst = now;
	cout << ans << '\n';
}
int main(){ sieve(1000000); int T; cin >> T; while (T--) cin >> n >> P,solve(); return 0; }