[uclibc-ng.git] / libm / float / floorf.c

/*							ceilf()
 *							floorf()
 *							frexpf()
 *							ldexpf()
 *							signbitf()
 *							isnanf()
 *							isfinitef()
 *
 *	Single precision floating point numeric utilities
 *
 *
 *
 * SYNOPSIS:
 *
 * float x, y;
 * float ceilf(), floorf(), frexpf(), ldexpf();
 * int signbit(), isnan(), isfinite();
 * int expnt, n;
 *
 * y = floorf(x);
 * y = ceilf(x);
 * y = frexpf( x, &expnt );
 * y = ldexpf( x, n );
 * n = signbit(x);
 * n = isnan(x);
 * n = isfinite(x);
 *
 *
 *
 * DESCRIPTION:
 *
 * All four routines return a single precision floating point
 * result.
 *
 * sfloor() returns the largest integer less than or equal to x.
 * It truncates toward minus infinity.
 *
 * sceil() returns the smallest integer greater than or equal
 * to x.  It truncates toward plus infinity.
 *
 * sfrexp() 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.
 *
 * ldexpf() 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.
 */
\f

/*
Cephes Math Library Release 2.1:  December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


#include <math.h>
#ifdef DEC
#undef DENORMAL
#define DENORMAL 0
#endif

#ifdef UNK
#undef UNK
#if BIGENDIAN
#define MIEEE 1
#else
#define IBMPC 1
#endif
/*
char *unkmsg = "ceil(), floor(), frexp(), ldexp() must be rewritten!\n";
*/
#endif

#define EXPMSK 0x807f
#define MEXP 255
#define NBITS 24


extern float MAXNUMF; /* (2^24 - 1) * 2^103 */
#ifdef ANSIC
float floorf(float);
#else
float floorf();
#endif

float ceilf( float x )
{
float y;

#ifdef UNK
printf( "%s\n", unkmsg );
return(0.0);
#endif

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


float floorf( float x )
{
unsigned short *p;
union
  {
    float y;
    unsigned short i[2];
  } u;
int e;

#ifdef UNK
printf( "%s\n", unkmsg );
return(0.0);
#endif

u.y = x;
/* find the exponent (power of 2) */
#ifdef DEC
p = &u.i[0];
e = (( *p  >> 7) & 0377) - 0201;
p += 3;
#endif

#ifdef IBMPC
p = &u.i[1];
e = (( *p >> 7) & 0xff) - 0x7f;
p -= 1;
#endif

#ifdef MIEEE
p = &u.i[0];
e = (( *p >> 7) & 0xff) - 0x7f;
p += 1;
#endif

if( e < 0 )
	{
	if( u.y < 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);
}


float frexpf( float x, int *pw2 )
{
union
  {
    float y;
    unsigned short i[2];
  } u;
int i, k;
short *q;

u.y = x;

#ifdef UNK
printf( "%s\n", unkmsg );
return(0.0);
#endif

#ifdef IBMPC
q = &u.i[1];
#endif

#ifdef DEC
q = &u.i[0];
#endif

#ifdef MIEEE
q = &u.i[0];
#endif

/* find the exponent (power of 2) */

i  = ( *q >> 7) & 0xff;
if( i == 0 )
	{
	if( u.y == 0.0 )
		{
		*pw2 = 0;
		return(0.0);
		}
/* Number is denormal or zero */
#if DENORMAL
/* Handle denormal number. */
	do
		{
		u.y *= 2.0;
		i -= 1;
		k  = ( *q >> 7) & 0xff;
		}
	while( k == 0 );
	i = i + k;
#else
	*pw2 = 0;
	return( 0.0 );
#endif /* DENORMAL */
	}
i -= 0x7e;
*pw2 = i;
*q &= 0x807f;	/* strip all exponent bits */
*q |= 0x3f00;	/* mantissa between 0.5 and 1 */
return( u.y );
}


float ldexpf( float x, int pw2 )
{
union
  {
    float y;
    unsigned short i[2];
  } u;
short *q;
int e;

#ifdef UNK
printf( "%s\n", unkmsg );
return(0.0);
#endif

u.y = x;
#ifdef DEC
q = &u.i[0];
#endif

#ifdef IBMPC
q = &u.i[1];
#endif
#ifdef MIEEE
q = &u.i[0];
#endif
while( (e = ( *q >> 7) & 0xff) == 0 )
	{
	if( u.y == (float )0.0 )
		{
		return( 0.0 );
		}
/* Input is denormal. */
	if( pw2 > 0 )
		{
		u.y *= 2.0;
		pw2 -= 1;
		}
	if( pw2 < 0 )
		{
		if( pw2 < -24 )
			return( 0.0 );
		u.y *= 0.5;
		pw2 += 1;
		}
	if( pw2 == 0 )
		return(u.y);
	}

e += pw2;

/* Handle overflow */
if( e > MEXP )
	{
	return( MAXNUMF );
	}

*q &= 0x807f;

/* Handle denormalized results */
if( e < 1 )
	{
#if DENORMAL
	if( e < -24 )
		return( 0.0 );
	*q |= 0x80; /* Set LSB of exponent. */
	/* 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 *= ldexpf(1.0f, e - 1);
	return(u.y);
#else
	return( 0.0 );
#endif
	}
*q |= (e & 0xff) << 7;
return(u.y);
}


/* Return 1 if the sign bit of x is 1, else 0.  */

int signbitf(x)
float x;
{
union
	{
	float f;
	short s[4];
	int i;
	} u;

u.f = x;

if( sizeof(int) == 4 )
	{
#ifdef IBMPC
	return( u.i < 0 );
#endif
#ifdef DEC
	return( u.s[1] < 0 );
#endif
#ifdef MIEEE
	return( u.i < 0 );
#endif
	}
else
	{
#ifdef IBMPC
	return( u.s[1] < 0 );
#endif
#ifdef DEC
	return( u.s[1] < 0 );
#endif
#ifdef MIEEE
	return( u.s[0] < 0 );
#endif
	}
}


/* Return 1 if x is a number that is Not a Number, else return 0.  */

int isnanf(x)
float x;
{
#ifdef NANS
union
	{
	float f;
	unsigned short s[2];
	unsigned int i;
	} u;

u.f = x;

if( sizeof(int) == 4 )
	{
#ifdef IBMPC
	if( ((u.i & 0x7f800000) == 0x7f800000)
	    && ((u.i & 0x007fffff) != 0) )
		return 1;
#endif
#ifdef DEC
	if( (u.s[1] & 0x7f80) == 0)
		{
		if( (u.s[1] | u.s[0]) != 0 )
			return(1);
		}
#endif
#ifdef MIEEE
	if( ((u.i & 0x7f800000) == 0x7f800000)
	    && ((u.i & 0x007fffff) != 0) )
		return 1;
#endif
	return(0);
	}
else
	{ /* size int not 4 */
#ifdef IBMPC
	if( (u.s[1] & 0x7f80) == 0x7f80)
		{
		if( ((u.s[1] & 0x007f) | u.s[0]) != 0 )
			return(1);
		}
#endif
#ifdef DEC
	if( (u.s[1] & 0x7f80) == 0)
		{
		if( (u.s[1] | u.s[0]) != 0 )
			return(1);
		}
#endif
#ifdef MIEEE
	if( (u.s[0] & 0x7f80) == 0x7f80)
		{
		if( ((u.s[0] & 0x000f) | u.s[1]) != 0 )
			return(1);
		}
#endif
	return(0);
	} /* size int not 4 */

#else
/* No NANS.  */
return(0);
#endif
}


/* Return 1 if x is not infinite and is not a NaN.  */

int isfinitef(x)
float x;
{
#ifdef INFINITIES
union
	{
	float f;
	unsigned short s[2];
	unsigned int i;
	} u;

u.f = x;

if( sizeof(int) == 4 )
	{
#ifdef IBMPC
	if( (u.i & 0x7f800000) != 0x7f800000)
		return 1;
#endif
#ifdef DEC
	if( (u.s[1] & 0x7f80) == 0)
		{
		if( (u.s[1] | u.s[0]) != 0 )
			return(1);
		}
#endif
#ifdef MIEEE
	if( (u.i & 0x7f800000) != 0x7f800000)
		return 1;
#endif
	return(0);
	}
else
	{
#ifdef IBMPC
	if( (u.s[1] & 0x7f80) != 0x7f80)
		return 1;
#endif
#ifdef DEC
	if( (u.s[1] & 0x7f80) == 0)
		{
		if( (u.s[1] | u.s[0]) != 0 )
			return(1);
		}
#endif
#ifdef MIEEE
	if( (u.s[0] & 0x7f80) != 0x7f80)
		return 1;
#endif
	return(0);
	}
#else
/* No INFINITY.  */
return(1);
#endif
}
Commit	Line	Data
1077fa4d EA	1	/* ceilf()
	2	* floorf()
	3	* frexpf()
	4	* ldexpf()
	5	* signbitf()
	6	* isnanf()
	7	* isfinitef()
	8	*
	9	* Single precision floating point numeric utilities
	10	*
	11	*
	12	*
	13	* SYNOPSIS:
	14	*
	15	* float x, y;
	16	* float ceilf(), floorf(), frexpf(), ldexpf();
	17	* int signbit(), isnan(), isfinite();
	18	* int expnt, n;
	19	*
	20	* y = floorf(x);
	21	* y = ceilf(x);
	22	* y = frexpf( x, &expnt );
	23	* y = ldexpf( x, n );
	24	* n = signbit(x);
	25	* n = isnan(x);
	26	* n = isfinite(x);
	27	*
	28	*
	29	*
	30	* DESCRIPTION:
	31	*
	32	* All four routines return a single precision floating point
	33	* result.
	34	*
	35	* sfloor() returns the largest integer less than or equal to x.
	36	* It truncates toward minus infinity.
	37	*
	38	* sceil() returns the smallest integer greater than or equal
	39	* to x. It truncates toward plus infinity.
	40	*
	41	* sfrexp() extracts the exponent from x. It returns an integer
	42	* power of two to expnt and the significand between 0.5 and 1
	43	* to y. Thus x = y * 2**expn.
	44	*
	45	* ldexpf() multiplies x by 2**n.
	46	*
	47	* signbit(x) returns 1 if the sign bit of x is 1, else 0.
	48	*
	49	* These functions are part of the standard C run time library
	50	* for many but not all C compilers. The ones supplied are
	51	* written in C for either DEC or IEEE arithmetic. They should
	52	* be used only if your compiler library does not already have
	53	* them.
	54	*
	55	* The IEEE versions assume that denormal numbers are implemented
	56	* in the arithmetic. Some modifications will be required if
	57	* the arithmetic has abrupt rather than gradual underflow.
	58	*/
	59	\f
	60
	61	/*
	62	Cephes Math Library Release 2.1: December, 1988
	63	Copyright 1984, 1987, 1988 by Stephen L. Moshier
	64	Direct inquiries to 30 Frost Street, Cambridge, MA 02140
65	*/
66
67
68	#include <math.h>
69	#ifdef DEC
70	#undef DENORMAL
71	#define DENORMAL 0
72	#endif
73
74	#ifdef UNK
75	#undef UNK
76	#if BIGENDIAN
77	#define MIEEE 1
78	#else
79	#define IBMPC 1
80	#endif
81	/*
82	char *unkmsg = "ceil(), floor(), frexp(), ldexp() must be rewritten!\n";
83	*/
84	#endif
85
86	#define EXPMSK 0x807f
87	#define MEXP 255
88	#define NBITS 24
89
90
91	extern float MAXNUMF; /* (2^24 - 1) * 2^103 */
92	#ifdef ANSIC
93	float floorf(float);
94	#else
95	float floorf();
96	#endif
97
98	float ceilf( float x )
99	{
100	float y;
101
102	#ifdef UNK
103	printf( "%s\n", unkmsg );
104	return(0.0);
105	#endif
106
107	y = floorf( (float )x );
108	if( y < x )
109	y += 1.0;
110	return(y);
111	}
112
113
114
115
116	/* Bit clearing masks: */
117
118	static unsigned short bmask[] = {
119	0xffff,
120	0xfffe,
121	0xfffc,
122	0xfff8,
123	0xfff0,
124	0xffe0,
125	0xffc0,
126	0xff80,
127	0xff00,
128	0xfe00,
129	0xfc00,
130	0xf800,
131	0xf000,
132	0xe000,
133	0xc000,
134	0x8000,
135	0x0000,
136	};
137
138
139
140	float floorf( float x )
141	{
142	unsigned short *p;
143	union
144	{
145	float y;
146	unsigned short i[2];
147	} u;
148	int e;
149
150	#ifdef UNK
151	printf( "%s\n", unkmsg );
152	return(0.0);
153	#endif
154
155	u.y = x;
156	/* find the exponent (power of 2) */
157	#ifdef DEC
158	p = &u.i[0];
159	e = (( *p >> 7) & 0377) - 0201;
160	p += 3;
161	#endif
162
163	#ifdef IBMPC
164	p = &u.i[1];
165	e = (( *p >> 7) & 0xff) - 0x7f;
166	p -= 1;
167	#endif
168
169	#ifdef MIEEE
170	p = &u.i[0];
171	e = (( *p >> 7) & 0xff) - 0x7f;
172	p += 1;
173	#endif
174
175	if( e < 0 )
176	{
177	if( u.y < 0 )
178	return( -1.0 );
179	else
180	return( 0.0 );
181	}
182
183	e = (NBITS -1) - e;
184	/* clean out 16 bits at a time */
185	while( e >= 16 )
186	{
187	#ifdef IBMPC
188	*p++ = 0;
189	#endif
190
191	#ifdef DEC
192	*p-- = 0;
193	#endif
194
195	#ifdef MIEEE
196	*p-- = 0;
197	#endif
198	e -= 16;
199	}
200
201	/* clear the remaining bits */
202	if( e > 0 )
203	*p &= bmask[e];
204
205	if( (x < 0) && (u.y != x) )
206	u.y -= 1.0;
207
208	return(u.y);
209	}
210
211
212
213	float frexpf( float x, int *pw2 )
214	{
215	union
216	{
217	float y;
218	unsigned short i[2];
219	} u;
220	int i, k;
221	short *q;
222
223	u.y = x;
224
225	#ifdef UNK
226	printf( "%s\n", unkmsg );
227	return(0.0);
228	#endif
229
230	#ifdef IBMPC
231	q = &u.i[1];
232	#endif
233
234	#ifdef DEC
235	q = &u.i[0];
236	#endif
237
238	#ifdef MIEEE
239	q = &u.i[0];
240	#endif
241
242	/* find the exponent (power of 2) */
243
244	i = ( *q >> 7) & 0xff;
245	if( i == 0 )
246	{
247	if( u.y == 0.0 )
248	{
249	*pw2 = 0;
250	return(0.0);
251	}
252	/* Number is denormal or zero */
253	#if DENORMAL
254	/* Handle denormal number. */
255	do
256	{
257	u.y *= 2.0;
258	i -= 1;
259	k = ( *q >> 7) & 0xff;
260	}
261	while( k == 0 );
262	i = i + k;
263	#else
264	*pw2 = 0;
265	return( 0.0 );
266	#endif /* DENORMAL */
267	}
268	i -= 0x7e;
269	*pw2 = i;
270	q &= 0x807f; / strip all exponent bits */
271	q \|= 0x3f00; / mantissa between 0.5 and 1 */
272	return( u.y );
273	}
274
275
276
277
278
279	float ldexpf( float x, int pw2 )
280	{
281	union
282	{
283	float y;
284	unsigned short i[2];
285	} u;
286	short *q;
287	int e;
288
289	#ifdef UNK
290	printf( "%s\n", unkmsg );
291	return(0.0);
292	#endif
293
294	u.y = x;
295	#ifdef DEC
296	q = &u.i[0];
297	#endif
298
299	#ifdef IBMPC
300	q = &u.i[1];
301	#endif
302	#ifdef MIEEE
303	q = &u.i[0];
304	#endif
305	while( (e = ( *q >> 7) & 0xff) == 0 )
306	{
307	if( u.y == (float )0.0 )
308	{
309	return( 0.0 );
310	}
311	/* Input is denormal. */
312	if( pw2 > 0 )
313	{
314	u.y *= 2.0;
315	pw2 -= 1;
316	}
317	if( pw2 < 0 )
318	{
319	if( pw2 < -24 )
320	return( 0.0 );
321	u.y *= 0.5;
322	pw2 += 1;
323	}
324	if( pw2 == 0 )
325	return(u.y);
326	}
327
328	e += pw2;
329
330	/* Handle overflow */
331	if( e > MEXP )
332	{
333	return( MAXNUMF );
334	}
335
336	*q &= 0x807f;
337
338	/* Handle denormalized results */
339	if( e < 1 )
340	{
341	#if DENORMAL
342	if( e < -24 )
343	return( 0.0 );
344	q \|= 0x80; / Set LSB of exponent. */
345	/* For denormals, significant bits may be lost even
346	when dividing by 2. Construct 2^-(1-e) so the result
347	is obtained with only one multiplication. */
348	u.y *= ldexpf(1.0f, e - 1);
349	return(u.y);
350	#else
351	return( 0.0 );
352	#endif
353	}
354	*q \|= (e & 0xff) << 7;
355	return(u.y);
356	}
357
358
359	/* Return 1 if the sign bit of x is 1, else 0. */
360
361	int signbitf(x)
362	float x;
363	{
364	union
365	{
366	float f;
367	short s[4];
368	int i;
369	} u;
370
371	u.f = x;
372
373	if( sizeof(int) == 4 )
374	{
375	#ifdef IBMPC
376	return( u.i < 0 );
377	#endif
378	#ifdef DEC
379	return( u.s[1] < 0 );
380	#endif
381	#ifdef MIEEE
382	return( u.i < 0 );
383	#endif
384	}
385	else
386	{
387	#ifdef IBMPC
388	return( u.s[1] < 0 );
389	#endif
390	#ifdef DEC
391	return( u.s[1] < 0 );
392	#endif
393	#ifdef MIEEE
394	return( u.s[0] < 0 );
395	#endif
396	}
397	}
398
399
400	/* Return 1 if x is a number that is Not a Number, else return 0. */
401
402	int isnanf(x)
403	float x;
404	{
405	#ifdef NANS
406	union
407	{
408	float f;
409	unsigned short s[2];
410	unsigned int i;
411	} u;
412
413	u.f = x;
414
415	if( sizeof(int) == 4 )
416	{
417	#ifdef IBMPC
418	if( ((u.i & 0x7f800000) == 0x7f800000)
419	&& ((u.i & 0x007fffff) != 0) )
420	return 1;
421	#endif
422	#ifdef DEC
423	if( (u.s[1] & 0x7f80) == 0)
424	{
425	if( (u.s[1] \| u.s[0]) != 0 )
426	return(1);
427	}
428	#endif
429	#ifdef MIEEE
430	if( ((u.i & 0x7f800000) == 0x7f800000)
431	&& ((u.i & 0x007fffff) != 0) )
432	return 1;
433	#endif
434	return(0);
435	}
436	else
437	{ /* size int not 4 */
438	#ifdef IBMPC
439	if( (u.s[1] & 0x7f80) == 0x7f80)
440	{
441	if( ((u.s[1] & 0x007f) \| u.s[0]) != 0 )
442	return(1);
443	}
444	#endif
445	#ifdef DEC
446	if( (u.s[1] & 0x7f80) == 0)
447	{
448	if( (u.s[1] \| u.s[0]) != 0 )
449	return(1);
450	}
451	#endif
452	#ifdef MIEEE
453	if( (u.s[0] & 0x7f80) == 0x7f80)
454	{
455	if( ((u.s[0] & 0x000f) \| u.s[1]) != 0 )
456	return(1);
457	}
458	#endif
459	return(0);
460	} /* size int not 4 */
461
462	#else
463	/* No NANS. */
464	return(0);
465	#endif
466	}
467
468
469	/* Return 1 if x is not infinite and is not a NaN. */
470
471	int isfinitef(x)
472	float x;
473	{
474	#ifdef INFINITIES
475	union
476	{
477	float f;
478	unsigned short s[2];
479	unsigned int i;
480	} u;
481
482	u.f = x;
483
484	if( sizeof(int) == 4 )
485	{
486	#ifdef IBMPC
487	if( (u.i & 0x7f800000) != 0x7f800000)
488	return 1;
489	#endif
490	#ifdef DEC
491	if( (u.s[1] & 0x7f80) == 0)
492	{
493	if( (u.s[1] \| u.s[0]) != 0 )
494	return(1);
495	}
496	#endif
497	#ifdef MIEEE
498	if( (u.i & 0x7f800000) != 0x7f800000)
499	return 1;
500	#endif
501	return(0);
502	}
503	else
504	{
505	#ifdef IBMPC
506	if( (u.s[1] & 0x7f80) != 0x7f80)
507	return 1;
508	#endif
509	#ifdef DEC
510	if( (u.s[1] & 0x7f80) == 0)
511	{
512	if( (u.s[1] \| u.s[0]) != 0 )
513	return(1);
514	}
515	#endif
516	#ifdef MIEEE
517	if( (u.s[0] & 0x7f80) != 0x7f80)
518	return 1;
519	#endif
520	return(0);
521	}
522	#else
523	/* No INFINITY. */
524	return(1);
525	#endif
526	}