import sys MOD = 10**9 + 9 def main(): N = int(sys.stdin.readline().strip()) if N == 0: return # Precompute factorial and inverse factorial modulo MOD fact = [1] * (N + 1) for i in range(1, N + 1): fact[i] = fact[i-1] * i % MOD # Compute g(x) = sum_{a=1 to N} (a+1)^2 x^a g = [0] * (N+1) for a in range(1, N+1): if a > N: continue g[a] = (a + 1) ** 2 % MOD # Now compute cos(g) + sin(g) - 1 mod x^(N+1) # Initialize C (cos) and S (sin) C = [0] * (N+1) S = [0] * (N+1) C[0] = 1 # cos(0) = 1 # S[0] = 0 # sin(0) = 0 # Precompute the derivative g'(x) gp = [0] * (N) for a in range(0, N): gp[a] = ( (a+1) * ( (a+2) ** 2 ) ) % MOD # Compute C and S using dynamic programming based on their differential equations for k in range(1, N+1): # Compute C[k] using C' = -g' * S # The coefficient of x^{k-1} in C' is given by the convolution of gp and S up to x^{k-1} # C' = -gp * S conv = 0 for i in range(max(0, k-1 - (N-1)), k): if i >= len(gp): continue j = (k-1) - i if j < 0 or j > N: continue conv += gp[i] * S[j] conv %= MOD Ck = (-conv) % MOD Ck = Ck * pow(k, MOD-2, MOD) % MOD C[k] = Ck # Compute S[k] using S' = gp * C conv = 0 for i in range(max(0, k-1 - (N-1)), k): if i >= len(gp): continue j = (k-1) - i if j < 0 or j > N: continue conv += gp[i] * C[j] conv %= MOD Sk = conv * pow(k, MOD-2, MOD) % MOD S[k] = Sk # Now compute S_total = C + S - 1 S_total = [ (C[i] + S[i]) % MOD for i in range(N+1) ] S_total[0] = (S_total[0] - 1) % MOD if N >= 0: S_total[0] %= MOD for k in range(1, N+1): res = S_total[k] * fact[N] % MOD print(res) if __name__ == "__main__": main()