/*
 * Floating point complex number routines specifically for Mathomatic.
 *
 * Copyright (C) 1987-2007 George Gesslein II.
 */

#include "includes.h"

/*
 * Convert doubles x and y from rectangular coordinates to polar coordinates.
 *
 * The amplitude is stored in *radiusp and the angle in radians is stored in *thetap.
 */
void
rect_to_polar(x, y, radiusp, thetap)
double	x, y, *radiusp, *thetap;
{
	*radiusp = sqrt(x * x + y * y);
	*thetap = atan2(y, x);
}

#if	!LIBRARY
/*
 * The roots command.
 */
int
roots_cmd(cp)
char	*cp;
{
#define	MAX_ROOT	10000.0	/* root limit needed because more roots become more inaccurate and take longer to calculate */

	complexs	c, c2, check;
	double		d, k;
	double		root;
	double		radius, theta;
	double		radius_root = 0.0;
	char		buf[MAX_CMD_LEN];

	if (*cp == '\0') {
		my_strlcpy(prompt_str, _("Enter root (positive integer): "), sizeof(prompt_str));
		if ((cp = get_string(buf, sizeof(buf))) == NULL)
			return false;
	}
	root = strtod(cp, &cp);
	if ((*cp && !isspace(*cp)) || root < 0.0 || root > MAX_ROOT || fmod(root, 1.0) != 0.0) {
		printf(_("Root must be a positive integer less than or equal to %.0f.\n"), MAX_ROOT);
		return false;
	}
	cp = skip_space(cp);
	if (*cp == '\0') {
		my_strlcpy(prompt_str, _("Enter real part (X): "), sizeof(prompt_str));
		if ((cp = get_string(buf, sizeof(buf))) == NULL)
			return false;
	}
	c.re = strtod(cp, &cp);
	if (*cp && !isspace(*cp)) {
		printf(_("Invalid number.\n"));
		return false;
	}
	cp = skip_space(cp);
	if (*cp == '\0') {
		my_strlcpy(prompt_str, _("Enter imaginary part (Y): "), sizeof(prompt_str));
		if ((cp = get_string(buf, sizeof(buf))) == NULL)
			return false;
	}
	c.im = strtod(cp, &cp);
	if (*cp) {
		printf(_("Invalid number.\n"));
		return false;
	}
	if (c.re == 0.0 && c.im == 0.0) {
		return false;
	}
/* convert to polar coordinates */
	errno = 0;
	rect_to_polar(c.re, c.im, &radius, &theta);
	if (root) {
		radius_root = pow(radius, 1.0 / root);
	}
	check_err();
	fprintf(gfp, _("\nThe polar coordinates are:\n%.12g amplitude and %.12g radians (%.12g degrees).\n\n"),
	    radius, theta, theta * 180.0 / M_PI);
	if (root) {
		if (c.im == 0.0) {
			fprintf(gfp, _("The %.12g roots of (%.12g)^(1/%.12g) are:\n\n"), root, c.re, root);
		} else {
			fprintf(gfp, _("The %.12g roots of (%.12g%+.12g*i)^(1/%.12g) are:\n\n"), root, c.re, c.im, root);
		}
		for (k = 0.0; k < root; k += 1.0) {
/* add constants to theta and convert back to rectangular coordinates */
			c2.re = radius_root * cos((theta + 2.0 * k * M_PI) / root);
			c2.im = radius_root * sin((theta + 2.0 * k * M_PI) / root);
			complex_fixup(&c2);
			if (c2.im == 0.0) {
				fprintf(gfp, "%.12g\n", c2.re);
			} else {
				fprintf(gfp, "%.12g %+.12g*i\n", c2.re, c2.im);
			}
			check = c2;
			for (d = 1.0; d < root; d += 1.0) {
				check = complex_mult(check, c2);
			}
			complex_fixup(&check);
			if (check.im == 0.0) {
				printf(_("Inverse check: %.12g\n\n"), check.re);
			} else {
				printf(_("Inverse check: %.12g %+.12g*i\n\n"), check.re, check.im);
			}
		}
	}
	return true;
}
#endif

/*
 * Approximate roots of complex numbers:
 * (complex^real) and (real^complex) and (complex^complex).
 * Only gives one root, when there may be many.
 *
 * Return true if expression was modified.
 */
int
complex_root_simp(equation, np)
token_type	*equation;	/* equation side pointer */
int		*np;		/* pointer to length of equation side */
{
	int		i, j;
	int		level;
	int		len;
	complexs	c, p;
	int		modified = false;

start_over:
	for (i = 1; i < *np; i += 2) {
		if (equation[i].token.operatr != POWER)
			continue;
		level = equation[i].level;
		for (j = i + 2; j < *np && equation[j].level >= level; j += 2)
			;
		len = j - (i + 1);
		if (!parse_complex(&equation[i+1], len, &p))
			continue;
		for (j = i - 1; j >= 0 && equation[j].level >= level; j--)
			;
		j++;
		if (!parse_complex(&equation[j], i - j, &c))
			continue;
		if (c.im == 0.0 && p.im == 0.0)
			continue;
		i += len + 1;
		c = complex_pow(c, p);
		if (*np + 5 - (i - j) > n_tokens) {
			error_huge();
		}
		blt(&equation[j+5], &equation[i], (*np - i) * sizeof(token_type));
		*np += 5 - (i - j);
		equation[j].level = level;
		equation[j].kind = CONSTANT;
		equation[j].token.constant = c.re;
		j++;
		equation[j].level = level;
		equation[j].kind = OPERATOR;
		equation[j].token.operatr = PLUS;
		j++;
		level++;
		equation[j].level = level;
		equation[j].kind = CONSTANT;
		equation[j].token.constant = c.im;
		j++;
		equation[j].level = level;
		equation[j].kind = OPERATOR;
		equation[j].token.operatr = TIMES;
		j++;
		equation[j].level = level;
		equation[j].kind = VARIABLE;
		equation[j].token.variable = IMAGINARY;
		modified = true;
		goto start_over;
	}
	if (modified) {
		debug_string(0, "Warning: complex number root approximated, solutions may be lost or inaccurate.");
	}
	return modified;
}

/*
 * Get a constant, if the passed expression is a constant.
 * This code only handles the case of a single constant, where "n" == 1.
 * I don't know how to easily fix this, should at least handle "n" == 3,
 * which would allow 2 constants separated by any operator.
 *
 * Return true if successful, with number in "*dp".
 */
int
get_constant(p1, n, dp)
token_type	*p1;	/* expression pointer */
int		n;	/* length of expression */
double		*dp;	/* pointer to returned double */
{
	if (n != 1)
		return false;
	switch (p1[0].kind) {
	case CONSTANT:
		*dp = p1[0].token.constant;
		return true;
	case VARIABLE:
		if (var_is_const(p1[0].token.variable, dp)) {
			return true;
		}
	}
	return false;
}

/*
 * Parse a complex number expression.
 *
 * If successful return true with complex number in "*cp".
 */
int
parse_complex(p1, n, cp)
token_type	*p1;	/* expression pointer */
int		n;	/* length of expression */
complexs	*cp;	/* pointer to returned complex number */
{
	int		j;
	int		imag_cnt = 0, plus_cnt = 0, times_cnt = 0;
	complexs	c;
	int		level2;
	double		junk;

	if (get_constant(p1, n, &c.re)) {
		c.im = 0.0;
		*cp = c;
		return true;
	}
	c.re = 0.0;
	c.im = 1.0;
	for (j = n - 1; j >= 0; j--) {
		switch (p1[j].kind) {
		case CONSTANT:
			break;
		case VARIABLE:
			if (var_is_const(p1[j].token.variable, &junk))
				break;
			if (p1[j].token.variable != IMAGINARY)
				return false;
			imag_cnt++;
			break;
		case OPERATOR:
			level2 = p1[j].level;
			switch (p1[j].token.operatr) {
			case TIMES:
				if (++times_cnt > 1)
					return false;
				if (p1[j-1].level != level2 || p1[j+1].level != level2)
					return false;
				if (p1[j-1].kind == VARIABLE && p1[j-1].token.variable == IMAGINARY) {
					if (!get_constant(&p1[j+1], 1, &c.im))
						return false;
					continue;
				}
				if (p1[j+1].kind == VARIABLE && p1[j+1].token.variable == IMAGINARY) {
					if (!get_constant(&p1[j-1], 1, &c.im))
						return false;
					continue;
				}
				return false;
			case DIVIDE:
				if (++times_cnt > 1)
					return false;
				if (p1[j-1].level != level2 || p1[j+1].level != level2)
					return false;
				if (p1[j-1].kind == VARIABLE && p1[j-1].token.variable == IMAGINARY) {
					if (!get_constant(&p1[j+1], 1, &c.im))
						return false;
					c.im = 1.0 / c.im;
					continue;
				}
				return false;
			case MINUS:
				if (imag_cnt) {
					c.im = -c.im;
				}
			case PLUS:
				if (++plus_cnt > 1)
					return false;
				if (p1[j-1].level == level2 && get_constant(&p1[j-1], 1, &c.re)) {
					continue;
				}
				if (p1[j+1].level == level2 && get_constant(&p1[j+1], 1, &c.re)) {
					if (p1[j].token.operatr == MINUS)
						c.re = -c.re;
					continue;
				}
			}
		default:
			return false;
		}
	}
	if (imag_cnt != 1)
		return false;
	*cp = c;
	return true;
}