結果
問題 | No.502 階乗を計算するだけ |
ユーザー |
|
提出日時 | 2021-02-28 19:38:39 |
言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
結果 |
AC
|
実行時間 | 112 ms / 1,000 ms |
コード長 | 4,488 bytes |
コンパイル時間 | 2,290 ms |
コンパイル使用メモリ | 205,508 KB |
最終ジャッジ日時 | 2025-01-19 08:16:19 |
ジャッジサーバーID (参考情報) |
judge3 / judge1 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
other | AC * 52 |
ソースコード
#include<bits/stdc++.h>using namespace std;using i64 = int64_t;using u64 = uint64_t;using base = complex<double>;void fft(vector<base> &a, bool inv){int n = a.size(), j = 0;for(int i=1; i<n; i++){int bit = (n >> 1);while(j >= bit){j -= bit;bit >>= 1;}j += bit;if(i < j) swap(a[i], a[j]);}double ang = 2 * acos(-1) / n * (inv ? -1 : 1);vector<base> roots(n/2);for(int i=0; i<n/2; i++)roots[i] = base(cos(ang * i), sin(ang * i));for(int i=2; i<=n; i<<=1){int step = n / i;for(int j=0; j<n; j+=i){for(int k=0; k<i/2; k++){base u = a[j+k], v = a[j+k+i/2] * roots[step * k];a[j+k] = u+v;a[j+k+i/2] = u-v;}}}if(inv) for(int i=0; i<n; i++) a[i] /= n;}vector<i64> multiply(vector<i64> &v, vector<i64> &w, i64 mod){int n = 2; while(n < v.size() + w.size()) n <<= 1;vector<base> v1(n), v2(n), r1(n), r2(n);for(int i=0; i<v.size(); i++)v1[i] = base(v[i] >> 15, v[i] & 32767);for(int i=0; i<w.size(); i++)v2[i] = base(w[i] >> 15, w[i] & 32767);fft(v1, 0);fft(v2, 0);for(int i=0; i<n; i++){int j = (i ? (n - i) : i);base ans1 = (v1[i] + conj(v1[j])) * base(0.5, 0);base ans2 = (v1[i] - conj(v1[j])) * base(0, -0.5);base ans3 = (v2[i] + conj(v2[j])) * base(0.5, 0);base ans4 = (v2[i] - conj(v2[j])) * base(0, -0.5);r1[i] = (ans1 * ans3) + (ans1 * ans4) * base(0, 1);r2[i] = (ans2 * ans3) + (ans2 * ans4) * base(0, 1);}fft(r1, 1);fft(r2, 1);vector<i64> ret(n);for(int i=0; i<n; i++){i64 av = (i64)round(r1[i].real());i64 bv = (i64)round(r1[i].imag()) + (i64)round(r2[i].real());i64 cv = (i64)round(r2[i].imag());av %= mod, bv %= mod, cv %= mod;ret[i] = (av << 30) + (bv << 15) + cv;ret[i] %= mod;ret[i] += mod;ret[i] %= mod;}return ret;}/*Input h: vector of length d+1: h(0), h(1), h(2), ..., h(d) where h is d-degree polynomialInput P: prime, where h degree polynomial is represented as moduloOutput: vector of length 4d+2: h(0), h(1), h(2), ..., h(4d+1)Lagrange interpolation is given as follows:h(x) = sum_(i=0 to d)[ h(i) prod(j=0 to d, i != j){(x-j)/(i-j)} ]= {prod(j=0 to d)(x-j)} *sum(i=0 to d){h(i)/(i!(d-i)!(-1)^(d-i)) * 1/(x-i) }Letting f(i) = h(i)/(i!(d-i)!(-1)^(d-i)), g(x) = prod(j=0 to d)(x-j)h(x) = g(x) * sum(i=0 to d){f(i) * 1/(x-i)} which is easily calculated by convolution.f(0)...f(d) and 1/1, 1/2, ..., 1/(4d+1) will be fed to multiply logic*/vector<i64> lagrange(vector<i64> h, i64 P) {int d = (int)h.size()-1; assert(d >= 0);auto mul = [P](i64 a, i64 b){return a*b%P;};auto ipow = [&mul](i64 a, i64 b) {i64 res = 1;while(b) {if(b&1) res = mul(res, a);a = mul(a, a);b >>= 1;}return res;};auto modInv = [&ipow, P](i64 a) {return ipow(a, P-2);};vector<i64> fact(4*d+2), invfact(4*d+2);fact[0] = 1;for(int i=1; i<=4*d+1; ++i)fact[i] = mul(i,fact[i-1]);invfact[4*d+1] = modInv(fact[4*d+1]);for(int i=4*d; i>=0; --i)invfact[i] = mul(invfact[i+1], i+1);vector<i64> f(d+1);for(int i=0; i<=d; ++i) {f[i] = h[i];f[i] = mul(invfact[i], f[i]);f[i] = mul(invfact[d-i], f[i]);if((d-i)%2==1) f[i] = P-f[i];if(f[i] == P) f[i] = 0;}vector<i64> inv(4*d+2);for(int i=1; i<4*d+2; ++i)inv[i] = mul(fact[i-1], invfact[i]);vector<i64> g(4*d+2);g[d+1] = 1;for(int j=0; j<=d; ++j) g[d+1] = mul(g[d+1], d+1-j);for(int i=d+2; i<4*d+2; ++i)g[i] = mul(g[i-1], mul(i, inv[i-d-1]));vector<i64> conv = multiply(f, inv, P);vector<i64> ret(4*d+2);for(int i=0; i<=d; ++i) ret[i] = h[i];for(int i=d+1; i<4*d+2; ++i) ret[i] = mul(g[i], conv[i]);return ret;}vector<i64> squarepoly(vector<i64> poly, i64 P) {vector<i64> ss = lagrange(poly, P);vector<i64> ret(ss.size()/2);auto mul = [P](i64 a, i64 b){return a*b%P;};for(int i=0; i<(int)ss.size()/2; ++i)ret[i] = mul(ss[2*i], ss[2*i+1]);return ret;}i64 factorial(i64 N,i64 P){i64 d = 1;vector<i64> fact_part = {1%P, 2%P};while(N > d*(d+1)) {fact_part = squarepoly(fact_part, P);d *= 2;}auto mul = [P](i64 a, i64 b) {return a*b%P;};i64 ans = 1, bucket = N/d;for(int i=0; i<bucket; ++i) ans = mul(ans, fact_part[i]);for(i64 i=bucket*d+1; i<=N; ++i) ans = mul(ans, i);return ans;}int main(){cin.tie(0);ios::sync_with_stdio(false);i64 N, P;cin >> N;P = 1e9 + 7;if(N>=P){cout << 0 << '\n';return 0;}cout << factorial(N,P) << '\n';return 0;}