Gram-Schmidt process with sympy using a custom scalar product.

gram_schmidt.py

import sympy as sp

x = sp.Symbol("x")

def gs(vecs, scalar_product):
    new_base = [
        None for _ in range(len(vecs))
    ]
    new_base[0] = vecs[0] / norm(vecs[0], scalar_product = scalar_product)

    for i in range(1, len(vecs)):
        b_tilde = vecs[i] - sum(
            scalar_product(new_base[k], vecs[i]) * new_base[k]
            for k in range(i)
        )
        new_base[i] = b_tilde / norm(b_tilde, scalar_product = scalar_product)

    return new_base

def norm(a, scalar_product):
    return sp.sqrt(scalar_product(a, a))

def scalar_product(a, b):
    l = [sp.Integer(1), sp.Integer(2), sp.Integer(3)]
    total = sp.Integer(0)
    for k in range(len(l)):
        total += a.subs(x, (l[k])) * b.subs(x, (l[k]))

    return total

old_base = [sp.Integer(1), x, x**2]
print(gs(old_base, scalar_product))