/* 
    complex - functions and definitions for complex types

    Copyright (C) 1995, Adam Fedor
*/

#include <math.h>
#include <stdlib.h>
#include "MathArray/complex.h"

/* Addition, multiplication of complex numbers */
complex_double 
c_add (complex_double x, complex_double y)
{
    x.real += y.real;
    x.imag += y.imag;
    return x;
}

complex_double 
c_sub (complex_double x, complex_double y)
{
    x.real -= y.real;
    x.imag -= y.imag;
    return x;
}

complex_double 
c_mult (complex_double x, complex_double y)
{
    complex_double t;

    t.real = (x.real*y.real - x.imag*y.imag);
    t.imag = (x.real*y.imag + x.imag*y.real);
    return t;
}

complex_double 
c_div (complex_double x, complex_double y)
{
    complex_double t;
    double r;
    double denom;

    /* FIXME: raise exception if denom == 0? */
    if (fabs(y.real) >= fabs(y.imag)) {
	r = y.imag / y.real;
	denom = y.real + r * y.imag;
	t.real = (x.real + r * x.imag) / denom;
	t.imag = (x.imag - r * x.real) / denom;
    } else {
	r = y.real / y.imag;
	denom = y.imag + r * y.real;
	t.real = (x.real * r + x.imag) / denom;
	t.imag = (x.imag * r - x.real) / denom;
    }
    return t;
}


/* Functions you typically find with complex numbers */
double 
c_abs (complex_double z)
{
    double x, y, t, temp;
    x = fabs(z.real);
    y = fabs(z.imag);
    if (x == 0.0)
	t = y;
    else if (y == 0.0)
	t = x;
    else if (x > y) {
	temp = y / x;
	t = x * sqrt(1.0+temp*temp);
    } else {
	temp = x / y;
	t = y * sqrt(1.0+temp*temp);
    }
    return t;
}

double 
c_imag (complex_double z)
{
    return z.imag;
}

double 
c_real (complex_double z)
{
    return z.real;
}

complex_double 
c_set (double real, double imag)
{
    complex_double d;
    d.real = real;
    d.imag = imag;
    return d;
}

complex_double 
c_f2d (complex_float f)
{
    complex_double d;

    d.real = f.real;
    d.imag = f.imag;
    return d;
}

complex_float  
c_d2f (complex_double d)
{
    complex_float f;

    f.real = d.real;
    f.imag = d.imag;
    return f;
}

complex_double 
c_cos (complex_double z)
{
    complex_double t;

    t.real = cos(z.real) * cosh(z.imag);
    t.imag = -sin(z.real) * sinh(z.imag);
    return t;
}

complex_double 
c_exp (complex_double z)
{
    complex_double t;

    t.real = exp(z.real) * cos(z.imag);
    t.imag = exp(z.real) * sin(z.imag);
    return t;
}

complex_double 
c_log (complex_double z)
{
    complex_double t;

    t.real = log(c_abs(z));
    t.imag = atan2(z.imag, z.real);
    return t;
}

complex_double 
c_log10 (complex_double z)
{
    complex_double t;

    /* FIXME: correct implementation of c_log10? */
    t.real = log10(c_abs(z));
    t.imag = atan2(z.imag, z.real);
    return t;
}

complex_double 
c_conj (complex_double z)
{
    z.imag = -z.imag;
    return z;
}

complex_double 
c_pow (complex_double x, complex_double y)
{
    return c_exp(c_mult(y, c_log(x)));
}

complex_double 
c_sin (complex_double z)
{
    complex_double t;

    t.real = sin(z.real) * cosh(z.imag);
    t.imag = cos(z.real) * sinh(z.imag);
    return t;
}

complex_double 
c_sqrt (complex_double z)
{
    complex_double t;
    if ((z.real == 0.0) && (z.imag == 0.0)) {
	t.real = 0.0;
	t.imag = 0.0;
    } else {
	double x, y, w, r;
	x = fabs(z.real);
	y = fabs(z.imag);
	if (x >=  y) {
	    r = y/x;
	    w = sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
	} else {
	    r = x/y;
	    w = sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
	}
	if (z.real >= 0.0) {
	    t.real = w;
	    t.imag = z.imag/(2.0*w);
	} else {
	    t.imag = (z.imag >= 0) ? w : -w;
	    t.real = z.imag/(2.0*t.imag);
	}
    }
    return t;
}


/* Not so common functions of complex numbers */
complex_double	
c_acos( complex_double x )
{
    complex_double one;
    complex_double nimag;
    one.real = 1.0;
    one.imag = 0.0;
    nimag.real = 0.0;
    nimag.imag = -1.0;

    return c_mult(nimag, c_log( c_add(x, c_sqrt( c_sub( c_mult(x, x), one)))));
}

complex_double	
c_asin( complex_double x )
{
    complex_double one;
    complex_double nimag;
    complex_double ix;
    one.real = 1.0;
    one.imag = 0.0;
    nimag.real = 0.0;
    nimag.imag = 1.0;
    ix = c_mult(nimag, x);
    nimag.real = 0.0;
    nimag.imag = -1.0;

    return c_mult(nimag, c_log( c_add(ix, c_sqrt( c_sub( c_mult(x, x), one)))));
}

complex_double	
c_atan( complex_double x )
{
    complex_double imag;
    complex_double ihalf;
    imag.real = 0.0;
    imag.imag = 1.0;
    ihalf.real = 0.0;
    ihalf.imag = 0.5;

    return c_mult( ihalf, c_log( c_div( c_add(imag, x), c_sub(imag,x ))));
}

complex_double	
c_atan2( complex_double x, complex_double __y )
{
    complex_double t;
    abort();
    return t;
}

complex_double	
c_cosh( complex_double x )
{
    complex_double two;
    complex_double nx;
    two.real = 2.0;
    two.imag = 0.0;
    nx.real = -x.real;
    nx.imag = -x.imag;

    return c_div( c_add( c_exp(x), c_exp(nx)), two);
}

complex_double	
c_mod( complex_double x, complex_double __y )
{
    complex_double t;
    abort();
    return t;
}

complex_double	
c_sinh( complex_double x )
{
    complex_double two;
    complex_double nx;
    two.real = 2.0;
    two.imag = 0.0;
    nx.real = -x.real;
    nx.imag = -x.imag;

    return c_div( c_sub( c_exp(x), c_exp(nx)), two);
}

complex_double	
c_tan( complex_double x )
{
    return c_div(c_sin(x), c_cos(x));
}

complex_double	
c_tanh( complex_double x )
{
    return c_div(c_sinh(x), c_cosh(x));
}