/*							qcbrt.c		*/

/* cube root */
/* 21 December 83, SLM */

#include "qhead.h"

#ifndef ANSIPROT
int qfrexp(), qldexp();
#endif

extern QELT qzero[], qone[], qtwo[], qthree[];

#define CBRT2 1.25992104989487316477
#define CBRT4 1.58740105196819947475
#define CBRT2I 0.793700525984099737376
#define CBRT4I 0.629960524947436582384

/* ??? For some reason, VAX C miscompiles this file if these
   variables are made local.  */
static QELT a[NQ], b[NQ], c[NQ];

int qcbrt( xx, y )
QELT *xx, *y;
{
int sign;
long e, rem;
double x, z;


if( xx[1] == 0 )
	{
	qclear( y );
	return( 0 );
	}

sign = xx[0];
qmov( xx, c );
c[0] = 0;

/* extract power of 2, leaving
 * mantissa between 0.5 and 1
 */
qfrexp( c, &e, a );

qtoe( a, (unsigned short *) &z );
/* Approximate cube root of number between .5 and 1,
 * peak relative error = 9.2e-6
 */
x = (((-1.3466110473359520655053e-1  * z
      + 5.4664601366395524503440e-1) * z
      - 9.5438224771509446525043e-1) * z
      + 1.1399983354717293273738e0 ) * z
      + 4.0238979564544752126924e-1;
/* Newton iterations. */
x -= ( x - (z/(x*x)) )*0.33333333333333333333;
x -= ( x - (z/(x*x)) )*0.33333333333333333333;

/* exponent divided by 3 */
if( e >= 0 )
	{
	rem = e;
	e = e / 3;
	rem = rem - 3 * e;
	if( rem == 1 )
		x = x * CBRT2;
	else if( rem == 2 )
		x = x * CBRT4;
	}
else
	{
	/* argument less than 1 */
	e = -e;
	rem = e;
	e /= 3;
	rem -= 3*e;
	if( rem == 1 )
		x *= CBRT2I;
	else if( rem == 2 )
		x *= CBRT4I;
	e = -e;
	}
etoq( (unsigned short *) &x, a );

/* multiply by power of 2 */
qldexp( a, e, a );

/* More Newton iterations. */
qmul( a, a, b );	/* 	x * x		*/
qdiv( b, c, b );	/*	z / x*x	*/
qsub( b, a, b );	/*	x - z/x*x	*/
qdiv( qthree, b, b );	/*	.../3		*/
qsub( b, a, a );	/*	x - ...		*/

qmul( a, a, b );
qdiv( b, c, b );
qsub( b, a, b );
qdiv( qthree, b, b );
qsub( b, a, a );

qmul( a, a, b );
qdiv( b, c, b );
qsub( b, a, b );
qdiv( qthree, b, b );
qsub( b, a, a );


if( sign != 0 )
	a[0] = -1;
qmov( a, y );
return( 0 );
}