/* fmod.c - fmod math routine *//* Copyright 1992-2002 Wind River Systems, Inc. *//*modification history--------------------03b,01apr02,pfl  clean up fix for SPR 20486 plus                  ISO/IEC 9899:1999 conformance02a,31oct01,md3  +gkc +ycl rewrite: fixed SPR #20486.01f,05feb93,jdi  doc changes based on kdl review.01e,02dec92,jdi  doc tweaks.01d,28oct92,jdi  documentation cleanup.01c,20sep92,smb  documentation additions01b,30jul92,kdl  changed _d_type() calls to fpTypeGet().01a,08jul92,smb  documentation.*//*DESCRIPTIONSEE ALSO: American National Standard X3.159-1989NOMANUAL*/#include "vxWorks.h"#include "math.h"#include "stddef.h"#include "errno.h"#include "private/mathP.h"#include "errnoLib.h"/* mask to extract exponent part of a double */#define  MASK_DB_EXP  0x7ff00000/* mask to extract high 20 significant bits */#define  MASK_DB_HI20  0x000fffff/* mask to filter exponent part of double */#define  MASK_DB_MANT  0x800fffff/* for big endian machine */#if (_BYTE_ORDER==_BIG_ENDIAN)   /* high 32 bits of Db - 1-bit sign + 11-bits exponent + 20-bit significant */#define DB_HI(d) ((unsigned int) *((int*)&d))/* hi 20 bits of a double's significant */#define DB_HI_SIG(d) ((unsigned int) *((int*)&d) & MASK_DB_HI20)  /* lower 32 bits of double's significant */#define DB_LO_SIG(d) ((unsigned int) *(1+(int*)&d))      /* for little endian machine */#else /* (_BYTE_ORDER != _BIG_ENDIAN) */  /* high 32 bits of Db - 1-bit sign + 11-bits exponent + 20-bit significant */#define DB_HI(d) ((unsigned int) *(1+(int*)&d))      /* hi 20 bits of a double's significant */#define DB_HI_SIG(d) ((unsigned int) *(1+(int*)&d) & MASK_DB_HI20)  /* lower 32 bits of double's significant */#define DB_LO_SIG(d) ((unsigned int) *((int*)&d))                           #endif /* (_BYTE_ORDER==_BIG_ENDIAN) *//* return (significant of d1) > (significant of d2) */#define DB_GT_SIG(d1,d2) \        ((DB_HI_SIG(d1) > DB_HI_SIG(d2)) || \         ((DB_HI_SIG(d1) == DB_HI_SIG(d2)) && (DB_LO_SIG(d1) > DB_LO_SIG(d2))))/* NAN value */#define NAN_VAL(d) {*(int *)(&d) = 0xffffffff; *(1+(int *)&d)=0xffffffff;}/********************************************************************************* fmod - compute the remainder of x/y (ANSI)** This routine returns the remainder of <x>/<y> with the sign of <x>,* in double precision.* * INCLUDE FILES: math.h** RETURNS: The value <x> - <i> * <y>, for some integer <i>.  If <y> is* non-zero, the result has the same sign as <x> and magnitude less than the* magnitude of <y>.  If <y> is zero, fmod() returns zero.** ERRNO: EDOM** SEE ALSO: mathALib*/double fmod    (    double x,    /* numerator   */    double y     /* denominator */    )    {    double       negative = 1.0;                /* used for sign of return   */    double       yMult2n;   	                /* computed next decrement   */    unsigned int zeroExp_y;            		/* factor for next decrement */    int          errx = fpTypeGet (x, NULL);    /* determine number type */    int          erry = fpTypeGet (y, NULL);    /* determine number type */    /*      * Per ISO/IEC 9899:1999, ANSI/IEEE 754-1985     * fmod ( +-0,     y) = +-0 (if y != 0)     * fmod (   x,     y) = NAN (if x  = INF or y = 0)      * fmod (   x, +-INF) = x   (if x != INF)     */       /* check for boundary conditions returning NAN */        if (errx == NAN || erry == NAN || errx == INF || erry == ZERO)        {        errnoSet(EDOM);         /* Domain error is set to maintain backward				   compatibility with ANSI X3.159-1989 */	NAN_VAL(x);        return (x);		/* return NAN */        }    /* check for boundary conditions returning X */    if (errx == ZERO || erry == INF)        return (x);    /* make Y absolute */    if (y < 0.0)        y = -y;    /* make X absolute, record sign for return. */    if (x < 0.0)        {        x = -x;        negative = -1.0;        }    /* initial decrementing conditions */    yMult2n = y;    zeroExp_y= MASK_DB_MANT & (DB_HI(y));    /* compute fmod by subtracting estimated decrementations */    while (x>=y)         {        if (DB_GT_SIG(y,x))            {             (DB_HI(yMult2n)) = ( (DB_HI(x)-0x100000) & MASK_DB_EXP ) | zeroExp_y;             }         else             {             (DB_HI(yMult2n)) = ( DB_HI(x) & MASK_DB_EXP ) | zeroExp_y;             }         x-= yMult2n ;         }    return (negative*x);    }