/* some hyper complex functions */
  /* see Fractint.c for a description of the "include"  hierarchy */
#include "port.h"
#include "prototyp.h"

void HComplexMult(_HCMPLX *arg1, _HCMPLX *arg2, _HCMPLX *out)
{
    /* it is possible to reoganize this code and reduce the multiplies
        from 16 to 10, but on my 486 it is SLOWER !!! so I left it
        like this - Tim Wegner */
    out->x = arg1->x * arg2->x - arg1->y * arg2->y
           - arg1->z * arg2->z + arg1->t * arg2->t;
    out->y = arg1->y * arg2->x + arg1->x * arg2->y
           - arg1->t * arg2->z - arg1->z * arg2->t;
    out->z = arg1->z * arg2->x - arg1->t * arg2->y
           + arg1->x * arg2->z - arg1->y * arg2->t;
    out->t = arg1->t * arg2->x + arg1->z * arg2->y
           + arg1->y * arg2->z + arg1->x * arg2->t;
}

void HComplexSqr(_HCMPLX *arg, _HCMPLX *out)
{
    out->x = arg->x * arg->x - arg->y * arg->y
           - arg->z * arg->z + arg->t * arg->t;
    out->y = 2 * arg->x * arg->y - 2 * arg->z * arg->t;
    out->z = 2 * arg->z * arg->x - 2 * arg->t * arg->y;
    out->t = 2 * arg->t * arg->x + 2 * arg->z * arg->y;
}

int HComplexInv(_HCMPLX *arg, _HCMPLX *out)
{
   double det, mod, xt_minus_yz;

   det = (sqr(arg->x - arg->t) + sqr(arg->y + arg->z))*
           (sqr(arg->x + arg->t) + sqr(arg->y - arg->z));

   if(det == 0.0)
      return(-1);
   mod = sqr(arg->x) + sqr(arg->y) + sqr(arg->z) + sqr(arg->t);
   xt_minus_yz = arg->x * arg->t - arg->y * arg->z;

   out->x = ( arg->x * mod - 2 * arg->t * xt_minus_yz)/det;
   out->y = (-arg->y * mod - 2 * arg->z * xt_minus_yz)/det;
   out->z = (-arg->z * mod - 2 * arg->y * xt_minus_yz)/det;
   out->t = ( arg->t * mod - 2 * arg->x * xt_minus_yz)/det;
   return(0);
}

void HComplexAdd(_HCMPLX *arg1, _HCMPLX *arg2, _HCMPLX *out)
{
    out->x = arg1->x + arg2->x;
    out->y = arg1->y + arg2->y;
    out->z = arg1->z + arg2->z;
    out->t = arg1->t + arg2->t;
}

void HComplexSub(_HCMPLX *arg1, _HCMPLX *arg2, _HCMPLX *out)
{
    out->x = arg1->x - arg2->x;
    out->y = arg1->y - arg2->y;
    out->z = arg1->z - arg2->z;
    out->t = arg1->t - arg2->t;
}

void HComplexMinus(_HCMPLX *arg1, _HCMPLX *out)
{
    out->x = -arg1->x;
    out->y = -arg1->y;
    out->z = -arg1->z;
    out->t = -arg1->t;
}

/* extends the unary function f to *h1 */
void HComplexTrig0(_HCMPLX *h, _HCMPLX *out)
{
    /* This is the whole beauty of Hypercomplex numbers - *ANY* unary
       complex valued function of a complex variable can easily
       be generalized to hypercomplex numbers */

    _CMPLX a,b, resulta,resultb;

    /* convert to duplex form */
    a.x = h->x - h->t;
    a.y = h->y + h->z;
    b.x = h->x + h->t;
    b.y = h->y - h->z;

    /* apply function to each part */
    CMPLXtrig0(a,resulta);
    CMPLXtrig0(b,resultb);

    /* convert back */
    out->x =  (resulta.x + resultb.x)/2;
    out->y =  (resulta.y + resultb.y)/2;
    out->z =  (resulta.y - resultb.y)/2;
    out->t =  (resultb.x - resulta.x)/2;
}