Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : FLOORL.C

/* ceill()
* floorl()
* frexpl()
* ldexpl()
* fabsl()
*
* Floating point numeric utilities
*
*
*
* SYNOPSIS:
*
* long double x, y;
* long double ceill(), floorl(), frexpl(), ldexpl(), fabsl();
* int expnt, n;
*
* y = floorl(x);
* y = ceill(x);
* y = frexpl( x, &expnt );
* y = ldexpl( x, n );
* y = fabsl( x );
*
*
*
* DESCRIPTION:
*
* All four routines return a long double precision floating point
* result.
*
* floorl() returns the largest integer less than or equal to x.
* It truncates toward minus infinity.
*
* ceill() returns the smallest integer greater than or equal
* to x. It truncates toward plus infinity.
*
* frexpl() 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.
*
* ldexpl() multiplies x by 2**n.
*
* fabsl() returns the absolute value of its argument.
*
* These functions are part of the standard C run time library
* for some but not all C compilers. The ones supplied are
* written in C for 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.2: July, 1992
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


#include "mconf.h"
#define DENORMAL 0

#ifdef UNK
char *unkmsg = "ceill(), floorl(), frexpl(), ldexpl() must be rewritten!\n";
#endif

#ifdef IBMPC
#define EXPMSK 0x800f
#define MEXP 0x7ff
#define NBITS 64
#endif

#ifdef MIEEE
#define EXPMSK 0x800f
#define MEXP 0x7ff
#define NBITS 64
#endif

extern double MAXNUML;


long double fabsl(x)
long double x;
{
if( x < 0 )
return( -x );
else
return( x );
}



long double ceill(x)
long double x;
{
long double y;
long double floorl();

#ifdef UNK
mtherr( "ceill", DOMAIN );
return(0.0L);
#endif

y = floorl(x);
if( y < x )
y += 1.0L;
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,
};




long double floorl(x)
long double x;
{
unsigned short *p;
long double y;
int e, i;

#ifdef UNK
mtherr( "floor", DOMAIN );
return(0.0L);
#endif

y = x;
/* find the exponent (power of 2) */
#ifdef IBMPC
p = (unsigned short *)&y + 4;
e = (*p & 0x7fff) - 0x3fff;
p -= 4;
#endif

#ifdef MIEEE
p = (unsigned short *)&y;
e = (*p & 0x7fff) - 0x3fff;
p += 5;
#endif

if( e < 0 )
{
if( y < 0.0L )
return( -1.0L );
else
return( 0.0L );
}

e = (NBITS -1) - e;
/* clean out 16 bits at a time */
while( e >= 16 )
{
#ifdef IBMPC
*p++ = 0;
#endif

#ifdef MIEEE
*p-- = 0;
#endif
e -= 16;
}

/* clear the remaining bits */
if( e > 0 )
*p &= bmask[e];

if( (x < 0) && (y != x) )
y -= 1.0L;

return(y);
}



long double frexpl( x, pw2 )
long double x;
int *pw2;
{
long double y;
int i, k;
short *q;

y = x;

#ifdef UNK
mtherr( "frexp", DOMAIN );
return(0.0L);
#endif

/* find the exponent (power of 2) */
#ifdef IBMPC
q = (short *)&y + 4;
i = *q & 0x7fff;
#endif

#ifdef MIEEE
q = (short *)&y;
i = *q & 0x7fff;
#endif


if( i == 0 )
{
if( y == 0.0L )
{
*pw2 = 0;
return(0.0L);
}
/* Number is denormal or zero */
#if DENORMAL
/* Handle denormal number. */
do
{
y *= 2.0L;
i -= 1;
k = *q & 0x7fff;
}
while( (k == 0) && (i > -66) );
i = i + k;
#else
*pw2 = 0;
return(0.0L);
#endif /* DENORMAL */
}

*pw2 = i - 0x3ffe;
*q = 0x3ffe;
return( y );
}






long double ldexpl( x, pw2 )
long double x;
int pw2;
{
long double y;
unsigned short *q;
long e;

#ifdef UNK
mtherr( "ldexp", DOMAIN );
return(0.0L);
#endif

y = x;
#ifdef IBMPC
q = (unsigned short *)&y + 4;
#endif
#ifdef MIEEE
q = (unsigned short *)&y;
#endif
while( (e = (*q & 0x7fffL)) == 0 )
{
#if DENORMAL
if( y == 0.0L )
{
return( 0.0L );
}
/* Input is denormal. */
if( pw2 > 0 )
{
y *= 2.0L;
pw2 -= 1;
}
if( pw2 < 0 )
{
if( pw2 < -64 )
return(0.0L);
y *= 0.5;
pw2 += 1;
}
if( pw2 == 0 )
return(y);
#else
return( 0.0L );
#endif
}

e = e + pw2;

/* Handle overflow */
if( e > 0x7fffL )
{
return( MAXNUML );
}
*q &= 0x8000;
/* Handle denormalized results */
if( e < 1 )
{
#if DENORMAL
if( e < -64 )
return(0.0L);

#ifdef IBMPC
*(q-1) |= 0x8000;
#endif
#ifdef MIEEE
*(q+2) |= 0x8000;
#endif

while( e < 1 )
{
y *= 0.5;
e += 1;
}
e = 0;
#else
return(0.0L);
#endif
}

*q |= e & 0x7fff;
return(y);
}