def pascal_matrix(size=100, zp=2):
    '''
    Generates pascal_matrix, i.e. the matrix whose diagonals are binomial coeficients.
    Arguments:
        @size: Matrix size (pascal_matrix is square)
        @zp: pascal_matrix entries are divided modulus zp

    Returns:
        The size x size matrix whose ith antidiagonal corresponds to the binomial coefficients of 
    (1 + x)^i mod zp
    '''
    M = ones((size, size), dtype=int)
    for i in range(1, size):
        for j in range(1, size):
            M[i, j] = int(M[i - 1, j] + M[i, j - 1])
        if zp != 1 and i > 1:
            for j in range(1, size):
                M[i - 1, j] %= zp
    if zp != 1:
        for j in range(1, size):
            M[-1, j] %= zp
    return M

#TODO: Check why pascal_matrix is much more faster in sage
M = pascal_matrix(zp=2)

imshow(M)