--- /dev/null
+/* ef_acos.c -- float version of e_acos.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
+one = 1.0000000000e+00, /* 0x3F800000 */
+pi = 3.1415925026e+00, /* 0x40490fda */
+pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */
+pio2_lo = 7.5497894159e-08, /* 0x33a22168 */
+pS0 = 1.6666667163e-01, /* 0x3e2aaaab */
+pS1 = -3.2556581497e-01, /* 0xbea6b090 */
+pS2 = 2.0121252537e-01, /* 0x3e4e0aa8 */
+pS3 = -4.0055535734e-02, /* 0xbd241146 */
+pS4 = 7.9153501429e-04, /* 0x3a4f7f04 */
+pS5 = 3.4793309169e-05, /* 0x3811ef08 */
+qS1 = -2.4033949375e+00, /* 0xc019d139 */
+qS2 = 2.0209457874e+00, /* 0x4001572d */
+qS3 = -6.8828397989e-01, /* 0xbf303361 */
+qS4 = 7.7038154006e-02; /* 0x3d9dc62e */
+
+#ifdef __STDC__
+ float __ieee754_acosf(float x)
+#else
+ float __ieee754_acosf(x)
+ float x;
+#endif
+{
+ float z,p,q,r,w,s,c,df;
+ __int32_t hx,ix;
+ GET_FLOAT_WORD(hx,x);
+ ix = hx&0x7fffffff;
+ if(ix==0x3f800000) { /* |x|==1 */
+ if(hx>0) return 0.0; /* acos(1) = 0 */
+ else return pi+(float)2.0*pio2_lo; /* acos(-1)= pi */
+ } else if(ix>0x3f800000) { /* |x| >= 1 */
+ return (x-x)/(x-x); /* acos(|x|>1) is NaN */
+ }
+ if(ix<0x3f000000) { /* |x| < 0.5 */
+ if(ix<=0x23000000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
+ z = x*x;
+ p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+ q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+ r = p/q;
+ return pio2_hi - (x - (pio2_lo-x*r));
+ } else if (hx<0) { /* x < -0.5 */
+ z = (one+x)*(float)0.5;
+ p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+ q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+ s = __ieee754_sqrtf(z);
+ r = p/q;
+ w = r*s-pio2_lo;
+ return pi - (float)2.0*(s+w);
+ } else { /* x > 0.5 */
+ __int32_t idf;
+ z = (one-x)*(float)0.5;
+ s = __ieee754_sqrtf(z);
+ df = s;
+ GET_FLOAT_WORD(idf,df);
+ SET_FLOAT_WORD(df,idf&0xfffff000);
+ c = (z-df*df)/(s+df);
+ p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+ q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+ r = p/q;
+ w = r*s+c;
+ return (float)2.0*(df+w);
+ }
+}
--- /dev/null
+/* ef_rem_pio2.c -- float version of e_rem_pio2.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.
+ * ====================================================
+ *
+ */
+
+/* __ieee754_rem_pio2f(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __kernel_rem_pio2f()
+ */
+
+#include "fdlibm.h"
+
+/*
+ * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
+ */
+#ifdef __STDC__
+static const __int32_t two_over_pi[] = {
+#else
+static __int32_t two_over_pi[] = {
+#endif
+0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
+0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
+0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
+0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
+0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
+0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
+0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
+0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
+0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
+0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
+0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
+0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
+0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
+0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
+0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
+0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
+0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
+0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
+0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
+0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
+0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
+0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
+};
+
+/* This array is like the one in e_rem_pio2.c, but the numbers are
+ single precision and the last 8 bits are forced to 0. */
+#ifdef __STDC__
+static const __int32_t npio2_hw[] = {
+#else
+static __int32_t npio2_hw[] = {
+#endif
+0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
+0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
+0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
+0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
+0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
+0x4242c700, 0x42490f00
+};
+
+/*
+ * invpio2: 24 bits of 2/pi
+ * pio2_1: first 17 bit of pi/2
+ * pio2_1t: pi/2 - pio2_1
+ * pio2_2: second 17 bit of pi/2
+ * pio2_2t: pi/2 - (pio2_1+pio2_2)
+ * pio2_3: third 17 bit of pi/2
+ * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+
+#ifdef __STDC__
+static const float
+#else
+static float
+#endif
+zero = 0.0000000000e+00, /* 0x00000000 */
+half = 5.0000000000e-01, /* 0x3f000000 */
+two8 = 2.5600000000e+02, /* 0x43800000 */
+invpio2 = 6.3661980629e-01, /* 0x3f22f984 */
+pio2_1 = 1.5707855225e+00, /* 0x3fc90f80 */
+pio2_1t = 1.0804334124e-05, /* 0x37354443 */
+pio2_2 = 1.0804273188e-05, /* 0x37354400 */
+pio2_2t = 6.0770999344e-11, /* 0x2e85a308 */
+pio2_3 = 6.0770943833e-11, /* 0x2e85a300 */
+pio2_3t = 6.1232342629e-17; /* 0x248d3132 */
+
+#ifdef __STDC__
+ __int32_t __ieee754_rem_pio2f(float x, float *y)
+#else
+ __int32_t __ieee754_rem_pio2f(x,y)
+ float x,y[];
+#endif
+{
+ float z,w,t,r,fn;
+ float tx[3];
+ __int32_t i,j,n,ix,hx;
+ int e0,nx;
+
+ GET_FLOAT_WORD(hx,x);
+ ix = hx&0x7fffffff;
+ if(ix<=0x3f490fd8) /* |x| ~<= pi/4 , no need for reduction */
+ {y[0] = x; y[1] = 0; return 0;}
+ if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */
+ if(hx>0) {
+ z = x - pio2_1;
+ if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
+ y[0] = z - pio2_1t;
+ y[1] = (z-y[0])-pio2_1t;
+ } else { /* near pi/2, use 24+24+24 bit pi */
+ z -= pio2_2;
+ y[0] = z - pio2_2t;
+ y[1] = (z-y[0])-pio2_2t;
+ }
+ return 1;
+ } else { /* negative x */
+ z = x + pio2_1;
+ if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
+ y[0] = z + pio2_1t;
+ y[1] = (z-y[0])+pio2_1t;
+ } else { /* near pi/2, use 24+24+24 bit pi */
+ z += pio2_2;
+ y[0] = z + pio2_2t;
+ y[1] = (z-y[0])+pio2_2t;
+ }
+ return -1;
+ }
+ }
+ if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
+ t = fabsf(x);
+ n = (__int32_t) (t*invpio2+half);
+ fn = (float)n;
+ r = t-fn*pio2_1;
+ w = fn*pio2_1t; /* 1st round good to 40 bit */
+ if(n<32&&(ix&0xffffff00)!=npio2_hw[n-1]) {
+ y[0] = r-w; /* quick check no cancellation */
+ } else {
+ __uint32_t high;
+ j = ix>>23;
+ y[0] = r-w;
+ GET_FLOAT_WORD(high,y[0]);
+ i = j-((high>>23)&0xff);
+ if(i>8) { /* 2nd iteration needed, good to 57 */
+ t = r;
+ w = fn*pio2_2;
+ r = t-w;
+ w = fn*pio2_2t-((t-r)-w);
+ y[0] = r-w;
+ GET_FLOAT_WORD(high,y[0]);
+ i = j-((high>>23)&0xff);
+ if(i>25) { /* 3rd iteration need, 74 bits acc */
+ t = r; /* will cover all possible cases */
+ w = fn*pio2_3;
+ r = t-w;
+ w = fn*pio2_3t-((t-r)-w);
+ y[0] = r-w;
+ }
+ }
+ }
+ y[1] = (r-y[0])-w;
+ if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
+ else return n;
+ }
+ /*
+ * all other (large) arguments
+ */
+ if(!FLT_UWORD_IS_FINITE(ix)) {
+ y[0]=y[1]=x-x; return 0;
+ }
+ /* set z = scalbn(|x|,ilogb(x)-7) */
+ e0 = (int)((ix>>23)-134); /* e0 = ilogb(z)-7; */
+ SET_FLOAT_WORD(z, ix - ((__int32_t)e0<<23));
+ for(i=0;i<2;i++) {
+ tx[i] = (float)((__int32_t)(z));
+ z = (z-tx[i])*two8;
+ }
+ tx[2] = z;
+ nx = 3;
+ while(tx[nx-1]==zero) nx--; /* skip zero term */
+ n = __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi);
+ if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
+ return n;
+}
--- /dev/null
+/* ef_sqrtf.c -- float version of e_sqrt.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, tiny=1.0e-30;
+#else
+static float one = 1.0, tiny=1.0e-30;
+#endif
+
+#ifdef __STDC__
+ float __ieee754_sqrtf(float x)
+#else
+ float __ieee754_sqrtf(x)
+ float x;
+#endif
+{
+ float z;
+ __uint32_t r,hx;
+ __int32_t ix,s,q,m,t,i;
+
+ GET_FLOAT_WORD(ix,x);
+ hx = ix&0x7fffffff;
+
+ /* take care of Inf and NaN */
+ if(!FLT_UWORD_IS_FINITE(hx))
+ return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
+ sqrt(-inf)=sNaN */
+ /* take care of zero and -ves */
+ if(FLT_UWORD_IS_ZERO(hx)) return x;/* sqrt(+-0) = +-0 */
+ if(ix<0) return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
+
+ /* normalize x */
+ m = (ix>>23);
+ if(FLT_UWORD_IS_SUBNORMAL(hx)) { /* subnormal x */
+ for(i=0;(ix&0x00800000L)==0;i++) ix<<=1;
+ m -= i-1;
+ }
+ m -= 127; /* unbias exponent */
+ ix = (ix&0x007fffffL)|0x00800000L;
+ if(m&1) /* odd m, double x to make it even */
+ ix += ix;
+ m >>= 1; /* m = [m/2] */
+
+ /* generate sqrt(x) bit by bit */
+ ix += ix;
+ q = s = 0; /* q = sqrt(x) */
+ r = 0x01000000L; /* r = moving bit from right to left */
+
+ while(r!=0) {
+ t = s+r;
+ if(t<=ix) {
+ s = t+r;
+ ix -= t;
+ q += r;
+ }
+ ix += ix;
+ r>>=1;
+ }
+
+ /* use floating add to find out rounding direction */
+ if(ix!=0) {
+ z = one-tiny; /* trigger inexact flag */
+ if (z>=one) {
+ z = one+tiny;
+ if (z>one)
+ q += 2;
+ else
+ q += (q&1);
+ }
+ }
+ ix = (q>>1)+0x3f000000L;
+ ix += (m <<23);
+ SET_FLOAT_WORD(z,ix);
+ return z;
+}
--- /dev/null
+
+/* @(#)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.
+ * ====================================================
+ */
+
+/* AltOS local */
+#include <math.h>
+#include <stdint.h>
+#define __int32_t int32_t
+#define __uint32_t uint32_t
+
+#define __ieee754_acosf acosf
+#define __ieee754_sqrtf sqrtf
+
+/* REDHAT LOCAL: Include files. */
+#include <math.h>
+/* #include <sys/types.h> */
+#include <machine/ieeefp.h>
+
+/* REDHAT LOCAL: Default to XOPEN_MODE. */
+#define _XOPEN_MODE
+
+/* Most routines need to check whether a float is finite, infinite, or not a
+ number, and many need to know whether the result of an operation will
+ overflow. These conditions depend on whether the largest exponent is
+ used for NaNs & infinities, or whether it's used for finite numbers. The
+ macros below wrap up that kind of information:
+
+ FLT_UWORD_IS_FINITE(X)
+ True if a positive float with bitmask X is finite.
+
+ FLT_UWORD_IS_NAN(X)
+ True if a positive float with bitmask X is not a number.
+
+ FLT_UWORD_IS_INFINITE(X)
+ True if a positive float with bitmask X is +infinity.
+
+ FLT_UWORD_MAX
+ The bitmask of FLT_MAX.
+
+ FLT_UWORD_HALF_MAX
+ The bitmask of FLT_MAX/2.
+
+ FLT_UWORD_EXP_MAX
+ The bitmask of the largest finite exponent (129 if the largest
+ exponent is used for finite numbers, 128 otherwise).
+
+ FLT_UWORD_LOG_MAX
+ The bitmask of log(FLT_MAX), rounded down. This value is the largest
+ input that can be passed to exp() without producing overflow.
+
+ FLT_UWORD_LOG_2MAX
+ The bitmask of log(2*FLT_MAX), rounded down. This value is the
+ largest input than can be passed to cosh() without producing
+ overflow.
+
+ FLT_LARGEST_EXP
+ The largest biased exponent that can be used for finite numbers
+ (255 if the largest exponent is used for finite numbers, 254
+ otherwise) */
+
+#ifdef _FLT_LARGEST_EXPONENT_IS_NORMAL
+#define FLT_UWORD_IS_FINITE(x) 1
+#define FLT_UWORD_IS_NAN(x) 0
+#define FLT_UWORD_IS_INFINITE(x) 0
+#define FLT_UWORD_MAX 0x7fffffff
+#define FLT_UWORD_EXP_MAX 0x43010000
+#define FLT_UWORD_LOG_MAX 0x42b2d4fc
+#define FLT_UWORD_LOG_2MAX 0x42b437e0
+#define HUGE ((float)0X1.FFFFFEP128)
+#else
+#define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L)
+#define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L)
+#define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L)
+#define FLT_UWORD_MAX 0x7f7fffffL
+#define FLT_UWORD_EXP_MAX 0x43000000
+#define FLT_UWORD_LOG_MAX 0x42b17217
+#define FLT_UWORD_LOG_2MAX 0x42b2d4fc
+#define HUGE ((float)3.40282346638528860e+38)
+#endif
+#define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1L<<23))
+#define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23)
+
+/* Many routines check for zero and subnormal numbers. Such things depend
+ on whether the target supports denormals or not:
+
+ FLT_UWORD_IS_ZERO(X)
+ True if a positive float with bitmask X is +0. Without denormals,
+ any float with a zero exponent is a +0 representation. With
+ denormals, the only +0 representation is a 0 bitmask.
+
+ FLT_UWORD_IS_SUBNORMAL(X)
+ True if a non-zero positive float with bitmask X is subnormal.
+ (Routines should check for zeros first.)
+
+ FLT_UWORD_MIN
+ The bitmask of the smallest float above +0. Call this number
+ REAL_FLT_MIN...
+
+ FLT_UWORD_EXP_MIN
+ The bitmask of the float representation of REAL_FLT_MIN's exponent.
+
+ FLT_UWORD_LOG_MIN
+ The bitmask of |log(REAL_FLT_MIN)|, rounding down.
+
+ FLT_SMALLEST_EXP
+ REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported,
+ -22 if they are).
+*/
+
+#ifdef _FLT_NO_DENORMALS
+#define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L)
+#define FLT_UWORD_IS_SUBNORMAL(x) 0
+#define FLT_UWORD_MIN 0x00800000
+#define FLT_UWORD_EXP_MIN 0x42fc0000
+#define FLT_UWORD_LOG_MIN 0x42aeac50
+#define FLT_SMALLEST_EXP 1
+#else
+#define FLT_UWORD_IS_ZERO(x) ((x)==0)
+#define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L)
+#define FLT_UWORD_MIN 0x00000001
+#define FLT_UWORD_EXP_MIN 0x43160000
+#define FLT_UWORD_LOG_MIN 0x42cff1b5
+#define FLT_SMALLEST_EXP -22
+#endif
+
+#ifdef __STDC__
+#undef __P
+#define __P(p) p
+#else
+#define __P(p) ()
+#endif
+
+/*
+ * 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>. */
+
+#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. */
+#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)
+
+/* Macros to avoid undefined behaviour that can arise if the amount
+ of a shift is exactly equal to the size of the shifted operand. */
+
+#define SAFE_LEFT_SHIFT(op,amt) \
+ (((amt) < 8 * sizeof(op)) ? ((op) << (amt)) : 0)
+
+#define SAFE_RIGHT_SHIFT(op,amt) \
+ (((amt) < 8 * sizeof(op)) ? ((op) >> (amt)) : 0)
+
+#ifdef _COMPLEX_H
+
+/*
+ * Quoting from ISO/IEC 9899:TC2:
+ *
+ * 6.2.5.13 Types
+ * Each complex type has the same representation and alignment requirements as
+ * an array type containing exactly two elements of the corresponding real type;
+ * the first element is equal to the real part, and the second element to the
+ * imaginary part, of the complex number.
+ */
+typedef union {
+ float complex z;
+ float parts[2];
+} float_complex;
+
+typedef union {
+ double complex z;
+ double parts[2];
+} double_complex;
+
+typedef union {
+ long double complex z;
+ long double parts[2];
+} long_double_complex;
+
+#define REAL_PART(z) ((z).parts[0])
+#define IMAG_PART(z) ((z).parts[1])
+
+#endif /* _COMPLEX_H */
+
--- /dev/null
+#ifndef _IEEE_FP_H_
+#define _IEEE_FP_H_
+
+#include "_ansi.h"
+
+#include <machine/ieeefp.h>
+
+_BEGIN_STD_C
+
+/* FIXME FIXME FIXME:
+ Neither of __ieee_{float,double}_shape_tape seem to be used anywhere
+ except in libm/test. If that is the case, please delete these from here.
+ If that is not the case, please insert documentation here describing why
+ they're needed. */
+
+#ifdef __IEEE_BIG_ENDIAN
+
+typedef union
+{
+ double value;
+ struct
+ {
+ unsigned int sign : 1;
+ unsigned int exponent: 11;
+ unsigned int fraction0:4;
+ unsigned int fraction1:16;
+ unsigned int fraction2:16;
+ unsigned int fraction3:16;
+
+ } number;
+ struct
+ {
+ unsigned int sign : 1;
+ unsigned int exponent: 11;
+ unsigned int quiet:1;
+ unsigned int function0:3;
+ unsigned int function1:16;
+ unsigned int function2:16;
+ unsigned int function3:16;
+ } nan;
+ struct
+ {
+ unsigned long msw;
+ unsigned long lsw;
+ } parts;
+ long aslong[2];
+} __ieee_double_shape_type;
+
+#endif
+
+#ifdef __IEEE_LITTLE_ENDIAN
+
+typedef union
+{
+ double value;
+ struct
+ {
+#ifdef __SMALL_BITFIELDS
+ unsigned int fraction3:16;
+ unsigned int fraction2:16;
+ unsigned int fraction1:16;
+ unsigned int fraction0: 4;
+#else
+ unsigned int fraction1:32;
+ unsigned int fraction0:20;
+#endif
+ unsigned int exponent :11;
+ unsigned int sign : 1;
+ } number;
+ struct
+ {
+#ifdef __SMALL_BITFIELDS
+ unsigned int function3:16;
+ unsigned int function2:16;
+ unsigned int function1:16;
+ unsigned int function0:3;
+#else
+ unsigned int function1:32;
+ unsigned int function0:19;
+#endif
+ unsigned int quiet:1;
+ unsigned int exponent: 11;
+ unsigned int sign : 1;
+ } nan;
+ struct
+ {
+ unsigned long lsw;
+ unsigned long msw;
+ } parts;
+
+ long aslong[2];
+
+} __ieee_double_shape_type;
+
+#endif
+
+#ifdef __IEEE_BIG_ENDIAN
+
+typedef union
+{
+ float value;
+ struct
+ {
+ unsigned int sign : 1;
+ unsigned int exponent: 8;
+ unsigned int fraction0: 7;
+ unsigned int fraction1: 16;
+ } number;
+ struct
+ {
+ unsigned int sign:1;
+ unsigned int exponent:8;
+ unsigned int quiet:1;
+ unsigned int function0:6;
+ unsigned int function1:16;
+ } nan;
+ long p1;
+
+} __ieee_float_shape_type;
+
+#endif
+
+#ifdef __IEEE_LITTLE_ENDIAN
+
+typedef union
+{
+ float value;
+ struct
+ {
+ unsigned int fraction0: 7;
+ unsigned int fraction1: 16;
+ unsigned int exponent: 8;
+ unsigned int sign : 1;
+ } number;
+ struct
+ {
+ unsigned int function1:16;
+ unsigned int function0:6;
+ unsigned int quiet:1;
+ unsigned int exponent:8;
+ unsigned int sign:1;
+ } nan;
+ long p1;
+
+} __ieee_float_shape_type;
+
+#endif
+
+
+
+
+
+/* FLOATING ROUNDING */
+
+typedef int fp_rnd;
+#define FP_RN 0 /* Round to nearest */
+#define FP_RM 1 /* Round down */
+#define FP_RP 2 /* Round up */
+#define FP_RZ 3 /* Round to zero (trunate) */
+
+fp_rnd _EXFUN(fpgetround,(void));
+fp_rnd _EXFUN(fpsetround, (fp_rnd));
+
+/* EXCEPTIONS */
+
+typedef int fp_except;
+#define FP_X_INV 0x10 /* Invalid operation */
+#define FP_X_DX 0x80 /* Divide by zero */
+#define FP_X_OFL 0x04 /* Overflow exception */
+#define FP_X_UFL 0x02 /* Underflow exception */
+#define FP_X_IMP 0x01 /* imprecise exception */
+
+fp_except _EXFUN(fpgetmask,(void));
+fp_except _EXFUN(fpsetmask,(fp_except));
+fp_except _EXFUN(fpgetsticky,(void));
+fp_except _EXFUN(fpsetsticky, (fp_except));
+
+/* INTEGER ROUNDING */
+
+typedef int fp_rdi;
+#define FP_RDI_TOZ 0 /* Round to Zero */
+#define FP_RDI_RD 1 /* Follow float mode */
+
+fp_rdi _EXFUN(fpgetroundtoi,(void));
+fp_rdi _EXFUN(fpsetroundtoi,(fp_rdi));
+
+#undef isnan
+#undef isinf
+
+int _EXFUN(isnan, (double));
+int _EXFUN(isinf, (double));
+int _EXFUN(finite, (double));
+
+
+
+int _EXFUN(isnanf, (float));
+int _EXFUN(isinff, (float));
+int _EXFUN(finitef, (float));
+
+#define __IEEE_DBL_EXPBIAS 1023
+#define __IEEE_FLT_EXPBIAS 127
+
+#define __IEEE_DBL_EXPLEN 11
+#define __IEEE_FLT_EXPLEN 8
+
+
+#define __IEEE_DBL_FRACLEN (64 - (__IEEE_DBL_EXPLEN + 1))
+#define __IEEE_FLT_FRACLEN (32 - (__IEEE_FLT_EXPLEN + 1))
+
+#define __IEEE_DBL_MAXPOWTWO ((double)(1L << 32 - 2) * (1L << (32-11) - 32 + 1))
+#define __IEEE_FLT_MAXPOWTWO ((float)(1L << (32-8) - 1))
+
+#define __IEEE_DBL_NAN_EXP 0x7ff
+#define __IEEE_FLT_NAN_EXP 0xff
+
+#ifndef __ieeefp_isnanf
+#define __ieeefp_isnanf(x) (((*(long *)&(x) & 0x7f800000L)==0x7f800000L) && \
+ ((*(long *)&(x) & 0x007fffffL)!=0000000000L))
+#endif
+#define isnanf(x) __ieeefp_isnanf(x)
+
+#ifndef __ieeefp_isinff
+#define __ieeefp_isinff(x) (((*(long *)&(x) & 0x7f800000L)==0x7f800000L) && \
+ ((*(long *)&(x) & 0x007fffffL)==0000000000L))
+#endif
+#define isinff(x) __ieeefp_isinff(x)
+
+#ifndef __ieeefp_finitef
+#define __ieeefp_finitef(x) (((*(long *)&(x) & 0x7f800000L)!=0x7f800000L))
+#endif
+#define finitef(x) __ieeefp_finitef(x)
+
+#ifdef _DOUBLE_IS_32BITS
+#undef __IEEE_DBL_EXPBIAS
+#define __IEEE_DBL_EXPBIAS __IEEE_FLT_EXPBIAS
+
+#undef __IEEE_DBL_EXPLEN
+#define __IEEE_DBL_EXPLEN __IEEE_FLT_EXPLEN
+
+#undef __IEEE_DBL_FRACLEN
+#define __IEEE_DBL_FRACLEN __IEEE_FLT_FRACLEN
+
+#undef __IEEE_DBL_MAXPOWTWO
+#define __IEEE_DBL_MAXPOWTWO __IEEE_FLT_MAXPOWTWO
+
+#undef __IEEE_DBL_NAN_EXP
+#define __IEEE_DBL_NAN_EXP __IEEE_FLT_NAN_EXP
+
+#undef __ieee_double_shape_type
+#define __ieee_double_shape_type __ieee_float_shape_type
+
+#endif /* _DOUBLE_IS_32BITS */
+
+_END_STD_C
+
+#endif /* _IEEE_FP_H_ */
--- /dev/null
+/* kf_cos.c -- float version of k_cos.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
+one = 1.0000000000e+00, /* 0x3f800000 */
+C1 = 4.1666667908e-02, /* 0x3d2aaaab */
+C2 = -1.3888889225e-03, /* 0xbab60b61 */
+C3 = 2.4801587642e-05, /* 0x37d00d01 */
+C4 = -2.7557314297e-07, /* 0xb493f27c */
+C5 = 2.0875723372e-09, /* 0x310f74f6 */
+C6 = -1.1359647598e-11; /* 0xad47d74e */
+
+#ifdef __STDC__
+ float __kernel_cosf(float x, float y)
+#else
+ float __kernel_cosf(x, y)
+ float x,y;
+#endif
+{
+ float a,hz,z,r,qx;
+ __int32_t ix;
+ GET_FLOAT_WORD(ix,x);
+ ix &= 0x7fffffff; /* ix = |x|'s high word*/
+ if(ix<0x32000000) { /* if x < 2**27 */
+ if(((int)x)==0) return one; /* generate inexact */
+ }
+ z = x*x;
+ r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
+ if(ix < 0x3e99999a) /* if |x| < 0.3 */
+ return one - ((float)0.5*z - (z*r - x*y));
+ else {
+ if(ix > 0x3f480000) { /* x > 0.78125 */
+ qx = (float)0.28125;
+ } else {
+ SET_FLOAT_WORD(qx,ix-0x01000000); /* x/4 */
+ }
+ hz = (float)0.5*z-qx;
+ a = one-qx;
+ return a - (hz - (z*r-x*y));
+ }
+}
--- /dev/null
+/* kf_rem_pio2.c -- float version of k_rem_pio2.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"
+
+/* In the float version, the input parameter x contains 8 bit
+ integers, not 24 bit integers. 113 bit precision is not supported. */
+
+#ifdef __STDC__
+static const int init_jk[] = {4,7,9}; /* initial value for jk */
+#else
+static int init_jk[] = {4,7,9};
+#endif
+
+#ifdef __STDC__
+static const float PIo2[] = {
+#else
+static float PIo2[] = {
+#endif
+ 1.5703125000e+00, /* 0x3fc90000 */
+ 4.5776367188e-04, /* 0x39f00000 */
+ 2.5987625122e-05, /* 0x37da0000 */
+ 7.5437128544e-08, /* 0x33a20000 */
+ 6.0026650317e-11, /* 0x2e840000 */
+ 7.3896444519e-13, /* 0x2b500000 */
+ 5.3845816694e-15, /* 0x27c20000 */
+ 5.6378512969e-18, /* 0x22d00000 */
+ 8.3009228831e-20, /* 0x1fc40000 */
+ 3.2756352257e-22, /* 0x1bc60000 */
+ 6.3331015649e-25, /* 0x17440000 */
+};
+
+#ifdef __STDC__
+static const float
+#else
+static float
+#endif
+zero = 0.0,
+one = 1.0,
+two8 = 2.5600000000e+02, /* 0x43800000 */
+twon8 = 3.9062500000e-03; /* 0x3b800000 */
+
+#ifdef __STDC__
+ int __kernel_rem_pio2f(float *x, float *y, int e0, int nx, int prec, const __int32_t *ipio2)
+#else
+ int __kernel_rem_pio2f(x,y,e0,nx,prec,ipio2)
+ float x[], y[]; int e0,nx,prec; __int32_t ipio2[];
+#endif
+{
+ __int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
+ float z,fw,f[20],fq[20],q[20];
+
+ /* initialize jk*/
+ jk = init_jk[prec];
+ jp = jk;
+
+ /* determine jx,jv,q0, note that 3>q0 */
+ jx = nx-1;
+ jv = (e0-3)/8; if(jv<0) jv=0;
+ q0 = e0-8*(jv+1);
+
+ /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
+ j = jv-jx; m = jx+jk;
+ for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (float) ipio2[j];
+
+ /* compute q[0],q[1],...q[jk] */
+ for (i=0;i<=jk;i++) {
+ for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
+ }
+
+ jz = jk;
+recompute:
+ /* distill q[] into iq[] reversingly */
+ for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
+ fw = (float)((__int32_t)(twon8* z));
+ iq[i] = (__int32_t)(z-two8*fw);
+ z = q[j-1]+fw;
+ }
+
+ /* compute n */
+ z = scalbnf(z,(int)q0); /* actual value of z */
+ z -= (float)8.0*floorf(z*(float)0.125); /* trim off integer >= 8 */
+ n = (__int32_t) z;
+ z -= (float)n;
+ ih = 0;
+ if(q0>0) { /* need iq[jz-1] to determine n */
+ i = (iq[jz-1]>>(8-q0)); n += i;
+ iq[jz-1] -= i<<(8-q0);
+ ih = iq[jz-1]>>(7-q0);
+ }
+ else if(q0==0) ih = iq[jz-1]>>8;
+ else if(z>=(float)0.5) ih=2;
+
+ if(ih>0) { /* q > 0.5 */
+ n += 1; carry = 0;
+ for(i=0;i<jz ;i++) { /* compute 1-q */
+ j = iq[i];
+ if(carry==0) {
+ if(j!=0) {
+ carry = 1; iq[i] = 0x100- j;
+ }
+ } else iq[i] = 0xff - j;
+ }
+ if(q0>0) { /* rare case: chance is 1 in 12 */
+ switch(q0) {
+ case 1:
+ iq[jz-1] &= 0x7f; break;
+ case 2:
+ iq[jz-1] &= 0x3f; break;
+ }
+ }
+ if(ih==2) {
+ z = one - z;
+ if(carry!=0) z -= scalbnf(one,(int)q0);
+ }
+ }
+
+ /* check if recomputation is needed */
+ if(z==zero) {
+ j = 0;
+ for (i=jz-1;i>=jk;i--) j |= iq[i];
+ if(j==0) { /* need recomputation */
+ for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */
+
+ for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */
+ f[jx+i] = (float) ipio2[jv+i];
+ for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
+ q[i] = fw;
+ }
+ jz += k;
+ goto recompute;
+ }
+ }
+
+ /* chop off zero terms */
+ if(z==(float)0.0) {
+ jz -= 1; q0 -= 8;
+ while(iq[jz]==0) { jz--; q0-=8;}
+ } else { /* break z into 8-bit if necessary */
+ z = scalbnf(z,-(int)q0);
+ if(z>=two8) {
+ fw = (float)((__int32_t)(twon8*z));
+ iq[jz] = (__int32_t)(z-two8*fw);
+ jz += 1; q0 += 8;
+ iq[jz] = (__int32_t) fw;
+ } else iq[jz] = (__int32_t) z ;
+ }
+
+ /* convert integer "bit" chunk to floating-point value */
+ fw = scalbnf(one,(int)q0);
+ for(i=jz;i>=0;i--) {
+ q[i] = fw*(float)iq[i]; fw*=twon8;
+ }
+
+ /* compute PIo2[0,...,jp]*q[jz,...,0] */
+ for(i=jz;i>=0;i--) {
+ for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
+ fq[jz-i] = fw;
+ }
+
+ /* compress fq[] into y[] */
+ switch(prec) {
+ case 0:
+ fw = 0.0;
+ for (i=jz;i>=0;i--) fw += fq[i];
+ y[0] = (ih==0)? fw: -fw;
+ break;
+ case 1:
+ case 2:
+ fw = 0.0;
+ for (i=jz;i>=0;i--) fw += fq[i];
+ y[0] = (ih==0)? fw: -fw;
+ fw = fq[0]-fw;
+ for (i=1;i<=jz;i++) fw += fq[i];
+ y[1] = (ih==0)? fw: -fw;
+ break;
+ case 3: /* painful */
+ for (i=jz;i>0;i--) {
+ fw = fq[i-1]+fq[i];
+ fq[i] += fq[i-1]-fw;
+ fq[i-1] = fw;
+ }
+ for (i=jz;i>1;i--) {
+ fw = fq[i-1]+fq[i];
+ fq[i] += fq[i-1]-fw;
+ fq[i-1] = fw;
+ }
+ for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
+ if(ih==0) {
+ y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
+ } else {
+ y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
+ }
+ }
+ return n&7;
+}
--- /dev/null
+/* kf_sin.c -- float version of k_sin.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
+half = 5.0000000000e-01,/* 0x3f000000 */
+S1 = -1.6666667163e-01, /* 0xbe2aaaab */
+S2 = 8.3333337680e-03, /* 0x3c088889 */
+S3 = -1.9841270114e-04, /* 0xb9500d01 */
+S4 = 2.7557314297e-06, /* 0x3638ef1b */
+S5 = -2.5050759689e-08, /* 0xb2d72f34 */
+S6 = 1.5896910177e-10; /* 0x2f2ec9d3 */
+
+#ifdef __STDC__
+ float __kernel_sinf(float x, float y, int iy)
+#else
+ float __kernel_sinf(x, y, iy)
+ float x,y; int iy; /* iy=0 if y is zero */
+#endif
+{
+ float z,r,v;
+ __int32_t ix;
+ GET_FLOAT_WORD(ix,x);
+ ix &= 0x7fffffff; /* high word of x */
+ if(ix<0x32000000) /* |x| < 2**-27 */
+ {if((int)x==0) return x;} /* generate inexact */
+ z = x*x;
+ v = z*x;
+ r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
+ if(iy==0) return x+v*(S1+z*r);
+ else return x-((z*(half*y-v*r)-y)-v*S1);
+}
--- /dev/null
+#ifndef __IEEE_BIG_ENDIAN
+#ifndef __IEEE_LITTLE_ENDIAN
+
+/* This file can define macros to choose variations of the IEEE float
+ format:
+
+ _FLT_LARGEST_EXPONENT_IS_NORMAL
+
+ Defined if the float format uses the largest exponent for finite
+ numbers rather than NaN and infinity representations. Such a
+ format cannot represent NaNs or infinities at all, but it's FLT_MAX
+ is twice the IEEE value.
+
+ _FLT_NO_DENORMALS
+
+ Defined if the float format does not support IEEE denormals. Every
+ float with a zero exponent is taken to be a zero representation.
+
+ ??? At the moment, there are no equivalent macros above for doubles and
+ the macros are not fully supported by --enable-newlib-hw-fp.
+
+ __IEEE_BIG_ENDIAN
+
+ Defined if the float format is big endian. This is mutually exclusive
+ with __IEEE_LITTLE_ENDIAN.
+
+ __IEEE_LITTLE_ENDIAN
+
+ Defined if the float format is little endian. This is mutually exclusive
+ with __IEEE_BIG_ENDIAN.
+
+ Note that one of __IEEE_BIG_ENDIAN or __IEEE_LITTLE_ENDIAN must be specified for a
+ platform or error will occur.
+
+ __IEEE_BYTES_LITTLE_ENDIAN
+
+ This flag is used in conjunction with __IEEE_BIG_ENDIAN to describe a situation
+ whereby multiple words of an IEEE floating point are in big endian order, but the
+ words themselves are little endian with respect to the bytes.
+
+ _DOUBLE_IS_32BITS
+
+ This is used on platforms that support double by using the 32-bit IEEE
+ float type.
+
+ _FLOAT_ARG
+
+ This represents what type a float arg is passed as. It is used when the type is
+ not promoted to double.
+
+*/
+
+#if (defined(__arm__) || defined(__thumb__)) && !defined(__MAVERICK__)
+/* ARM traditionally used big-endian words; and within those words the
+ byte ordering was big or little endian depending upon the target.
+ Modern floating-point formats are naturally ordered; in this case
+ __VFP_FP__ will be defined, even if soft-float. */
+#ifdef __VFP_FP__
+# ifdef __ARMEL__
+# define __IEEE_LITTLE_ENDIAN
+# else
+# define __IEEE_BIG_ENDIAN
+# endif
+#else
+# define __IEEE_BIG_ENDIAN
+# ifdef __ARMEL__
+# define __IEEE_BYTES_LITTLE_ENDIAN
+# endif
+#endif
+#endif
+
+#ifdef __hppa__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __SPU__
+#define __IEEE_BIG_ENDIAN
+
+#define isfinite(__y) \
+ (__extension__ ({int __cy; \
+ (sizeof (__y) == sizeof (float)) ? (1) : \
+ (__cy = fpclassify(__y)) != FP_INFINITE && __cy != FP_NAN;}))
+
+#define isinf(__x) ((sizeof (__x) == sizeof (float)) ? (0) : __isinfd(__x))
+#define isnan(__x) ((sizeof (__x) == sizeof (float)) ? (0) : __isnand(__x))
+
+/*
+ * Macros for use in ieeefp.h. We can't just define the real ones here
+ * (like those above) as we have name space issues when this is *not*
+ * included via generic the ieeefp.h.
+ */
+#define __ieeefp_isnanf(x) 0
+#define __ieeefp_isinff(x) 0
+#define __ieeefp_finitef(x) 1
+#endif
+
+#ifdef __sparc__
+#ifdef __LITTLE_ENDIAN_DATA__
+#define __IEEE_LITTLE_ENDIAN
+#else
+#define __IEEE_BIG_ENDIAN
+#endif
+#endif
+
+#if defined(__m68k__) || defined(__mc68000__)
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#if defined(__mc68hc11__) || defined(__mc68hc12__) || defined(__mc68hc1x__)
+#define __IEEE_BIG_ENDIAN
+#ifdef __HAVE_SHORT_DOUBLE__
+# define _DOUBLE_IS_32BITS
+#endif
+#endif
+
+#if defined (__H8300__) || defined (__H8300H__) || defined (__H8300S__) || defined (__H8500__) || defined (__H8300SX__)
+#define __IEEE_BIG_ENDIAN
+#define _FLOAT_ARG float
+#define _DOUBLE_IS_32BITS
+#endif
+
+#if defined (__xc16x__) || defined (__xc16xL__) || defined (__xc16xS__)
+#define __IEEE_LITTLE_ENDIAN
+#define _FLOAT_ARG float
+#define _DOUBLE_IS_32BITS
+#endif
+
+
+#ifdef __sh__
+#ifdef __LITTLE_ENDIAN__
+#define __IEEE_LITTLE_ENDIAN
+#else
+#define __IEEE_BIG_ENDIAN
+#endif
+#if defined(__SH2E__) || defined(__SH3E__) || defined(__SH4_SINGLE_ONLY__) || defined(__SH2A_SINGLE_ONLY__)
+#define _DOUBLE_IS_32BITS
+#endif
+#endif
+
+#ifdef _AM29K
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef _WIN32
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __i386__
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __i960__
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __lm32__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __M32R__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#if defined(_C4x) || defined(_C3x)
+#define __IEEE_BIG_ENDIAN
+#define _DOUBLE_IS_32BITS
+#endif
+
+#ifdef __TMS320C6X__
+#ifdef _BIG_ENDIAN
+#define __IEEE_BIG_ENDIAN
+#else
+#define __IEEE_LITTLE_ENDIAN
+#endif
+#endif
+
+#ifdef __TIC80__
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __MIPSEL__
+#define __IEEE_LITTLE_ENDIAN
+#endif
+#ifdef __MIPSEB__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __MMIX__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __D30V__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+/* necv70 was __IEEE_LITTLE_ENDIAN. */
+
+#ifdef __W65__
+#define __IEEE_LITTLE_ENDIAN
+#define _DOUBLE_IS_32BITS
+#endif
+
+#if defined(__Z8001__) || defined(__Z8002__)
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __m88k__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __mn10300__
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __mn10200__
+#define __IEEE_LITTLE_ENDIAN
+#define _DOUBLE_IS_32BITS
+#endif
+
+#ifdef __v800
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __v850
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __D10V__
+#define __IEEE_BIG_ENDIAN
+#if __DOUBLE__ == 32
+#define _DOUBLE_IS_32BITS
+#endif
+#endif
+
+#ifdef __PPC__
+#if (defined(_BIG_ENDIAN) && _BIG_ENDIAN) || (defined(_AIX) && _AIX)
+#define __IEEE_BIG_ENDIAN
+#else
+#if (defined(_LITTLE_ENDIAN) && _LITTLE_ENDIAN) || (defined(__sun__) && __sun__) || (defined(_WIN32) && _WIN32)
+#define __IEEE_LITTLE_ENDIAN
+#endif
+#endif
+#endif
+
+#ifdef __xstormy16__
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __arc__
+#ifdef __big_endian__
+#define __IEEE_BIG_ENDIAN
+#else
+#define __IEEE_LITTLE_ENDIAN
+#endif
+#endif
+
+#ifdef __CRX__
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __fr30__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __mcore__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __mt__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __frv__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __moxie__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __ia64__
+#ifdef __BIG_ENDIAN__
+#define __IEEE_BIG_ENDIAN
+#else
+#define __IEEE_LITTLE_ENDIAN
+#endif
+#endif
+
+#ifdef __AVR__
+#define __IEEE_LITTLE_ENDIAN
+#define _DOUBLE_IS_32BITS
+#endif
+
+#if defined(__or32__) || defined(__or1k__) || defined(__or16__)
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __IP2K__
+#define __IEEE_BIG_ENDIAN
+#define __SMALL_BITFIELDS
+#define _DOUBLE_IS_32BITS
+#endif
+
+#ifdef __iq2000__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __MAVERICK__
+#ifdef __ARMEL__
+# define __IEEE_LITTLE_ENDIAN
+#else /* must be __ARMEB__ */
+# define __IEEE_BIG_ENDIAN
+#endif /* __ARMEL__ */
+#endif /* __MAVERICK__ */
+
+#ifdef __m32c__
+#define __IEEE_LITTLE_ENDIAN
+#define __SMALL_BITFIELDS
+#endif
+
+#ifdef __CRIS__
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __BFIN__
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __x86_64__
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifdef __mep__
+#ifdef __LITTLE_ENDIAN__
+#define __IEEE_LITTLE_ENDIAN
+#else
+#define __IEEE_BIG_ENDIAN
+#endif
+#endif
+
+#ifdef __MICROBLAZE__
+#define __IEEE_BIG_ENDIAN
+#endif
+
+#ifdef __RL78__
+#define __IEEE_LITTLE_ENDIAN
+#define __SMALL_BITFIELDS /* 16 Bit INT */
+#define _DOUBLE_IS_32BITS
+#endif
+
+#ifdef __RX__
+
+#ifdef __RX_BIG_ENDIAN__
+#define __IEEE_BIG_ENDIAN
+#else
+#define __IEEE_LITTLE_ENDIAN
+#endif
+
+#ifndef __RX_64BIT_DOUBLES__
+#define _DOUBLE_IS_32BITS
+#endif
+
+#ifdef __RX_16BIT_INTS__
+#define __SMALL_BITFIELDS
+#endif
+
+#endif
+
+#if (defined(__CR16__) || defined(__CR16C__) ||defined(__CR16CP__))
+#define __IEEE_LITTLE_ENDIAN
+#define __SMALL_BITFIELDS /* 16 Bit INT */
+#endif
+
+#ifndef __IEEE_BIG_ENDIAN
+#ifndef __IEEE_LITTLE_ENDIAN
+#error Endianess not declared!!
+#endif /* not __IEEE_LITTLE_ENDIAN */
+#endif /* not __IEEE_BIG_ENDIAN */
+
+#endif /* not __IEEE_LITTLE_ENDIAN */
+#endif /* not __IEEE_BIG_ENDIAN */
+
--- /dev/null
+/*
+ * Copyright © 2013 Keith Packard <keithp@keithp.com>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; version 2 of the License.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
+ */
+
+#ifndef _MATH_H_
+#define _MATH_H_
+
+float acosf(float x);
+
+float cosf(float x);
+
+float sqrtf(float x);
+
+float fabsf(float x);
+
+float floorf(float x);
+
+float scalbnf(float x, int n);
+
+float copysignf(float x, float y);
+
+#endif
--- /dev/null
+/* 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) */
--- /dev/null
+/* sf_cos.c -- float version of s_cos.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 cosf(float x)
+#else
+ float cosf(x)
+ float x;
+#endif
+{
+ float y[2],z=0.0;
+ __int32_t n,ix;
+
+ GET_FLOAT_WORD(ix,x);
+
+ /* |x| ~< pi/4 */
+ ix &= 0x7fffffff;
+ if(ix <= 0x3f490fd8) return __kernel_cosf(x,z);
+
+ /* cos(Inf or NaN) is NaN */
+ else if (!FLT_UWORD_IS_FINITE(ix)) return x-x;
+
+ /* argument reduction needed */
+ else {
+ n = __ieee754_rem_pio2f(x,y);
+ switch(n&3) {
+ case 0: return __kernel_cosf(y[0],y[1]);
+ case 1: return -__kernel_sinf(y[0],y[1],1);
+ case 2: return -__kernel_cosf(y[0],y[1]);
+ default:
+ return __kernel_sinf(y[0],y[1],1);
+ }
+ }
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+ double cos(double x)
+#else
+ double cos(x)
+ double x;
+#endif
+{
+ return (double) cosf((float) x);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */
--- /dev/null
+/* sf_fabs.c -- float version of s_fabs.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.
+ * ====================================================
+ */
+
+/*
+ * fabsf(x) returns the absolute value of x.
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+ float fabsf(float x)
+#else
+ float fabsf(x)
+ float x;
+#endif
+{
+ __uint32_t ix;
+ GET_FLOAT_WORD(ix,x);
+ SET_FLOAT_WORD(x,ix&0x7fffffff);
+ return x;
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+ double fabs(double x)
+#else
+ double fabs(x)
+ double x;
+#endif
+{
+ return (double) fabsf((float) x);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */
--- /dev/null
+/* sf_floor.c -- float version of s_floor.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.
+ * ====================================================
+ */
+
+/*
+ * floorf(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ * Bit twiddling.
+ * Exception:
+ * Inexact flag raised if x not equal to floorf(x).
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+static const float huge = 1.0e30;
+#else
+static float huge = 1.0e30;
+#endif
+
+#ifdef __STDC__
+ float floorf(float x)
+#else
+ float floorf(x)
+ float x;
+#endif
+{
+ __int32_t i0,j0;
+ __uint32_t i,ix;
+ GET_FLOAT_WORD(i0,x);
+ ix = (i0&0x7fffffff);
+ j0 = (ix>>23)-0x7f;
+ if(j0<23) {
+ if(j0<0) { /* raise inexact if x != 0 */
+ if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */
+ if(i0>=0) {i0=0;}
+ else if(!FLT_UWORD_IS_ZERO(ix))
+ { i0=0xbf800000;}
+ }
+ } else {
+ i = (0x007fffff)>>j0;
+ if((i0&i)==0) return x; /* x is integral */
+ if(huge+x>(float)0.0) { /* raise inexact flag */
+ if(i0<0) i0 += (0x00800000)>>j0;
+ i0 &= (~i);
+ }
+ }
+ } else {
+ if(!FLT_UWORD_IS_FINITE(ix)) return x+x; /* inf or NaN */
+ else return x; /* x is integral */
+ }
+ SET_FLOAT_WORD(x,i0);
+ return x;
+}
+
+#ifdef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+ double floor(double x)
+#else
+ double floor(x)
+ double x;
+#endif
+{
+ return (double) floorf((float) x);
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */
--- /dev/null
+/* 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>
+#include <float.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;
+ __uint32_t hx;
+
+ GET_FLOAT_WORD(ix,x);
+ hx = ix&0x7fffffff;
+ k = hx>>23; /* extract exponent */
+ if (FLT_UWORD_IS_ZERO(hx))
+ return x;
+ if (!FLT_UWORD_IS_FINITE(hx))
+ return x+x; /* NaN or Inf */
+ if (FLT_UWORD_IS_SUBNORMAL(hx)) {
+ x *= two25;
+ GET_FLOAT_WORD(ix,x);
+ k = ((ix&0x7f800000)>>23) - 25;
+ if (n< -50000) return tiny*x; /*underflow*/
+ }
+ k = k+n;
+ if (k > FLT_LARGEST_EXP) return huge*copysignf(huge,x); /* overflow */
+ if (k > 0) /* normal result */
+ {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
+ if (k < FLT_SMALLEST_EXP) {
+ 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) */
projects / fw / altos / commitdiff