454 lines
6.1 KiB
C
454 lines
6.1 KiB
C
|
/* ceil()
|
|||
|
|
* floor()
|
|||
|
|
* frexp()
|
|||
|
|
* ldexp()
|
|||
|
|
* signbit()
|
|||
|
|
* isnan()
|
|||
|
|
* isfinite()
|
|||
|
|
*
|
|||
|
|
* Floating point numeric utilities
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* SYNOPSIS:
|
|||
|
|
*
|
|||
|
|
* double ceil(), floor(), frexp(), ldexp();
|
|||
|
|
* int signbit(), isnan(), isfinite();
|
|||
|
|
* double x, y;
|
|||
|
|
* int expnt, n;
|
|||
|
|
*
|
|||
|
|
* y = floor(x);
|
|||
|
|
* y = ceil(x);
|
|||
|
|
* y = frexp( x, &expnt );
|
|||
|
|
* y = ldexp( x, n );
|
|||
|
|
* n = signbit(x);
|
|||
|
|
* n = isnan(x);
|
|||
|
|
* n = isfinite(x);
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* DESCRIPTION:
|
|||
|
|
*
|
|||
|
|
* All four routines return a double precision floating point
|
|||
|
|
* result.
|
|||
|
|
*
|
|||
|
|
* floor() returns the largest integer less than or equal to x.
|
|||
|
|
* It truncates toward minus infinity.
|
|||
|
|
*
|
|||
|
|
* ceil() returns the smallest integer greater than or equal
|
|||
|
|
* to x. It truncates toward plus infinity.
|
|||
|
|
*
|
|||
|
|
* frexp() extracts the exponent from x. It returns an integer
|
|||
|
|
* power of two to expnt and the significand between 0.5 and 1
|
|||
|
|
* to y. Thus x = y * 2**expn.
|
|||
|
|
*
|
|||
|
|
* ldexp() multiplies x by 2**n.
|
|||
|
|
*
|
|||
|
|
* signbit(x) returns 1 if the sign bit of x is 1, else 0.
|
|||
|
|
*
|
|||
|
|
* These functions are part of the standard C run time library
|
|||
|
|
* for many but not all C compilers. The ones supplied are
|
|||
|
|
* written in C for either DEC or IEEE arithmetic. They should
|
|||
|
|
* be used only if your compiler library does not already have
|
|||
|
|
* them.
|
|||
|
|
*
|
|||
|
|
* The IEEE versions assume that denormal numbers are implemented
|
|||
|
|
* in the arithmetic. Some modifications will be required if
|
|||
|
|
* the arithmetic has abrupt rather than gradual underflow.
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
|
|||
|
|
/*
|
|||
|
|
Cephes Math Library Release 2.8: June, 2000
|
|||
|
|
Copyright 1984, 1995, 2000 by Stephen L. Moshier
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
|
|||
|
|
#include "mconf.h"
|
|||
|
|
|
|||
|
|
#ifdef UNK
|
|||
|
|
/* ceil(), floor(), frexp(), ldexp() may need to be rewritten. */
|
|||
|
|
#undef UNK
|
|||
|
|
#if BIGENDIAN
|
|||
|
|
#define MIEEE 1
|
|||
|
|
#else
|
|||
|
|
#define IBMPC 1
|
|||
|
|
#endif
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef DEC
|
|||
|
|
#define EXPMSK 0x807f
|
|||
|
|
#define MEXP 255
|
|||
|
|
#define NBITS 56
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef IBMPC
|
|||
|
|
#define EXPMSK 0x800f
|
|||
|
|
#define MEXP 0x7ff
|
|||
|
|
#define NBITS 53
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef MIEEE
|
|||
|
|
#define EXPMSK 0x800f
|
|||
|
|
#define MEXP 0x7ff
|
|||
|
|
#define NBITS 53
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
extern double MAXNUM, NEGZERO;
|
|||
|
|
#ifdef ANSIPROT
|
|||
|
|
double floor ( double );
|
|||
|
|
int isnan ( double );
|
|||
|
|
int isfinite ( double );
|
|||
|
|
double ldexp ( double, int );
|
|||
|
|
#else
|
|||
|
|
double floor();
|
|||
|
|
int isnan(), isfinite();
|
|||
|
|
double ldexp();
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
double ceil(x)
|
|||
|
|
double x;
|
|||
|
|
{
|
|||
|
|
double y;
|
|||
|
|
|
|||
|
|
#ifdef UNK
|
|||
|
|
mtherr( "ceil", DOMAIN );
|
|||
|
|
return(0.0);
|
|||
|
|
#endif
|
|||
|
|
#ifdef NANS
|
|||
|
|
if( isnan(x) )
|
|||
|
|
return( x );
|
|||
|
|
#endif
|
|||
|
|
#ifdef INFINITIES
|
|||
|
|
if(!isfinite(x))
|
|||
|
|
return(x);
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
y = floor(x);
|
|||
|
|
if( y < x )
|
|||
|
|
y += 1.0;
|
|||
|
|
#ifdef MINUSZERO
|
|||
|
|
if( y == 0.0 && x < 0.0 )
|
|||
|
|
return( NEGZERO );
|
|||
|
|
#endif
|
|||
|
|
return(y);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Bit clearing masks: */
|
|||
|
|
|
|||
|
|
static unsigned short bmask[] = {
|
|||
|
|
0xffff,
|
|||
|
|
0xfffe,
|
|||
|
|
0xfffc,
|
|||
|
|
0xfff8,
|
|||
|
|
0xfff0,
|
|||
|
|
0xffe0,
|
|||
|
|
0xffc0,
|
|||
|
|
0xff80,
|
|||
|
|
0xff00,
|
|||
|
|
0xfe00,
|
|||
|
|
0xfc00,
|
|||
|
|
0xf800,
|
|||
|
|
0xf000,
|
|||
|
|
0xe000,
|
|||
|
|
0xc000,
|
|||
|
|
0x8000,
|
|||
|
|
0x0000,
|
|||
|
|
};
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
double floor(x)
|
|||
|
|
double x;
|
|||
|
|
{
|
|||
|
|
union
|
|||
|
|
{
|
|||
|
|
double y;
|
|||
|
|
unsigned short sh[4];
|
|||
|
|
} u;
|
|||
|
|
unsigned short *p;
|
|||
|
|
int e;
|
|||
|
|
|
|||
|
|
#ifdef UNK
|
|||
|
|
mtherr( "floor", DOMAIN );
|
|||
|
|
return(0.0);
|
|||
|
|
#endif
|
|||
|
|
#ifdef NANS
|
|||
|
|
if( isnan(x) )
|
|||
|
|
return( x );
|
|||
|
|
#endif
|
|||
|
|
#ifdef INFINITIES
|
|||
|
|
if(!isfinite(x))
|
|||
|
|
return(x);
|
|||
|
|
#endif
|
|||
|
|
#ifdef MINUSZERO
|
|||
|
|
if(x == 0.0L)
|
|||
|
|
return(x);
|
|||
|
|
#endif
|
|||
|
|
u.y = x;
|
|||
|
|
/* find the exponent (power of 2) */
|
|||
|
|
#ifdef DEC
|
|||
|
|
p = (unsigned short *)&u.sh[0];
|
|||
|
|
e = (( *p >> 7) & 0377) - 0201;
|
|||
|
|
p += 3;
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef IBMPC
|
|||
|
|
p = (unsigned short *)&u.sh[3];
|
|||
|
|
e = (( *p >> 4) & 0x7ff) - 0x3ff;
|
|||
|
|
p -= 3;
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef MIEEE
|
|||
|
|
p = (unsigned short *)&u.sh[0];
|
|||
|
|
e = (( *p >> 4) & 0x7ff) - 0x3ff;
|
|||
|
|
p += 3;
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
if( e < 0 )
|
|||
|
|
{
|
|||
|
|
if( u.y < 0.0 )
|
|||
|
|
return( -1.0 );
|
|||
|
|
else
|
|||
|
|
return( 0.0 );
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
e = (NBITS -1) - e;
|
|||
|
|
/* clean out 16 bits at a time */
|
|||
|
|
while( e >= 16 )
|
|||
|
|
{
|
|||
|
|
#ifdef IBMPC
|
|||
|
|
*p++ = 0;
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef DEC
|
|||
|
|
*p-- = 0;
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef MIEEE
|
|||
|
|
*p-- = 0;
|
|||
|
|
#endif
|
|||
|
|
e -= 16;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* clear the remaining bits */
|
|||
|
|
if( e > 0 )
|
|||
|
|
*p &= bmask[e];
|
|||
|
|
|
|||
|
|
if( (x < 0) && (u.y != x) )
|
|||
|
|
u.y -= 1.0;
|
|||
|
|
|
|||
|
|
return(u.y);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
double frexp( x, pw2 )
|
|||
|
|
double x;
|
|||
|
|
int *pw2;
|
|||
|
|
{
|
|||
|
|
union
|
|||
|
|
{
|
|||
|
|
double y;
|
|||
|
|
unsigned short sh[4];
|
|||
|
|
} u;
|
|||
|
|
int i;
|
|||
|
|
#ifdef DENORMAL
|
|||
|
|
int k;
|
|||
|
|
#endif
|
|||
|
|
short *q;
|
|||
|
|
|
|||
|
|
u.y = x;
|
|||
|
|
|
|||
|
|
#ifdef UNK
|
|||
|
|
mtherr( "frexp", DOMAIN );
|
|||
|
|
return(0.0);
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef IBMPC
|
|||
|
|
q = (short *)&u.sh[3];
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef DEC
|
|||
|
|
q = (short *)&u.sh[0];
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef MIEEE
|
|||
|
|
q = (short *)&u.sh[0];
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
/* find the exponent (power of 2) */
|
|||
|
|
#ifdef DEC
|
|||
|
|
i = ( *q >> 7) & 0377;
|
|||
|
|
if( i == 0 )
|
|||
|
|
{
|
|||
|
|
*pw2 = 0;
|
|||
|
|
return(0.0);
|
|||
|
|
}
|
|||
|
|
i -= 0200;
|
|||
|
|
*pw2 = i;
|
|||
|
|
*q &= 0x807f; /* strip all exponent bits */
|
|||
|
|
*q |= 040000; /* mantissa between 0.5 and 1 */
|
|||
|
|
return(u.y);
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef IBMPC
|
|||
|
|
i = ( *q >> 4) & 0x7ff;
|
|||
|
|
if( i != 0 )
|
|||
|
|
goto ieeedon;
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef MIEEE
|
|||
|
|
i = *q >> 4;
|
|||
|
|
i &= 0x7ff;
|
|||
|
|
if( i != 0 )
|
|||
|
|
goto ieeedon;
|
|||
|
|
#ifdef DENORMAL
|
|||
|
|
|
|||
|
|
#else
|
|||
|
|
*pw2 = 0;
|
|||
|
|
return(0.0);
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
|
|||
|
|
#ifndef DEC
|
|||
|
|
/* Number is denormal or zero */
|
|||
|
|
#ifdef DENORMAL
|
|||
|
|
if( u.y == 0.0 )
|
|||
|
|
{
|
|||
|
|
*pw2 = 0;
|
|||
|
|
return( 0.0 );
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Handle denormal number. */
|
|||
|
|
do
|
|||
|
|
{
|
|||
|
|
u.y *= 2.0;
|
|||
|
|
i -= 1;
|
|||
|
|
k = ( *q >> 4) & 0x7ff;
|
|||
|
|
}
|
|||
|
|
while( k == 0 );
|
|||
|
|
i = i + k;
|
|||
|
|
#endif /* DENORMAL */
|
|||
|
|
|
|||
|
|
ieeedon:
|
|||
|
|
|
|||
|
|
i -= 0x3fe;
|
|||
|
|
*pw2 = i;
|
|||
|
|
*q &= 0x800f;
|
|||
|
|
*q |= 0x3fe0;
|
|||
|
|
return( u.y );
|
|||
|
|
#endif
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
double ldexp( x, pw2 )
|
|||
|
|
double x;
|
|||
|
|
int pw2;
|
|||
|
|
{
|
|||
|
|
union
|
|||
|
|
{
|
|||
|
|
double y;
|
|||
|
|
unsigned short sh[4];
|
|||
|
|
} u;
|
|||
|
|
short *q;
|
|||
|
|
int e;
|
|||
|
|
|
|||
|
|
#ifdef UNK
|
|||
|
|
mtherr( "ldexp", DOMAIN );
|
|||
|
|
return(0.0);
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
u.y = x;
|
|||
|
|
#ifdef DEC
|
|||
|
|
q = (short *)&u.sh[0];
|
|||
|
|
e = ( *q >> 7) & 0377;
|
|||
|
|
if( e == 0 )
|
|||
|
|
return(0.0);
|
|||
|
|
#else
|
|||
|
|
|
|||
|
|
#ifdef IBMPC
|
|||
|
|
q = (short *)&u.sh[3];
|
|||
|
|
#endif
|
|||
|
|
#ifdef MIEEE
|
|||
|
|
q = (short *)&u.sh[0];
|
|||
|
|
#endif
|
|||
|
|
while( (e = (*q & 0x7ff0) >> 4) == 0 )
|
|||
|
|
{
|
|||
|
|
if( u.y == 0.0 )
|
|||
|
|
{
|
|||
|
|
return( 0.0 );
|
|||
|
|
}
|
|||
|
|
/* Input is denormal. */
|
|||
|
|
if( pw2 > 0 )
|
|||
|
|
{
|
|||
|
|
u.y *= 2.0;
|
|||
|
|
pw2 -= 1;
|
|||
|
|
}
|
|||
|
|
if( pw2 < 0 )
|
|||
|
|
{
|
|||
|
|
if( pw2 < -53 )
|
|||
|
|
return(0.0);
|
|||
|
|
u.y /= 2.0;
|
|||
|
|
pw2 += 1;
|
|||
|
|
}
|
|||
|
|
if( pw2 == 0 )
|
|||
|
|
return(u.y);
|
|||
|
|
}
|
|||
|
|
#endif /* not DEC */
|
|||
|
|
|
|||
|
|
e += pw2;
|
|||
|
|
|
|||
|
|
/* Handle overflow */
|
|||
|
|
#ifdef DEC
|
|||
|
|
if( e > MEXP )
|
|||
|
|
return( MAXNUM );
|
|||
|
|
#else
|
|||
|
|
if( e >= MEXP )
|
|||
|
|
return( 2.0*MAXNUM );
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
/* Handle denormalized results */
|
|||
|
|
if( e < 1 )
|
|||
|
|
{
|
|||
|
|
#ifdef DENORMAL
|
|||
|
|
if( e < -53 )
|
|||
|
|
return(0.0);
|
|||
|
|
*q &= 0x800f;
|
|||
|
|
*q |= 0x10;
|
|||
|
|
/* For denormals, significant bits may be lost even
|
|||
|
|
when dividing by 2. Construct 2^-(1-e) so the result
|
|||
|
|
is obtained with only one multiplication. */
|
|||
|
|
u.y *= ldexp(1.0, e-1);
|
|||
|
|
return(u.y);
|
|||
|
|
#else
|
|||
|
|
return(0.0);
|
|||
|
|
#endif
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
#ifdef DEC
|
|||
|
|
*q &= 0x807f; /* strip all exponent bits */
|
|||
|
|
*q |= (e & 0xff) << 7;
|
|||
|
|
#else
|
|||
|
|
*q &= 0x800f;
|
|||
|
|
*q |= (e & 0x7ff) << 4;
|
|||
|
|
#endif
|
|||
|
|
return(u.y);
|
|||
|
|
}
|
|||
|
|
}
|