Given are a prime number p and a sequence of p integers a_0, \ldots, a_{p-1} consisting of zeros and ones.
Find a polynomial of degree at most p-1, f(x) = b_{p-1} x^{p-1} + b_{p-2} x^{p-2} + \ldots + b_0, satisfying the following conditions:
 - For each i (0 \leq i \leq p-1), b_i is an integer such that 0 \leq b_i \leq p-1.
 - For each i (0 \leq i \leq p-1), f(i) \equiv a_i \pmod p.

-----Constraints-----
 - 2 \leq p \leq 2999
 - p is a prime number.
 - 0 \leq a_i \leq 1

-----Input-----
Input is given from Standard Input in the following format:
p
a_0 a_1 \ldots a_{p-1}

-----Output-----
Print b_0, b_1, \ldots, b_{p-1} of a polynomial f(x) satisfying the conditions, in this order, with spaces in between.
It can be proved that a solution always exists. If multiple solutions exist, any of them will be accepted.

-----Sample Input-----
2
1 0

-----Sample Output-----
1 1

f(x) = x + 1 satisfies the conditions, as follows:
 - f(0) = 0 + 1 = 1 \equiv 1 \pmod 2
 - f(1) = 1 + 1 = 2 \equiv 0 \pmod 2

Difficulty: interview

Starter code:

Solutions:
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10 ** 7)

import numpy as np
import itertools

P = int(input())
A = [int(x) for x in input().split()]

# k,i -> i^k
power = np.ones((P,P),dtype = np.int64)
for k in range(1,P):
    power[k] = power[k-1] * np.arange(P,dtype=np.int64) % P

f = np.zeros(P,dtype=np.int64)
for i,a in enumerate(A):
    if a == 1:
        f[0] += 1
        f -= power[:,i][::-1]
f %= P

f %= P
print((*f))



import numpy as np

p = int(input())
aaa = list(map(int, input().split()))

bins = [1]
for i in range(1, p):
    bins.append(bins[-1] * (p - i) * pow(i, p - 2, p) % p)
bins = np.array(bins, dtype=np.int32)

pows_base = -np.arange(p, dtype=np.int32)
pows = np.zeros((p, p), dtype=np.int32)
pows[:, p - 1] = 1
for i in range(p - 2, -1, -1):
    pows[:, i] = pows[:, i + 1] * pows_base % p

coefs = np.zeros(p, dtype=np.int32)

for i, a in enumerate(aaa):
    if a == 0:
        continue

    coefs[0] += 1
    coefs = (coefs - bins * pows[i]) % p

print((*coefs))


def solve():
    p = int(input())
    *A, = map(int, input().split())

    L = p-1
    fact = [1]*(L+1)
    rfact = [1]*(L+1)
    r = 1
    for i in range(1, L+1):
        fact[i] = r = r * i % p
    rfact[L] = r = pow(r, p-2, p)
    for i in range(L, 0, -1):
        rfact[i-1] = r = r * i % p
    V = [(-fact[p-1] * rfact[k] * rfact[p-1-k]) % p for k in range(p)]

    R = [0]*p
    I = [i for i, a in enumerate(A) if a]
    C = [a for a in A]

    for k in range(p):
        R[k] = V[k] * sum(C) % p

        for i in I:
            C[i] = C[i] * -i % p
    R[-1] = (R[-1] + len(I)) % p
    R.reverse()
    print(*R)
solve()

p = int(input())
b_list = [0 for _ in range(p-1)]
a_list = list(map(int,input().split()))
def frac(n):
    if a_list[n] == 1:
        t = 1
        for j in range(p-1):
            b_list[j] -= t
            t = (t*n)%p
if a_list[0] == 1:
    b_list[0] -= 1
for i in range(p-1):
    frac(i+1)
b_list.reverse()
print(a_list[0],end = ' ')
for j in b_list:
    print(j % p,end = ' ')


# -*- coding: utf-8 -*-
from sys import stdin
import numpy as np
# import sys
# sys.setrecursionlimit(10**4)
# import time

def _li(): return list(map(int, stdin.readline().split()))
def _li_(): return list([int(x)-1 for x in stdin.readline().split()])
def _lf(): return list(map(float, stdin.readline().split()))
def _ls(): return stdin.readline().split()
def _i(): return int(stdin.readline())
def _f(): return float(stdin.readline())
def _s(): return stdin.readline()[:-1]


p = _i()
a_list = np.array(_li())

factorial_invs = {0: 1, 1: 1}
for i in range(2, p):
    factorial_invs[i] = (factorial_invs[i-1] * pow(i, p-2, p)) % p
np_fact_invs = np.tile(np.array(list(factorial_invs.values()))
                       .reshape(1, -1), (p, 1))

pow_table = np.ones((p, p), dtype=int)
pow_table[1, :] = np.arange(p)
for i in range(2, p):
    pow_table[i] = (pow_table[i-1] * pow_table[1]) % p
pow_table = pow_table.T

fact_pminus1 = pow(factorial_invs[p-1], p-2, p)

coefs = np.zeros(p, dtype=int)
# for i, a in enumerate(a_list):
#     if a == 1:
#         for k in range(1, p):
#             coefs[k] = coefs[k] - (-1)**(k % 2) * pow(i, p-1-k, p)\
#                     * fact_pminus1\
#                     * factorial_invs[k]\
#                     * factorial_invs[p-1-k]
#             coefs[k] %= p

coefs[0] = a_list[0]
sign_array = np.tile(((-1) ** (np.arange(1, p) % 2)).reshape(1, -1), (p, 1))
coefs[1:] = - ((((sign_array[a_list == 1, :]
                  * pow_table[a_list == 1, :-1][:, ::-1]
                  * fact_pminus1) % p
                 * np_fact_invs[a_list == 1, 1:]) % p
                * np_fact_invs[a_list == 1, :-1][:, ::-1]
                ) % p).sum(axis=0)
coefs %= p

print((*coefs))


import sys
input = sys.stdin.readline
sys.setrecursionlimit(10 ** 7)

import numpy as np
from numpy.fft import rfft, irfft
import itertools

P = int(input())
A = [int(x) for x in input().split()]

fft_len = 6144
def fft_convolve(A,B):
    C = irfft(rfft(A,fft_len) * rfft(B,fft_len))
    C = (C + .5).astype(np.int64)
    C[:P] += C[P:P+P]
    return C[:P] % P
  
one = np.array([1], dtype=np.int64)
arr = [np.array([(-x)%P,1],dtype=np.int64) for x in range(P)]
L = [one] + list(itertools.accumulate(arr[:-1],lambda a,b: np.convolve(a,b) % P))
R = ([one] + list(itertools.accumulate(arr[1:][::-1],lambda a,b: np.convolve(a,b) % P)))[::-1]
X = [fft_convolve(A,B) % P for A,B in zip(L,R)]

f = np.zeros(P,dtype=np.int64)
for a,x in zip(A,X):
    if a == 1:
        f -= x

f %= P
print((*f))


def solve():
    p = int(input())
    *A, = map(int, input().split())

    L = p-1
    fact = [1]*(L+1)
    rfact = [1]*(L+1)
    r = 1
    for i in range(1, L+1):
        fact[i] = r = r * i % p
    rfact[L] = r = pow(r, p-2, p)
    for i in range(L, 0, -1):
        rfact[i-1] = r = r * i % p
    V = [(-fact[p-1] * rfact[k] * rfact[p-1-k]) % p for k in range(p)]

    B = [-i for i in range(p)]
    R = [0]*p
    I = [i for i, a in enumerate(A) if a]
    C = [a for a in A]

    for k in range(p):
        r = 0
        R[k] = V[k] * sum(C) % p

        for i in I:
            C[i] = C[i] * B[i] % p
    R[-1] += len(I)
    for k in range(p):
        R[k] %= p
    R.reverse()
    print(*R)
solve()

import os
import sys

import numpy as np

if os.getenv("LOCAL"):
    sys.stdin = open("_in.txt", "r")

sys.setrecursionlimit(10 ** 9)
INF = float("inf")
IINF = 10 ** 18


# MOD = 10 ** 9 + 7


# MOD = 998244353


def get_factorials(max, mod=None):
    """
    階乗 0!, 1!, 2!, ..., max!
    :param int max:
    :param int mod:
    :return:
    """
    ret = [1]
    n = 1
    if mod:
        for i in range(1, max + 1):
            n *= i
            n %= mod
            ret.append(n)
    else:
        for i in range(1, max + 1):
            n *= i
            ret.append(n)
    return ret


def mod_invs(max, mod):
    """
    逆元のリスト 0 から max まで
    :param int max:
    :param int mod:
    """
    invs = [1] * (max + 1)
    for x in range(2, max + 1):
        invs[x] = (-(mod // x) * invs[mod % x]) % mod
    return invs


def factorial_invs(max, mod):
    """
    階乗 0!, 1!, 2!, ..., max! の逆元
    :param int max:
    :param int mod:
    """
    ret = []
    r = 1
    for inv in mod_invs(max, mod):
        r = r * inv % mod
        ret.append(r)
    return ret


class Combination:
    def __init__(self, max, mod):
        """
        :param int max:
        :param int mod: 3 以上の素数であること
        """
        self._factorials = get_factorials(max, mod)
        self._finvs = factorial_invs(max, mod)
        self._mod = mod

    def ncr(self, n, r):
        """
        :param int n:
        :param int r:
        :rtype: int
        """
        if n < r:
            return 0
        return (
                self._factorials[n]
                * self._finvs[r]
                % self._mod
                * self._finvs[n - r]
                % self._mod
        )


P = int(sys.stdin.readline())
A = list(map(int, sys.stdin.readline().split()))


# A = np.array(A, dtype=int)
# vpow = np.vectorize(lambda a, b: pow(a, b, P))
# mat = []
# for i in range(P):
#     mat.append(vpow(i, np.arange(P)))
# print(np.array(mat).sum(axis=1) % P)
# print(np.array(mat).sum(axis=0) % P)
# mat = np.hstack((mat, [[a] for a in A]))
#
# print(mat)
#


def power_table():
    ret = np.ones((P, P), dtype=int)
    row = -np.arange(P, dtype=int)
    for i in range(1, P):
        ret[i] = ret[i - 1] * row % P
    return ret


# 解説AC
# f(x) = 1−(x - j)^(P - 1) は、x == j のとき1、それ以外のとき0
# A[j] == 1 である j について上記の式を足し合わせる
comb = Combination(max=P - 1, mod=P)
ncr = np.array([comb.ncr(P - 1, i) for i in range(P)])
# power[i][j]: (-j)^i
power = power_table()

B = np.zeros(P, dtype=int)
for j, a in enumerate(A):
    if a == 0:
        continue
    B -= power[:, j] * ncr % P
    B[-1] += 1
    B %= P
print((*B[::-1]))


import numpy as np
mod = 998244353

def poly_mul(f, g):
    # 参考: https://judge.yosupo.jp/submission/2380
    Lf = len(f); Lg = len(g); L = Lf + Lg - 1
    if Lf <= 16 or Lg <= 16:
        if Lf == 0 or Lg == 0:
            return np.zeros((0,), dtype=np.int64)
        return (np.convolve(f.astype(np.uint64), g.astype(np.uint64)) % mod).astype(np.int64)
    fft = np.fft.rfft; ifft = np.fft.irfft
    fft_len = 1 << L.bit_length()
    fl = f & (1 << 15) - 1; fh = f >> 15
    gl = g & (1 << 15) - 1; gh = g >> 15
    x = (ifft(fft(fl, fft_len) * fft(gl, fft_len))[:L] + 0.5).astype(np.int64) % mod
    y = (ifft(fft(fl+fh, fft_len) * fft(gl+gh, fft_len))[:L] + 0.5).astype(np.int64) % mod
    z = (ifft(fft(fh, fft_len) * fft(gh, fft_len))[:L] + 0.5).astype(np.int64) % mod
    return (x + ((y - x - z) << 15) + (z << 30)) % mod

def poly_inv(fps, n=None):
    assert fps[0] != 0
    if n is None:
        n = len(fps)
    res = np.zeros(1<<(n-1).bit_length(), dtype=np.int64)
    res[0] = pow(int(fps[0]), mod-2, mod)
    i = 1
    while i < n:
        i <<= 1
        res[:i] = ((res[:i]<<1) - poly_mul(poly_mul(res[:i>>1], res[:i>>1]), fps[:i])[:i]) % mod
    return res[:n]

def poly_div(fps1, fps2):
    n1, n2 = len(fps1), len(fps2)
    if n1 < n2:
        return np.zeros((0,), dtype=np.int64)
    n = n1 - n2 + 1
    res = poly_mul(fps1[-1:-n-1:-1], poly_inv(fps2[::-1], n))[n-1::-1]
    return res

def poly_mod(fps1, fps2):
    n1, n2 = len(fps1), len(fps2)
    if n1 < n2:
        return fps1
    res = fps1[:n2-1] - poly_mul(poly_div(fps1, fps2), fps2)[:n2-1]
    return res % mod

def multipoint_evaluation(fps, xs):
    threshold = 8
    n_xs = len(xs)
    bit = (n_xs-1).bit_length()
    if bit <= threshold:
        res = np.zeros_like(xs)
        xs_cumprod = np.ones_like(xs)
        for coef in fps:
            res += xs_cumprod * coef
            xs_cumprod *= xs
            xs_cumprod %= mod
        return res
    k = 1<<bit
    fpss = np.zeros((bit+1, k+1), dtype=fps.dtype)
    fpss[0, :n_xs] = -xs % mod
    fpss[1, :k:2] = fpss[0, :k:2] * fpss[0, 1::2] % mod
    fpss[1, 1::2] = (fpss[0, :k:2] + fpss[0, 1::2]) % mod
    for i in range(1, bit):
        step = 2<<i
        half = step>>1
        for j in range(0, k, step):
            f1 = fpss[i, j:j+half+1].copy()
            f2 = fpss[i, j+half:j+step+1].copy()
            f1[-1] = f2[-1] = 1
            f = poly_mul(f1, f2)
            fpss[i+1, j:j+step] = f[:-1]
    f = poly_mod(fps, f)
    fpss[-1, :len(f)] = f
    fpss[-1, len(f):] = 0
    for i in range(bit-1, threshold-1, -1):
        step = 2<<i
        half = step>>1
        for j in range(0, k, step):
            f = fpss[i+1, j:j+step]
            f1 = fpss[i, j:j+half+1].copy()
            f2 = fpss[i, j+half:j+step+1].copy()
            f1[-1] = f2[-1] = 1
            fpss[i, j:j+half] = poly_mod(f, f1)
            fpss[i, j+half:j+step] = poly_mod(f, f2)
    xs = (-fpss[0, :k] % mod).reshape(-1, 1<<threshold)
    xs_cumprod = np.ones_like(xs)
    res = np.zeros_like(xs)
    for i in range(1<<threshold):
        res += fpss[threshold, i:k:1<<threshold, None] * xs_cumprod % mod
        xs_cumprod *= xs
        xs_cumprod %= mod
    return res.reshape(-1)[:n_xs] % mod

def poly_differential(fps):
    return fps[1:] * np.arange(1, len(fps)) % mod

def lagrange_interpolation(X, Y, mod):
    # old
    n = len(X)
    g = [0]*(n+1)
    g[0] = 1
    for i, x in enumerate(X):
        for j in range(i, -1, -1):
            g[j+1] += g[j] * (-x) % mod
    res = [0]*n
    for x, y in zip(X, Y):
        f = g[:]
        denom = 0
        v = 1
        pow_x = [1]  # x の idx 乗
        for _ in range(n-1):
            v = v * x % mod
            pow_x.append(v)
        pow_x.reverse()  # n-1 乗 ~ 0 乗
        for i, po in enumerate(pow_x):
            f_i = f[i]
            f[i+1] += f_i * x % mod  # f = g / (x - x_i) を組立除法で求める
            denom = (denom + f_i * po) % mod
        denom_inv = pow(denom, mod-2, mod)
        for i, f_i in enumerate(f[:n]):
            res[i] += (f_i * y * denom_inv)# % mod  # mod が大きいと 64bit に収まらなくなるのでひとつずつ mod 取った方がいいか？
    return [v % mod for v in res]

def polynomial_interpolation(xs, ys):
    # 参考: https://rsk0315.hatenablog.com/entry/2020/04/05/203210
    assert len(xs) == len(ys)
    threshold = 8
    as_strided = np.lib.stride_tricks.as_strided
    n = len(xs)
    if n==1:
        return ys.copy()
    bit = (n-1).bit_length()
    if bit <= threshold:
        res = lagrange_interpolation(xs.tolist(), ys.tolist(), mod)
        return np.array(res[::-1], dtype=np.int64)
    k = 1<<bit
    fpss = np.zeros((bit+1, n+1), dtype=np.int64)
    fpss[0, :n] = -xs % mod
    for i in range(bit):
        step = 2 << i
        half = step >> 1
        for j in range(0, n, step):
            if j+half >= n:
                fpss[i+1, j:n] = fpss[i, j:n]
                continue
            f1 = fpss[i, j:j+half+1].copy()
            f2 = fpss[i, j+half:j+step+1].copy()
            f1[-1] = f2[-1] = 1
            f = poly_mul(f1, f2)
            fpss[i+1, j:j+len(f)-1] = f[:-1]
    fpss2 = np.zeros((bit+1, k+1), dtype=np.int64)
    fpss2[bit, :n] = poly_differential(f)
    for i in range(bit-1, threshold-1, -1):
        step = 2<<i
        half = step>>1
        for j in range(0, n, step):
            if j+half >= n:
                fpss2[i, j:n] = fpss2[i+1, j:n]
                continue
            f = fpss2[i+1, j:min(j+step, n)]
            f1 = fpss[i, j:j+half+1].copy()
            f2 = fpss[i, j+half:min(j+step, n)+1].copy()
            f1[-1] = f2[-1] = 1
            fpss2[i, j:j+half] = poly_mod(f, f1)
            fpss2[i, j+half:min(j+step, n)] = poly_mod(f, f2)
    xs = as_strided(xs, (k>>threshold, 1<<threshold), (8<<threshold, 8))
    xs_cumprod = np.ones_like(xs)
    f = np.zeros_like(xs)
    for i in range(1<<threshold):
        f += fpss2[threshold, i:k:1<<threshold, None] * xs_cumprod % mod
        xs_cumprod *= xs
        xs_cumprod %= mod
    f = f.ravel()
    for j in range(n):
        fpss2[0, j] = ys[j] * pow(int(f[j]), mod-2, mod) % mod
    for i in range(bit):
        step = 2 << i
        half = step >> 1
        for j in range(0, k, step):
            if j+half >= n:
                fpss2[i+1, j:n] = fpss2[i, j:n]
                continue
            f1 = fpss[i, j:j+half+1].copy()
            f2 = fpss[i, j+half:j+step+1].copy()
            f1[-1] = f2[-1] = 1
            fpss2[i+1, j:min(j+step, n)] = (
                poly_mul(fpss2[i, j:j+half], f2)
                + poly_mul(fpss2[i, j+half:min(j+step, n)], f1)
            ) % mod
    return fpss2[bit, :n]

mod = int(input())
A = np.array(input().split(), dtype=np.int64)
X = np.arange(mod, dtype=np.int64)
Ans = polynomial_interpolation(X, A)
print((" ".join(map(str, Ans.tolist()))))