#define __normdf norm	/* for now */
/*
 * double modf(value, iptr)
 * double value;
 * double *iptr;
 *
 * returns fractional part of value
 *	   in *iptr returns the integral part
 *	   such that (*iptr + fractional) == value
 *
 * ++jrb	bammi@dsrgsun.ces.cwru.edu
 *
 * special thanks to peter housel  housel@ecn.purdue.edu
 * this (improved) version is based on his X86 code
 *
 */
#include "flonum.h"

#define BIAS 1023
#define B1   1022	/* BIAS - 1 */
#define B2   1075	/* (BIAS - 1) + 53 == max bias for integ part    */
#define B3   1011	/* (BIAS - 1) - 11 == the max bias for a frac. # */

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

double modf(value, iptr)
double value, *iptr;
{
    struct bitdouble *bd = (struct bitdouble *)&value;
    unsigned int expo = bd->exp;
#ifdef __STDC__
    double __normdf( double, int, int, int );
#else
    extern double __normdf();
#endif
    
    if(expo <= B1)	/* all frac, no integ */
    {
	*iptr = 0.0;
	return value;
    }
    if(expo >= B2)	/* all integ, no frac */
    {
	*iptr = value;
	return 0.0;
    }

    /* both integ and frac */
    {
	register unsigned long i0, i1;	/* integral   part accumulator */
	register unsigned long f0, f1;	/* fractional part accumulator */
	unsigned int sign  = bd->sign;
	
	i0 = bd->mant1 | 0x00100000L;
	i1 = bd->mant2;
	f0 = f1 = 0L;

	do
	{
	    /* shift R integ/frac, with bit falling into frac acc */
	    asm volatile(" lsrl  #1,%1;
 			   roxrl #1,%0"
		: "=d" (i1), "=d" (i0)
		: "0"  (i1), "1"  (i0) );

	    asm volatile(" roxrl #1,%1;
 			   roxrl #1,%0"
		: "=d" (f1), "=d" (f0)
		: "0"  (f1), "1"  (f0) );

	} while(++expo < B2);

	/* stuff the results, and normalize */
	((struct ldouble *)&value)->hi = i0;	/* integral part */
	((struct ldouble *)&value)->lo = i1;
	*iptr = __normdf(value, expo, sign , 0);

	/* dont call norm if frac is all zero, as B3 exp will screw it up */
	if((f0 == 0) && (f1 == 0))
	    return 0.0;

	((struct ldouble *)&value)->hi = f0;	/* fractional part */
	((struct ldouble *)&value)->lo = f1;
	return __normdf(value, B3, sign, 0);
    }
}


#ifdef TEST
#include <stdio.h>
#ifdef __MSHORT__
#error "please run this test in 32 bit int mode"
#endif

#define NTEST	100000L
#define MAXRAND 0x7fffffffL	/* assuming 32 bit ints */
#define ABS(X)	( (X < 0) ? (-X) : X )
extern long rand();

int main()
{
    double frac, integ, e, expected, result, max_abs_err;
    register long i;
    register long errs = 0;
    register int s;

    max_abs_err = 0.0;
    for(i = 0; i < NTEST; i++)
    {
	expected = ((double)(s = rand()) + (double)rand())/(double)(MAXRAND);
	if(s > (MAXRAND >> 1)) expected = -expected;
	frac = modf(expected, &integ);
	if(ABS(frac) >= 1.0)
	{
	    printf("|frac| >= 1, %.6e:  integ %.6e  frac %.6e\n",
		   expected, integ, frac);
	    errs++;
	    continue;
	}
	if( (integ != 0) && (ABS(integ) < 1.0) )
	{
	    printf("|integ| < 1, %.6e:  integ %.6e  frac %.6e\n",
		   expected, integ, frac);
	    errs++;
	    continue;
	}

	result = integ + frac;
	e = (expected == 0.0) ? result : (result - expected)/expected;
	if(e < 0) e = (-e);
	if(e > 1.0e-6)
	{
	    printf("%.4e: I %.4e F %.4e R %.4e E %.8e\n",
		   expected, integ, frac, result, e);
	    errs++;
	}
	if (e > max_abs_err) max_abs_err = e;
    }
    
    printf("%ld Error(s), Max abs err %.14e\n", errs, max_abs_err);
    return errs;
}
#endif /* TEST */