Proof of Euler's four square identity using Python and SymPy

Euler_four_square_identity.py

"""
Euler's four square theorem states that
  (x1^2 + x2^2 + x3^2 + x4^2) * 
  (y1^2 + y2^2 + y3^2 + y4^2) =
  z1^2 + z2^2 + z3^2 + z4^2, 

where

  z1 = x1*y1 + x2*y2 + x3*y3 + x4*y4
  z2 = x1*y2 - x2*y1 - x3*y4 + x4*y3
  z3 = x1*y3 - x3*y1 + x2*y4 - x4*y2
  z4 = x1*y4 - x4*y1 - x2*y3 + x3*y2

PROOF:
"""

from sympy import *

x1, x2, x3, x4 = symbols("x1 x2 x3 x4")
y1, y2, y3, y4 = symbols("y1 y2 y3 y4")

z1 = x1*y1 + x2*y2 + x3*y3 + x4*y4
z2 = x1*y2 - x2*y1 - x3*y4 + x4*y3
z3 = x1*y3 - x3*y1 + x2*y4 - x4*y2
z4 = x1*y4 - x4*y1 - x2*y3 + x3*y2

xx = x1**2 + x2**2 + x3**2 + x4**2
yy = y1**2 + y2**2 + y3**2 + y4**2
zz = z1**2 + z2**2 + z3**2 + z4**2

print expand(xx*yy - zz)