SIMD-friendly complex numbers with redundancy

complex3.c

/* example implementation -- not recommended for actual use, but it would work!
 * (SIMD-friendly in principle; not shown)
 * (probably not ideal for memory-bound workloads; uses 50% more space)
 * based on two insights:
 * - complex multiplication (a+bi)(c+di) with i^2 = -1 is equivalent to 
 *   split-complex multiplication (a+bk)(c+dk) - 2*b*d with k^2 = +1
 * - the split-complex numbers are isomorphic to the direct product R*R with the
 *   basis (1+k)/sqrt(2), (1-k)/sqrt(2), both idempotents, thus split-complex
 *   multiplication (a+bk)(c+dk) can be done by (a+b)(c+d) and (a-b)(c-d)
 */

/* we use long, but could replace with double, char etc */
typedef struct {long a, b, c;} complex3;

complex3 encode3(long x, long y) {
  complex3 ret = {.a = x + y, .b =  x - y, .c = y << 1};
  return ret;
}

/* addition, negation, subtraction all componentwise */
complex3 add3(complex3 l, complex3 r) {
  return (complex3) {.a = l.a + r.a, .b = l.b + r.b, .c = l.c + r.c};
}

complex3 neg3(complex3 z) {
  return (complex3) {.a = -z.a, .b = -z.b, .c = -z.c};
}

complex3 sub3(complex3 l, complex3 r) {
  return add3(l, neg3(r));
}

/* the good part: multiplication via
 * - one round of 3 multiplications (parallelizeable)
 * - one one-bit bitshift (fast, I guess?)
 * - one round of 3 subtractions (parallelizeable)
 */
complex3 multiply3(complex3 l, complex3 r) {
  long p[3] = {l.a * r.a, l.b * r.b, l.c * r.c >> 1};
  return (complex3) {.a = p[0] - p[2], .b =  p[1] - p[2], .c =  p[0] - p[1]};
}

/* even norms and conjugates are pretty easy! */
long norm3(complex3 z) {
  return (z.a*z.a + z.b*z.b) >> 1;
}

complex3 conj3(complex3 z) {
  return (complex3) {z.b, z.a, -z.c};
}

complex3 scalardiv3(long k, complex3 z) {
  return (complex3) {z.a / k, z.b / k, z.c / k};
}

complex3 divide3(complex3 l, complex3 r) {
  return scalardiv3(norm3(r), multiply3(l, conj3(r)));
}

long real3(complex3 z) {
  return (z.a + z.b) >> 1;
}

long imag3(complex3 z) {
  return z.c >> 1;
}