Hello,

I'd like to create some classes to deal with Lorentz vectors (the "natural" language for special relativity). They're basically vectors on R⁴, with a so-called minkowskian metric η = diag(1, -1, -1, -1). One can "transpose" a vector using the metric, vᵗ = ηv. Then, the "square" of a vector is defined as v² = vᵗv = vηv (1). Component-wise, the square looks like (a b c d)² = a² - b² - c² - d² (2).

I wonder how to compute transposed vectors (and square) efficiently (using Eigen of course). I can find two easy ways: computing the product as in (1) or using the components directly as in (2).

Thanks
Louis

There is currently no easy way to avoid the product by 1,-1,-1,-1, and both (1) and (2) will produce the exact same assembly, as verified using the following code:

Code: Select all

static const Vector4d eta(1,-1,-1,-1);

void t(const Vector4d &v, Vector4d &vt) {
  vt = eta.array() * v.array();
}
double sqNorm_v1(const Vector4d &v) {
  return (v.transpose() * eta.asDiagonal() * v)(0,0);
}
double sqNorm_v2(const Vector4d &v) {
  return (eta.array() * v.array().square()).sum();
}


sqNorm_v* produces: (using clang+AVX2)

Code: Select all

        vmovupd   (%rdi), %ymm0
   vmulpd   LCPI2_0(%rip), %ymm0, %ymm1
   vmulpd   %ymm1, %ymm0, %ymm0
   vperm2f128   $1, %ymm0, %ymm0, %ymm1 ## ymm1 = ymm0[2,3,0,1]
   vhaddpd   %ymm1, %ymm0, %ymm0
   vhaddpd   %ymm0, %ymm0, %ymm0


In theory, the following should be more efficient:

Code: Select all

double sqNorm_v3(const Vector4d &v) {
  Array4d w = v;
  w.x() = -w.x();
  return -(v.array() * w).sum();
}


but the compiler seems to fail to implement w.x() = -w.x(); as a simple vxorpd without register spilling.

Thank you ggael, playing with Eigen is gonna be fun.