libsimple_geom/unit_tests/sym.cpp

#include <iostream>
#include <string>
#include <cassert>
#include <sstream>
#include <algorithm>
#include "simple/geom/vector.hpp"
#include "simple/geom/algorithm.hpp"

using namespace std::literals;
using namespace simple::geom;

std::stringstream ss;

struct sym
{
	std::string val;

#define DEFINE_SYM_OPERATOR(op) \
	sym operator op(sym s) const \
	{ \
		if(!val.empty()) \
			s.val = val + #op + s.val; \
		return s; \
	} \
	sym& operator op##=(sym s) \
	{ \
		(*this) = (*this) op s; \
		return (*this); \
	}
	DEFINE_SYM_OPERATOR(+)
	DEFINE_SYM_OPERATOR(-)
	DEFINE_SYM_OPERATOR(*)
	DEFINE_SYM_OPERATOR(/)
	DEFINE_SYM_OPERATOR(%)
#undef DEFINE_SYM_OPERATOR

	friend std::ostream& operator<<(std::ostream& os, const sym& s)
	{
		return os << s.val;
	}
};

template <typename... Cs>
auto sym_vector(Cs... cs)
{
	return vector(sym{cs}...);
}

template<typename T, size_t D, typename O>
const auto& vector_only_geometric_product = dot_wedge<T, D, O>;

template<typename T, size_t D, typename O>
const auto& make_rotor = dot_wedge<T, D, O>;


void Explode()
{

	const auto wrap_in_vector = [](const auto &x) {return vector(x);};

	auto a = vector( sym_vector("x1", "x2", "x3") ); // 3x1
	auto b = sym_vector("y1", "y2", "y3").transformed(wrap_in_vector); // 1x3

	ss << a(b); // 3x1 * 1x3 = 3x3

	assert( ss.str() ==
R"(^^^
(x1*y1, x2*y1, x3*y1)
(x1*y2, x2*y2, x3*y2)
(x1*y3, x2*y3, x3*y3)
vvv
)"
	); ss.str("");

	ss << b(a); // 1x3 * 3x1 = 1x1
	assert( ss.str() ==
R"(^
(y1*x1+y2*x2+y3*x3)
v
)"
	); ss.str("");


	auto ab = vector(a(b)); // 3x3x1

	ss << ab(b); // 3x3x1 * 1x3 = 3x3x3
	assert( ss.str() ==
R"(^^^
^^^
(x1*y1*y1, x2*y1*y1, x3*y1*y1)
(x1*y2*y1, x2*y2*y1, x3*y2*y1)
(x1*y3*y1, x2*y3*y1, x3*y3*y1)
vvv

^^^
(x1*y1*y2, x2*y1*y2, x3*y1*y2)
(x1*y2*y2, x2*y2*y2, x3*y2*y2)
(x1*y3*y2, x2*y3*y2, x3*y3*y2)
vvv

^^^
(x1*y1*y3, x2*y1*y3, x3*y1*y3)
(x1*y2*y3, x2*y2*y3, x3*y2*y3)
(x1*y3*y3, x2*y3*y3, x3*y3*y3)
vvv

vvv
)"
	); ss.str("");

	ss << b(ab); // 1x3 * 3x3x1 = 1x3x1
	assert( ss.str() ==
R"(^^^
^
(y1*x1*y1+y2*x2*y1+y3*x3*y1)
(y1*x1*y2+y2*x2*y2+y3*x3*y2)
(y1*x1*y3+y2*x2*y3+y3*x3*y3)
v

vvv
)"
	); ss.str("");

	ss << (ab(b))(ab); // brace yourselves: 3x3x3 * 3x3x1 = 3x3x3x1
	assert( ss.str() ==
R"(^^^
^^^
^^^
(x1*y1*y1*x1*y1+x1*y1*y2*x2*y1+x1*y1*y3*x3*y1, x2*y1*y1*x1*y1+x2*y1*y2*x2*y1+x2*y1*y3*x3*y1, x3*y1*y1*x1*y1+x3*y1*y2*x2*y1+x3*y1*y3*x3*y1)
(x1*y2*y1*x1*y1+x1*y2*y2*x2*y1+x1*y2*y3*x3*y1, x2*y2*y1*x1*y1+x2*y2*y2*x2*y1+x2*y2*y3*x3*y1, x3*y2*y1*x1*y1+x3*y2*y2*x2*y1+x3*y2*y3*x3*y1)
(x1*y3*y1*x1*y1+x1*y3*y2*x2*y1+x1*y3*y3*x3*y1, x2*y3*y1*x1*y1+x2*y3*y2*x2*y1+x2*y3*y3*x3*y1, x3*y3*y1*x1*y1+x3*y3*y2*x2*y1+x3*y3*y3*x3*y1)
vvv

^^^
(x1*y1*y1*x1*y2+x1*y1*y2*x2*y2+x1*y1*y3*x3*y2, x2*y1*y1*x1*y2+x2*y1*y2*x2*y2+x2*y1*y3*x3*y2, x3*y1*y1*x1*y2+x3*y1*y2*x2*y2+x3*y1*y3*x3*y2)
(x1*y2*y1*x1*y2+x1*y2*y2*x2*y2+x1*y2*y3*x3*y2, x2*y2*y1*x1*y2+x2*y2*y2*x2*y2+x2*y2*y3*x3*y2, x3*y2*y1*x1*y2+x3*y2*y2*x2*y2+x3*y2*y3*x3*y2)
(x1*y3*y1*x1*y2+x1*y3*y2*x2*y2+x1*y3*y3*x3*y2, x2*y3*y1*x1*y2+x2*y3*y2*x2*y2+x2*y3*y3*x3*y2, x3*y3*y1*x1*y2+x3*y3*y2*x2*y2+x3*y3*y3*x3*y2)
vvv

^^^
(x1*y1*y1*x1*y3+x1*y1*y2*x2*y3+x1*y1*y3*x3*y3, x2*y1*y1*x1*y3+x2*y1*y2*x2*y3+x2*y1*y3*x3*y3, x3*y1*y1*x1*y3+x3*y1*y2*x2*y3+x3*y1*y3*x3*y3)
(x1*y2*y1*x1*y3+x1*y2*y2*x2*y3+x1*y2*y3*x3*y3, x2*y2*y1*x1*y3+x2*y2*y2*x2*y3+x2*y2*y3*x3*y3, x3*y2*y1*x1*y3+x3*y2*y2*x2*y3+x3*y2*y3*x3*y3)
(x1*y3*y1*x1*y3+x1*y3*y2*x2*y3+x1*y3*y3*x3*y3, x2*y3*y1*x1*y3+x2*y3*y2*x2*y3+x2*y3*y3*x3*y3, x3*y3*y1*x1*y3+x3*y3*y2*x2*y3+x3*y3*y3*x3*y3)
vvv

vvv

vvv
)"
	); ss.str("");


	auto a2 = vector( sym_vector("a", "b") ); // 2x1
	auto b2 = sym_vector("c", "d").transformed(wrap_in_vector); // 1x2

	ss << a2(b2) << '\n';

	ss << dot_wedge(
		sym_vector("a", "b"),
		sym_vector("c", "d")
	) << '\n';

	ss << dot_wedge(
		sym_vector("u1", "u2", "u3"),
		sym_vector("v1", "v2", "v3")
	) << '\n';

	// (a + ib)(c + id) = ac + aid + ibc + iibd = (ac - bd) + (ad + bc)i
	// (a + ib)(c - id) = ac - aid + ibc - iibd = (ac + bd) + (bc - ad)i
	// ((ac - bd) + (ad + bc)i)(a-ib) = (ac - bd)(a-ib) + (ad+bc)i(a-ib) =
	// = aca - acib - bda + bdib + adia + bcia - adiib - bciib =
	// = aca - bda + bdib + adia + adb + bcb

	assert(ss.str() ==
R"(^^
(a*c, b*c)
(a*d, b*d)
vv

(b*c-a*d, a*c+b*d)
(u2*v1-u1*v2, u3*v1-u1*v3, u3*v2-u2*v3, u1*v1+u2*v2+u3*v3)
)"
	); ss.str("");


	auto aa = vector(
		sym_vector("x1", "x2"),
		sym_vector("x3", "x4")
	);
	auto bb = vector(
		sym_vector("y1", "y2"),
		sym_vector("y3", "y4")
	);
	auto A = vector( aa ); // 2x2x1
	auto B = bb.transformed([&](auto x){ return x.transformed(wrap_in_vector);}); // 1x2x2
	ss << A(B); // 2x2x1 * 1x2x2 = 2x2x2x2
	assert(ss.str() ==
R"(^^
^^
^^
(x1*y1, x2*y1)
(x3*y1, x4*y1)
vv

^^
(x1*y2, x2*y2)
(x3*y2, x4*y2)
vv

vv

^^
^^
(x1*y3, x2*y3)
(x3*y3, x4*y3)
vv

^^
(x1*y4, x2*y4)
(x3*y4, x4*y4)
vv

vv

vv
)"
	); ss.str("");

}

int main()
{
	Explode();

	{
		auto a = sym_vector("a", "b", "c");
		auto m = vector(
			sym_vector("x1", "y1", "z1"),
			sym_vector("x2", "y2", "z2"),
			sym_vector("x3", "y3", "z3")
		);

		std::cout << m(a) << '\n';
		std::cout << a(m) << '\n';
	}

	{
		auto a = sym_vector("x1", "y1", "z1");
		auto b = sym_vector("x2", "y2", "z2");
		std::cout << "=============== EXPAND ============" << '\n';
		std::cout << dot_wedge(b,a) << '\n';
	}


	// return {y1*z2 - y2*z1, z1*x2 - z2*x1, x1*y2 - x2*y1};

	return 0;
}