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

/* qlog.c */
/* natural logarithm */

#include "qhead.h"

extern short qone[], qtwo[], qlog2[], qsqrt2[];

qlog( x, y )
short *x, *y;
{
short xx[NQ], z[NQ], a[NQ], b[NQ], t[NQ], qj[NQ];
long i, j, n, nsq;
long ex;

if( x[0] < 0 )
{
qclear(y);
printf( "qlog domain error\n" );
return;
}
if( x[1] == 0 )
{
qinfin( y );
y[0] = -1;
printf( "qlog singularity\n" );
return(0);
}

/* range reduction: log x = log( 2**ex * m ) = ex * log2 + log m */
qmov(x, xx );
ex = *(xx+1);
ex -= 040000;
xx[1] = 040000;

/* Adjust range to 1/sqrt(2), sqrt(2) */
qsqrt2[1] -= 1;
if( qcmp( xx, qsqrt2 ) < 0 )
{
ex -= 1;
xx[1] += 1;
}
qsqrt2[1] += 1;

qadd( qone, xx, b );
qsub( qone, xx, a );
if( a[1] == 0 )
{
qclear(y);
goto bdone;
}
qdiv( b, a, y ); /* store (x-1)/(x+1) in y */

qmul( y, y, z );

qmov( qone, a );
qmov( qone, b );
qmov( qone, qj );
do
{
qadd( qtwo, qj, qj ); /* 2 * i + 1 */
qmul( z, a, a );
qdiv( qj, a, t );
qadd( t, b, b );
}
while( (b[1] - t[1]) < NBITS );


qmul( b, y, y );
y[1] += 1;

bdone:

/* now add log of 2**ex */
if( ex != 0 )
{
ltoq( &ex, b );
qmul( qlog2, b, b );
qadd( b, y, y );
}
}