summaryrefslogtreecommitdiff
path: root/libc/math
diff options
context:
space:
mode:
Diffstat (limited to 'libc/math')
-rw-r--r--libc/math/fdlibm.h261
-rw-r--r--libc/math/s_copysign.c82
-rw-r--r--libc/math/s_finite.c35
-rw-r--r--libc/math/s_frexp.c114
-rw-r--r--libc/math/s_infconst.c15
-rw-r--r--libc/math/s_isinf.c26
-rw-r--r--libc/math/s_isnan.c122
-rw-r--r--libc/math/s_ldexp.c81
-rw-r--r--libc/math/s_modf.c131
-rw-r--r--libc/math/s_scalbn.c103
-rw-r--r--libc/math/sf_copysign.c50
-rw-r--r--libc/math/sf_finite.c47
-rw-r--r--libc/math/sf_frexp.c61
-rw-r--r--libc/math/sf_isinf.c35
-rw-r--r--libc/math/sf_isnan.c49
-rw-r--r--libc/math/sf_ldexp.c44
-rw-r--r--libc/math/sf_modf.c73
-rw-r--r--libc/math/sf_scalbn.c79
18 files changed, 1408 insertions, 0 deletions
diff --git a/libc/math/fdlibm.h b/libc/math/fdlibm.h
new file mode 100644
index 0000000..752d242
--- /dev/null
+++ b/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)
diff --git a/libc/math/s_copysign.c b/libc/math/s_copysign.c
new file mode 100644
index 0000000..bfc546d
--- /dev/null
+++ b/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 */
diff --git a/libc/math/s_finite.c b/libc/math/s_finite.c
new file mode 100644
index 0000000..17f4e84
--- /dev/null
+++ b/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 */
diff --git a/libc/math/s_frexp.c b/libc/math/s_frexp.c
new file mode 100644
index 0000000..aaa3606
--- /dev/null
+++ b/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 */
diff --git a/libc/math/s_infconst.c b/libc/math/s_infconst.c
new file mode 100644
index 0000000..85b3b68
--- /dev/null
+++ b/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) */
diff --git a/libc/math/s_isinf.c b/libc/math/s_isinf.c
new file mode 100644
index 0000000..87f0995
--- /dev/null
+++ b/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 */
diff --git a/libc/math/s_isnan.c b/libc/math/s_isnan.c
new file mode 100644
index 0000000..5d83fc0
--- /dev/null
+++ b/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 */
diff --git a/libc/math/s_ldexp.c b/libc/math/s_ldexp.c
new file mode 100644
index 0000000..ccf7171
--- /dev/null
+++ b/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 */
diff --git a/libc/math/s_modf.c b/libc/math/s_modf.c
new file mode 100644
index 0000000..0115139
--- /dev/null
+++ b/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 */
diff --git a/libc/math/s_scalbn.c b/libc/math/s_scalbn.c
new file mode 100644
index 0000000..245888f
--- /dev/null
+++ b/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 */
diff --git a/libc/math/sf_copysign.c b/libc/math/sf_copysign.c
new file mode 100644
index 0000000..f547c82
--- /dev/null
+++ b/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) */
diff --git a/libc/math/sf_finite.c b/libc/math/sf_finite.c
new file mode 100644
index 0000000..4c48f40
--- /dev/null
+++ b/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) */
diff --git a/libc/math/sf_frexp.c b/libc/math/sf_frexp.c
new file mode 100644
index 0000000..271fb9d
--- /dev/null
+++ b/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) */
diff --git a/libc/math/sf_isinf.c b/libc/math/sf_isinf.c
new file mode 100644
index 0000000..1af4aab
--- /dev/null
+++ b/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) */
diff --git a/libc/math/sf_isnan.c b/libc/math/sf_isnan.c
new file mode 100644
index 0000000..8fdb6ad
--- /dev/null
+++ b/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) */
diff --git a/libc/math/sf_ldexp.c b/libc/math/sf_ldexp.c
new file mode 100644
index 0000000..2781304
--- /dev/null
+++ b/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) */
diff --git a/libc/math/sf_modf.c b/libc/math/sf_modf.c
new file mode 100644
index 0000000..6c64e3f
--- /dev/null
+++ b/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) */
diff --git a/libc/math/sf_scalbn.c b/libc/math/sf_scalbn.c
new file mode 100644
index 0000000..ee65f40
--- /dev/null
+++ b/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) */
generated by cgit v1.2.3 (git 2.25.1) at 2025-05-06 17:28:03 +0000