#!/usr/bin/env python3

# Deriving the expression of the cross product given only the dot product, a
# standard basis, and known identities.

# https://www.youtube.com/watch?v=L68VR5VuwnI&t=1h6m15s
# https://guide.handmadehero.org/code/day376/#3975

from sympy import *

vec = lambda v: [Symbol(f'{v}_{i}') for i in range(1, 3+1)]
dot = lambda a, b: sum(a_i * b_i for a_i, b_i in zip(a, b))

a = vec('a')
b = vec('b')
c = vec('c')

# `c = cross(a, b)`, solve for `c`.
# |c|   = |a|   |b|   sin(theta)
# |c|^2 = |a|^2 |b|^2 sin(theta)^2
# |c|^2 = |a|^2 |b|^2 (1 - cos(theta)^2)
# |c|^2 = |a|^2 |b|^2 (1 - dot(a, b)^2 / |a|^2 |b|^2)
# |c|^2 = |a|^2 |b|^2    - dot(a, b)^2
# We loose the sign when we square, recovered later.
eqs = [
    Eq(dot(a, c), 0),
    Eq(dot(b, c), 0),
    Eq(dot(c, c), dot(a, a) * dot(b, b) - dot(a, b) * dot(a, b)),
]

sols = solve(eqs, *c, dict=True)
for sol in sols:
    # We check the sign here.
    # dot(cross([1, 0, 0], [0, 1, 0]), [0, 0, 1]) == 1
    c_3 = dot(sol.values(), [0, 0, 1]).subs([*zip(a, [1, 0, 0]), *zip(b, [0, 1, 0])])
    if c_3 == 1:
        for kv in sol.items():
            k, v = map(pretty, kv)
            pprint(f"{k}: {v}")

# c₁: a₂⋅b₃ - a₃⋅b₂
# c₂: -a₁⋅b₃ + a₃⋅b₁
# c₃: a₁⋅b₂ - a₂⋅b₁