Problem: 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.

Submission attempts (in order):
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):
    
    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):
    
    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):
    
    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):
        
        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 = []

#     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)
#





comb = Combination(max=P - 1, mod=P)
ncr = np.vectorize(lambda r: comb.ncr(P - 1, r))
vpow = np.vectorize(lambda a, b: pow(a, b, P))
B = np.zeros(P, dtype=int)
for j, a in enumerate(A):
    if a == 0:
        continue
    B -= vpow(-j, np.arange(P)) * ncr(np.arange(P)) % P
    B[-1] += 1
    B %= P
print(*B[::-1])


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):
    
    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):
    
    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):
    
    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):
        
        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 = []

#     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





comb = Combination(max=P - 1, mod=P)
ncr = np.array([comb.ncr(P - 1, i) for i in range(P)])

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])


Verdicts: Time Limit Exceeded

Accepted