#include "config.h"
#include <complex.h>
#include <string.h>
#include "vmath.h"

#define la_int int

extern void
F77_FUNC(zaxpy) (const la_int *n, const void *za, const void *zx,
	         const la_int *incx, void *zy, const la_int *incy);


static double complex
complex_mul (double complex x, double complex y)
{
	double xr = creal(x), xi = cimag(x);
	double yr = creal(y), yi = cimag(y);
	double zr = xr * yr - xi * yi;
	double zi = xr * yi + xi * yr;
        double complex z = zr + I * zi;	
	return z;
}


static void
vzClear (int n, double complex *z)
{
        memset(z, 0, n * 2 * sizeof(double));
}

static void
vzCopy (int n, const double complex *x, double complex *z)
{
        if (x != z) {
                memcpy(z, x, n * 2 * sizeof(double));
        }
}

void
vzConj (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = conj(x[i]);
        }
}

void
vzScale (int n, const double complex *palpha, const double complex *x, double complex *z)
{
        double complex alpha = *palpha;
        
        if (alpha == 1) {
                vzCopy(n, x, z);
        } else if (alpha == 0) {
                vzClear(n, z);
        } else {
                int i;
                for (i = 0; i < n; i++) {
                        z[i] = complex_mul(alpha, x[i]);
                }
        }
        
        /* Note: Using zscal sometimes gives differences in the least
         * significant bit of the answer */
}

void
vzShift (int n, const double complex *palpha, const double complex *x, double complex *z)
{
        double complex alpha = *palpha;
        
        if (alpha == 0) {
                vzCopy(n, x, z);
        } else {
                int i;
                for (i = 0; i < n; i++) {
                        z[i] = alpha + x[i];
                }
        }
}

void
vzNeg (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = -x[i];
        }
}

void
vzAbs (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = cabs(x[i]) + 0*I;
        }
}

static void
zsgn (const double complex x, double complex *z)
{
        if (x == 0) {
                *z = x;
        } else {
                double r = cabs(x);
                *z = creal(x)/r + cimag(x)/r*I;
        }
}

void
vzSgn (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                zsgn(x[i], &(z[i]));
        }
}

void
vzInv (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = 1.0 / x[i];
        }
}

static void
vzAxpy (int n, const double complex *palpha, const double complex *x, const double complex *y, double complex *z)
{
        double complex alpha = *palpha;
        la_int one = 1;
        
        if (alpha == 0) {
                vzCopy(n, y, z);
        } else if (y == z) {
                F77_FUNC(zaxpy) (&n, palpha, x, &one, z, &one);
        } else if (alpha == 1 && x == z) {
                F77_FUNC(zaxpy) (&n, palpha, y, &one, z, &one);
        } else {
                int i;
                for (i = 0; i < n; i++) {
                        z[i] = alpha * x[i] + y[i];
                }
        }
}

void
vzAxpby (int n, const double complex *palpha, const double complex *x, const double complex *pbeta, const double complex *y, double complex *z)
{
        double complex alpha = *palpha;
        double complex beta = *pbeta;        
                
        if (alpha == 0) {
                vzScale(n, pbeta, y, z);
        } else if (alpha == 1) {
                vzAxpy(n, pbeta, y, x, z);
        } else if (beta == 0) {
                vzScale(n, palpha, x, z);
        } else if (beta == 1) {
                vzAxpy(n, palpha, x, y, z);
        } else {
                int i;
                for (i = 0; i < n; i++) {
                        z[i] = alpha * x[i] + beta * y[i];
                }
        }
}

void
vzAdd (int n, const double complex *x, const double complex *y, double complex *z)
{
        vdAdd(2 * n, (const double *)x, (const double *)y, (double *)z);
}

void vzSub (int n, const double complex *x, const double complex *y, double complex *z)
{
        vdSub(2 * n, (const double *)x, (const double *)y, (double *)z);
}

void
vzMul (int n, const double complex *x, const double complex *y, double complex *z)
{
        /* Note: using ztbmv sometimes gives different answers in the lower
         * bits of precision */
        /*
        if (y == z) {
                blas_ztbmv(BlasUpper, BlasNoTrans, BlasNonUnit, n, 0, x, 1,
                           z, 1);
        } else if (x == z) {
                blas_ztbmv(BlasUpper, BlasNoTrans, BlasNonUnit, n, 0, y, 1,
                           z, 1);
        } */
        int i;
        for (i = 0; i < n; i++) {
                z[i] = complex_mul(x[i], y[i]);
        }
}

void
vzDiv (int n, const double complex *x, const double complex *y, double complex *z)
{
        /* Note: using dtbsv gives slightly different answers
         */
        /*
        if (y == z) {
                blas_ztbsv(BlasUpper, BlasNoTrans, BlasNonUnit, n, 0, x, 1,
                           z, 1);
        } else */
        if (x == z) {
                int i;
                for (i = 0; i < n; i++) {
                        z[i] = z[i] / y[i];
                }
        } else {
                int i;
                for (i = 0; i < n; i++) {
                        z[i] = x[i] / y[i];
                }
        }
}

void
vzExp (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = cexp(x[i]);
        }
}

void
vzSqrt (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = csqrt(x[i]);
        }
}

void
vzLog (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = clog(x[i]);
        }        
}

void
vzPow (int n, const double complex *x, const double complex *y, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = cpow(x[i], y[i]);
        }        
}

void
vzSin (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = csin(x[i]);
        }                
}

void
vzCos (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = ccos(x[i]);
        }                        
}

void
vzTan (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = ctan(x[i]);
        }                                
}

void
vzASin (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = casin(x[i]);
        }                                        
}

void
vzACos (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = cacos(x[i]);
        }                                                
}

void
vzATan (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = catan(x[i]);
        }
}

void
vzSinh (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = csinh(x[i]);
        }                
}

void
vzCosh (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = ccosh(x[i]);
        }                        
}

void
vzTanh (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = ctanh(x[i]);
        }                                
}

void
vzASinh (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = casinh(x[i]);
        }                                        
}

void
vzACosh (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = cacosh(x[i]);
        }                                                
}

void
vzATanh (int n, const double complex *x, double complex *z)
{
        int i;
        for (i = 0; i < n; i++) {
                z[i] = catanh(x[i]);
        }
}