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

 
Output of file : QSQRTA.C contained in archive : CEPHES22.ZIP
/* qsqrta.c */
/* Square root check routine, done by long division. */
/* Copyright (C) 1984-1988 by Stephen L. Moshier. */


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

extern unsigned short qzero[], qone[];
static unsigned short rndbit[NQ+1] = {0};
static short rndset = 0;

int qsqrt( x, y )
short *x, *y;
{
unsigned short temp[NQ+1], num[NQ+1], sq[NQ+1], xx[NQ+1];
int i, j, k, n;
long m, exp;

if( rndset == 0 )
{
qclear( rndbit );
rndbit[NQ] = 0;
rndbit[NQ-1] = 1;
rndset = 1;
}
/* Check for arg <= 0 */
/*i = qcmp( x, qzero );*/
if( x[1] == 0 )
goto iszero;
if( x[0] != 0 )
{
mtherr( "qsqrt", DOMAIN );
iszero:
qclear(y);
return;
}

/* Bring in the arg and renormalized if it is denormal. */
qmovz( x, xx );
m = xx[1]; /* local long word exponent */
/*
if( m == 0 )
m -= normlz( xx );
*/
/* Divide exponent by 2 */
m -= 0x4000;
exp = (unsigned short )( (m / 2) + 0x4000 );

/* Adjust if exponent odd */
if( (m & 1) != 0 )
{
if( m > 0 )
exp += 1;
shdn1( xx );
}

qclear( sq );
sq[NQ] = 0;
qclear( num );
num[NQ] = 0;
n = 8;

for( j=0; j<2*NQ-6; j++ )
{
/* bring in next word of arg */
if( j < NQ-3 )
num[NQ] = xx[j+3];
/* Do additional bit on last outer loop, for roundoff. */
if( j == 2*NQ-7 )
n = 9;
for( i=0; i {
/* Next 2 bits of arg */
shup1( num );
shup1( num );
/* Shift up answer */
shup1( sq );
/* Make trial divisor */
for( k=0; k<=NQ; k++ )
temp[k] = sq[k];
shup1( temp );
addm( rndbit, temp );
/* Subtract and insert answer bit if it goes in */
if( cmpm( temp, num ) <= 0 )
{
subm( temp, num );
sq[NQ-1] |= 1;
}
}
}

exp -= 1; /* Adjust for extra, roundoff loop done. */
sq[1] = exp;

/* Sticky bit = 1 if the remainder is nonzero. */
k = 0;
for( i=3; i<=NQ; i++ )
k |= num[i];

/* Renormalize and round off. */
mdnorm( sq, k );
qmov( sq, y );
}


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : QSQRTA.C

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/