void q_exp(point c, const point q)
{
static double n, f, ex;
n = sqrt(q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
ex = exp(q[0]);
f = ex*sin(n);
if (n != 0.0) f /= n;
c[0] = ex*cos(n);
c[1] = f*q[1];
c[2] = f*q[2];
c[3] = f*q[3];
}
void q_log(point c, const point q)
{
static double n, f;
n = sqrt(q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
if (n != 0.0)
{
f = atan2(n, q[0]);
f /= n;
c[0] = 0.5*log(q[0]*q[0] + n*n);
c[1] = f*q[1];
c[2] = f*q[2];
c[3] = f*q[3];
return;
}
else /* special case: real number */
{
c[0] = 0.5*log(q[0]*q[0]);
c[1] = atan2(0.0,q[0]);
c[2] = 0.0;
c[3] = 0.0;
return;
/* comment on asymmetry in i, j and k:
e^(i*pi)=-1 && e^(j*pi)=-1 && e^(k*pi)=-1
[ and even linear combinations of i, j and k give:
e^(a*i+b*j+c*k) = -1 for every a, b, c E R with:
sqrt(a^2+b^2+c^2) = pi ]
the above means that there is an infinite number of solutions to ln(-1).
Practically, one must choose one of them and thus break symmetry.
[ Even then there still is an infinite number of solutions due to
periodicity of sin and cos... ] */
}
}
void q_pow(point c, const point a, const point b)
{
/* if a == zero: exp(-inf*b) = 0, if b>0;
exp(-inf*b) = inf, if b<0
if b isn't real, exp(-inf*b) isn't defined, because
lim sin(x) for x->-inf (same with cos) doesn't exists */
static double an, bnp;
static point p;
an = a[0]*a[0] + a[1]*a[1] + a[2]*a[2] + a[3]*a[3];
bnp = b[1]*b[1] + b[2]*b[2] + b[3]*b[3];
if (fabs(an) < 1E-200 && (b[0] > 0.0 || bnp != 0.0) ) {
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0; c[3] = 0.0;
return;
}
q_log(c, a);
q_mul(p, c, b);
q_exp(c, p);
return;
}
I must write both functions directly in assembly, really hard stuff!