/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:22 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "spolz2.h"
#include <float.h>
#include <stdlib.h>
		/* PARAMETER translations */
#define	C16	16.e0
#define	HALF	.5e0
#define	ONE	1.e0
#define	ZERO	0.e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ spolz2(
float a[],
float z[])
{
	long int _l0;
	float aq, au, d, f, p, q, r1, r2, r3, r4, u, v, w;
	static float c1, c2;
	static LOGICAL32 first = TRUE;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const A = &a[0] - 1;
	float *const Z = &z[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1996-04-27 SPOLZ2 Krogh   Changes to use .C. and C%%.
	 *>> 1996-03-30 SPOLZ2 Krogh   Added external statement.
	 *>> 1996-01-18 SPOLZ2 Krogh   Added code for M77CON for conversion to C.
	 *>> 1992-03-13 SPOLZ2 FTK     Removed implicit statements.
	 *>> 1987-02-25 SPOLZ2 Lawson  Initial code.
	 *--S Replaces "?": ?POLZ2
	 *++ Default NO_COMPLEX = .C. | (.N. == 'D')
	 *++ Default COMPLEX = ~NO_COMPLEX
	 *
	 *     ------------------------------------------------------------------
	 *
	 *     Find the two roots of the quadratic polynomial
	 *               A(1)*X*X + A(2)*X + A(3)
	 *     Return the roots as the complex numbers (Z(1), Z(2)) and
	 *     (Z(3), Z(4)) if Z is not complex, else return the two complex
	 *     numbers, Z(1) and Z(2).
	 *     Require A(1) .ne. 0.
	 *
	 *     Method:
	 *     Divide through by A(1). New polynomial is
	 *            X*X + P*X + Q
	 *     Let U = -P/2
	 *     Roots are U + SQRT(U*U-Q) and U - SQRT(U*U-Q).
	 *     Avoid computing U*U explicity if it might overflow or
	 *     underflow. In case of real roots, U + Z and U - Z, the first
	 *     root is  U + SIGN(U)*ABS(Z) and the second is Q / Z(1).
	 *
	 *     C. L. Lawson & S. Chiu, JPL, 1987 Feb 16.
	 *     ------------------------------------------------------------------
	 * */
	/*++ CODE for COMPLEX is inactive
	 *      COMPLEX Z(2), R1, R2
	 *++ CODE for NO_COMPLEX is active */
	/*++ END
	 * */
 
	/*     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
	if (first)
	{
		first = FALSE;
 
		/*     R1MACH(1) is the underflow limit.
		 * */
		c1 = sqrtf( FLT_MIN*C16 );
 
		/*     R1MACH(2) is the overflow limit.
		 * */
		c2 = sqrtf( FLT_MAX/C16 );
	}
 
	if (A[1] == ZERO)
	{
		ermsg( "SPOLZ2", 1, 0, "A(1) .EQ. 0.", '.' );
		r1 = ZERO;
		r2 = ZERO;
		/*++ CODE for NO_COMPLEX is active */
		r3 = ZERO;
		r4 = ZERO;
		/*++ END */
		goto L_99;
	}
	else
	{
		p = A[2]/A[1];
		q = A[3]/A[1];
	}
 
	if (q == ZERO)
	{
		r1 = ZERO;
		/*++ CODE for COMPLEX is inactive
		 *        R2 = -P
		 *++ CODE for NO_COMPLEX is active */
		r2 = ZERO;
		r3 = -p;
		r4 = ZERO;
		/*++ END */
		goto L_99;
	}
 
	if (p == ZERO)
	{
		w = sqrtf( fabsf( q ) );
		if (q > ZERO)
		{
			/*++ CODE for COMPLEX is inactive
			 *          R1 = CMPLX(ZERO,W)
			 *          R2 = CMPLX(ZERO,-W)
			 *++ CODE for NO_COMPLEX is active */
			r1 = ZERO;
			r2 = w;
			r3 = ZERO;
			r4 = -w;
			/*++ END */
		}
		else
		{
			/*++ CODE for COMPLEX is inactive
			 *          R1 = CMPLX(W,ZERO)
			 *          R2 = CMPLX(-W,ZERO)
			 *++ CODE for NO_COMPLEX is active */
			r1 = w;
			r2 = ZERO;
			r3 = -w;
			r4 = ZERO;
			/*++ END */
		}
		goto L_99;
	}
	u = -p*HALF;
 
	/*     Compute D having the sign of U*U-Q
	 *         and F = SQRT(ABS(U*U-Q))
	 * */
	au = fabsf( u );
	if (au > c2)
	{
		d = ONE - (q/u)/u;
		f = au*sqrtf( fabsf( d ) );
	}
	else if (au < c1)
	{
		aq = fabsf( q );
		d = u*(u/aq) - signf( ONE, q );
		f = sqrtf( aq )*sqrtf( fabsf( d ) );
	}
	else
	{
		d = u*u - q;
		f = sqrtf( fabsf( d ) );
	}
 
	if (d == ZERO)
	{
		/*++ CODE for COMPLEX is inactive
		 *        R1 = CMPLX(U,ZERO)
		 *        R2 = R1
		 *++ CODE for NO_COMPLEX is active */
		r1 = u;
		r2 = ZERO;
		r3 = r1;
		r4 = ZERO;
		/*++ END */
	}
	else
	{
		if (d > ZERO)
		{
			if (u > ZERO)
			{
				v = u + f;
			}
			else
			{
				v = u - f;
			}
			/*++ CODE for COMPLEX is inactive
			 *          R1 = CMPLX(V,ZERO)
			 *          R2 = CMPLX(Q/V,ZERO)
			 *++ CODE for NO_COMPLEX is active */
			r1 = v;
			r2 = ZERO;
			r3 = q/v;
			r4 = ZERO;
			/*++ END */
		}
		else
		{
			/*++ CODE for COMPLEX is inactive
			 *          R1 = CMPLX(U,F)
			 *          R2 = CMPLX(U,-F)
			 *++ CODE for NO_COMPLEX is active */
			r1 = u;
			r2 = f;
			r3 = u;
			r4 = -f;
			/*++ END */
		}
	}
 
L_99:
	;
	Z[1] = r1;
	Z[2] = r2;
	/*++ CODE for NO_COMPLEX is active */
	Z[3] = r3;
	Z[4] = r4;
	/*++ END */
	return;
} /* end of function */