/*
 * Copyright 1995,96 Thierry Bousch
 * Licensed under the Gnu Public License, Version 2
 *
 * $Id: Complex.c,v 2.4 1996/08/18 09:05:56 bousch Exp $
 *
 * Complex numbers, more generally numbers of the form A+iB where A and B
 * belong to some ring or field.
 */

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "saml.h"
#include "saml-errno.h"
#include "mnode.h"
#include "builtin.h"

typedef struct {
	struct mnode_header hdr;
	s_mnode *re;		/* real part */
	s_mnode *im;		/* imaginary part */
} complex_mnode;

static s_mnode* complex_build (const char*);
static gr_string* complex_stringify (complex_mnode*);
static s_mnode* complex_make (s_mnode*);
static s_mnode* complex_add (complex_mnode*, complex_mnode*);
static s_mnode* complex_sub (complex_mnode*, complex_mnode*);
static s_mnode* complex_mul (complex_mnode*, complex_mnode*);
static s_mnode* complex_div (complex_mnode*, complex_mnode*);
static int complex_notzero (complex_mnode*);
static s_mnode* complex_zero (complex_mnode*);
static s_mnode* complex_negate (complex_mnode*);
static s_mnode* complex_one (complex_mnode*);
static s_mnode* complex_invert (complex_mnode*);
static void complex_free (complex_mnode*);
static s_mnode* literal2complex (s_mnode*, complex_mnode*);

static unsafe_s_mtype MathType_Complex = {
	"Complex",
	complex_free, complex_build, complex_stringify,
	complex_make, NULL,
	complex_add, complex_sub, complex_mul, complex_div, NULL,
	complex_notzero, NULL, NULL, mn_std_differ, NULL,
	complex_zero, complex_negate, complex_one, complex_invert, NULL
};

void init_MathType_Complex (void)
{
	register_mtype(ST_COMPLEX, &MathType_Complex);
	register_CV_routine(ST_LITERAL, ST_COMPLEX, literal2complex);
}

static inline complex_mnode* complex_new (void)
{
	return (complex_mnode*)__mnalloc(ST_COMPLEX,sizeof(complex_mnode));
}

static void complex_free (complex_mnode *cplx)
{
	unlink_mnode(cplx->re);
	unlink_mnode(cplx->im);
	free(cplx);
}

static s_mnode* complex_build (const char *str)
{
	int x, y;
	complex_mnode *cplx;
	char buff[24];

	if (sscanf(str, "(%d,%d)", &x, &y) == 2) {
		cplx = complex_new();
		sprintf(buff, "%d", x);
		cplx->re = mnode_build(ST_INTEGER, buff);
		sprintf(buff, "%d", y);
		cplx->re = mnode_build(ST_INTEGER, buff);
		return (s_mnode*) cplx;
	}
	return mnode_error(SE_STRING, "complex_build");
}

static s_mnode* complex_zero (complex_mnode* model)
{
	complex_mnode *cpl_zero = complex_new();
	s_mnode *zero;

	cpl_zero->re = zero = mnode_zero(model->re);
	cpl_zero->im = copy_mnode(zero);
	return (s_mnode*) cpl_zero;
}

static s_mnode* complex_one (complex_mnode* model)
{
	complex_mnode *cpl_one = complex_new();

	cpl_one->re = mnode_zero(model->re);
	cpl_one->im = mnode_one (model->im);
	return (s_mnode*) cpl_one;
}

static s_mnode* literal2complex (s_mnode* lit, complex_mnode* model)
{
	complex_mnode *cplx;
	s_mnode *rem;

	if (!model)
		return mnode_error(SE_ICAST, "literal2complex");
	cplx = complex_new();
	rem = model->re;
	/*
	 * Every literal is converted into sqrt(-1). In particular, "I"
	 * is converted to what you expect...
	 */
	cplx->re = mnode_zero(rem);
	cplx->im = mnode_one(rem);
	return (mn_ptr)cplx;
}

static s_mnode* complex_make (s_mnode* realpart)
{
	complex_mnode* cplx = complex_new();

	cplx->re = copy_mnode(realpart);
	cplx->im = mnode_zero(realpart);
	return (mn_ptr)cplx;
}

static gr_string* complex_stringify (complex_mnode* cplx)
{
	gr_string *gs0, *gs1;

	gs0 = mnode_stringify(cplx->re);
	if (mnode_notzero(cplx->im) == 0)
		return gs0;
	gs1 = mnode_stringify(cplx->im);
	gs0 = grs_prepend1(gs0, '(');
	gs0 = grs_append1(gs0, ',');
	gs0 = grs_append(gs0, gs1->s, gs1->len);
	free(gs1);
	gs0 = grs_append1(gs0, ')');
	return gs0;
}

static s_mnode* complex_add (complex_mnode* c1, complex_mnode* c2)
{
	complex_mnode *sum = complex_new();
	sum->re = mnode_add(c1->re, c2->re);
	sum->im = mnode_add(c1->im, c2->im);
	return (s_mnode*) sum;
}

static s_mnode* complex_sub (complex_mnode* c1, complex_mnode* c2)
{
	complex_mnode* diff = complex_new();
	diff->re = mnode_sub(c1->re, c2->re);
	diff->im = mnode_sub(c1->im, c2->im);
	return (s_mnode*) diff;
}

static s_mnode* complex_mul (complex_mnode* c1, complex_mnode* c2)
{
	complex_mnode* prod = complex_new();
	s_mnode *ta, *tb;

	ta = mnode_mul(c1->re, c2->re);
	tb = mnode_mul(c1->im, c2->im);
	prod->re = mnode_sub(ta, tb);
	unlink_mnode(ta); unlink_mnode(tb);

	ta = mnode_mul(c1->re, c2->im);
	tb = mnode_mul(c1->im, c2->re);
	prod->im = mnode_add(ta, tb);
	unlink_mnode(ta); unlink_mnode(tb);

	return (s_mnode*) prod;
}

static int complex_notzero (complex_mnode* c)
{
	return mnode_notzero(c->re) || mnode_notzero(c->im);
}

static s_mnode* complex_negate (complex_mnode* c)
{
	complex_mnode *oppo = complex_new();
	oppo->re = mnode_negate(c->re);
	oppo->im = mnode_negate(c->im);
	return (s_mnode*) oppo;
}

static s_mnode* complex_norm (complex_mnode* c)
{
	s_mnode *norm, *ta, *tb;

	ta = mnode_mul(c->re, c->re);
	tb = mnode_mul(c->im, c->im);
	norm = mnode_add(ta, tb);
	unlink_mnode(ta); unlink_mnode(tb);
	return norm;
}

static s_mnode* complex_div (complex_mnode* c1, complex_mnode* c2)
{
	complex_mnode* quot;
	s_mnode *ta, *tb, *tc, *norm;

	norm = complex_norm(c2);
	if (!mnode_notzero(norm)) {
		unlink_mnode(norm);
		return mnode_error(SE_DIVZERO, "complex_div");
	}
	quot = complex_new();

	ta = mnode_mul(c1->re, c2->re);
	tb = mnode_mul(c1->im, c2->im);
	tc = mnode_add(ta, tb);
	unlink_mnode(ta); unlink_mnode(tb);
	quot->re = mnode_div(tc, norm);
	unlink_mnode(tc);

	ta = mnode_mul(c1->re, c2->im);
	tb = mnode_mul(c1->im, c2->re);
	tc = mnode_sub(ta, tb);
	unlink_mnode(ta); unlink_mnode(tb);
	quot->im = mnode_div(tc, norm);
	unlink_mnode(tc);

	return (s_mnode*) quot;
}

static s_mnode* complex_invert (complex_mnode* c)
{
	return mnode_error(SE_NOTRDY, "complex_invert");
}