import numpy.random as rd
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial


def frequence(n,p,repetitions):
    L=[0 for i in range(0, n*p+1)]
    for nb in range(repetitions):
        X = rd.randint(1, n + 1, p)
        s = sum(X)
        L[s]+=1
    L=[L[i]/repetitions for i in range(0, n*p+1)]
    return L

def proba(n,p):
    L=[0]+n*[1/n]
    P=Polynomial(L)
    Q=P**p
    return Q.coef



n = 20
p = 10
repetitions = 5000

plt.plot(range(0, n*p+1), frequence(n,p,repetitions))
plt.plot(range(0, n*p+1), proba(n,p))
plt.show()