175 lines
4.4 KiB
C
175 lines
4.4 KiB
C
/* 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 conformance
|
|
02a,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 additions
|
|
01b,30jul92,kdl changed _d_type() calls to fpTypeGet().
|
|
01a,08jul92,smb documentation.
|
|
*/
|
|
|
|
/*
|
|
DESCRIPTION
|
|
SEE ALSO: American National Standard X3.159-1989
|
|
NOMANUAL
|
|
*/
|
|
|
|
#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);
|
|
}
|