gistfile1.txt

import sympy as sp
import numpy as np

terms, max_pow = 2, 2
q, p = sp.Indexed('q', i), sp.Indexed('p', i)

H = sp.Sum(q**2 + p**2, (i, 0, terms-1)).doit()

def pbrackets(f, g):
    return sp.Sum(f.diff(q)*g.diff(p) - f.diff(p)*g.diff(q), (i, 0, terms-1)).doit().simplify()

poly_terms = [
    np.prod([
        q.base[n]**q_pow[n] * p.base[n]**p_pow[n]
        for n in range(terms)
    ])
    for q_pow in product(list(range(max_pow)), repeat=terms)
    for p_pow in product(list(range(max_pow)), repeat=terms)
]

Op = np.array([
    [sp.Poly(result, *poly_terms).coeff_monomial(term) for term in poly_terms]
    for result in map(lambda t: pbrackets(H, t), poly_terms)
], dtype=np.int)

for v in sp.Matrix(*Op.shape, Op.flatten()).nullspace():
    print(v.dot(poly_terms))

gistfile2.txt

1
p[1]*q[1]
-p[0]*q[1] + p[1]*q[0]
p[0]*q[0]
p[0]*p[1] + q[0]*q[1]
p[0]*p[1]*q[0]*q[1]