import from github

agbcc/libc/math/fdlibm.h

@@ -0,0 +1,261 @@
+
+/* @(#)fdlibm.h 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* CYGNUS LOCAL: Include files.  */
+#include <math.h>
+#include <machine/ieeefp.h>
+
+/* CYGNUS LOCAL: Default to XOPEN_MODE.  */
+#define _XOPEN_MODE
+
+#ifdef __STDC__
+#define	__P(p)	p
+#else
+#define	__P(p)	()
+#endif
+
+#define	HUGE	((float)3.40282346638528860e+38)
+
+/* 
+ * set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
+ * (one may replace the following line by "#include <values.h>")
+ */
+
+#define X_TLOSS		1.41484755040568800000e+16 
+
+/* Functions that are not documented, and are not in <math.h>.  */
+
+extern double logb __P((double));
+#ifdef _SCALB_INT
+extern double scalb __P((double, int));
+#else
+extern double scalb __P((double, double));
+#endif
+extern double significand __P((double));
+
+/* ieee style elementary functions */
+extern double __ieee754_sqrt __P((double));			
+extern double __ieee754_acos __P((double));			
+extern double __ieee754_acosh __P((double));			
+extern double __ieee754_log __P((double));			
+extern double __ieee754_atanh __P((double));			
+extern double __ieee754_asin __P((double));			
+extern double __ieee754_atan2 __P((double,double));			
+extern double __ieee754_exp __P((double));
+extern double __ieee754_cosh __P((double));
+extern double __ieee754_fmod __P((double,double));
+extern double __ieee754_pow __P((double,double));
+extern double __ieee754_lgamma_r __P((double,int *));
+extern double __ieee754_gamma_r __P((double,int *));
+extern double __ieee754_log10 __P((double));
+extern double __ieee754_sinh __P((double));
+extern double __ieee754_hypot __P((double,double));
+extern double __ieee754_j0 __P((double));
+extern double __ieee754_j1 __P((double));
+extern double __ieee754_y0 __P((double));
+extern double __ieee754_y1 __P((double));
+extern double __ieee754_jn __P((int,double));
+extern double __ieee754_yn __P((int,double));
+extern double __ieee754_remainder __P((double,double));
+extern __int32_t __ieee754_rem_pio2 __P((double,double*));
+#ifdef _SCALB_INT
+extern double __ieee754_scalb __P((double,int));
+#else
+extern double __ieee754_scalb __P((double,double));
+#endif
+
+/* fdlibm kernel function */
+extern double __kernel_standard __P((double,double,int));
+extern double __kernel_sin __P((double,double,int));
+extern double __kernel_cos __P((double,double));
+extern double __kernel_tan __P((double,double,int));
+extern int    __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*));
+
+/* Undocumented float functions.  */
+extern float logbf __P((float));
+#ifdef _SCALB_INT
+extern float scalbf __P((float, int));
+#else
+extern float scalbf __P((float, float));
+#endif
+extern float significandf __P((float));
+
+/* ieee style elementary float functions */
+extern float __ieee754_sqrtf __P((float));			
+extern float __ieee754_acosf __P((float));			
+extern float __ieee754_acoshf __P((float));			
+extern float __ieee754_logf __P((float));			
+extern float __ieee754_atanhf __P((float));			
+extern float __ieee754_asinf __P((float));			
+extern float __ieee754_atan2f __P((float,float));			
+extern float __ieee754_expf __P((float));
+extern float __ieee754_coshf __P((float));
+extern float __ieee754_fmodf __P((float,float));
+extern float __ieee754_powf __P((float,float));
+extern float __ieee754_lgammaf_r __P((float,int *));
+extern float __ieee754_gammaf_r __P((float,int *));
+extern float __ieee754_log10f __P((float));
+extern float __ieee754_sinhf __P((float));
+extern float __ieee754_hypotf __P((float,float));
+extern float __ieee754_j0f __P((float));
+extern float __ieee754_j1f __P((float));
+extern float __ieee754_y0f __P((float));
+extern float __ieee754_y1f __P((float));
+extern float __ieee754_jnf __P((int,float));
+extern float __ieee754_ynf __P((int,float));
+extern float __ieee754_remainderf __P((float,float));
+extern __int32_t __ieee754_rem_pio2f __P((float,float*));
+#ifdef _SCALB_INT
+extern float __ieee754_scalbf __P((float,int));
+#else
+extern float __ieee754_scalbf __P((float,float));
+#endif
+
+/* float versions of fdlibm kernel functions */
+extern float __kernel_sinf __P((float,float,int));
+extern float __kernel_cosf __P((float,float));
+extern float __kernel_tanf __P((float,float,int));
+extern int   __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*));
+
+/* The original code used statements like
+	n0 = ((*(int*)&one)>>29)^1;		* index of high word *
+	ix0 = *(n0+(int*)&x);			* high word of x *
+	ix1 = *((1-n0)+(int*)&x);		* low word of x *
+   to dig two 32 bit words out of the 64 bit IEEE floating point
+   value.  That is non-ANSI, and, moreover, the gcc instruction
+   scheduler gets it wrong.  We instead use the following macros.
+   Unlike the original code, we determine the endianness at compile
+   time, not at run time; I don't see much benefit to selecting
+   endianness at run time.  */
+
+#ifndef __IEEE_BIG_ENDIAN
+#ifndef __IEEE_LITTLE_ENDIAN
+ #error Must define endianness
+#endif
+#endif
+
+/* A union which permits us to convert between a double and two 32 bit
+   ints.  */
+
+#ifdef __IEEE_BIG_ENDIAN
+
+typedef union 
+{
+  double value;
+  struct 
+  {
+    __uint32_t msw;
+    __uint32_t lsw;
+  } parts;
+} ieee_double_shape_type;
+
+#endif
+
+#ifdef __IEEE_LITTLE_ENDIAN
+
+typedef union 
+{
+  double value;
+  struct 
+  {
+    __uint32_t lsw;
+    __uint32_t msw;
+  } parts;
+} ieee_double_shape_type;
+
+#endif
+
+/* Get two 32 bit ints from a double.  */
+
+#define EXTRACT_WORDS(ix0,ix1,d)				\
+do {								\
+  ieee_double_shape_type ew_u;					\
+  ew_u.value = (d);						\
+  (ix0) = ew_u.parts.msw;					\
+  (ix1) = ew_u.parts.lsw;					\
+} while (0)
+
+/* Get the more significant 32 bit int from a double.  */
+
+#define GET_HIGH_WORD(i,d)					\
+do {								\
+  ieee_double_shape_type gh_u;					\
+  gh_u.value = (d);						\
+  (i) = gh_u.parts.msw;						\
+} while (0)
+
+/* Get the less significant 32 bit int from a double.  */
+
+#define GET_LOW_WORD(i,d)					\
+do {								\
+  ieee_double_shape_type gl_u;					\
+  gl_u.value = (d);						\
+  (i) = gl_u.parts.lsw;						\
+} while (0)
+
+/* Set a double from two 32 bit ints.  */
+
+#define INSERT_WORDS(d,ix0,ix1)					\
+do {								\
+  ieee_double_shape_type iw_u;					\
+  iw_u.parts.msw = (ix0);					\
+  iw_u.parts.lsw = (ix1);					\
+  (d) = iw_u.value;						\
+} while (0)
+
+/* Set the more significant 32 bits of a double from an int.  */
+
+#define SET_HIGH_WORD(d,v)					\
+do {								\
+  ieee_double_shape_type sh_u;					\
+  sh_u.value = (d);						\
+  sh_u.parts.msw = (v);						\
+  (d) = sh_u.value;						\
+} while (0)
+
+/* Set the less significant 32 bits of a double from an int.  */
+
+#define SET_LOW_WORD(d,v)					\
+do {								\
+  ieee_double_shape_type sl_u;					\
+  sl_u.value = (d);						\
+  sl_u.parts.lsw = (v);						\
+  (d) = sl_u.value;						\
+} while (0)
+
+/* A union which permits us to convert between a float and a 32 bit
+   int.  */
+
+typedef union
+{
+  float value;
+  __uint32_t word;
+} ieee_float_shape_type;
+
+/* Get a 32 bit int from a float.  */
+
+#define GET_FLOAT_WORD(i,d)					\
+do {								\
+  ieee_float_shape_type gf_u;					\
+  gf_u.value = (d);						\
+  (i) = gf_u.word;						\
+} while (0)
+
+/* Set a float from a 32 bit int.  */
+
+#define SET_FLOAT_WORD(d,i)					\
+do {								\
+  ieee_float_shape_type sf_u;					\
+  sf_u.word = (i);						\
+  (d) = sf_u.value;						\
+} while (0)

agbcc/libc/math/s_copysign.c

@@ -0,0 +1,82 @@
+
+/* @(#)s_copysign.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+FUNCTION
+<<copysign>>, <<copysignf>>---sign of <[y]>, magnitude of <[x]>
+
+INDEX
+	copysign
+INDEX
+	copysignf
+
+ANSI_SYNOPSIS
+	#include <math.h>
+	double copysign (double <[x]>, double <[y]>);
+	float copysignf (float <[x]>, float <[y]>);
+
+TRAD_SYNOPSIS
+	#include <math.h>
+	double copysign (<[x]>, <[y]>)
+	double <[x]>;
+	double <[y]>;
+
+	float copysignf (<[x]>, <[y]>)
+	float <[x]>;
+	float <[y]>;
+
+DESCRIPTION
+<<copysign>> constructs a number with the magnitude (absolute value)
+of its first argument, <[x]>, and the sign of its second argument,
+<[y]>.
+
+<<copysignf>> does the same thing; the two functions differ only in
+the type of their arguments and result.
+
+RETURNS
+<<copysign>> returns a <<double>> with the magnitude of
+<[x]> and the sign of <[y]>.
+<<copysignf>> returns a <<float>> with the magnitude of
+<[x]> and the sign of <[y]>.
+
+PORTABILITY
+<<copysign>> is not required by either ANSI C or the System V Interface
+Definition (Issue 2).
+
+*/
+
+/*
+ * copysign(double x, double y)
+ * copysign(x,y) returns a value with the magnitude of x and
+ * with the sign bit of y.
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	double copysign(double x, double y)
+#else
+	double copysign(x,y)
+	double x,y;
+#endif
+{
+	__uint32_t hx,hy;
+	GET_HIGH_WORD(hx,x);
+	GET_HIGH_WORD(hy,y);
+	SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
+        return x;
+}
+
+#endif /* _DOUBLE_IS_32BITS */

agbcc/libc/math/s_finite.c

@@ -0,0 +1,35 @@
+
+/* @(#)s_finite.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * finite(x) returns 1 is x is finite, else 0;
+ * no branching!
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	int finite(double x)
+#else
+	int finite(x)
+	double x;
+#endif
+{
+	__int32_t hx;
+	GET_HIGH_WORD(hx,x);
+	return  (int)((__uint32_t)((hx&0x7fffffff)-0x7ff00000)>>31);
+}
+
+#endif /* _DOUBLE_IS_32BITS */

agbcc/libc/math/s_frexp.c

@@ -0,0 +1,114 @@
+
+/* @(#)s_frexp.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+FUNCTION
+       <<frexp>>, <<frexpf>>---split floating-point number
+INDEX
+	frexp
+INDEX
+	frexpf
+
+ANSI_SYNOPSIS
+	#include <math.h>
+        double frexp(double <[val]>, int *<[exp]>);
+        float frexpf(float <[val]>, int *<[exp]>);
+
+TRAD_SYNOPSIS
+	#include <math.h>
+        double frexp(<[val]>, <[exp]>)
+        double <[val]>;
+        int *<[exp]>;
+
+        float frexpf(<[val]>, <[exp]>)
+        float <[val]>;
+        int *<[exp]>;
+
+
+DESCRIPTION
+	All non zero, normal numbers can be described as <[m]> * 2**<[p]>.
+	<<frexp>> represents the double <[val]> as a mantissa <[m]>
+	and a power of two <[p]>. The resulting mantissa will always
+	be greater than or equal to <<0.5>>, and less than <<1.0>> (as
+	long as <[val]> is nonzero). The power of two will be stored
+	in <<*>><[exp]>. 
+
+@ifinfo
+<[m]> and <[p]> are calculated so that
+<[val]> is <[m]> times <<2>> to the power <[p]>.
+@end ifinfo
+@tex
+<[m]> and <[p]> are calculated so that
+$ val = m \times 2^p $.
+@end tex
+
+<<frexpf>> is identical, other than taking and returning
+floats rather than doubles.
+
+RETURNS
+<<frexp>> returns the mantissa <[m]>. If <[val]> is <<0>>, infinity,
+or Nan, <<frexp>> will set <<*>><[exp]> to <<0>> and return <[val]>.
+
+PORTABILITY
+<<frexp>> is ANSI.
+<<frexpf>> is an extension.
+
+
+*/
+
+/*
+ * for non-zero x 
+ *	x = frexp(arg,&exp);
+ * return a double fp quantity x such that 0.5 <= |x| <1.0
+ * and the corresponding binary exponent "exp". That is
+ *	arg = x*2^exp.
+ * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg 
+ * with *exp=0. 
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+static const double
+#else
+static double
+#endif
+two54 =  1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
+
+#ifdef __STDC__
+	double frexp(double x, int *eptr)
+#else
+	double frexp(x, eptr)
+	double x; int *eptr;
+#endif
+{
+	__int32_t hx, ix, lx;
+	EXTRACT_WORDS(hx,lx,x);
+	ix = 0x7fffffff&hx;
+	*eptr = 0;
+	if(ix>=0x7ff00000||((ix|lx)==0)) return x;	/* 0,inf,nan */
+	if (ix<0x00100000) {		/* subnormal */
+	    x *= two54;
+	    GET_HIGH_WORD(hx,x);
+	    ix = hx&0x7fffffff;
+	    *eptr = -54;
+	}
+	*eptr += (ix>>20)-1022;
+	hx = (hx&0x800fffff)|0x3fe00000;
+	SET_HIGH_WORD(x,hx);
+	return x;
+}
+
+#endif /* _DOUBLE_IS_32BITS */

agbcc/libc/math/s_infconst.c

@@ -0,0 +1,15 @@
+/* Infinity as a constant value.   This is used for HUGE_VAL.
+ * Added by Cygnus Support.
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+#ifdef __IEEE_BIG_ENDIAN
+const union __dmath __infinity = { 0x7ff00000, 0 };
+#else
+const union __dmath __infinity = { 0, 0x7ff00000 };
+#endif
+#else /* defined (_DOUBLE_IS_32BITS) */
+const union __dmath __infinity = { 0x7f800000, 0 };
+#endif /* defined (_DOUBLE_IS_32BITS) */

agbcc/libc/math/s_isinf.c

@@ -0,0 +1,26 @@
+/*
+ * isinf(x) returns 1 if x is infinity, else 0;
+ * no branching!
+ * Added by Cygnus Support.
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	int isinf(double x)
+#else
+	int isinf(x)
+	double x;
+#endif
+{
+	__int32_t hx,lx;
+	EXTRACT_WORDS(hx,lx,x);
+	hx &= 0x7fffffff;
+	hx |= (__uint32_t)(lx|(-lx))>>31;	
+	hx = 0x7ff00000 - hx;
+	return 1 - (int)((__uint32_t)(hx|(-hx))>>31);
+}
+
+#endif /* _DOUBLE_IS_32BITS */

agbcc/libc/math/s_isnan.c

@@ -0,0 +1,122 @@
+
+/* @(#)s_isnan.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+FUNCTION
+	<<isnan>>,<<isnanf>>,<<isinf>>,<<isinff>>,<<finite>>,<<finitef>>---test for exceptional numbers
+
+INDEX
+	isnan
+INDEX
+	isinf
+INDEX
+	finite
+
+INDEX
+	isnanf
+INDEX
+	isinff
+INDEX
+	finitef
+
+ANSI_SYNOPSIS
+	#include <ieeefp.h>
+	int isnan(double <[arg]>);
+	int isinf(double <[arg]>);
+	int finite(double <[arg]>);
+	int isnanf(float <[arg]>);
+	int isinff(float <[arg]>);
+	int finitef(float <[arg]>);
+
+TRAD_SYNOPSIS
+	#include <ieeefp.h>
+	int isnan(<[arg]>)
+	double <[arg]>;
+	int isinf(<[arg]>)
+	double <[arg]>;
+	int finite(<[arg]>);
+	double <[arg]>;
+	int isnanf(<[arg]>);
+	float <[arg]>;
+	int isinff(<[arg]>);
+	float <[arg]>;
+	int finitef(<[arg]>);
+	float <[arg]>;
+
+
+DESCRIPTION
+	These functions provide information on the floating point
+	argument supplied.
+
+	There are five major number formats - 
+	o+
+	o zero
+	 a number which contains all zero bits.
+	o subnormal
+	 Is used to represent  number with a zero exponent, but a non zero fraction.
+	 o normal
+	  A number with an exponent, and a fraction
+     	o infinity
+	  A number with an all 1's exponent and a zero fraction.
+	o NAN
+	  A number with an all 1's exponent and a non zero fraction.
+
+	o-
+
+	<<isnan>> returns 1 if the argument is a nan. <<isinf>>
+	returns 1 if the argument is infinity.  <<finite>> returns 1 if the
+	argument is zero, subnormal or normal.
+	
+	The <<isnanf>>, <<isinff>> and <<finitef>> perform the same
+	operations as their <<isnan>>, <<isinf>> and <<finite>>
+	counterparts, but on single precision floating point numbers.
+
+QUICKREF
+	isnan - pure
+QUICKREF
+	isinf - pure
+QUICKREF
+	finite - pure
+QUICKREF
+	isnan - pure
+QUICKREF
+	isinf - pure
+QUICKREF
+	finite - pure
+*/
+
+/*
+ * isnan(x) returns 1 is x is nan, else 0;
+ * no branching!
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	int isnan(double x)
+#else
+	int isnan(x)
+	double x;
+#endif
+{
+	__int32_t hx,lx;
+	EXTRACT_WORDS(hx,lx,x);
+	hx &= 0x7fffffff;
+	hx |= (__uint32_t)(lx|(-lx))>>31;	
+	hx = 0x7ff00000 - hx;
+	return (int)(((__uint32_t)(hx))>>31);
+}
+
+#endif /* _DOUBLE_IS_32BITS */

agbcc/libc/math/s_ldexp.c

@@ -0,0 +1,81 @@
+
+/* @(#)s_ldexp.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+FUNCTION
+       <<ldexp>>, <<ldexpf>>---load exponent
+
+INDEX
+	ldexp
+INDEX
+	ldexpf
+
+ANSI_SYNOPSIS
+       #include <math.h>
+       double ldexp(double <[val]>, int <[exp]>);
+       float ldexpf(float <[val]>, int <[exp]>);
+
+TRAD_SYNOPSIS
+       #include <math.h>
+
+       double ldexp(<[val]>, <[exp]>)
+              double <[val]>;
+              int <[exp]>;
+
+       float ldexpf(<[val]>, <[exp]>)
+              float <[val]>;
+              int <[exp]>;
+
+
+DESCRIPTION
+<<ldexp>> calculates the value 
+@ifinfo
+<[val]> times 2 to the power <[exp]>.
+@end ifinfo
+@tex
+$val\times 2^{exp}$.
+@end tex
+<<ldexpf>> is identical, save that it takes and returns <<float>>
+rather than <<double>> values.
+
+RETURNS
+<<ldexp>> returns the calculated value.
+
+Underflow and overflow both set <<errno>> to <<ERANGE>>.
+On underflow, <<ldexp>> and <<ldexpf>> return 0.0.
+On overflow, <<ldexp>> returns plus or minus <<HUGE_VAL>>.
+
+PORTABILITY
+<<ldexp>> is ANSI, <<ldexpf>> is an extension.
+              
+*/   
+
+#include "fdlibm.h"
+#include <errno.h>
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	double ldexp(double value, int exp)
+#else
+	double ldexp(value, exp)
+	double value; int exp;
+#endif
+{
+	if(!finite(value)||value==0.0) return value;
+	value = scalbn(value,exp);
+	if(!finite(value)||value==0.0) errno = ERANGE;
+	return value;
+}
+
+#endif /* _DOUBLE_IS_32BITS */

agbcc/libc/math/s_modf.c

@@ -0,0 +1,131 @@
+
+/* @(#)s_modf.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+FUNCTION
+       <<modf>>, <<modff>>---split fractional and integer parts
+
+INDEX
+	modf
+INDEX
+	modff
+
+ANSI_SYNOPSIS
+	#include <math.h>
+	double modf(double <[val]>, double *<[ipart]>);
+        float modff(float <[val]>, float *<[ipart]>);
+
+TRAD_SYNOPSIS
+	#include <math.h>
+	double modf(<[val]>, <[ipart]>)
+        double <[val]>;
+        double *<[ipart]>;
+
+	float modff(<[val]>, <[ipart]>)
+	float <[val]>;
+        float *<[ipart]>;
+
+DESCRIPTION
+	<<modf>> splits the double <[val]> apart into an integer part
+	and a fractional part, returning the fractional part and
+	storing the integer part in <<*<[ipart]>>>.  No rounding
+	whatsoever is done; the sum of the integer and fractional
+	parts is guaranteed to be exactly  equal to <[val]>.   That
+	is, if . <[realpart]> = modf(<[val]>, &<[intpart]>); then
+	`<<<[realpart]>+<[intpart]>>>' is the same as <[val]>.
+	<<modff>> is identical, save that it takes and returns
+	<<float>> rather than <<double>> values. 
+
+RETURNS
+	The fractional part is returned.  Each result has the same
+	sign as the supplied argument <[val]>.
+
+PORTABILITY
+	<<modf>> is ANSI C. <<modff>> is an extension.
+
+QUICKREF
+	modf  ansi pure 
+	modff - pure
+
+*/
+
+/*
+ * modf(double x, double *iptr) 
+ * return fraction part of x, and return x's integral part in *iptr.
+ * Method:
+ *	Bit twiddling.
+ *
+ * Exception:
+ *	No exception.
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+static const double one = 1.0;
+#else
+static double one = 1.0;
+#endif
+
+#ifdef __STDC__
+	double modf(double x, double *iptr)
+#else
+	double modf(x, iptr)
+	double x,*iptr;
+#endif
+{
+	__int32_t i0,i1,j0;
+	__uint32_t i;
+	EXTRACT_WORDS(i0,i1,x);
+	j0 = ((i0>>20)&0x7ff)-0x3ff;	/* exponent of x */
+	if(j0<20) {			/* integer part in high x */
+	    if(j0<0) {			/* |x|<1 */
+	        INSERT_WORDS(*iptr,i0&0x80000000,0);	/* *iptr = +-0 */
+		return x;
+	    } else {
+		i = (0x000fffff)>>j0;
+		if(((i0&i)|i1)==0) {		/* x is integral */
+		    __uint32_t high;
+		    *iptr = x;
+		    GET_HIGH_WORD(high,x);
+		    INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
+		    return x;
+		} else {
+		    INSERT_WORDS(*iptr,i0&(~i),0);
+		    return x - *iptr;
+		}
+	    }
+	} else if (j0>51) {		/* no fraction part */
+	    __uint32_t high;
+	    *iptr = x*one;
+	    GET_HIGH_WORD(high,x);
+	    INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
+	    return x;
+	} else {			/* fraction part in low x */
+	    i = ((__uint32_t)(0xffffffff))>>(j0-20);
+	    if((i1&i)==0) { 		/* x is integral */
+	        __uint32_t high;
+		*iptr = x;
+		GET_HIGH_WORD(high,x);
+		INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
+		return x;
+	    } else {
+	        INSERT_WORDS(*iptr,i0,i1&(~i));
+		return x - *iptr;
+	    }
+	}
+}
+
+#endif /* _DOUBLE_IS_32BITS */

agbcc/libc/math/s_scalbn.c

@@ -0,0 +1,103 @@
+
+/* @(#)s_scalbn.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+FUNCTION
+<<scalbn>>, <<scalbnf>>---scale by integer
+INDEX
+	scalbn
+INDEX
+	scalbnf
+
+ANSI_SYNOPSIS
+	#include <math.h>
+	double scalbn(double <[x]>, int <[y]>);
+	float scalbnf(float <[x]>, int <[y]>);
+
+TRAD_SYNOPSIS
+	#include <math.h>
+	double scalbn(<[x]>,<[y]>)
+	double <[x]>;
+	int <[y]>;
+	float scalbnf(<[x]>,<[y]>)
+	float <[x]>;
+	int <[y]>;
+
+DESCRIPTION
+<<scalbn>> and <<scalbnf>> scale <[x]> by <[n]>, returning <[x]> times
+2 to the power <[n]>.  The result is computed by manipulating the
+exponent, rather than by actually performing an exponentiation or
+multiplication.
+
+RETURNS
+<[x]> times 2 to the power <[n]>.
+
+PORTABILITY
+Neither <<scalbn>> nor <<scalbnf>> is required by ANSI C or by the System V
+Interface Definition (Issue 2).
+
+*/
+
+/* 
+ * scalbn (double x, int n)
+ * scalbn(x,n) returns x* 2**n  computed by  exponent  
+ * manipulation rather than by actually performing an 
+ * exponentiation or a multiplication.
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+static const double
+#else
+static double
+#endif
+two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
+huge   = 1.0e+300,
+tiny   = 1.0e-300;
+
+#ifdef __STDC__
+	double scalbn (double x, int n)
+#else
+	double scalbn (x,n)
+	double x; int n;
+#endif
+{
+	__int32_t  k,hx,lx;
+	EXTRACT_WORDS(hx,lx,x);
+        k = (hx&0x7ff00000)>>20;		/* extract exponent */
+        if (k==0) {				/* 0 or subnormal x */
+            if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
+	    x *= two54; 
+	    GET_HIGH_WORD(hx,x);
+	    k = ((hx&0x7ff00000)>>20) - 54; 
+            if (n< -50000) return tiny*x; 	/*underflow*/
+	    }
+        if (k==0x7ff) return x+x;		/* NaN or Inf */
+        k = k+n; 
+        if (k >  0x7fe) return huge*copysign(huge,x); /* overflow  */
+        if (k > 0) 				/* normal result */
+	    {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
+        if (k <= -54)
+            if (n > 50000) 	/* in case integer overflow in n+k */
+		return huge*copysign(huge,x);	/*overflow*/
+	    else return tiny*copysign(tiny,x); 	/*underflow*/
+        k += 54;				/* subnormal result */
+	SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
+        return x*twom54;
+}
+
+#endif /* _DOUBLE_IS_32BITS */

agbcc/libc/math/sf_copysign.c

@@ -0,0 +1,50 @@
+/* sf_copysign.c -- float version of s_copysign.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * copysignf(float x, float y)
+ * copysignf(x,y) returns a value with the magnitude of x and
+ * with the sign bit of y.
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+	float copysignf(float x, float y)
+#else
+	float copysignf(x,y)
+	float x,y;
+#endif
+{
+	__uint32_t ix,iy;
+	GET_FLOAT_WORD(ix,x);
+	GET_FLOAT_WORD(iy,y);
+	SET_FLOAT_WORD(x,(ix&0x7fffffff)|(iy&0x80000000));
+        return x;
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	double copysign(double x, double y)
+#else
+	double copysign(x,y)
+	double x,y;
+#endif
+{
+	return (double) copysignf((float) x, (float) y);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */

agbcc/libc/math/sf_finite.c

@@ -0,0 +1,47 @@
+/* sf_finite.c -- float version of s_finite.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * finitef(x) returns 1 is x is finite, else 0;
+ * no branching!
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+	int finitef(float x)
+#else
+	int finitef(x)
+	float x;
+#endif
+{
+	__int32_t ix;
+	GET_FLOAT_WORD(ix,x);
+	return  (int)((__uint32_t)((ix&0x7fffffff)-0x7f800000)>>31);
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	int finite(double x)
+#else
+	int finite(x)
+	double x;
+#endif
+{
+	return finitef((float) x);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */

agbcc/libc/math/sf_frexp.c

@@ -0,0 +1,61 @@
+/* sf_frexp.c -- float version of s_frexp.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+static const float
+#else
+static float
+#endif
+two25 =  3.3554432000e+07; /* 0x4c000000 */
+
+#ifdef __STDC__
+	float frexpf(float x, int *eptr)
+#else
+	float frexpf(x, eptr)
+	float x; int *eptr;
+#endif
+{
+	__int32_t hx, ix;
+	GET_FLOAT_WORD(hx,x);
+	ix = 0x7fffffff&hx;
+	*eptr = 0;
+	if(ix>=0x7f800000||(ix==0)) return x;	/* 0,inf,nan */
+	if (ix<0x00800000) {		/* subnormal */
+	    x *= two25;
+	    GET_FLOAT_WORD(hx,x);
+	    ix = hx&0x7fffffff;
+	    *eptr = -25;
+	}
+	*eptr += (ix>>23)-126;
+	hx = (hx&0x807fffff)|0x3f000000;
+	SET_FLOAT_WORD(x,hx);
+	return x;
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	double frexp(double x, int *eptr)
+#else
+	double frexp(x, eptr)
+	double x; int *eptr;
+#endif
+{
+	return (double) frexpf((float) x, eptr);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */

agbcc/libc/math/sf_isinf.c

@@ -0,0 +1,35 @@
+/*
+ * isinff(x) returns 1 if x is infinity, else 0;
+ * no branching!
+ * Added by Cygnus Support.
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+	int isinff(float x)
+#else
+	int isinff(x)
+	float x;
+#endif
+{
+	__int32_t ix;
+	GET_FLOAT_WORD(ix,x);
+	ix &= 0x7fffffff;
+	ix = 0x7f800000 - ix;
+	return 1 - (int)((__uint32_t)(ix|(-ix))>>31);
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	int isinf(double x)
+#else
+	int isinf(x)
+	double x;
+#endif
+{
+	return isinff((float) x);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */

agbcc/libc/math/sf_isnan.c

@@ -0,0 +1,49 @@
+/* sf_isnan.c -- float version of s_isnan.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * isnanf(x) returns 1 is x is nan, else 0;
+ * no branching!
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+	int isnanf(float x)
+#else
+	int isnanf(x)
+	float x;
+#endif
+{
+	__int32_t ix;
+	GET_FLOAT_WORD(ix,x);
+	ix &= 0x7fffffff;
+	ix = 0x7f800000 - ix;
+	return (int)(((__uint32_t)(ix))>>31);
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	int isnan(double x)
+#else
+	int isnan(x)
+	double x;
+#endif
+{
+	return isnanf((float) x);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */

agbcc/libc/math/sf_ldexp.c

@@ -0,0 +1,44 @@
+/* sf_ldexp.c -- float version of s_ldexp.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include "fdlibm.h"
+#include <errno.h>
+
+#ifdef __STDC__
+	float ldexpf(float value, int exp)
+#else
+	float ldexpf(value, exp)
+	float value; int exp;
+#endif
+{
+	if(!finitef(value)||value==(float)0.0) return value;
+	value = scalbnf(value,exp);
+	if(!finitef(value)||value==(float)0.0) errno = ERANGE;
+	return value;
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	double ldexp(double value, int exp)
+#else
+	double ldexp(value, exp)
+	double value; int exp;
+#endif
+{
+	return (double) ldexpf((float) value, exp);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */

agbcc/libc/math/sf_modf.c

@@ -0,0 +1,73 @@
+/* sf_modf.c -- float version of s_modf.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+static const float one = 1.0;
+#else
+static float one = 1.0;
+#endif
+
+#ifdef __STDC__
+	float modff(float x, float *iptr)
+#else
+	float modff(x, iptr)
+	float x,*iptr;
+#endif
+{
+	__int32_t i0,j0;
+	__uint32_t i;
+	GET_FLOAT_WORD(i0,x);
+	j0 = ((i0>>23)&0xff)-0x7f;	/* exponent of x */
+	if(j0<23) {			/* integer part in x */
+	    if(j0<0) {			/* |x|<1 */
+	        SET_FLOAT_WORD(*iptr,i0&0x80000000);	/* *iptr = +-0 */
+		return x;
+	    } else {
+		i = (0x007fffff)>>j0;
+		if((i0&i)==0) {			/* x is integral */
+		    __uint32_t ix;
+		    *iptr = x;
+		    GET_FLOAT_WORD(ix,x);
+		    SET_FLOAT_WORD(x,ix&0x80000000);	/* return +-0 */
+		    return x;
+		} else {
+		    SET_FLOAT_WORD(*iptr,i0&(~i));
+		    return x - *iptr;
+		}
+	    }
+	} else {			/* no fraction part */
+	    __uint32_t ix;
+	    *iptr = x*one;
+	    GET_FLOAT_WORD(ix,x);
+	    SET_FLOAT_WORD(x,ix&0x80000000);	/* return +-0 */
+	    return x;
+	}
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	double modf(double x, double *iptr)
+#else
+	double modf(x, iptr)
+	double x,*iptr;
+#endif
+{
+	return (double) modff((float) x, (float *) iptr);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */

agbcc/libc/math/sf_scalbn.c

@@ -0,0 +1,79 @@
+/* sf_scalbn.c -- float version of s_scalbn.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include "fdlibm.h"
+#include <limits.h>
+
+#if INT_MAX > 50000
+#define OVERFLOW_INT 50000
+#else
+#define OVERFLOW_INT 30000
+#endif
+
+#ifdef __STDC__
+static const float
+#else
+static float
+#endif
+two25   =  3.355443200e+07,	/* 0x4c000000 */
+twom25  =  2.9802322388e-08,	/* 0x33000000 */
+huge   = 1.0e+30,
+tiny   = 1.0e-30;
+
+#ifdef __STDC__
+	float scalbnf (float x, int n)
+#else
+	float scalbnf (x,n)
+	float x; int n;
+#endif
+{
+	__int32_t  k,ix;
+	GET_FLOAT_WORD(ix,x);
+        k = (ix&0x7f800000)>>23;		/* extract exponent */
+        if (k==0) {				/* 0 or subnormal x */
+            if ((ix&0x7fffffff)==0) return x; /* +-0 */
+	    x *= two25;
+	    GET_FLOAT_WORD(ix,x);
+	    k = ((ix&0x7f800000)>>23) - 25; 
+            if (n< -50000) return tiny*x; 	/*underflow*/
+	    }
+        if (k==0xff) return x+x;		/* NaN or Inf */
+        k = k+n; 
+        if (k >  0xfe) return huge*copysignf(huge,x); /* overflow  */
+        if (k > 0) 				/* normal result */
+	    {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
+        if (k <= -25)
+            if (n > OVERFLOW_INT) 	/* in case integer overflow in n+k */
+		return huge*copysignf(huge,x);	/*overflow*/
+	    else return tiny*copysignf(tiny,x);	/*underflow*/
+        k += 25;				/* subnormal result */
+	SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
+        return x*twom25;
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	double scalbn(double x, int n)
+#else
+	double scalbn(x,n)
+	double x;
+	int n;
+#endif
+{
+	return (double) scalbnf((float) x, n);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */