vec xmul(vec a, vec b) {
  // thanks fractalforums.com
  vec res = void;
  auto ab = a * b;
  auto
    rsq = abs(a.dot(b)),
    r = sqrt(rsq),
    yang = atan2(ab.xy().length, a.z),
    zang = atan2(a.y, b.x);
  res.x = rsq * sin(yang*2 + PI/2) * cos(zang*2 + PI);
  res.y = rsq * sin(yang*2 + PI/2) * sin(zang*2 + PI);
  res.z = rsq * cos(yang*2 + PI/2);
  return res;
}

vec xsqr(vec v) {
  vec res = void;
  auto vv = v*v;
  auto
    r = v.length,
    yang = atan2(vv.xy().length, v.z),
    zang = atan2(v.y, v.x);
  res.x = (r*r) * sin(yang*2 + PI/2) * cos(zang*2 + PI);
  res.y = (r*r) * sin(yang*2 + PI/2) * sin(zang*2 + PI);
  res.z = (r*r) * cos(yang*2 + PI/2);
  return res;
}

vec xpow(int i)(vec v) {
  static if (i == 1) return v;
  else static if (i == 2) return v*v;
  else static if ((i%2) == 1) return xmul(xsqr(xpow!(i/2)(v)), v);
  else return xsqr(xpow!(i/2)(v));
}