/*
 * double modf(value, iptr)
 * double value;
 * double *iptr;
 *
 * returns fractional part of value
 *	   in *iptr rets. the integral part such that (*iptr + modf) == value
 *
 * ++jrb	bammi@dsrgsun.ces.cwru.edu
 */
#include "flonum.h"

/*
 * extract integral portion of an Ieee double -- return it as a double
 * ++jrb
 */

#define BIAS  1023
#define B1    1022
#define MANT1 21
#define MANT2 32
#define MASK  0xFFFFFFFFL
struct ldouble {
    unsigned long hi, lo;
};
static double extract_int(val)
double val;
{
    struct bitdouble *bd = (struct bitdouble *)&val;
    struct ldouble *l    = (struct ldouble *)&val;
    register int ex = bd->exp - BIAS, sign = bd->sign;
#ifdef __STDC__
    double norm( double, int, int, int );
#else
    extern double norm();
#endif
    
    /* trivial case */
    if(ex < 0) return 0.0;
    
    bd->sign = 0;
    bd->exp  = 1;  /* hidden bit */
    /* exponent <= # of bits in mant1 */
    if(ex <= MANT1)
    {
	l->hi &= (MASK << (MANT1 - ex));
	l->lo = 0L;
    }
    else /* if (ex <= (MANT1 + MANT2)) */
    {
	l->lo &= (MASK << ((MANT1 + MANT2) - ex));
    }
    /* else the value is too big */

    return norm(val, ex+BIAS, sign, 0);
}


double modf(value, iptr)
double value, *iptr;
{
    double integral;
    
    if((-1.0 < value) && (value < 1.0))
    {
	*iptr = 0.0;
	return value;
    }

    integral = extract_int(value);
    *iptr = integral;
    return value - integral;
}

#ifdef TEST1
#ifndef TEST
#define TEST
#endif
#endif

#ifdef TEST

#define  MAXRAND 0x7fffffffL

#ifndef TEST1
#define TESTSIZE 10240
double testin[TESTSIZE];

static void init()
{
    register int i;
    extern long rand();
    
    for(i = 0; i < TESTSIZE; i++)
	testin[i]  = ((double)rand()) + (((double)rand()) / ((double)MAXRAND));
}

#else

#define TESTSIZE 23
double testin[TESTSIZE];
static void init()
{
    register int i;
    double v;
    double inc = 0.1;
    
    for(v = -1.1, i = 0; i < TESTSIZE; v += inc, i++)
	testin[i] = v;
}
#endif /* TEST1 */

#define ABS(x) (((x) < 0.0)? (-x) : (x))
#define MARGIN 1.0e-6		/* 6 is arb */

#include <stdio.h>
int main()
{
    double frac, integ, e, r;
    register int i;
    register int errors = 0;
    extern double modf();

    init();
    for(i = 0; i < TESTSIZE; i++)
    {
	frac = modf(testin[i], &integ);
	if(frac >= 1.0)
	{
	    printf("Error, frac >= 1, testin %f integ %f frac %f\n",
		   testin[i], integ, frac);
	    errors++;
	    continue;
	}

	r = integ + frac;
	e = testin[i] - r;
	if(ABS(e) > MARGIN)
	{
	    printf("In %f\tInt %f Frac %f Res %f\tError %f\n",
		   testin[i], integ, frac, r, e);
	    errors++;
	}
#ifdef TEST1
	printf("%f\t%f  %f\n", testin[i], integ, frac);
#endif
    }
    
    printf("%d error(s)\n", errors);
    return errors;
}
#endif /* TEST */