/*
*	icalc - complex-expression parser
*
*	Math routines for complex-number expression parser.
*	In the routines, rv is variable holding 'return value'.
*
*	(C) Martin W Scott, 1991.
*/
#include <stdio.h>
#include <math.h>
#include "complex.h"
#include "constant.h"

#define NOERRCHECK	/* domain checking doesn't work when using */
			/* IEEE libraries... */
#ifdef NOERRCHECK
#define Log(x) log(x)
#endif

extern Symbol *ans;

static int prec = 8;		/* number of decimal places to print */
				/* (when they exist) */

void cprin(fp, prefix, suffix, z)	/* print a complex number to file fp */
	FILE *fp;
	char *prefix, *suffix;
	Complex z;
{
	fprintf(fp, prefix);

	if (z.imag == 0.0)
		fprintf(fp, "%10.*lg", prec, z.real);
	else if (z.real == 0.0)
		fprintf(fp, "%10.*lg i", prec, z.imag);
	else
		fprintf(fp, "%.*g %c %.*g i",
			prec, z.real, sign(z.imag), prec, abs(z.imag));

	fprintf(fp, suffix);
}

Complex printres(z)	/* prints and returns its argument */
	Complex z;
{
	cprin(stdout, "\t", "\n", z);
	return z;
}

Complex precision(z)	/* adjust decimal places shown */
	Complex z;
{
	prec = (int)z.real;
	return ans->u.val;
}

static double integer(x)	/* round x to NEAREST integer */
	double x;		/* if fractional part is 0.5, round UP */
{
	double f,i;

	f = modf(x,&i);
	if (abs(f) >= 0.5)
		if (x >= 0.0)
			return i + 1.0;
		else
			return i - 1.0;
	else return i;
}

Complex cinteger(z)	/* round both real and imag parts to nearest integer */
	Complex z;
{
	z.real = integer(z.real);
	z.imag = integer(z.imag);

	return z;
}
	
Complex cceil(z)	/* ceiling of real and imaginary parts of z */
	Complex z;
{
	z.real = ceil(z.real);
	z.imag = ceil(z.imag);
	return z;
}
	
Complex cfloor(z)	/* ceiling of real and imaginary parts of z */
	Complex z;
{
	z.real = floor(z.real);
	z.imag = floor(z.imag);
	return z;
}

Complex Re(z)		/* real part of complex number */
	Complex z;
{
	z.imag = 0.0;
	return z;
}

Complex Im(z)		/* imaginary part of complex number */
	Complex z;
{
	z.real = z.imag;
	z.imag = 0.0;
	return z;
}

#define marg(z) atan2((z).imag,(z).real)	/* macro arg */
#define rnorm(z) (z.real*z.real+z.imag*z.imag)	/* real norm */

Complex arg(z)		/* argument of complex number, in range (-PI,PI] */
	Complex z;
{
	z.real = atan2(z.imag, z.real);
	z.imag = 0.0;
	return z;
}

Complex norm(z)		/* norm of a complex number */
	Complex z;
{
	z.real = sqr(z.real) + sqr(z.imag);
	z.imag = 0.0;
	return z;
}

Complex cabs(z)		/* absolute value of a complex number */
	Complex z;
{
	z.real = sqrt(rnorm(z));
	z.imag = 0.0;
	return z;
}

Complex cadd(w,z)	/* complex addition */
	Complex w,z;
{
	w.real = w.real + z.real;
	w.imag = w.imag + z.imag;

	return w;
}

Complex csub(w,z)	/* complex subtraction */
	Complex w,z;
{
	w.real = w.real - z.real;
	w.imag = w.imag - z.imag;

	return w;
}

Complex cmul(w,z)	/* complex multiplication */
	Complex w,z;
{
	Complex rv;

	rv.real = w.real*z.real - w.imag*z.imag;
	rv.imag = w.real*z.imag + w.imag*z.real;

	return rv;
}

Complex cdiv(w,z)	/* complex division */
	Complex w,z;
{
	Complex rv;
	double temp = sqr(z.real)+sqr(z.imag);

	if (temp == 0.0)
		execerror("division by zero", NULL);	

	rv.real = (w.real*z.real + w.imag*z.imag)/temp;
	rv.imag = (w.imag*z.real - w.real*z.imag)/temp;

	return rv;
}

Complex cneg(z)		/* complex negation */
	Complex z;
{
	z.real = -z.real;
	z.imag = -z.imag;

	return z;
}

Complex csqr(z)		/* complex square, w^2 */
	Complex z;
{
	Complex rv;

	if (z.imag == 0.0)	/* small optimisation for reals */
	{
		z.real *= z.real;
		return z;
	}
	rv.real = sqr(z.real) - sqr(z.imag);
	rv.imag = 2*z.real*z.imag;

	return rv;
}

Complex csqrt(z)	/* complex square-root */
	Complex z;
{
	Complex rv;
	double temp = sqrt(rnorm(z));

	rv.real = sqrt((temp + z.real)*0.5);
	rv.imag = sqrt((temp - z.real)*0.5);

	return rv;
}

Complex conj(z)		/* conjugate of w */
	Complex z;
{
	z.imag = -z.imag;

	return z;
}

Complex cexp(z)		/* complex exponential function e^z */
	Complex z;
{
	double temp = exp(z.real);

	if (z.imag == 0.0)	/* small optimisation for reals */
		z.real = temp;
	else
	{
		z.real = temp*cos(z.imag);
		z.imag = temp*sin(z.imag);
	}
	return z;
}

Complex clog(z)		/* complex natural logarithm */
	Complex z;
{
	Complex rv;

	rv.real = Log(rnorm(z))*0.5;
	rv.imag = marg(z);

	return rv;
}

Complex cpow(w,z)	/* complex exponential, w^z */
	Complex w,z;
{
	return cexp(cmul(z,clog(w)));
}

Complex csin(z)		/* complex sine */
	Complex z;
{
	if (z.imag == 0.0)	/* small optimisation for reals */
	{
		z.real = sin(z.real);
		return z;
	}
	else
	{
		Complex rv;

		rv.real = sin(z.real)*cosh(z.imag);
		rv.imag = cos(z.real)*sinh(z.imag);

		return rv;
	}
}

Complex ccos(z)		/* complex cosine */
	Complex z;
{
	if (z.imag == 0.0)	/* small optimisation for reals */
	{
		z.real = cos(z.real);
		return z;
	}
	else
	{
		Complex rv;

		rv.real = cos(z.real)*cosh(z.imag);
		rv.imag = -(sin(z.real)*sinh(z.imag));

		return rv;
	}
}

Complex ctan(z)		/* complex tangent */
	Complex z;
{
	if (z.imag == 0.0)	/* small optimisation for reals */
	{
		z.real = tan(z.real);
		return z;
	}
	else
		return cdiv(csin(z),ccos(z));
}

Complex casin(z)	/* complex arcsine - lazy version */
	Complex z;
{
	/* asin z = -ilog(iz + sqrt(1-z^2)) */

	if (abs(z.real) <= 1.0 && z.imag == 0.0)
	{
		z.real = asin(z.real);
		return z;
	}
	else
		return cmul(minuseye,clog(cadd(cmul(eye,z),csqrt(csub(one,csqr(z))))));
}

Complex cacos(z)	/* complex arccosine - lazy version */
	Complex z;
{
	/* acos z = -ilog(z + sqrt(z^2-1)) */

	if (abs(z.real) <= 1.0 && z.imag == 0.0)
	{
		z.real = acos(z.real);
		return z;
	}
	else
		return cmul(minuseye,clog(cadd(z,csqrt(csub(csqr(z),one)))));
}

Complex catan(z)	/* complex arctangent - not so lazy version */
	Complex z;
{
	if (z.imag == 0.0)	/* small optimisation for reals */
		z.real = atan(z.real);
	else
	{
		Complex ctemp;
		double temp = rnorm(z);

		ctemp.real = ctemp.imag = 1.0 / (1.0 + temp - 2.0 * z.imag);
		ctemp.real *= 1.0 - temp;
		ctemp.imag *= -2.0*z.real;
		ctemp = clog(ctemp);
		z.real = -0.5*ctemp.imag;
		z.imag = 0.5*ctemp.real;
	}

	return z;
}

Complex csinh(z)	/* complex hyperbolic sine */
	Complex z;
{
	Complex rv;

	rv.real = cos(z.imag)*sinh(z.real);
	rv.imag = sin(z.imag)*cosh(z.real);

	return rv;
}

Complex ccosh(z)		/* complex hyperbolic cosine */
	Complex z;
{
	Complex rv;

	rv.real = cos(z.imag)*cosh(z.real);
	rv.imag = sin(z.imag)*sinh(z.real);

	return rv;
}

Complex ctanh(z)		/* complex tangent */
	Complex z;
{
	return cdiv(csinh(z),ccosh(z));
}