/*
 *  Ieee double normalization function
 *
 *  double norm( double d, int exp, int sign, int rbits )
 *
 *  inputs:
 *	m	a 63 bit mantissa (passed as 8 bytes in a double's storage)
 *	exp	BIASED exponent
 *	sign	sign of the result (1 == -ve  0 == +ve)
 *	rbits	8 bits of rounding info
 *	  on input radix point is at extreme right (exp should reflect this)
 * output:
 *	the normalized double in ieee double format
 *
 * algorithm:
 *	see posting by p. housel housel@ecn.purdue.edu (thanks peter)
 *	    and  knuth v2 section 4.2.1	algorithm N
 *
 *	++jrb	bammi@dsrgsun.ces.cwru.edu
 * bugs:
 *	strictly mc68X coding
 *	(can be coded using gnu c `long long' type very easily)
 */
#include "flonum.h"

struct ldouble {
    unsigned long hi, lo;
};		/* yet another alias for a double */

#if 0
double norm( double m, int exp, int sign, int rbits )
#else
double __normdf( double m, int exp, int sign, int rbits )
#endif
{
    /* we'll use constraints later for values to end up in D reggies */
    register unsigned long d0, d1;	/*  64 bit accumulator     */
    register char rounding; /* rounding reg-signed char for convenience */
    double result;			/* the normalized result    */
    struct bitdouble *r = (struct bitdouble *)&result; /* alias    */

    /* init */
    d0 = ((struct ldouble *)&m)->hi;
    d1 = ((struct ldouble *)&m)->lo;
    rounding = rbits;
    
    if( (d0 == 0) && (d1 == 0) && (rounding == 0) )
	goto stuff_result; /* nothing to do, but stuff bits into result */

    /* keep shifting R & accumulating rounding bits until we have 53 bits */
    /* lsb of rounding register is sticky for 1 bits shifted out          */

    asm volatile("ReNormalize:");	/* ok, i admit it :-) */
    while( (d0 & 0xffe00000L) != 0)
    {
#if 0 /* only 5 operands allowed in asm() - why?? */
	asm volatile (" lsrl	#1,%2 ; 
		        roxrl	#1,%1 ;
 			roxrb	#1,%0 ;
 			jcc	NoCarry;
 			orb	#1,%0 ;
 			NoCarry:      "
        : "=d" (rounding), "=d" (d1), "=d" (d0)	   /* outputs */
        : "0"  (rounding), "1"  (d1), "2"  (d0) ); /* inputs  */
#else
	asm volatile (" lsrl	#1,%1 ; 
		        roxrl	#1,%0 "
        : "=d" (d1), "=d" (d0)	   /* outputs */
        : "0"  (d1), "1"  (d0) ); /* inputs  */
        asm volatile ("	roxrb	#1,%0 ;
 			jcc	NoCarry;
 			orb	#1,%0 ;
 			NoCarry:      "
        : "=d" (rounding)	   /* outputs */
        : "0"  (rounding) ); /* inputs  */
#endif
	/* increment exponent to reflect the shift */
	exp++;
    }

    /* shift L until there is a 1 in the hidden bit pos, shifting
       the rounding reggie into the lsb */
    while( (d0 & 0x00100000L) == 0)
    {
#if 0 /* same here */
	asm volatile ("	lslb	#1,%2 ;
 			roxll	#1,%1 ;
 			roxll	#1,%0 "
        : "=d" (d0), "=d" (d1), "=d" (rounding)	   /* outputs */
        : "0"  (d0), "1"  (d1), "2"  (rounding) ); /* inputs  */
#else
	asm volatile ("	lslb	#1,%1 ;
 			roxll	#1,%0 "
        : "=d" (d1), "=d" (rounding)	   /* outputs */
        : "0"  (d1), "1"  (rounding) );	   /* inputs  */

	asm volatile ("	roxll	#1,%0 "
        : "=d" (d0)	   /* outputs */
        : "0"  (d0)); 	   /* inputs  */
#endif
	/* decrement exponent to reflect the shift */
	--exp;
    }

    /* check rounding register */
    if (rounding < 0)
    {	/* not round down */
	if( (((unsigned char)rounding) & 0x80) == 0)
	{	/* tie: round to even */
	    d1 &= 0xfffffffeL;	/* set lsb to 0 */
	}
	else
	{	/* round up */
	    rounding = 0;
	asm volatile ("
		addql  #1,%1 ;
		jcc    NoCarry1 ;
		addql  #1,%0 ;
    		bra   ReNormalize; |just in case of overflow into hidden bit\n
                NoCarry1:           "
        : "=d" (d0), "=d" (d1)	   /* outputs */
        : "0"  (d0), "1"  (d1) );  /* inputs  */
	}
    }
    
    /* exponent underflow ?? */
    if(exp <= 0)
    {
	printf("Underflow %lx %lx %d\n", d0, d1, exp);
	d0 = 0;
	d1 = 0;
	sign = 0;
	exp = 0;
	goto stuff_result;
    }
    /* exp overflow ? */
    if(exp >= 2047)
    {
	/* cause overflow trap, but how ?? */
	printf("Overflow %lx %lx %d\n", d0, d1, exp);
    }
    
    stuff_result: /* stuff bit in result and ret */
    r->sign = sign;
    r->exp  = exp;
    r->mant1 = d0;
    r->mant2 = d1;
    
    return result;
}

#ifdef TEST
#include <stdio.h>

main()
{
    register unsigned long d0, d1;
    double r1, r2, r3;
    struct bitdouble *pr1 = (struct bitdouble *)&r1;
    struct bitdouble *pr3 = (struct bitdouble *)&r3;
    unsigned long *l = (unsigned long *)&r2;
    
    r2 = r1 = 3.14159265358979323846265e23;

    *l &= 0x000FFFFFL;	/* wallup sign and exponent  in r2 */
    *l |= 0x00100000L;  /* add hidden bit	    	   */
    
    /* try the straight case first */
    r3 = norm(r2, (unsigned int)pr1->exp, (unsigned int)pr1->sign, 0);
    
    printf("%30.25e %05lx %08lx %4d %1d\n", r1, pr1->mant1, pr1->mant2,
	   pr1->exp, pr1->sign);
    printf("%30.25e %05lx %08lx %4d %1d\n\n", r3, pr3->mant1, pr3->mant2,
	   pr3->exp, pr3->sign);

    /* shift L and try */
    d0 = l[0];
    d1 = l[1];
    asm volatile ("	lsll	#1,%1 ;
			roxll	#1,%0 "
        : "=d" (d0), "=d" (d1)	   /* outputs */
        : "0"  (d0), "1"  (d1) );  /* inputs  */
    l[0] = d0;
    l[1] = d1;
    
    r3 = norm(r2, (unsigned int)pr1->exp - 1, (unsigned int)pr1->sign, 0);
    
    printf("%30.25e %05lx %08lx %4d %1d\n", r1, pr1->mant1, pr1->mant2,
	   pr1->exp, pr1->sign);
    printf("%30.25e %05lx %08lx %4d %1d\n\n", r3, pr3->mant1, pr3->mant2,
	   pr3->exp, pr3->sign);

    /* how about a shift R */
    r2 =r1;

    *l &= 0x000FFFFFL;	/* wallup sign and exponent  in r2 */
    *l |= 0x00100000L;  /* add hidden bit	    	   */

    d0 = l[0];
    d1 = l[1];
    asm volatile ("	lsrl	#1,%0 ;
			roxrl	#1,%1 "
        : "=d" (d0), "=d" (d1)	   /* outputs */
        : "0"  (d0), "1"  (d1) );  /* inputs  */
    l[0] = d0;
    l[1] = d1;
    
    r3 = norm(r2, (unsigned int)pr1->exp + 1, (unsigned int)pr1->sign, 0);
    
    printf("%30.25e %05lx %08lx %4d %1d\n", r1, pr1->mant1, pr1->mant2,
	   pr1->exp, pr1->sign);
    printf("%30.25e %05lx %08lx %4d %1d\n\n", r3, pr3->mant1, pr3->mant2,
	   pr3->exp, pr3->sign);
}
#endif