1 /* some hyper complex functions */
    2   /* see Fractint.c for a description of the "include"  hierarchy */
    3 #include "port.h"
    4 #include "prototyp.h"
    5 
    6 void HComplexMult(_HCMPLX *arg1, _HCMPLX *arg2, _HCMPLX *out)
    7 {
    8     /* it is possible to reoganize this code and reduce the multiplies
    9         from 16 to 10, but on my 486 it is SLOWER !!! so I left it
   10         like this - Tim Wegner */
   11     out->x = arg1->x * arg2->x - arg1->y * arg2->y
   12            - arg1->z * arg2->z + arg1->t * arg2->t;
   13     out->y = arg1->y * arg2->x + arg1->x * arg2->y
   14            - arg1->t * arg2->z - arg1->z * arg2->t;
   15     out->z = arg1->z * arg2->x - arg1->t * arg2->y
   16            + arg1->x * arg2->z - arg1->y * arg2->t;
   17     out->t = arg1->t * arg2->x + arg1->z * arg2->y
   18            + arg1->y * arg2->z + arg1->x * arg2->t;
   19 }
   20 
   21 void HComplexSqr(_HCMPLX *arg, _HCMPLX *out)
   22 {
   23     out->x = arg->x * arg->x - arg->y * arg->y
   24            - arg->z * arg->z + arg->t * arg->t;
   25     out->y = 2 * arg->x * arg->y - 2 * arg->z * arg->t;
   26     out->z = 2 * arg->z * arg->x - 2 * arg->t * arg->y;
   27     out->t = 2 * arg->t * arg->x + 2 * arg->z * arg->y;
   28 }
   29 
   30 int HComplexInv(_HCMPLX *arg, _HCMPLX *out)
   31 {
   32    double det, mod, xt_minus_yz;
   33 
   34    det = (sqr(arg->x - arg->t) + sqr(arg->y + arg->z))*
   35            (sqr(arg->x + arg->t) + sqr(arg->y - arg->z));
   36 
   37    if(det == 0.0)
   38       return(-1);
   39    mod = sqr(arg->x) + sqr(arg->y) + sqr(arg->z) + sqr(arg->t);
   40    xt_minus_yz = arg->x * arg->t - arg->y * arg->z;
   41 
   42    out->x = ( arg->x * mod - 2 * arg->t * xt_minus_yz)/det;
   43    out->y = (-arg->y * mod - 2 * arg->z * xt_minus_yz)/det;
   44    out->z = (-arg->z * mod - 2 * arg->y * xt_minus_yz)/det;
   45    out->t = ( arg->t * mod - 2 * arg->x * xt_minus_yz)/det;
   46    return(0);
   47 }
   48 
   49 void HComplexAdd(_HCMPLX *arg1, _HCMPLX *arg2, _HCMPLX *out)
   50 {
   51     out->x = arg1->x + arg2->x;
   52     out->y = arg1->y + arg2->y;
   53     out->z = arg1->z + arg2->z;
   54     out->t = arg1->t + arg2->t;
   55 }
   56 
   57 void HComplexSub(_HCMPLX *arg1, _HCMPLX *arg2, _HCMPLX *out)
   58 {
   59     out->x = arg1->x - arg2->x;
   60     out->y = arg1->y - arg2->y;
   61     out->z = arg1->z - arg2->z;
   62     out->t = arg1->t - arg2->t;
   63 }
   64 
   65 void HComplexMinus(_HCMPLX *arg1, _HCMPLX *out)
   66 {
   67     out->x = -arg1->x;
   68     out->y = -arg1->y;
   69     out->z = -arg1->z;
   70     out->t = -arg1->t;
   71 }
   72 
   73 /* extends the unary function f to *h1 */
   74 void HComplexTrig0(_HCMPLX *h, _HCMPLX *out)
   75 {
   76     /* This is the whole beauty of Hypercomplex numbers - *ANY* unary
   77        complex valued function of a complex variable can easily
   78        be generalized to hypercomplex numbers */
   79 
   80     _CMPLX a,b, resulta,resultb;
   81 
   82     /* convert to duplex form */
   83     a.x = h->x - h->t;
   84     a.y = h->y + h->z;
   85     b.x = h->x + h->t;
   86     b.y = h->y - h->z;
   87 
   88     /* apply function to each part */
   89     CMPLXtrig0(a,resulta);
   90     CMPLXtrig0(b,resultb);
   91 
   92     /* convert back */
   93     out->x =  (resulta.x + resultb.x)/2;
   94     out->y =  (resulta.y + resultb.y)/2;
   95     out->z =  (resulta.y - resultb.y)/2;
   96     out->t =  (resultb.x - resulta.x)/2;
   97 }
   98