ieee.c

win32program disassembler

p = a;
q = b + (NE-1); /* point to output exponent */
/* combine sign and exponent */
i = *p++;
if( i )
	*q-- = *p++ | 0x8000;
else
	*q-- = *p++;
#ifdef INFINITY
if( *(p-1) == 0x7fff )
	{
#ifdef NANS
	if( eiisnan(a) )
		{
		enan( b, NBITS );
		return;
		}
#endif
	einfin(b);
	return;
	}
#endif
/* skip over guard word */
++p;
/* move the significand */
for( i=0; i<NE-1; i++ )
	*q-- = *p++;
}




/* Clear out internal format number.
 */

void ecleaz( xi )
register unsigned short *xi;
{
register int i;

for( i=0; i<NI; i++ )
	*xi++ = 0;
}

/* same, but don't touch the sign. */

void ecleazs( xi )
register unsigned short *xi;
{
register int i;

++xi;
for(i=0; i<NI-1; i++)
	*xi++ = 0;
}




/* Move internal format number from a to b.
 */
void emovz( a, b )
register unsigned short *a, *b;
{
register int i;

for( i=0; i<NI-1; i++ )
	*b++ = *a++;
/* clear low guard word */
*b = 0;
}

/* Return nonzero if internal format number is a NaN.
 */

int eiisnan (x)
unsigned short x[];
{
int i;

if( (x[E] & 0x7fff) == 0x7fff )
	{
	for( i=M+1; i<NI; i++ )
		{
		if( x[i] != 0 )
			return(1);
		}
	}
return(0);
}

#ifdef INFINITY
/* Return nonzero if internal format number is infinite. */

static int 
eiisinf (x)
     unsigned short x[];
{

#ifdef NANS
  if (eiisnan (x))
    return (0);
#endif
  if ((x[E] & 0x7fff) == 0x7fff)
    return (1);
  return (0);
}
#endif

/*
;	Compare significands of numbers in internal format.
;	Guard words are included in the comparison.
;
;	unsigned short a[NI], b[NI];
;	cmpm( a, b );
;
;	for the significands:
;	returns	+1 if a > b
;		 0 if a == b
;		-1 if a < b
*/
int ecmpm( a, b )
register unsigned short *a, *b;
{
int i;

a += M; /* skip up to significand area */
b += M;
for( i=M; i<NI; i++ )
	{
	if( *a++ != *b++ )
		goto difrnt;
	}
return(0);

difrnt:
if( *(--a) > *(--b) )
	return(1);
else
	return(-1);
}


/*
;	Shift significand down by 1 bit
*/

void eshdn1(x)
register unsigned short *x;
{
register unsigned short bits;
int i;

x += M;	/* point to significand area */

bits = 0;
for( i=M; i<NI; i++ )
	{
	if( *x & 1 )
		bits |= 1;
	*x >>= 1;
	if( bits & 2 )
		*x |= 0x8000;
	bits <<= 1;
	++x;
	}	
}



/*
;	Shift significand up by 1 bit
*/

void eshup1(x)
register unsigned short *x;
{
register unsigned short bits;
int i;

x += NI-1;
bits = 0;

for( i=M; i<NI; i++ )
	{
	if( *x & 0x8000 )
		bits |= 1;
	*x <<= 1;
	if( bits & 2 )
		*x |= 1;
	bits <<= 1;
	--x;
	}
}



/*
;	Shift significand down by 8 bits
*/

void eshdn8(x)
register unsigned short *x;
{
register unsigned short newbyt, oldbyt;
int i;

x += M;
oldbyt = 0;
for( i=M; i<NI; i++ )
	{
	newbyt = *x << 8;
	*x >>= 8;
	*x |= oldbyt;
	oldbyt = newbyt;
	++x;
	}
}

/*
;	Shift significand up by 8 bits
*/

void eshup8(x)
register unsigned short *x;
{
int i;
register unsigned short newbyt, oldbyt;

x += NI-1;
oldbyt = 0;

for( i=M; i<NI; i++ )
	{
	newbyt = *x >> 8;
	*x <<= 8;
	*x |= oldbyt;
	oldbyt = newbyt;
	--x;
	}
}

/*
;	Shift significand up by 16 bits
*/

void eshup6(x)
register unsigned short *x;
{
int i;
register unsigned short *p;

p = x + M;
x += M + 1;

for( i=M; i<NI-1; i++ )
	*p++ = *x++;

*p = 0;
}

/*
;	Shift significand down by 16 bits
*/

void eshdn6(x)
register unsigned short *x;
{
int i;
register unsigned short *p;

x += NI-1;
p = x + 1;

for( i=M; i<NI-1; i++ )
	*(--p) = *(--x);

*(--p) = 0;
}

/*
;	Add significands
;	x + y replaces y
*/

void eaddm( x, y )
unsigned short *x, *y;
{
register unsigned long a;
int i;
unsigned int carry;

x += NI-1;
y += NI-1;
carry = 0;
for( i=M; i<NI; i++ )
	{
	a = (unsigned long )(*x) + (unsigned long )(*y) + carry;
	if( a & 0x10000 )
		carry = 1;
	else
		carry = 0;
	*y = (unsigned short )a;
	--x;
	--y;
	}
}

/*
;	Subtract significands
;	y - x replaces y
*/

void esubm( x, y )
unsigned short *x, *y;
{
unsigned long a;
int i;
unsigned int carry;

x += NI-1;
y += NI-1;
carry = 0;
for( i=M; i<NI; i++ )
	{
	a = (unsigned long )(*y) - (unsigned long )(*x) - carry;
	if( a & 0x10000 )
		carry = 1;
	else
		carry = 0;
	*y = (unsigned short )a;
	--x;
	--y;
	}
}


/* Divide significands */

static unsigned short equot[NI] = {0}; /* was static */

#if 0
int edivm( den, num )
unsigned short den[], num[];
{
int i;
register unsigned short *p, *q;
unsigned short j;

p = &equot[0];
*p++ = num[0];
*p++ = num[1];

for( i=M; i<NI; i++ )
	{
	*p++ = 0;
	}

/* Use faster compare and subtraction if denominator
 * has only 15 bits of significance.
 */
p = &den[M+2];
if( *p++ == 0 )
	{
	for( i=M+3; i<NI; i++ )
		{
		if( *p++ != 0 )
			goto fulldiv;
		}
	if( (den[M+1] & 1) != 0 )
		goto fulldiv;
	eshdn1(num);
	eshdn1(den);

	p = &den[M+1];
	q = &num[M+1];

	for( i=0; i<NBITS+2; i++ )
		{
		if( *p <= *q )
			{
			*q -= *p;
			j = 1;
			}
		else
			{
			j = 0;
			}
		eshup1(equot);
		equot[NI-2] |= j;
		eshup1(num);
		}
	goto divdon;
	}

/* The number of quotient bits to calculate is
 * NBITS + 1 scaling guard bit + 1 roundoff bit.
 */
fulldiv:

p = &equot[NI-2];
for( i=0; i<NBITS+2; i++ )
	{
	if( ecmpm(den,num) <= 0 )
		{
		esubm(den, num);
		j = 1;	/* quotient bit = 1 */
		}
	else
		j = 0;
	eshup1(equot);
	*p |= j;
	eshup1(num);
	}

divdon:

eshdn1( equot );
eshdn1( equot );

/* test for nonzero remainder after roundoff bit */
p = &num[M];
j = 0;
for( i=M; i<NI; i++ )
	{
	j |= *p++;
	}
if( j )
	j = 1;


for( i=0; i<NI; i++ )
	num[i] = equot[i];
return( (int )j );
}

/* Multiply significands */
int emulm( a, b )
unsigned short a[], b[];
{
unsigned short *p, *q;
int i, j, k;

equot[0] = b[0];
equot[1] = b[1];
for( i=M; i<NI; i++ )
	equot[i] = 0;

p = &a[NI-2];
k = NBITS;
while( *p == 0 ) /* significand is not supposed to be all zero */
	{
	eshdn6(a);
	k -= 16;
	}
if( (*p & 0xff) == 0 )
	{
	eshdn8(a);
	k -= 8;
	}

q = &equot[NI-1];
j = 0;
for( i=0; i<k; i++ )
	{
	if( *p & 1 )
		eaddm(b, equot);
/* remember if there were any nonzero bits shifted out */
	if( *q & 1 )
		j |= 1;
	eshdn1(a);
	eshdn1(equot);
	}

for( i=0; i<NI; i++ )
	b[i] = equot[i];

/* return flag for lost nonzero bits */
return(j);
}

#else

/* Multiply significand of e-type number b
by 16-bit quantity a, e-type result to c. */

void m16m( a, b, c )
unsigned short a;
unsigned short b[], c[];
{
register unsigned short *pp;
register unsigned long carry;
unsigned short *ps;
unsigned short p[NI];
unsigned long aa, m;
int i;

aa = a;
pp = &p[NI-2];
*pp++ = 0;
*pp = 0;
ps = &b[NI-1];

for( i=M+1; i<NI; i++ )
	{
	if( *ps == 0 )
		{
		--ps;
		--pp;
		*(pp-1) = 0;
		}
	else
		{
		m = (unsigned long) aa * *ps--;
		carry = (m & 0xffff) + *pp;
		*pp-- = (unsigned short )carry;
		carry = (carry >> 16) + (m >> 16) + *pp;
		*pp = (unsigned short )carry;
		*(pp-1) = carry >> 16;
		}
	}
for( i=M; i<NI; i++ )
	c[i] = p[i];
}


/* Divide significands. Neither the numerator nor the denominator
is permitted to have its high guard word nonzero.  */


int edivm( den, num )
unsigned short den[], num[];
{
int i;
register unsigned short *p;
unsigned long tnum;
unsigned short j, tdenm, tquot;
unsigned short tprod[NI+1];

p = &equot[0];
*p++ = num[0];
*p++ = num[1];

for( i=M; i<NI; i++ )
	{
	*p++ = 0;
	}
eshdn1( num );
tdenm = den[M+1];
for( i=M; i<NI; i++ )
	{
	/* Find trial quotient digit (the radix is 65536). */
	tnum = (((unsigned long) num[M]) << 16) + num[M+1];

	/* Do not execute the divide instruction if it will overflow. */
        if( (tdenm * 0xffffL) < tnum )
		tquot = 0xffff;
	else
		tquot = tnum / tdenm;

		/* Prove that the divide worked. */
/*
	tcheck = (unsigned long )tquot * tdenm;
	if( tnum - tcheck > tdenm )
		tquot = 0xffff;
*/
	/* Multiply denominator by trial quotient digit. */
	m16m( tquot, den, tprod );
	/* The quotient digit may have been overestimated. */
	if( ecmpm( tprod, num ) > 0 )
		{
		tquot -= 1;
		esubm( den, tprod );
		if( ecmpm( tprod, num ) > 0 )
			{
			tquot -= 1;
			esubm( den, tprod );
			}
		}
/*
	if( ecmpm( tprod, num ) > 0 )
		{
		eshow( "tprod", tprod );
		eshow( "num  ", num );
		printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
			 tnum, den[M+1], tquot );
		}
*/
	esubm( tprod, num );
/*
	if( ecmpm( num, den ) >= 0 )
		{
		eshow( "num  ", num );
		eshow( "den  ", den );
		printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
			 tnum, den[M+1], tquot );
		}
*/
	equot[i] = tquot;
	eshup6(num);
	}
/* test for nonzero remainder after roundoff bit */
p = &num[M];
j = 0;
for( i=M; i<NI; i++ )
	{
	j |= *p++;
	}
if( j )
	j = 1;

for( i=0; i<NI; i++ )
	num[i] = equot[i];

return( (int )j );
}



/* Multiply significands */
int emulm( a, b )
unsigned short a[], b[];
{
unsigned short *p, *q;
unsigned short pprod[NI];
unsigned short j;
int i;

equot[0] = b[0];
equot[1] = b[1];
for( i=M; i<NI; i++ )
	equot[i] = 0;

j = 0;
p = &a[NI-1];
q = &equot[NI-1];
for( i=M+1; i<NI; i++ )
	{
	if( *p == 0 )
		{
		--p;
		}
	else
		{
		m16m( *p--, b, pprod );
		eaddm(pprod, equot);
		}
	j |= *q;
	eshdn6(equot);
	}

for( i=0; i<NI; i++ )
	b[i] = equot[i];

/* return flag for lost nonzero bits */
return( (int)j );
}


/*
eshow(str, x)
char *str;
unsigned short *x;
{
int i;

printf( "%s ", str );
for( i=0; i<NI; i++ )
	printf( "%04x ", *x++ );
printf( "\n" );
}
*/
#endif



/*
 * Normalize and round off.
 *
 * The internal format number to be rounded is "s".
 * Input "lost" indicates whether the number is exact.
 * This is the so-called sticky bit.
 *
 * Input "subflg" indicates whether the number was obtained
 * by a subtraction operation.  In that case if lost is nonzero
 * then the number is slightly smaller than indicated.
 *
 * Input "exp" is the biased exponent, which may be negative.
 * the exponent field of "s" is ignored but is replaced by
 * "exp" as adjusted by normalization and rounding.
 *
 * Input "rcntrl" is the rounding control.
 */

static int rlast = -1;
static int rw = 0;
static unsigned short rmsk = 0;
static unsigned short rmbit = 0;
static unsigned short rebit = 0;
static int re = 0;
static unsigned short rbit[NI] = {0,0,0,0,0,0,0,0};

void emdnorm( s, lost, subflg, exp, rcntrl )
unsigned short s[];