Revision history [back]

You can use polynomial rings over polynomial rings:

A.<a1,b1,c1,d1> = PolynomialRing(QQ, order='lex')
B.<a2,b2,c2,d2> = PolynomialRing(A)
Q.<i,j,k> = QuaternionAlgebra(B.fraction_field(), -1, -1)
q1 = a1 + b1*i + c1*j + d1*k
q2 = a2 + b2*i + c2*j + d2*k
I = A.ideal(sum([B(f).coefficients() for f in (q1*q2 - q2*q1).coefficient_tuple()], []))
I.groebner_basis()

Output:

[b1, c1, d1]

meaning $b_1=c_1=d_1=0$ needs to hold for $a_1 + b_1i + c_1j + d_1k$ to be in the center, so the center is $\mathbb{R}$.