#include #define rep(i,n) for(int i=0;i<(int)(n);i++) #define rep1(i,n) for(int i=1;i<=(int)(n);i++) #define all(c) c.begin(),c.end() #define pb push_back #define fs first #define sc second #define chmin(x,y) x=min(x,y) #define chmax(x,y) x=max(x,y) using namespace std; template ostream& operator<<(ostream& o,const pair &p){ return o<<"("< ostream& operator<<(ostream& o,const vector &vc){ o<<"{"; for(const T& v:vc) o< using V = vector; template using VV = vector>; constexpr ll TEN(int n) { return (n == 0) ? 1 : 10 * TEN(n-1); } #ifdef LOCAL #define show(x) cerr << "LINE" << __LINE__ << " : " << #x << " = " << (x) << endl #else #define show(x) true #endif template struct ModInt{ using uint = unsigned int; using ll = long long; using ull = unsigned long long; constexpr static uint mod = mod_; uint v; ModInt():v(0){} ModInt(ll _v):v(normS(_v%mod+mod)){} explicit operator bool() const {return v!=0;} static uint normS(const uint &x){return (x [0 , mod-1] static ModInt make(const uint &x){ModInt m; m.v=x; return m;} ModInt operator+(const ModInt& b) const { return make(normS(v+b.v));} ModInt operator-(const ModInt& b) const { return make(normS(v+mod-b.v));} ModInt operator-() const { return make(normS(mod-v)); } ModInt operator*(const ModInt& b) const { return make((ull)v*b.v%mod);} ModInt operator/(const ModInt& b) const { return *this*b.inv();} ModInt& operator+=(const ModInt& b){ return *this=*this+b;} ModInt& operator-=(const ModInt& b){ return *this=*this-b;} ModInt& operator*=(const ModInt& b){ return *this=*this*b;} ModInt& operator/=(const ModInt& b){ return *this=*this/b;} ModInt& operator++(int){ return *this=*this+1;} ModInt& operator--(int){ return *this=*this-1;} ll extgcd(ll a,ll b,ll &x,ll &y) const{ ll p[]={a,1,0},q[]={b,0,1}; while(*q){ ll t=*p/ *q; rep(i,3) swap(p[i]-=t*q[i],q[i]); } if(p[0]<0) rep(i,3) p[i]=-p[i]; x=p[1],y=p[2]; return p[0]; } ModInt inv() const { ll x,y; extgcd(v,mod,x,y); return make(normS(x+mod)); } ModInt pow(ll p) const { ModInt a = 1; ModInt x = *this; while(p){ if(p&1) a *= x; x *= x; p >>= 1; } return a; } bool operator==(const ModInt& b) const { return v==b.v;} bool operator!=(const ModInt& b) const { return v!=b.v;} friend istream& operator>>(istream &o,ModInt& x){ ll tmp; o>>tmp; x=ModInt(tmp); return o; } friend ostream& operator<<(ostream &o,const ModInt& x){ return o<; int bsr(int x) { return 31 - __builtin_clz(x); } using D = double; const D pi = acos(-1); using Pc = complex; void fft(bool type, vector &c){ //multiply : false -> mult -> true static vector buf[30]; int N = c.size(); int s = bsr(N); assert(1<(N); rep(i,N) buf[s][i] = polar(1,i*2*pi/N); } vector a(N),b(N); copy(begin(c),end(c),begin(a)); rep1(i,s){ int W = 1<<(s-i); for(int y=0;y>1] = l-r; } } swap(a,b); } copy(begin(a),end(a),begin(c)); } template vector multiply_fft(const vector& x,const vector& y){ if(x.empty() || y.empty()) return {}; const int B = 15; const int K = 2; int S = x.size()+y.size()-1; int N = 1; while(N a[K],b[K]; rep(t,K){ a[t] = vector(N); b[t] = vector(N); rep(i,x.size()) a[t][i] = Pc( (x[i].v >> (t*B)) & ((1<> (t*B)) & ((1< z(S); vector c(N); Mint base = 1; rep(t,K+K-1){ fill_n(begin(c),N,Pc(0,0)); rep(xt,K){ int yt = t-xt; if(0<=yt && yt struct Poly{ vector v; int size() const{ return v.size();} //deg+1 Poly(){} Poly(vector _v) : v(_v){shrink();} Poly& shrink(){ while(!v.empty()&&v.back()==D(0)) v.pop_back(); return *this; } D at(int i) const{ return (i=size() && !x) return; while(i>=size()) v.push_back(D(0)); v[i]=x; shrink(); return; } D operator()(D x) const { D res = 0; int n = size(); D a = 1; rep(i,n){ res += a*v[i]; a *= x; } return res; } Poly operator+(const Poly &r) const{ int N=max(size(),r.size()); vector ret(N); rep(i,N) ret[i]=at(i)+r.at(i); return Poly(ret); } Poly operator-(const Poly &r) const{ int N=max(size(),r.size()); vector ret(N); rep(i,N) ret[i]=at(i)-r.at(i); return Poly(ret); } Poly operator-() const{ int N=size(); vector ret(N); rep(i,N) ret[i] = -at(i); return Poly(ret); } Poly operator*(const Poly &r) const{ if(size()==0||r.size()==0) return Poly(); return mul_fft(r); // FFT or NTT ? } Poly operator*(const D &r) const{ int N=size(); vector ret(N); rep(i,N) ret[i]=v[i]*r; return Poly(ret); } Poly operator/(const D &r) const{ return *this * r.inv(); } Poly operator/(const Poly &y) const{ return div_fast(y); } Poly operator%(const Poly &y) const{ return rem_fast(y); // return rem_naive(y); } Poly operator<<(const int &n) const{ // *=x^n assert(n>=0); int N=size(); vector ret(N+n); rep(i,N) ret[i+n]=v[i]; return Poly(ret); } Poly operator>>(const int &n) const{ // /=x^n assert(n>=0); int N=size(); if(N<=n) return Poly(); vector ret(N-n); rep(i,N-n) ret[i]=v[i+n]; return Poly(ret); } bool operator==(const Poly &y) const{ return v==y.v; } bool operator!=(const Poly &y) const{ return v!=y.v; } Poly& operator+=(const Poly &r) {return *this = *this+r;} Poly& operator-=(const Poly &r) {return *this = *this-r;} Poly& operator*=(const Poly &r) {return *this = *this*r;} Poly& operator*=(const D &r) {return *this = *this*r;} Poly& operator/=(const Poly &r) {return *this = *this/r;} Poly& operator/=(const D &r) {return *this = *this/r;} Poly& operator%=(const Poly &y) {return *this = *this%y;} Poly& operator<<=(const int &n) {return *this = *this<>=(const int &n) {return *this = *this>>n;} Poly diff() const { int n = size(); if(n == 0) return Poly(); V u(n-1); rep(i,n-1) u[i] = at(i+1) * (i+1); return Poly(u); } Poly intg() const { int n = size(); V u(n+1); rep(i,n) u[i+1] = at(i) / (i+1); return Poly(u); } Poly pow(long long n, int L) const { // f^n, ignoring x^L,x^{L+1},.. Poly a({1}); Poly x = *this; while(n){ if(n&1){ a *= x; a = a.strip(L); } x *= x; x = x.strip(L); n /= 2; } return a; } /* [x^0~n] exp(f) = 1 + f + f^2 / 2 + f^3 / 6 + .. f(0) should be 0 O((N+n) log n) (N = size()) NTT, -O3 - N = n = 100000 : 200 [ms] - N = n = 200000 : 400 [ms] - N = n = 500000 : 1000 [ms] */ Poly exp(int n) const { assert(at(0) == 0); Poly f({1}), g({1}); for(int i=1;i<=n;i*=2){ g = (g*2 - f*g*g).strip(i); Poly q = (this->diff()).strip(i-1); Poly w = (q + g * (f.diff() - f*q)) .strip(2*i-1); f = (f + f * (*this - w.intg()).strip(2*i)) .strip(2*i); } return f.strip(n+1); } /* [x^0~n] log(f) = log(1-(1-f)) = - (1-f) - (1-f)^2 / 2 - (1-f)^3 / 3 - ... f(0) should be 1 O(n log n) NTT, -O3 1e5 : 140 [ms] 2e5 : 296 [ms] 5e5 : 640 [ms] 1e6 : 1343 [ms] */ Poly log(int n) const { assert(at(0) == 1); auto f = strip(n+1); return (f.diff() * f.inv(n)).strip(n).intg(); } /* [x^0~n] sqrt(f) f(0) should be 1 いや平方剰余なら何でもいいと思うけど探すのがめんどくさいので +- 2通りだけど 定数項が 1 の方 O(n log n) NTT, -O3 1e5 : 234 [ms] 2e5 : 484 [ms] 5e5 : 1000 [ms] 1e6 : 2109 [ms] */ Poly sqrt(int n) const { assert(at(0) == 1); Poly f = strip(n+1); Poly g({1}); for(int i=1; i<=n; i*=2){ g = (g + f.strip(2*i)*g.inv(2*i-1)) / 2; } return g.strip(n+1); } /* [x^0~n] f^-1 = (1-(1-f))^-1 = (1-f) + (1-f)^2 + ... f * f.inv(n) = 1 + x^n * poly f(0) should be non0 O(n log n) */ Poly inv(int n) const { assert(at(0) != 0); Poly f = strip(n+1); Poly g({at(0).inv()}); for(int i=1; i<=n; i*=2){ //need to strip!! g *= (Poly({2}) - f.strip(2*i)*g).strip(2*i); } return g.strip(n+1); } Poly exp_naive(int n) const { assert(at(0) == 0); Poly res; Poly fk({1}); rep(k,n+1){ res += fk; fk *= *this; fk = fk.strip(n+1) / (k+1); } return res; } Poly log_naive(int n) const { assert(at(0) == 1); Poly res; Poly g({1}); rep1(k,n){ g *= (Poly({1}) - *this); g = g.strip(n+1); res -= g / k; } return res; } Poly mul_naive(const Poly &r) const{ int N=size(),M=r.size(); vector ret(N+M-1); rep(i,N) rep(j,M) ret[i+j]+=at(i)*r.at(j); return Poly(ret); } Poly mul_ntt(const Poly &r) const{ return Poly(multiply_ntt(v,r.v)); } Poly mul_fft(const Poly &r) const{ return Poly(multiply_fft(v,r.v)); } Poly div_fast_with_inv(const Poly &inv, int B) const { return (*this * inv)>>(B-1); } Poly div_fast(const Poly &y) const{ if(size() res = v; res.resize(min(n,size())); return Poly(res); } Poly rev(int n = -1) const { //ignore x^n ~ -> return x^(n-1) * f(1/x) vector res = v; if(n!=-1) res.resize(n); reverse(all(res)); return Poly(res); } /* f.inv_div(n) = x^n / f f should be non0 O((N+n) log n) for division */ Poly inv_div(int n) const { n++; int d = size() - 1; assert(d != -1); if(n < d) return Poly(); Poly a = rev(); Poly g({at(d).inv()}); for(int i=1; i+d<=n; i*=2){ //need to strip!! g *= (Poly({2})-a.strip(2*i)*g).strip(2*i); } return g.rev(n-d); } friend ostream& operator<<(ostream &o,const Poly& x){ if(x.size()==0) return o<<0; rep(i,x.size()) if(x.v[i]!=D(0)){ o< fact,ifact; mint Choose(int a,int b){ if(b<0 || a=0;i--) ifact[i] = ifact[i+1] * (i+1); } vector extended(int n, const vector< vector >& coeffs, const vector& terms) { vector ret(max(n + 1, terms.size())); copy(terms.begin(), terms.end(), ret.begin()); const int order = coeffs.size() - 1; const int deg = coeffs[0].size() - 1; assert((int) terms.size() >= order); for (int m = terms.size(); m <= n; ++m) { mint s = 0; for (int i = 1; i <= order; ++i) { int k = m - i; mint t = ret[k]; for (int d = 0; d <= deg; ++d) { s += t * coeffs[i][d]; t *= k; } } mint denom = 0, mpow = 1; for (int d = 0; d <= deg; ++d) { denom += mpow * coeffs[0][d]; mpow *= m; } ret[m] = -s/denom; } return ret; } vector< vector > find_recurrence_relation(vector terms, int deg, int ord = -1, bool verify=true) { if(ord != -1){ //given order int n = (deg+1)*(ord+1)+ord-1; while((int)terms.size()>n) terms.pop_back(); } const int n = terms.size(); const int B = (n + 2) / (deg + 2); // number of blocks const int C = B * (deg + 1); // number of columns const int R = n - (B - 1); // number of rows assert(B >= 2); assert(R >= C - 1); auto error = [] (int order, int deg) { fprintf(stderr, "Error: Could not find a recurrence relation " "of order <= %d and degree <= %d.\n\n", order, deg); assert(0); }; vector< vector > mat(R, vector(C)); for (int y = 0; y < R; ++y) { for (int b = 0; b < B; ++b) { mint v = terms[y+b]; for (int d = 0; d <= deg; ++d) { mat[y][b * (deg + 1) + d] = v; v *= y+b; } } } int rank = 0; for (int x = 0; x < C; ++x) { int pivot = -1; for (int y = rank; y < R; ++y) if (mat[y][x] != 0) { pivot = y; break; } if (pivot < 0) break; if (pivot != rank) swap(mat[rank], mat[pivot]); mint inv = mat[rank][x].inv(); for (int x2 = x; x2 < C; ++x2) mat[rank][x2] *= inv; for (int y = rank + 1; y < R; ++y) if (mat[y][x]) { mint c = -mat[y][x]; for (int x2 = x; x2 < C; ++x2) { mat[y][x2] += c * mat[rank][x2]; } } ++rank; } if (rank == C) error(B - 1, deg); for (int y = rank - 1; y >= 0; --y) if (mat[y][rank]) { assert(mat[y][y] == 1); mint c = -mat[y][rank]; for (int y2 = 0; y2 < y; ++y2) { mat[y2][rank] += c * mat[y2][y]; } } int order = rank / (deg + 1); vector< vector > ret(order + 1, vector(deg + 1)); ret[0][rank % (deg + 1)] = 1; for (int y = rank - 1; y >= 0; --y) { int k = order - y / (deg + 1), d = y % (deg + 1); ret[k][d] = -mat[y][rank]; } if (verify) { auto extended_terms = extended(n - 1, ret, vector(terms.begin(), terms.begin() + order)); for (int i = 0; i < (int) terms.size(); ++i) { if (terms[i] != extended_terms[i]) error(B - 1, deg); } } auto verbose = [&] { int last = verify ? n - 1 : order + R - 1; fprintf(stderr, "[ Found a recurrence relation ]\n" "- order %d\n" "- degree %d\n" "- verified up to a(%d) (number of non-trivial terms: %d)\n", order, deg, last, (last + 1) - ((deg + 2) * (order + 1) - 2) ); fprintf(stderr, "{\n"); for (int k = 0; k <= order; ++k) { fprintf(stderr, " {"); for (int d = 0; d <= deg; ++d) { if (d) fprintf(stderr, ", "); fprintf(stderr, "%d", ret[k][d].v); } fprintf(stderr, "}%s\n", k == order ? "" : ","); } fprintf(stderr, "}\n\n"); }; verbose(); return ret; } void show_extended_sequence(int n, const vector& terms, int degree, int order = -1) { auto coeffs = find_recurrence_relation(terms, degree, order); auto extended_terms = extended(n, coeffs, terms); for (int i = 0; i < (int) extended_terms.size(); ++i) { printf("%d %d\n", i, extended_terms[i].v); } puts(""); } V get_extended_sequence(int n, const vector& terms, int degree, int order = -1) { auto coeffs = find_recurrence_relation(terms, degree, order); return extended(n, coeffs, terms); } int main(){ cin.tie(0); ios::sync_with_stdio(false); //DON'T USE scanf/printf/puts !! cout << fixed << setprecision(20); InitFact(5000010); int X,Y,Z; cin >> X >> Y >> Z; int ZZ = Z; chmin(Z,100000); Poly y({1,-2,1}); Poly f({1}); while((int)f.size() < X+Y+Z+100){ f *= (y + Poly({1})); y *= y; } show(f.size()); mint ans = 0; if(X+Y+Z == 0){ cout << 1 << endl; return 0; } V vals; rep(d,20){ mint ans = 0; rep1(k,f.size()){ ans += Choose(X+k-1,k-1) * Choose(Y+k-1,k-1) * Choose(Z+d+k-1,k-1) * f.at(k-1); } vals.pb(ans); } auto v = get_extended_sequence(ZZ-Z,vals,2,3); cout << v[ZZ-Z] << endl; }