/*
 *			C O M P L E X . C
 *
 *  Functions -
 *	CxDiv		Complex Division
 *	CxSqrt		Complex Square Root
 *	rt_pr_roots	Print complex roots
 *
 *  Authors -
 *	Douglas A Gwyn		(Original Version)
 *	Michael John Muuss	(Macro Version)
 *  
 *  Source -
 *	SECAD/VLD Computing Consortium, Bldg 394
 *	The U. S. Army Ballistic Research Laboratory
 *	Aberdeen Proving Ground, Maryland  21005
 *  
 *  Copyright Notice -
 *	This software is Copyright (C) 1985 by the United States Army.
 *	All rights reserved.
 */
#ifndef lint
static char RCSid[] = "@(#)$Header: complex.c,v 3.0 86/06/10 01:32:49 mike Exp $ (BRL)";
#endif

#include <stdio.h>
#include <math.h>
#include "../h/machine.h"
#include "../h/vmath.h"
#include "../h/raytrace.h"
#include "complex.h"

/* arbitrary numerical arguments, integer value: */
#define	Sgn( x )	((x) == 0 ? 0 : (x) > 0 ? 1 : -1)
#define	Abs( a )	((a) >= 0 ? (a) : -(a))

/*
 *			C X D I V
 *
 *	Divide one complex by another
 *
 *	CxDiv( &a, &b )	divides  a  by  b  and returns  &a;
 *			zero divisor fails
 */
void
CxDiv( ap, bp )
register complex	*ap, *bp;	/* may coincide (?) */
{
	LOCAL fastf_t	r, s;
	LOCAL fastf_t	ap__re;

	/* Note: classical formula may cause unnecessary overflow */
	ap__re = ap->re;
	r = bp->re;
	s = bp->im;
	if ( Abs( r ) >= Abs( s ) )  {
		r = s / r;			/* <= 1 */
		s = 1.0 / (bp->re + r * s);
		ap->re = (ap->re + ap->im * r) * s;
		ap->im = (ap->im - ap__re * r) * s;
	}  else  {
		/* Abs( s ) > Abs( r ) */
		r = r / s;			/* < 1 */
		s = 1.0 / (s + r * bp->re);
		ap->re = (ap->re * r + ap->im) * s;
		ap->im = (ap->im * r - ap__re) * s;
	}
}

/*
 *			C X S Q R T
 *
 *  Compute square root of complex number
 *
 *	CxSqrt( &c )	replaces  c  by  sqrt(c)  and returns  &c
 *
 *	Note:	This is a double-valued function; the result of
 *		CxSqrt() always has nonnegative imaginary part.
 */
complex *
CxSqrt( cp )
register complex	*cp;
{
	/* record signs of original real & imaginary parts */
	LOCAL int	re_sign;
	LOCAL int	im_sign;

	/* special cases are not necessary; they are here for speed */
	im_sign = Sgn( cp->im );
	if( (re_sign = Sgn(cp->re))  == 0 )  {
		if ( im_sign == 0 )
			;		/* (0,0) already there */
		else if ( im_sign > 0 )
			cp->re = cp->im = sqrt( cp->im * 0.5 );
		else			/* im_sign < 0 */
			cp->re = -(cp->im = sqrt( cp->im * -0.5 ));
	} else if ( im_sign == 0 )  {
		if ( re_sign > 0 )
			cp->re = sqrt( cp->re );
/*			cp->im = 0.0;	/* 0 already there */
		else  {		/* re_sign < 0 */
			cp->im = sqrt( -cp->re );
			cp->re = 0.0;
		}
	}  else  {			/* no shortcuts */
		FAST fastf_t	ampl, temp;

		ampl = CxAmpl( cp );
		if( (temp = (ampl - cp->re) * 0.5) < 0.0 )  {
			/* This case happens rather often, when the
			 *  hypot() in CxAmpl() returns an ampl ever so
			 *  slightly smaller than cp->re.  This can easily
			 *  happen when cp->re ~= 10**20.
			 *  Just ignore the imaginary part.
			 */
			cp->im = 0;
		} else
			cp->im = sqrt( temp );

		if( (temp = (ampl + cp->re) * 0.5) < 0.0 )  {
			cp->re = 0.0;
		} else {
			if( im_sign > 0 )
				cp->re = sqrt(temp);
			else			/* im_sign < 0 */
				cp->re = -sqrt(temp);
		}
	}
	return cp;
}

/*
 *			R T _ P R _ R O O T S
 */
void
rt_pr_roots( nroots, roots )
register int 		nroots;
register complex	roots[];
{
	register int	m;

	rt_log("Roots are:\n");
	for ( m=0; m < nroots; ++m ) {
		rt_log("\t%g + i%g\n",roots[m].re,roots[m].im);
	}
}
