[NOI P simulation] large factorial (Stirling number)

Keywords: number theory Math

meaning of the title

seek 16 16 In hexadecimal, n ! n! n! Remove the tail 0 0 Take mold after 0 2 64 2^{64} 264.

n < 2 64 n<2^{64} n<264

altogether T ≤ 10 T\leq10 T ≤ 10 sets of data.

Problem solution

There is an at the end 0 0 0 means that it can be 2 4 2^4 Divide by 24. Let's take this number 2 2 Find the times of 2 and take the mold 4 4 4. Take it again. We can use the similar method of extending Lucas to continuously divide the target by 2 to make statistics, that is, the sequential statistics can be calculated 2 1 , 2 2 , 2 3 , . . . , 2 63 2^1,2^2,2^3,...,2^{63} 21, 22, 23,..., 263, add up. Eliminate all 2 2 After 2, the answer becomes
∏ i = 0 f ( n 2 i ) \prod_{i=0} f(\frac{n}{2^i}) i=0∏​f(2in​)

among f ( n ) = ∏ i = 1 ( n + 1 ) / 2 − 1 ( 2 i + 1 ) f(n)=\prod_{i=1}^{(n+1)/2-1}(2i+1) f(n)=∏i=1(n+1)/2−1​(2i+1) .

Now our problem becomes begging f ( n ) f(n) f(n), let's expand it,
oh   Yo ,   Collapse   Collapse   la _ {oh \, yo, \, collapse \, collapse \, La} Ouch, collapse

there n n n is too big. It seems that we have to change our mind. We make 2 = x 2=x 2=x, during mold taking 2 64 2^{64} 264, we'll find that only x x x 1 ∼ 63 1\sim63 1 ∼ 63 times. Let's expand:
∑ i = 0 63 g ( n + 1 2 − 1 , i ) ⋅ x i \sum_{i=0}^{63} g(\frac{n+1}{2}-1,i)\cdot x^i i=0∑63​g(2n+1​−1,i)⋅xi

there g ( n , m ) g(n,m) g(n,m) indicates from 1 ∼ n 1\sim n 1 ∼ n successful selection m m Sum of all schemes multiplied by m numbers.

Unexpected directions: g ( n , m ) g(n,m) g(n,m) is recursive, and it is not difficult to find that the recursive formula is
g ( n , m ) = g ( n − 1 , m ) + g ( n − 1 , m − 1 ) × n g(n,m)=g(n-1,m)+g(n-1,m-1)\times n g(n,m)=g(n−1,m)+g(n−1,m−1)×n

This formula looks a little familiar. Let's think of another recursive formula:
s ( n , m ) = s ( n − 1 , m − 1 ) + s ( n − 1 , m ) × ( n − 1 ) s(n,m)=s(n-1,m-1)+s(n-1,m)\times(n-1) s(n,m)=s(n−1,m−1)+s(n−1,m)×(n−1)

So it is not difficult to find the first kind of Stirling number s ( n , m ) s(n,m) s(n,m) and g ( n , m ) g(n,m) Relationship between g(n,m):
g ( n , m ) = s ( n + 1 , n − m + 1 ) g(n,m)=s(n+1,n-m+1) g(n,m)=s(n+1,n−m+1)

Great, we can directly calculate the first kind of Stirling number:
oh   Yo ,   Collapse   Collapse   la _ {oh \, yo, \, collapse \, collapse \, La} Ouch, collapse

Still can't, because here n + 1 n+1 n+1 is too big, m m m is small, but n − m + 1 n-m+1 n − m+1 is too large.

Considering that the definition of the first kind of Stirling number is that several arrays are arranged in a circle, we can use m m m. Because the size exceeds 1 1 The number of circle permutations of 1 is less than or equal to m m m . So:
s ( n + 1 , n − m + 1 ) = ∑ i = 0 m d p [ m ] [ i ] ⋅ ( n + 1 m + i ) s(n+1,n-m+1)=\sum_{i=0}^{m} dp[m][i]\cdot{n+1\choose m+i} s(n+1,n−m+1)=i=0∑m​dp[m][i]⋅(m+in+1​)

among d p [ i ] [ j ] dp[i][j] dp[i][j] indicates i + j i+j i+j number formation j j j not in size 1 1 The number of schemes arranged in circles of 1,

This can also be recursive. Discuss whether the arrangement size of the last circle is greater than 2 2 2. The recursive formula can be obtained
d p [ i ] [ j ] = ( d p [ i − 1 ] [ j ] + d p [ i − 1 ] [ j − 1 ] ) × ( i + j − 1 ) dp[i][j]=\big(dp[i-1][j] + dp[i-1][j-1]\big)\times\big(i+j-1\big) dp[i][j]=(dp[i−1][j]+dp[i−1][j−1])×(i+j−1)

The calculation of combinatorial number is a difficult problem because it needs division, but 2 64 2^{64} 264 is not a good module.

So, we're going to take all the elements and denominators 2 2 2 are put forward, then divided, and then put into the combinatorial number.

Overall, the time complexity O ( T log ⁡ 4 n ) O(T\log^4n) O(Tlog4n), passable.

CODE

#include<set>
#include<map>
#include<queue>
#include<stack>
#include<random>
#include<vector>
#include<bitset>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define MAXN 10000005
#define LL long long
#define ULL unsigned long long
#define DB double
#define lowbit(x) (-(x) & (x))
#define ENDL putchar('\n')
#define FI first
#define SE second
LL read() {
    LL f=1,x=0;int s = getchar(); 
    while(s < '0' || s > '9') {if(s<0)return -1;if(s=='-')f=-f;s = getchar();}
    while(s >= '0' && s <= '9') {x = (x<<3) + (x<<1) + (s^48); s = getchar();}
    return f*x;
}
const char sxt[16] = {'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'};
ULL readull() {
	ULL x=0;int s = getchar(); 
	while(s < '0' || s > '9') {if(s < 0) return 0;s = getchar();}
	while(s >= '0' && s <= '9') {x=(x<<3)+(x<<1)+(s^48); s = getchar();}
	return x;
}
void putpos(ULL x) {if(!x)return ;putpos(x>>4);putchar(sxt[x&15]);}
void putnum(ULL x) {
    if(!x) {putchar('0');return ;}
    return putpos(x);
}
void AIput(ULL x,int c) {putnum(x);putchar(c);}

int n,m,s,o,k;
ULL inv[128];
ULL facf(ULL x) {
	if(!x) return 1ull;
	ULL dp[128][128] = {};
	dp[0][0] = 1;
	for(int i = 1;i <= 126 && i <= x;i ++) {
		for(int j = 1;j <= i;j ++) {
			dp[i][j] = dp[i-1][j] * (i-1+j) + dp[i-1][j-1] * (i-1+j);
		}
	}
	if(!(x&1)) x --;
	ULL n = x/2 + 1;
	ULL ans = 0,pw = 1;
	for(int i = 0;i <= 63 && i <= n;i ++,pw <<= 1) {
		ULL stl = 0;
		ULL c1[128],c2 = 1,tlc2 = 1;
		for(int j = 1;j <= i*2 && j <= n;j ++) {
			c1[j] = n-j+1;
			int b = j,cn = 0;
			while(!(b&1)) b>>=1,cn++;
			for(int k = 1;k <= j && cn;k ++) {
				while(cn && !(c1[k]&1)) c1[k]>>=1,cn --;
			}
			c2 *= inv[b];
			tlc2 *= b;
			if(j > i) {
				ULL C = c2;
				for(int k = 1;k <= j;k ++) C *= c1[k];
				stl += C * dp[i][j-i];
			}
		}
		if(i == 0) stl = 1;
		ans += pw * stl;
	}
	return ans;
}
int main() {
	freopen("multiplication.in","r",stdin);
	freopen("multiplication.out","w",stdout);
	inv[0] = 1;
	for(int i = 1;i <= 128;i += 2) {
		ULL a = i;
		inv[i] = i;
		for(int j = 1;j <= 62;j ++) {
			a = a * a;
			inv[i] *= a;
		}
	}
	int T = read();
	while(T --) {
		ULL N = readull();
		int ct = 0;
		ULL b = 1,as = facf(N);
		for(int i = 1;i < 64;i ++) {
			b <<= 1;
			ct += (N/b) & 3;
			ct &= 3;
			as *= facf(N/b);
		}
		for(int i = 1;i <= ct;i ++) as *= 2ull;
		AIput(as,'\n');
	}
	return 0;
}

Posted by asterinex on Sun, 31 Oct 2021 08:21:42 -0700