Language: C++

#include "bits/stdc++.h"
#include <gmp.h>
using namespace std;
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i))
#define reu(i,l,u) for(int (i)=(int)(l);(i)<(int)(u);++(i))
static const int INF = 0x3f3f3f3f; static const long long INFL = 0x3f3f3f3f3f3f3f3fLL;
typedef vector<int> vi; typedef pair<int, int> pii; typedef vector<pair<int, int> > vpii; typedef long long ll;
template<typename T, typename U> static void amin(T &x, U y) { if(y < x) x = y; }
template<typename T, typename U> static void amax(T &x, U y) { if(x < y) x = y; }


template<int MOD>
struct ModInt {
	static const int Mod = MOD;
	unsigned x;
	ModInt() : x(0) {}
	ModInt(signed sig) { int sigt = sig % MOD; if(sigt < 0) sigt += MOD; x = sigt; }
	ModInt(unsigned unsig) { x = unsig % MOD; }
	ModInt(signed long long sig) { int sigt = sig % MOD; if(sigt < 0) sigt += MOD; x = sigt; }
	int get() const { return (int)x; }

	ModInt &operator+=(ModInt that) { if((x += that.x) >= MOD) x -= MOD; return *this; }
	ModInt &operator-=(ModInt that) { if((x += MOD - that.x) >= MOD) x -= MOD; return *this; }
	ModInt &operator*=(ModInt that) { x = (unsigned long long)x * that.x % MOD; return *this; }
	ModInt &operator/=(ModInt that) { return *this *= that.inverse(); }

	ModInt operator+(ModInt that) const { return ModInt(*this) += that; }
	ModInt operator-(ModInt that) const { return ModInt(*this) -= that; }
	ModInt operator*(ModInt that) const { return ModInt(*this) *= that; }
	ModInt operator/(ModInt that) const { return ModInt(*this) /= that; }

	ModInt inverse() const {
		signed a = x, b = MOD, u = 1, v = 0;
		while(b) {
			signed t = a / b;
			a -= t * b; std::swap(a, b);
			u -= t * v; std::swap(u, v);
		}
		if(u < 0) u += Mod;
		ModInt res; res.x = (unsigned)u;
		return res;
	}
};
template<int MOD> ModInt<MOD> operator^(ModInt<MOD> a, unsigned long long k) {
	ModInt<MOD> r = 1;
	while(k) {
		if(k & 1) r *= a;
		a *= a;
		k >>= 1;
	}
	return r;
}
typedef ModInt<1000000007> mint;

#if defined(_WIN32)
extern "C" {
	void* __stdcall LoadLibraryA(const char *name);
	void* __stdcall GetProcAddress(void *handle, const char *name);
}
struct DynamicLibrary {
	DynamicLibrary(const char *name) { _handle = LoadLibraryA(name); }
	void *get(const char *name) const { void *p = GetProcAddress(_handle, name); if(!p) std::cerr << "can't find: " << name << std::endl, std::abort(); return p; }
	void *_handle;
} gmp("mpir.dll");
#else
extern "C" {
	void* __libc_dlopen_mode(const char *x, int y);
	void* __libc_dlsym(void *x, const char *y);
}
struct DynamicLibrary {
	DynamicLibrary(const char *name) { _handle = __libc_dlopen_mode(name, 2); }
	void *get(const char *name) const { void *p = __libc_dlsym(_handle, name); if(!p) std::cerr << "can't find: " << name << std::endl, std::abort(); return p; }
	void *_handle;
} gmp("/usr/lib/i386-linux-gnu/libgmp.so.10");
#endif

#define EXPAND_MACRO_TO_STR_2(a) #a
#define EXPAND_MACRO_TO_STR(a) EXPAND_MACRO_TO_STR_2(a)
#define DECL_GMP(name) const auto my_##name = (decltype(::name)*)gmp.get(EXPAND_MACRO_TO_STR(name))
DECL_GMP(mpn_mul);
DECL_GMP(mpn_sqr);


inline int getBitSize(unsigned x) {
#ifndef __GNUC__
	unsigned long res;
	_BitScanReverse(&res, x);
	return static_cast<int>(res) + 1;
#else
	return 32 - __builtin_clz(x);
#endif
}

enum { WordBits = sizeof(mp_limb_t) * 8 };
typedef unsigned Num;
enum { NumBits = sizeof(Num) * 8 };

typedef int Size;
typedef const mint *SrcPoly;
typedef mint *DstPoly, SrcDstPoly;
typedef const mp_limb_t *SrcPtr;
typedef mp_limb_t *DstPtr, *SrcDstPtr;

void zeroPolynomial(DstPoly resp, Size resn) {
	for(Size i = 0; i < resn; ++ i)
		resp[i] = mint();
}

void addPolynomial(DstPoly resp, SrcPoly xp, SrcPoly yp, Size xyn) {
	for(Size i = 0; i < xyn; ++ i)
		resp[i] = xp[i] + yp[i];
}
void subtractPolynomial(DstPoly resp, SrcPoly xp, SrcPoly yp, Size xyn) {
	for(Size i = 0; i < xyn; ++ i)
		resp[i] = xp[i] - yp[i];
}

void multiplyConstantPolynomial(DstPoly resp, SrcPoly xp, Size xn, mint y) {
	for(Size i = 0; i < xn; ++ i)
		resp[i] = xp[i] * y;
}

void packToBits(DstPtr resp, int resn, SrcPoly xp, Size xn, int width, int stride) {
	if(width > WordBits) {
		std::cerr << "unimplemented" << std::endl;
		std::abort();
	}
	for(Size i = 0; i < resn; ++ i)
		resp[i] = 0;
	Size pos = 0;
	int sa = stride / WordBits, sb = stride % WordBits;
	int offset = 0;
	for(Size i = 0; i < xn; ++ i) {
		assert(pos < resn);
		mp_limb_t x = static_cast<mp_limb_t>(xp[i].get());
		int right = offset + width;
		resp[pos] |= x << offset;
		if(right > WordBits) {
			assert(pos + 1 < resn);
			resp[pos + 1] |= x >> (WordBits - offset);
		}
		pos += sa;
		if((offset += sb) >= WordBits)
			offset -= WordBits, ++ pos;
	}
}

void unpackFromBitsMod(DstPoly resp, Size resn, SrcPtr xp, Size xn, int stride) {
	assert(stride > 0);
	int K = (stride - 1) / NumBits + 1, rems = stride - (K - 1) * NumBits;
	Num remmask = ((Num)2 << (rems - 1)) - 1;
	const mint Base = mint(~0U) + 1;
	std::unique_ptr<Num[]> buf(new Num[K]);
	Num *nums = buf.get();
	Size pos = 0, offset = 0;
	for(Size i = 0; i < resn; ++ i) {
		for(int k = 0; k < K; ++ k) {
			int smallwidth = k < K - 1 ? NumBits : rems;
			if(pos >= xn) {
				nums[k] = 0;
				offset += smallwidth;
				continue;
			}
			Num num = static_cast<Num>(xp[pos] >> offset);
			offset += smallwidth;
			if(pos + 1 < xn && offset > WordBits)
				num |= static_cast<Num>(xp[pos + 1] << (WordBits - (offset - smallwidth)));
			nums[k] = num;
			if(offset >= WordBits)
				offset -= WordBits, ++ pos;
		}
		nums[K - 1] &= remmask;
		mint cur = mint(nums[K - 1]);
		for(int k = K - 2; k >= 0; -- k) {
			cur *= Base;
			cur += nums[k];
		}
		resp[i] = cur;
	}
}

void multiplyPolynomial(DstPoly resp, Size resn, SrcPoly xp, Size xn, SrcPoly yp, Size yn) {
	if(xn > resn) xn = resn;
	if(yn > resn) yn = resn;
	if(resn > xn + yn - 1) {
		zeroPolynomial(resp + (xn + yn - 1), resn - (xn + yn - 1));
		resn = xn + yn - 1;
	}
	if(xn < yn) std::swap(xp, yp), std::swap(xn, yn);
	if(yn == 0) {
		zeroPolynomial(resp, resn);
		return;
	}
	if(yn == 1) {
		multiplyConstantPolynomial(resp, xp, xn, yp[0]);
		zeroPolynomial(resp + xn, resn - xn);
		return;
	}
	int width = getBitSize(mint::Mod - 1);
	int stride = width * 2 + getBitSize(yn - 1);
	if(xp == yp && xn == yn) {
		Size tmpxn = static_cast<Size>(((long long)xn * stride - 1) / WordBits + 1);
		Size tmpzn = tmpxn + tmpxn;
		std::unique_ptr<mp_limb_t[]> tmpbuf(new mp_limb_t[tmpxn + tmpzn]);
		mp_limb_t *tmpx = tmpbuf.get(), *tmpz = tmpx + tmpxn;
		packToBits(tmpx, tmpxn, xp, xn, width, stride);
		my_mpn_sqr(tmpz, tmpx, tmpxn);
		unpackFromBitsMod(resp, resn, tmpz, tmpzn, stride);
	} else {
		Size tmpxn = static_cast<Size>(((long long)xn * stride - 1) / WordBits + 1);
		Size tmpyn = static_cast<Size>(((long long)yn * stride - 1) / WordBits + 1);
		Size tmpzn = tmpxn + tmpyn;
		std::unique_ptr<mp_limb_t[]> tmpbuf(new mp_limb_t[tmpxn + tmpyn + tmpzn]);
		mp_limb_t *tmpx = tmpbuf.get(), *tmpy = tmpx + tmpxn, *tmpz = tmpy + tmpyn;
		packToBits(tmpx, tmpxn, xp, xn, width, stride);
		packToBits(tmpy, tmpyn, yp, yn, width, stride);
		my_mpn_mul(tmpz, tmpx, tmpxn, tmpy, tmpyn);
		unpackFromBitsMod(resp, resn, tmpz, tmpzn, stride);
	}
}



vector<mint> fact, factinv;
void nCr_computeFactinv(int N) {
	N = min(N, mint::Mod - 1);
	fact.resize(N + 1); factinv.resize(N + 1);
	fact[0] = 1;
	rer(i, 1, N) fact[i] = fact[i - 1] * i;
	factinv[N] = fact[N].inverse();
	for(int i = N; i >= 1; i --) factinv[i - 1] = factinv[i] * i;
}
mint nCr(int n, int r) {
	if(n >= mint::Mod)
		return nCr(n % mint::Mod, r % mint::Mod) * nCr(n / mint::Mod, r / mint::Mod);
	return r > n ? 0 : fact[n] * factinv[n - r] * factinv[r];
}

int main() {
	int N; int M; int K;
	while(~scanf("%d%d%d", &N, &M, &K)) {
		nCr_computeFactinv(N + M + K);
		vector<mint> p(M + 1), q(K + 1);
		rep(i, M + 1) p[i] = factinv[i];
		rep(i, K + 1) q[i] = factinv[i];
		vector<mint> pq(M + K + 1);
		multiplyPolynomial(pq.data(), pq.size(), p.data(), p.size(), q.data(), q.size());
		mint ans;
		rer(mk, 0, M + K) {
			mint sum = pq[mk];
			sum *= fact[N - 1 + mk];
			sum *= factinv[N - 1];
			ans = ans * 3 + sum;
		}
		printf("%d\n", ans.get());
	}
	return 0;
}


Compiler output:
/usr/bin/ld: /tmp/cc1PXXUO.o: in function `DynamicLibrary::DynamicLibrary(char const*)':
a.cc:(.text._ZN14DynamicLibraryC2EPKc[_ZN14DynamicLibraryC5EPKc]+0x1d): undefined reference to `__libc_dlopen_mode'
/usr/bin/ld: /tmp/cc1PXXUO.o: in function `DynamicLibrary::get(char const*) const':
a.cc:(.text._ZNK14DynamicLibrary3getEPKc[_ZNK14DynamicLibrary3getEPKc]+0x22): undefined reference to `__libc_dlsym'
collect2: error: ld returned 1 exit status