Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : QPOW.C

/* qpow.c */
/* power function: z = x**y */

#include "qhead.h"
#include "mconf.h"

extern short qone[];


qpow( x, y, z )
short *x, *y, *z;
{
short w[NQ];
long li;

qfloor( y, w );
if( qcmp(y,w) == 0 )
{
qifrac( y, &li, w );
if( li < 0 )
li = -li;
if( li < 0x7fffffff )
{
qpowi( x, y, z );
return;
}
}
/* z = exp( y * log(x) ) */

qlog( x, w );
qmul( y, w, w );
qexp( w, z );
}


/* y is integer valued. */

qpowi( x, y, z )
short x[], y[], z[];
{
short w[NQ];
long li, lx;
int signx, signy;

qifrac( y, &li, w );
if( li < 0 )
lx = -li;
else
lx = li;

if( lx == 0x7fffffff )
{
qpow( x, y, z );
return;
}


if( x[1] == 0.0 )
{
if( li == 0 )
{
qmov( qone, z );
return;
}
else if( li < 0 )
{
qinfin( z );
return;
}
else
{
qclear( z );
return;
}
}

if( li == 0L )
{
qmov( qone, z );
return;
}

qmov( x, w );
signx = w[0];
w[0] = 0;


if( li < 0 )
{
li = -li;
signy = -1;
}
else
signy = 0;


/* First bit of the power */
if( li & 1 )
{
qmov( w, z );
}
else
{
qmov( qone, z );
signx = 0;
}

/* Overflow detection */

lx = signy * li * (w[1] - 0x4001);
if( lx > 16383L )
{
qinfin( z );
mtherr( "qpowi", OVERFLOW );
goto done;
}
if( lx < -16383L )
{
qclear( z );
return;
}


li >>= 1;
while( li != 0L )
{
qmul( w, w, w ); /* arg to the 2-to-the-kth power */
if( li & 1L ) /* if that bit is set, then include in product */
qmul( w, z, z );
li >>= 1;
}


done:

if( signx )
qneg( z ); /* odd power of negative number */

if( signy )
{
if( z[1] != 0 )
{
qdiv( z, qone, z );
}
else
{
qinfin( z );
mtherr( "qpowi", OVERFLOW );
}
}
}