sym.py

import re

from sympy import simplify, expand, sqrt, Sum
from sympy.matrices import Matrix
from sympy.physics.vector.frame import ReferenceFrame
from sympy.vector import CoordSys3D
from sympy.physics.mechanics import cross
from sympy.core.symbol import symbols
from sympy.vector.vector import Cross, Vector

C = CoordSys3D('C')

x1, y1, z1 = symbols("x_1 y_1 z_1", real=True)
x2, y2, z2 = symbols("x_2 y_2 z_2", real=True)
x3, y3, z3 = symbols("x_3 y_3 z_3", real=True)
P1 = [x1, y1, z1]
P2 = [x2, y2, z2]
P3 = [x3, y3, z3]

# https://ictsr4.com/py/l0120.html
def inner(a, b):
    v_o = []
    lng = len(a)
    for i in range(lng * 2):
        v_o.append(a[(i+1) % lng]*b[(i+2) % lng]-a[(i+2) % lng]*b[(i+1) % lng])
    return v_o

def minus(a, b):
    return [ae - be for ae, be in zip(a, b)]

def simpl(a):
    a = expand(a)
    a = simplify(a)
    return a

def convertify(a):
    a = repr(a)
    a = re.sub(r'([xyz])_(\d+)', r'{\g<1>} rsub {\g<2>}', a)
    a = a.replace('*', '')
    return a

def canonicalize(a):
    a_ = 0
    for xx in a:
        a_ += xx
    a_ = sqrt(a_)
    a_ = simpl(a_)
    return [x / a_ for x in a]

P21 = minus(P2, P1)
P32 = minus(P3, P2)

N = inner(P21, P32)
N = [simpl(z) for z in N]
NC = canonicalize(N)
N = [convertify(z) for z in N]

print(N[0])
print(N[1])
print(N[2])

print()

print(NC[0])
print(NC[1])
print(NC[2])