diff --git a/liba/Makefile b/liba/Makefile index 60a339563..592e92c62 100644 --- a/liba/Makefile +++ b/liba/Makefile @@ -4,7 +4,8 @@ liba/src/external/sqlite/mem5.o: CFLAGS += -w objs += $(addprefix liba/src/, assert.o errno.o malloc.o memcpy.o memset.o strcmp.o strlcpy.o strlen.o external/sqlite/mem5.o) -objs += $(addprefix liba/src/external/msun/, \ +objs += $(addprefix liba/src/external/openbsd/, \ + e_logf.o \ e_log10f.o \ e_powf.o \ e_sqrtf.o \ @@ -13,8 +14,7 @@ objs += $(addprefix liba/src/external/msun/, \ s_scalbnf.o \ ) -liba/src/external/msun/%.o: CFLAGS += -Iliba/src/external/msun/include -Dlint=1 -Wno-parentheses -liba/src/external/msun/%.o: SFLAGS:=$(filter-out -DDEBUG=1,$(SFLAGS)) +liba/src/external/openbsd/%.o: CFLAGS += -Iliba/src/external/openbsd/include -Wno-parentheses tests += $(addprefix liba/test/, stdint.c strlcpy.c) diff --git a/liba/include/math.h b/liba/include/math.h index 41271f1ff..0f9914c4b 100644 --- a/liba/include/math.h +++ b/liba/include/math.h @@ -5,6 +5,7 @@ LIBA_BEGIN_DECLS +float logf(float x); float log10f(float x); float powf(float x, float y); float fabsf(float x); diff --git a/liba/src/external/msun/e_log10f.c b/liba/src/external/msun/e_log10f.c deleted file mode 100644 index 9856df2e7..000000000 --- a/liba/src/external/msun/e_log10f.c +++ /dev/null @@ -1,72 +0,0 @@ -/* - * ==================================================== - * 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 -__FBSDID("$FreeBSD$"); - -/* - * Float version of e_log10.c. See the latter for most comments. - */ - -#include "math.h" -#include "math_private.h" -#include "k_logf.h" - -static const float -two25 = 3.3554432000e+07, /* 0x4c000000 */ -ivln10hi = 4.3432617188e-01, /* 0x3ede6000 */ -ivln10lo = -3.1689971365e-05, /* 0xb804ead9 */ -log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ -log10_2lo = 7.9034151668e-07; /* 0x355427db */ - -static const float zero = 0.0; -static volatile float vzero = 0.0; - -float -__ieee754_log10f(float x) -{ - float f,hfsq,hi,lo,r,y; - int32_t i,k,hx; - - GET_FLOAT_WORD(hx,x); - - k=0; - if (hx < 0x00800000) { /* x < 2**-126 */ - if ((hx&0x7fffffff)==0) - return -two25/vzero; /* log(+-0)=-inf */ - if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ - k -= 25; x *= two25; /* subnormal number, scale up x */ - GET_FLOAT_WORD(hx,x); - } - if (hx >= 0x7f800000) return x+x; - if (hx == 0x3f800000) - return zero; /* log(1) = +0 */ - k += (hx>>23)-127; - hx &= 0x007fffff; - i = (hx+(0x4afb0d))&0x800000; - SET_FLOAT_WORD(x,hx|(i^0x3f800000)); /* normalize x or x/2 */ - k += (i>>23); - y = (float)k; - f = x - (float)1.0; - hfsq = (float)0.5*f*f; - r = k_log1pf(f); - - /* See e_log2f.c and e_log2.c for details. */ - if (sizeof(float_t) > sizeof(float)) - return (r - hfsq + f) * ((float_t)ivln10lo + ivln10hi) + - y * ((float_t)log10_2lo + log10_2hi); - hi = f - hfsq; - GET_FLOAT_WORD(hx,hi); - SET_FLOAT_WORD(hi,hx&0xfffff000); - lo = (f - hi) - hfsq + r; - return y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi + hi*ivln10hi + - y*log10_2hi; -} diff --git a/liba/src/external/msun/include/machine/endian.h b/liba/src/external/msun/include/machine/endian.h deleted file mode 100644 index de32f53f6..000000000 --- a/liba/src/external/msun/include/machine/endian.h +++ /dev/null @@ -1,7 +0,0 @@ -#ifndef LIBA_MSUN_MACHINE_ENDIAN_H -#define LIBA_MSUN_MACHINE_ENDIAN_H - -#define LITTLE_ENDIAN 1234 -#define BIG_ENDIAN 4321 - -#endif diff --git a/liba/src/external/msun/include/sys/cdefs.h b/liba/src/external/msun/include/sys/cdefs.h deleted file mode 100644 index 37bb34d68..000000000 --- a/liba/src/external/msun/include/sys/cdefs.h +++ /dev/null @@ -1,9 +0,0 @@ -#ifndef LIBA_MSUN_SYS_CDEFS_H -#define LIBA_MSUN_SYS_CDEFS_H - -// This file permits standalone compilation of FreeBSD's msun -#define __FBSDID(x) - -#define __strong_reference(x,y) - -#endif diff --git a/liba/src/external/msun/include/sys/types.h b/liba/src/external/msun/include/sys/types.h deleted file mode 100644 index 8f70dfc23..000000000 --- a/liba/src/external/msun/include/sys/types.h +++ /dev/null @@ -1,15 +0,0 @@ -#ifndef LIBA_MSUN_SYS_TYPES_H -#define LIBA_MSUN_SYS_TYPES_H - -#include - -typedef uint32_t u_int32_t; -typedef uint64_t u_int64_t; - -typedef uint16_t __uint16_t; -typedef uint32_t __uint32_t; -typedef uint64_t __uint64_t; - -typedef float float_t; - -#endif diff --git a/liba/src/external/msun/k_logf.h b/liba/src/external/msun/k_logf.h deleted file mode 100644 index 71c547e88..000000000 --- a/liba/src/external/msun/k_logf.h +++ /dev/null @@ -1,39 +0,0 @@ -/* - * ==================================================== - * 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 -__FBSDID("$FreeBSD$"); - -/* - * Float version of k_log.h. See the latter for most comments. - */ - -static const float -/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ -Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */ -Lg2 = 0xccce13.0p-25, /* 0.40000972152 */ -Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */ -Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */ - -static inline float -k_log1pf(float f) -{ - float hfsq,s,z,R,w,t1,t2; - - s = f/((float)2.0+f); - z = s*s; - w = z*z; - t1= w*(Lg2+w*Lg4); - t2= z*(Lg1+w*Lg3); - R = t2+t1; - hfsq=(float)0.5*f*f; - return s*(hfsq+R); -} diff --git a/liba/src/external/msun/math_private.h b/liba/src/external/msun/math_private.h deleted file mode 100644 index afaf201ae..000000000 --- a/liba/src/external/msun/math_private.h +++ /dev/null @@ -1,776 +0,0 @@ -/* - * ==================================================== - * 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. - * ==================================================== - */ - -/* - * from: @(#)fdlibm.h 5.1 93/09/24 - * $FreeBSD$ - */ - -#ifndef _MATH_PRIVATE_H_ -#define _MATH_PRIVATE_H_ - -#include -#include - -/* - * The original fdlibm 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. - */ - -/* - * A union which permits us to convert between a double and two 32 bit - * ints. - */ - -#ifdef __arm__ -#if defined(__VFP_FP__) || defined(__ARM_EABI__) -#define IEEE_WORD_ORDER BYTE_ORDER -#else -#define IEEE_WORD_ORDER BIG_ENDIAN -#endif -#else /* __arm__ */ -#define IEEE_WORD_ORDER BYTE_ORDER -#endif - -#if IEEE_WORD_ORDER == BIG_ENDIAN - -typedef union -{ - double value; - struct - { - u_int32_t msw; - u_int32_t lsw; - } parts; - struct - { - u_int64_t w; - } xparts; -} ieee_double_shape_type; - -#endif - -#if IEEE_WORD_ORDER == LITTLE_ENDIAN - -typedef union -{ - double value; - struct - { - u_int32_t lsw; - u_int32_t msw; - } parts; - struct - { - u_int64_t w; - } xparts; -} 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 a 64-bit int from a double. */ -#define EXTRACT_WORD64(ix,d) \ -do { \ - ieee_double_shape_type ew_u; \ - ew_u.value = (d); \ - (ix) = ew_u.xparts.w; \ -} 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 a double from a 64-bit int. */ -#define INSERT_WORD64(d,ix) \ -do { \ - ieee_double_shape_type iw_u; \ - iw_u.xparts.w = (ix); \ - (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; - /* FIXME: Assumes 32 bit int. */ - unsigned int 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) - -/* - * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long - * double. - */ - -#define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \ -do { \ - union IEEEl2bits ew_u; \ - ew_u.e = (d); \ - (ix0) = ew_u.xbits.expsign; \ - (ix1) = ew_u.xbits.man; \ -} while (0) - -/* - * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit - * long double. - */ - -#define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \ -do { \ - union IEEEl2bits ew_u; \ - ew_u.e = (d); \ - (ix0) = ew_u.xbits.expsign; \ - (ix1) = ew_u.xbits.manh; \ - (ix2) = ew_u.xbits.manl; \ -} while (0) - -/* Get expsign as a 16 bit int from a long double. */ - -#define GET_LDBL_EXPSIGN(i,d) \ -do { \ - union IEEEl2bits ge_u; \ - ge_u.e = (d); \ - (i) = ge_u.xbits.expsign; \ -} while (0) - -/* - * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int - * mantissa. - */ - -#define INSERT_LDBL80_WORDS(d,ix0,ix1) \ -do { \ - union IEEEl2bits iw_u; \ - iw_u.xbits.expsign = (ix0); \ - iw_u.xbits.man = (ix1); \ - (d) = iw_u.e; \ -} while (0) - -/* - * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints - * comprising the mantissa. - */ - -#define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \ -do { \ - union IEEEl2bits iw_u; \ - iw_u.xbits.expsign = (ix0); \ - iw_u.xbits.manh = (ix1); \ - iw_u.xbits.manl = (ix2); \ - (d) = iw_u.e; \ -} while (0) - -/* Set expsign of a long double from a 16 bit int. */ - -#define SET_LDBL_EXPSIGN(d,v) \ -do { \ - union IEEEl2bits se_u; \ - se_u.e = (d); \ - se_u.xbits.expsign = (v); \ - (d) = se_u.e; \ -} while (0) - -#ifdef __i386__ -/* Long double constants are broken on i386. */ -#define LD80C(m, ex, v) { \ - .xbits.man = __CONCAT(m, ULL), \ - .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0), \ -} -#else -/* The above works on non-i386 too, but we use this to check v. */ -#define LD80C(m, ex, v) { .e = (v), } -#endif - -#ifdef FLT_EVAL_METHOD -/* - * Attempt to get strict C99 semantics for assignment with non-C99 compilers. - */ -#if FLT_EVAL_METHOD == 0 || __GNUC__ == 0 -#define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval)) -#else -#define STRICT_ASSIGN(type, lval, rval) do { \ - volatile type __lval; \ - \ - if (sizeof(type) >= sizeof(long double)) \ - (lval) = (rval); \ - else { \ - __lval = (rval); \ - (lval) = __lval; \ - } \ -} while (0) -#endif -#endif /* FLT_EVAL_METHOD */ - -/* Support switching the mode to FP_PE if necessary. */ -#if defined(__i386__) && !defined(NO_FPSETPREC) -#define ENTERI() \ - long double __retval; \ - fp_prec_t __oprec; \ - \ - if ((__oprec = fpgetprec()) != FP_PE) \ - fpsetprec(FP_PE) -#define RETURNI(x) do { \ - __retval = (x); \ - if (__oprec != FP_PE) \ - fpsetprec(__oprec); \ - RETURNF(__retval); \ -} while (0) -#else -#define ENTERI(x) -#define RETURNI(x) RETURNF(x) -#endif - -/* Default return statement if hack*_t() is not used. */ -#define RETURNF(v) return (v) - -/* - * 2sum gives the same result as 2sumF without requiring |a| >= |b| or - * a == 0, but is slower. - */ -#define _2sum(a, b) do { \ - __typeof(a) __s, __w; \ - \ - __w = (a) + (b); \ - __s = __w - (a); \ - (b) = ((a) - (__w - __s)) + ((b) - __s); \ - (a) = __w; \ -} while (0) - -/* - * 2sumF algorithm. - * - * "Normalize" the terms in the infinite-precision expression a + b for - * the sum of 2 floating point values so that b is as small as possible - * relative to 'a'. (The resulting 'a' is the value of the expression in - * the same precision as 'a' and the resulting b is the rounding error.) - * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and - * exponent overflow or underflow must not occur. This uses a Theorem of - * Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum" - * is apparently due to Skewchuk (1997). - * - * For this to always work, assignment of a + b to 'a' must not retain any - * extra precision in a + b. This is required by C standards but broken - * in many compilers. The brokenness cannot be worked around using - * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this - * algorithm would be destroyed by non-null strict assignments. (The - * compilers are correct to be broken -- the efficiency of all floating - * point code calculations would be destroyed similarly if they forced the - * conversions.) - * - * Fortunately, a case that works well can usually be arranged by building - * any extra precision into the type of 'a' -- 'a' should have type float_t, - * double_t or long double. b's type should be no larger than 'a's type. - * Callers should use these types with scopes as large as possible, to - * reduce their own extra-precision and efficiciency problems. In - * particular, they shouldn't convert back and forth just to call here. - */ -#ifdef DEBUG -#define _2sumF(a, b) do { \ - __typeof(a) __w; \ - volatile __typeof(a) __ia, __ib, __r, __vw; \ - \ - __ia = (a); \ - __ib = (b); \ - assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \ - \ - __w = (a) + (b); \ - (b) = ((a) - __w) + (b); \ - (a) = __w; \ - \ - /* The next 2 assertions are weak if (a) is already long double. */ \ - assert((long double)__ia + __ib == (long double)(a) + (b)); \ - __vw = __ia + __ib; \ - __r = __ia - __vw; \ - __r += __ib; \ - assert(__vw == (a) && __r == (b)); \ -} while (0) -#else /* !DEBUG */ -#define _2sumF(a, b) do { \ - __typeof(a) __w; \ - \ - __w = (a) + (b); \ - (b) = ((a) - __w) + (b); \ - (a) = __w; \ -} while (0) -#endif /* DEBUG */ - -/* - * Set x += c, where x is represented in extra precision as a + b. - * x must be sufficiently normalized and sufficiently larger than c, - * and the result is then sufficiently normalized. - * - * The details of ordering are that |a| must be >= |c| (so that (a, c) - * can be normalized without extra work to swap 'a' with c). The details of - * the normalization are that b must be small relative to the normalized 'a'. - * Normalization of (a, c) makes the normalized c tiny relative to the - * normalized a, so b remains small relative to 'a' in the result. However, - * b need not ever be tiny relative to 'a'. For example, b might be about - * 2**20 times smaller than 'a' to give about 20 extra bits of precision. - * That is usually enough, and adding c (which by normalization is about - * 2**53 times smaller than a) cannot change b significantly. However, - * cancellation of 'a' with c in normalization of (a, c) may reduce 'a' - * significantly relative to b. The caller must ensure that significant - * cancellation doesn't occur, either by having c of the same sign as 'a', - * or by having |c| a few percent smaller than |a|. Pre-normalization of - * (a, b) may help. - * - * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2 - * exercise 19). We gain considerable efficiency by requiring the terms to - * be sufficiently normalized and sufficiently increasing. - */ -#define _3sumF(a, b, c) do { \ - __typeof(a) __tmp; \ - \ - __tmp = (c); \ - _2sumF(__tmp, (a)); \ - (b) += (a); \ - (a) = __tmp; \ -} while (0) - -/* - * Common routine to process the arguments to nan(), nanf(), and nanl(). - */ -void _scan_nan(uint32_t *__words, int __num_words, const char *__s); - -#ifdef _COMPLEX_H - -/* - * C99 specifies that complex numbers have the same representation as - * an array of two elements, where the first element is the real part - * and the second element is the imaginary part. - */ -typedef union { - float complex f; - float a[2]; -} float_complex; -typedef union { - double complex f; - double a[2]; -} double_complex; -typedef union { - long double complex f; - long double a[2]; -} long_double_complex; -#define REALPART(z) ((z).a[0]) -#define IMAGPART(z) ((z).a[1]) - -/* - * Inline functions that can be used to construct complex values. - * - * The C99 standard intends x+I*y to be used for this, but x+I*y is - * currently unusable in general since gcc introduces many overflow, - * underflow, sign and efficiency bugs by rewriting I*y as - * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product. - * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted - * to -0.0+I*0.0. - * - * The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL() - * to construct complex values. Compilers that conform to the C99 - * standard require the following functions to avoid the above issues. - */ - -#ifndef CMPLXF -static __inline float complex -CMPLXF(float x, float y) -{ - float_complex z; - - REALPART(z) = x; - IMAGPART(z) = y; - return (z.f); -} -#endif - -#ifndef CMPLX -static __inline double complex -CMPLX(double x, double y) -{ - double_complex z; - - REALPART(z) = x; - IMAGPART(z) = y; - return (z.f); -} -#endif - -#ifndef CMPLXL -static __inline long double complex -CMPLXL(long double x, long double y) -{ - long_double_complex z; - - REALPART(z) = x; - IMAGPART(z) = y; - return (z.f); -} -#endif - -#endif /* _COMPLEX_H */ - -#ifdef __GNUCLIKE_ASM - -/* Asm versions of some functions. */ - -#ifdef __amd64__ -static __inline int -irint(double x) -{ - int n; - - asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x)); - return (n); -} -#define HAVE_EFFICIENT_IRINT -#endif - -#ifdef __i386__ -static __inline int -irint(double x) -{ - int n; - - asm("fistl %0" : "=m" (n) : "t" (x)); - return (n); -} -#define HAVE_EFFICIENT_IRINT -#endif - -#if defined(__amd64__) || defined(__i386__) -static __inline int -irintl(long double x) -{ - int n; - - asm("fistl %0" : "=m" (n) : "t" (x)); - return (n); -} -#define HAVE_EFFICIENT_IRINTL -#endif - -#endif /* __GNUCLIKE_ASM */ - -#ifdef DEBUG -#if defined(__amd64__) || defined(__i386__) -#define breakpoint() asm("int $3") -#else -#include - -#define breakpoint() raise(SIGTRAP) -#endif -#endif - -/* Write a pari script to test things externally. */ -#ifdef DOPRINT -#include - -#ifndef DOPRINT_SWIZZLE -#define DOPRINT_SWIZZLE 0 -#endif - -#ifdef DOPRINT_LD80 - -#define DOPRINT_START(xp) do { \ - uint64_t __lx; \ - uint16_t __hx; \ - \ - /* Hack to give more-problematic args. */ \ - EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \ - __lx ^= DOPRINT_SWIZZLE; \ - INSERT_LDBL80_WORDS(*xp, __hx, __lx); \ - printf("x = %.21Lg; ", (long double)*xp); \ -} while (0) -#define DOPRINT_END1(v) \ - printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) -#define DOPRINT_END2(hi, lo) \ - printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ - (long double)(hi), (long double)(lo)) - -#elif defined(DOPRINT_D64) - -#define DOPRINT_START(xp) do { \ - uint32_t __hx, __lx; \ - \ - EXTRACT_WORDS(__hx, __lx, *xp); \ - __lx ^= DOPRINT_SWIZZLE; \ - INSERT_WORDS(*xp, __hx, __lx); \ - printf("x = %.21Lg; ", (long double)*xp); \ -} while (0) -#define DOPRINT_END1(v) \ - printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) -#define DOPRINT_END2(hi, lo) \ - printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ - (long double)(hi), (long double)(lo)) - -#elif defined(DOPRINT_F32) - -#define DOPRINT_START(xp) do { \ - uint32_t __hx; \ - \ - GET_FLOAT_WORD(__hx, *xp); \ - __hx ^= DOPRINT_SWIZZLE; \ - SET_FLOAT_WORD(*xp, __hx); \ - printf("x = %.21Lg; ", (long double)*xp); \ -} while (0) -#define DOPRINT_END1(v) \ - printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) -#define DOPRINT_END2(hi, lo) \ - printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ - (long double)(hi), (long double)(lo)) - -#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */ - -#ifndef DOPRINT_SWIZZLE_HIGH -#define DOPRINT_SWIZZLE_HIGH 0 -#endif - -#define DOPRINT_START(xp) do { \ - uint64_t __lx, __llx; \ - uint16_t __hx; \ - \ - EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \ - __llx ^= DOPRINT_SWIZZLE; \ - __lx ^= DOPRINT_SWIZZLE_HIGH; \ - INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \ - printf("x = %.36Lg; ", (long double)*xp); \ -} while (0) -#define DOPRINT_END1(v) \ - printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v)) -#define DOPRINT_END2(hi, lo) \ - printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \ - (long double)(hi), (long double)(lo)) - -#endif /* DOPRINT_LD80 */ - -#else /* !DOPRINT */ -#define DOPRINT_START(xp) -#define DOPRINT_END1(v) -#define DOPRINT_END2(hi, lo) -#endif /* DOPRINT */ - -#define RETURNP(x) do { \ - DOPRINT_END1(x); \ - RETURNF(x); \ -} while (0) -#define RETURNPI(x) do { \ - DOPRINT_END1(x); \ - RETURNI(x); \ -} while (0) -#define RETURN2P(x, y) do { \ - DOPRINT_END2((x), (y)); \ - RETURNF((x) + (y)); \ -} while (0) -#define RETURN2PI(x, y) do { \ - DOPRINT_END2((x), (y)); \ - RETURNI((x) + (y)); \ -} while (0) -#ifdef STRUCT_RETURN -#define RETURNSP(rp) do { \ - if (!(rp)->lo_set) \ - RETURNP((rp)->hi); \ - RETURN2P((rp)->hi, (rp)->lo); \ -} while (0) -#define RETURNSPI(rp) do { \ - if (!(rp)->lo_set) \ - RETURNPI((rp)->hi); \ - RETURN2PI((rp)->hi, (rp)->lo); \ -} while (0) -#endif -#define SUM2P(x, y) ({ \ - const __typeof (x) __x = (x); \ - const __typeof (y) __y = (y); \ - \ - DOPRINT_END2(__x, __y); \ - __x + __y; \ -}) - -/* - * ieee style elementary functions - * - * We rename functions here to improve other sources' diffability - * against fdlibm. - */ -#define __ieee754_sqrt sqrt -#define __ieee754_acos acos -#define __ieee754_acosh acosh -#define __ieee754_log log -#define __ieee754_log2 log2 -#define __ieee754_atanh atanh -#define __ieee754_asin asin -#define __ieee754_atan2 atan2 -#define __ieee754_exp exp -#define __ieee754_cosh cosh -#define __ieee754_fmod fmod -#define __ieee754_pow pow -#define __ieee754_lgamma lgamma -#define __ieee754_gamma gamma -#define __ieee754_lgamma_r lgamma_r -#define __ieee754_gamma_r gamma_r -#define __ieee754_log10 log10 -#define __ieee754_sinh sinh -#define __ieee754_hypot hypot -#define __ieee754_j0 j0 -#define __ieee754_j1 j1 -#define __ieee754_y0 y0 -#define __ieee754_y1 y1 -#define __ieee754_jn jn -#define __ieee754_yn yn -#define __ieee754_remainder remainder -#define __ieee754_scalb scalb -#define __ieee754_sqrtf sqrtf -#define __ieee754_acosf acosf -#define __ieee754_acoshf acoshf -#define __ieee754_logf logf -#define __ieee754_atanhf atanhf -#define __ieee754_asinf asinf -#define __ieee754_atan2f atan2f -#define __ieee754_expf expf -#define __ieee754_coshf coshf -#define __ieee754_fmodf fmodf -#define __ieee754_powf powf -#define __ieee754_lgammaf lgammaf -#define __ieee754_gammaf gammaf -#define __ieee754_lgammaf_r lgammaf_r -#define __ieee754_gammaf_r gammaf_r -#define __ieee754_log10f log10f -#define __ieee754_log2f log2f -#define __ieee754_sinhf sinhf -#define __ieee754_hypotf hypotf -#define __ieee754_j0f j0f -#define __ieee754_j1f j1f -#define __ieee754_y0f y0f -#define __ieee754_y1f y1f -#define __ieee754_jnf jnf -#define __ieee754_ynf ynf -#define __ieee754_remainderf remainderf -#define __ieee754_scalbf scalbf - -/* fdlibm kernel function */ -int __kernel_rem_pio2(double*,double*,int,int,int); - -/* double precision kernel functions */ -#ifndef INLINE_REM_PIO2 -int __ieee754_rem_pio2(double,double*); -#endif -double __kernel_sin(double,double,int); -double __kernel_cos(double,double); -double __kernel_tan(double,double,int); -double __ldexp_exp(double,int); -#ifdef _COMPLEX_H -double complex __ldexp_cexp(double complex,int); -#endif - -/* float precision kernel functions */ -#ifndef INLINE_REM_PIO2F -int __ieee754_rem_pio2f(float,double*); -#endif -#ifndef INLINE_KERNEL_SINDF -float __kernel_sindf(double); -#endif -#ifndef INLINE_KERNEL_COSDF -float __kernel_cosdf(double); -#endif -#ifndef INLINE_KERNEL_TANDF -float __kernel_tandf(double,int); -#endif -float __ldexp_expf(float,int); -#ifdef _COMPLEX_H -float complex __ldexp_cexpf(float complex,int); -#endif - -/* long double precision kernel functions */ -long double __kernel_sinl(long double, long double, int); -long double __kernel_cosl(long double, long double); -long double __kernel_tanl(long double, long double, int); - -#endif /* !_MATH_PRIVATE_H_ */ diff --git a/liba/src/external/openbsd/e_log10f.c b/liba/src/external/openbsd/e_log10f.c new file mode 100644 index 000000000..c9e7c0990 --- /dev/null +++ b/liba/src/external/openbsd/e_log10f.c @@ -0,0 +1,51 @@ +/* e_log10f.c -- float version of e_log10.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 "math.h" +#include "math_private.h" + +static const float +two25 = 3.3554432000e+07, /* 0x4c000000 */ +ivln10 = 4.3429449201e-01, /* 0x3ede5bd9 */ +log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */ +log10_2lo = 7.9034151668e-07; /* 0x355427db */ + +static const float zero = 0.0; + +float +log10f(float x) +{ + float y,z; + int32_t i,k,hx; + + GET_FLOAT_WORD(hx,x); + + k=0; + if (hx < 0x00800000) { /* x < 2**-126 */ + if ((hx&0x7fffffff)==0) + return -two25/zero; /* log(+-0)=-inf */ + if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ + k -= 25; x *= two25; /* subnormal number, scale up x */ + GET_FLOAT_WORD(hx,x); + } + if (hx >= 0x7f800000) return x+x; + k += (hx>>23)-127; + i = ((u_int32_t)k&0x80000000)>>31; + hx = (hx&0x007fffff)|((0x7f-i)<<23); + y = (float)(k+i); + SET_FLOAT_WORD(x,hx); + z = y*log10_2lo + ivln10*logf(x); + return z+y*log10_2hi; +} diff --git a/liba/src/external/openbsd/e_logf.c b/liba/src/external/openbsd/e_logf.c new file mode 100644 index 000000000..543f33a0d --- /dev/null +++ b/liba/src/external/openbsd/e_logf.c @@ -0,0 +1,81 @@ +/* e_logf.c -- float version of e_log.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 "math.h" +#include "math_private.h" + +static const float +ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ +ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ +two25 = 3.355443200e+07, /* 0x4c000000 */ +Lg1 = 6.6666668653e-01, /* 3F2AAAAB */ +Lg2 = 4.0000000596e-01, /* 3ECCCCCD */ +Lg3 = 2.8571429849e-01, /* 3E924925 */ +Lg4 = 2.2222198546e-01, /* 3E638E29 */ +Lg5 = 1.8183572590e-01, /* 3E3A3325 */ +Lg6 = 1.5313838422e-01, /* 3E1CD04F */ +Lg7 = 1.4798198640e-01; /* 3E178897 */ + +static const float zero = 0.0; + +float +logf(float x) +{ + float hfsq,f,s,z,R,w,t1,t2,dk; + int32_t k,ix,i,j; + + GET_FLOAT_WORD(ix,x); + + k=0; + if (ix < 0x00800000) { /* x < 2**-126 */ + if ((ix&0x7fffffff)==0) + return -two25/zero; /* log(+-0)=-inf */ + if (ix<0) return (x-x)/zero; /* log(-#) = NaN */ + k -= 25; x *= two25; /* subnormal number, scale up x */ + GET_FLOAT_WORD(ix,x); + } + if (ix >= 0x7f800000) return x+x; + k += (ix>>23)-127; + ix &= 0x007fffff; + i = (ix+(0x95f64<<3))&0x800000; + SET_FLOAT_WORD(x,ix|(i^0x3f800000)); /* normalize x or x/2 */ + k += (i>>23); + f = x-(float)1.0; + if((0x007fffff&(15+ix))<16) { /* |f| < 2**-20 */ + if(f==zero) {if(k==0) return zero; else {dk=(float)k; + return dk*ln2_hi+dk*ln2_lo;}} + R = f*f*((float)0.5-(float)0.33333333333333333*f); + if(k==0) return f-R; else {dk=(float)k; + return dk*ln2_hi-((R-dk*ln2_lo)-f);} + } + s = f/((float)2.0+f); + dk = (float)k; + z = s*s; + i = ix-(0x6147a<<3); + w = z*z; + j = (0x6b851<<3)-ix; + t1= w*(Lg2+w*(Lg4+w*Lg6)); + t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); + i |= j; + R = t2+t1; + if(i>0) { + hfsq=(float)0.5*f*f; + if(k==0) return f-(hfsq-s*(hfsq+R)); else + return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); + } else { + if(k==0) return f-s*(f-R); else + return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); + } +} diff --git a/liba/src/external/msun/e_powf.c b/liba/src/external/openbsd/e_powf.c similarity index 81% rename from liba/src/external/msun/e_powf.c rename to liba/src/external/openbsd/e_powf.c index 5c4647887..3c4269d71 100644 --- a/liba/src/external/msun/e_powf.c +++ b/liba/src/external/openbsd/e_powf.c @@ -8,17 +8,16 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#include -__FBSDID("$FreeBSD$"); - #include "math.h" #include "math_private.h" +static const volatile float huge = 1.0e+30, tiny = 1.0e-30; + static const float bp[] = {1.0, 1.5,}, dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */ @@ -27,8 +26,6 @@ zero = 0.0, one = 1.0, two = 2.0, two24 = 16777216.0, /* 0x4b800000 */ -huge = 1.0e30, -tiny = 1.0e-30, /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ L1 = 6.0000002384e-01, /* 0x3f19999a */ L2 = 4.2857143283e-01, /* 0x3edb6db7 */ @@ -46,17 +43,17 @@ lg2_h = 6.93145752e-01, /* 0x3f317200 */ lg2_l = 1.42860654e-06, /* 0x35bfbe8c */ ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */ cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */ -cp_h = 9.6191406250e-01, /* 0x3f764000 =12b cp */ -cp_l = -1.1736857402e-04, /* 0xb8f623c6 =tail of cp_h */ +cp_h = 9.6179199219e-01, /* 0x3f763800 =head of cp */ +cp_l = 4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */ ivln2 = 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */ ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ float -__ieee754_powf(float x, float y) +powf(float x, float y) { float z,ax,z_h,z_l,p_h,p_l; - float y1,t1,t2,r,s,sn,t,u,v,w; + float yy1,t1,t2,r,s,t,u,v,w; int32_t i,j,k,yisint,n; int32_t hx,hy,ix,iy,is; @@ -65,15 +62,15 @@ __ieee754_powf(float x, float y) ix = hx&0x7fffffff; iy = hy&0x7fffffff; /* y==zero: x**0 = 1 */ - if(iy==0) return one; + if(iy==0) return one; /* x==1: 1**y = 1, even if y is NaN */ if (hx==0x3f800000) return one; - /* y!=zero: result is NaN if either arg is NaN */ + /* +-NaN return x+y */ if(ix > 0x7f800000 || iy > 0x7f800000) - return (x+0.0F)+(y+0.0F); + return x+y; /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer @@ -81,14 +78,14 @@ __ieee754_powf(float x, float y) * yisint = 2 ... y is an even int */ yisint = 0; - if(hx<0) { + if(hx<0) { if(iy>=0x4b800000) yisint = 2; /* even integer y */ else if(iy>=0x3f800000) { k = (iy>>23)-0x7f; /* exponent */ j = iy>>(23-k); if((j<<(23-k))==iy) yisint = 2-(j&1); - } - } + } + } /* special value of y */ if (iy==0x7f800000) { /* y is +-inf */ @@ -98,14 +95,14 @@ __ieee754_powf(float x, float y) return (hy>=0)? y: zero; else /* (|x|<1)**-,+inf = inf,0 */ return (hy<0)?-y: zero; - } + } if(iy==0x3f800000) { /* y is +-1 */ if(hy<0) return one/x; else return x; } if(hy==0x40000000) return x*x; /* y is 2 */ if(hy==0x3f000000) { /* y is 0.5 */ if(hx>=0) /* x >= +0 */ - return __ieee754_sqrtf(x); + return sqrtf(x); } ax = fabsf(x); @@ -116,28 +113,23 @@ __ieee754_powf(float x, float y) if(hx<0) { if(((ix-0x3f800000)|yisint)==0) { z = (z-z)/(z-z); /* (-1)**non-int is NaN */ - } else if(yisint==1) + } else if(yisint==1) z = -z; /* (x<0)**odd = -(|x|**odd) */ } return z; } - - n = ((u_int32_t)hx>>31)-1; - + /* (x<0)**(non-int) is NaN */ - if((n|yisint)==0) return (x-x)/(x-x); - - sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */ - if((n|(yisint-1))==0) sn = -one;/* (-ve)**(odd int) */ + if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x); /* |y| is huge */ if(iy>0x4d000000) { /* if |y| > 2**27 */ /* over/underflow if x is not close to one */ - if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny; - if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny; - /* now |1-x| is tiny <= 2**-20, suffice to compute + if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny; + if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny; + /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ - t = ax-1; /* t has 20 trailing zeros */ + t = ax-one; /* t has 20 trailing zeros */ w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25)); u = ivln2_h*t; /* ivln2_h has 16 sig. bits */ v = t*ivln2_l-w*ivln2; @@ -168,8 +160,7 @@ __ieee754_powf(float x, float y) GET_FLOAT_WORD(is,s_h); SET_FLOAT_WORD(s_h,is&0xfffff000); /* t_h=ax+bp[k] High */ - is = ((ix>>1)&0xfffff000)|0x20000000; - SET_FLOAT_WORD(t_h,is+0x00400000+(k<<21)); + SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x0040000+(k<<21)); t_l = ax - (t_h-bp[k]); s_l = v*((u-s_h*t_h)-s_h*t_l); /* compute log(ax) */ @@ -199,22 +190,26 @@ __ieee754_powf(float x, float y) t2 = z_l-(((t1-t)-dp_h[k])-z_h); } - /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ + s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ + if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0) + s = -one; /* (-ve)**(odd int) */ + + /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */ GET_FLOAT_WORD(is,y); - SET_FLOAT_WORD(y1,is&0xfffff000); - p_l = (y-y1)*t1+y*t2; - p_h = y1*t1; + SET_FLOAT_WORD(yy1,is&0xfffff000); + p_l = (y-yy1)*t1+y*t2; + p_h = yy1*t1; z = p_l+p_h; GET_FLOAT_WORD(j,z); if (j>0x43000000) /* if z > 128 */ - return sn*huge*huge; /* overflow */ + return s*huge*huge; /* overflow */ else if (j==0x43000000) { /* if z == 128 */ - if(p_l+ovt>z-p_h) return sn*huge*huge; /* overflow */ + if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ } else if ((j&0x7fffffff)>0x43160000) /* z <= -150 */ - return sn*tiny*tiny; /* underflow */ + return s*tiny*tiny; /* underflow */ else if (j==0xc3160000){ /* z == -150 */ - if(p_l<=z-p_h) return sn*tiny*tiny; /* underflow */ + if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ } /* * compute 2**(p_h+p_l) @@ -229,10 +224,10 @@ __ieee754_powf(float x, float y) n = ((n&0x007fffff)|0x00800000)>>(23-k); if(j<0) n = -n; p_h -= t; - } + } t = p_l+p_h; GET_FLOAT_WORD(is,t); - SET_FLOAT_WORD(t,is&0xffff8000); + SET_FLOAT_WORD(t,is&0xfffff000); u = t*lg2_h; v = (p_l-(t-p_h))*lg2+t*lg2_l; z = u+v; @@ -245,5 +240,5 @@ __ieee754_powf(float x, float y) j += (n<<23); if((j>>23)<=0) z = scalbnf(z,n); /* subnormal output */ else SET_FLOAT_WORD(z,j); - return sn*z; + return s*z; } diff --git a/liba/src/external/msun/e_sqrtf.c b/liba/src/external/openbsd/e_sqrtf.c similarity index 85% rename from liba/src/external/msun/e_sqrtf.c rename to liba/src/external/openbsd/e_sqrtf.c index 7eba4d07f..80d35278c 100644 --- a/liba/src/external/msun/e_sqrtf.c +++ b/liba/src/external/openbsd/e_sqrtf.c @@ -8,35 +8,31 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#ifndef lint -static char rcsid[] = "$FreeBSD$"; -#endif - #include "math.h" #include "math_private.h" static const float one = 1.0, tiny=1.0e-30; float -__ieee754_sqrtf(float x) +sqrtf(float x) { float z; - int32_t sign = (int)0x80000000; + int32_t sign = (int)0x80000000; int32_t ix,s,q,m,t,i; u_int32_t r; GET_FLOAT_WORD(ix,x); /* take care of Inf and NaN */ - if((ix&0x7f800000)==0x7f800000) { + if((ix&0x7f800000)==0x7f800000) { return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf sqrt(-inf)=sNaN */ - } + } /* take care of zero */ if(ix<=0) { if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */ @@ -61,12 +57,12 @@ __ieee754_sqrtf(float x) r = 0x01000000; /* r = moving bit from right to left */ while(r!=0) { - t = s+r; - if(t<=ix) { - s = t+r; - ix -= t; - q += r; - } + t = s+r; + if(t<=ix) { + s = t+r; + ix -= t; + q += r; + } ix += ix; r>>=1; } diff --git a/liba/src/external/openbsd/include/sys/types.h b/liba/src/external/openbsd/include/sys/types.h new file mode 100644 index 000000000..0ca13ccbf --- /dev/null +++ b/liba/src/external/openbsd/include/sys/types.h @@ -0,0 +1,13 @@ +#ifndef LIBA_OPENBSD_SYS_TYPES_H +#define LIBA_OPENBSD_SYS_TYPES_H + +#include + +#define LITTLE_ENDIAN 0x1234 +#define BIG_ENDIAN 0x4321 +#define BYTE_ORDER LITTLE_ENDIAN + +typedef uint32_t u_int32_t; +typedef uint64_t u_int64_t; + +#endif diff --git a/liba/src/external/openbsd/math_private.h b/liba/src/external/openbsd/math_private.h new file mode 100644 index 000000000..5bece0cbe --- /dev/null +++ b/liba/src/external/openbsd/math_private.h @@ -0,0 +1,417 @@ +/* $OpenBSD: math_private.h,v 1.17 2014/06/02 19:31:17 kettenis Exp $ */ +/* + * ==================================================== + * 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. + * ==================================================== + */ + +/* + * from: @(#)fdlibm.h 5.1 93/09/24 + */ + +#ifndef _MATH_PRIVATE_H_ +#define _MATH_PRIVATE_H_ + +#include + +/* The original fdlibm 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. */ + +/* A union which permits us to convert between a long double and + four 32 bit ints. */ + +#if BYTE_ORDER == BIG_ENDIAN + +typedef union +{ + long double value; + struct { + u_int32_t mswhi; + u_int32_t mswlo; + u_int32_t lswhi; + u_int32_t lswlo; + } parts32; + struct { + u_int64_t msw; + u_int64_t lsw; + } parts64; +} ieee_quad_shape_type; + +#endif + +#if BYTE_ORDER == LITTLE_ENDIAN + +typedef union +{ + long double value; + struct { + u_int32_t lswlo; + u_int32_t lswhi; + u_int32_t mswlo; + u_int32_t mswhi; + } parts32; + struct { + u_int64_t lsw; + u_int64_t msw; + } parts64; +} ieee_quad_shape_type; + +#endif + +/* Get two 64 bit ints from a long double. */ + +#define GET_LDOUBLE_WORDS64(ix0,ix1,d) \ +do { \ + ieee_quad_shape_type qw_u; \ + qw_u.value = (d); \ + (ix0) = qw_u.parts64.msw; \ + (ix1) = qw_u.parts64.lsw; \ +} while (0) + +/* Set a long double from two 64 bit ints. */ + +#define SET_LDOUBLE_WORDS64(d,ix0,ix1) \ +do { \ + ieee_quad_shape_type qw_u; \ + qw_u.parts64.msw = (ix0); \ + qw_u.parts64.lsw = (ix1); \ + (d) = qw_u.value; \ +} while (0) + +/* Get the more significant 64 bits of a long double mantissa. */ + +#define GET_LDOUBLE_MSW64(v,d) \ +do { \ + ieee_quad_shape_type sh_u; \ + sh_u.value = (d); \ + (v) = sh_u.parts64.msw; \ +} while (0) + +/* Set the more significant 64 bits of a long double mantissa from an int. */ + +#define SET_LDOUBLE_MSW64(d,v) \ +do { \ + ieee_quad_shape_type sh_u; \ + sh_u.value = (d); \ + sh_u.parts64.msw = (v); \ + (d) = sh_u.value; \ +} while (0) + +/* Get the least significant 64 bits of a long double mantissa. */ + +#define GET_LDOUBLE_LSW64(v,d) \ +do { \ + ieee_quad_shape_type sh_u; \ + sh_u.value = (d); \ + (v) = sh_u.parts64.lsw; \ +} while (0) + +/* A union which permits us to convert between a long double and + three 32 bit ints. */ + +#if BYTE_ORDER == BIG_ENDIAN + +typedef union +{ + long double value; + struct { +#ifdef __LP64__ + int padh:32; +#endif + int exp:16; + int padl:16; + u_int32_t msw; + u_int32_t lsw; + } parts; +} ieee_extended_shape_type; + +#endif + +#if BYTE_ORDER == LITTLE_ENDIAN + +typedef union +{ + long double value; + struct { + u_int32_t lsw; + u_int32_t msw; + int exp:16; + int padl:16; +#ifdef __LP64__ + int padh:32; +#endif + } parts; +} ieee_extended_shape_type; + +#endif + +/* Get three 32 bit ints from a double. */ + +#define GET_LDOUBLE_WORDS(se,ix0,ix1,d) \ +do { \ + ieee_extended_shape_type ew_u; \ + ew_u.value = (d); \ + (se) = ew_u.parts.exp; \ + (ix0) = ew_u.parts.msw; \ + (ix1) = ew_u.parts.lsw; \ +} while (0) + +/* Set a double from two 32 bit ints. */ + +#define SET_LDOUBLE_WORDS(d,se,ix0,ix1) \ +do { \ + ieee_extended_shape_type iw_u; \ + iw_u.parts.exp = (se); \ + iw_u.parts.msw = (ix0); \ + iw_u.parts.lsw = (ix1); \ + (d) = iw_u.value; \ +} while (0) + +/* Get the more significant 32 bits of a long double mantissa. */ + +#define GET_LDOUBLE_MSW(v,d) \ +do { \ + ieee_extended_shape_type sh_u; \ + sh_u.value = (d); \ + (v) = sh_u.parts.msw; \ +} while (0) + +/* Set the more significant 32 bits of a long double mantissa from an int. */ + +#define SET_LDOUBLE_MSW(d,v) \ +do { \ + ieee_extended_shape_type sh_u; \ + sh_u.value = (d); \ + sh_u.parts.msw = (v); \ + (d) = sh_u.value; \ +} while (0) + +/* Get int from the exponent of a long double. */ + +#define GET_LDOUBLE_EXP(se,d) \ +do { \ + ieee_extended_shape_type ge_u; \ + ge_u.value = (d); \ + (se) = ge_u.parts.exp; \ +} while (0) + +/* Set exponent of a long double from an int. */ + +#define SET_LDOUBLE_EXP(d,se) \ +do { \ + ieee_extended_shape_type se_u; \ + se_u.value = (d); \ + se_u.parts.exp = (se); \ + (d) = se_u.value; \ +} while (0) + +/* A union which permits us to convert between a double and two 32 bit + ints. */ + +/* + * The arm port is little endian except for the FP word order which is + * big endian. + */ + +#if (BYTE_ORDER == BIG_ENDIAN) || (defined(__arm__) && !defined(__VFP_FP__)) + +typedef union +{ + double value; + struct + { + u_int32_t msw; + u_int32_t lsw; + } parts; +} ieee_double_shape_type; + +#endif + +#if (BYTE_ORDER == LITTLE_ENDIAN) && !(defined(__arm__) && !defined(__VFP_FP__)) + +typedef union +{ + double value; + struct + { + u_int32_t lsw; + u_int32_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; + u_int32_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) + +#ifdef FLT_EVAL_METHOD +/* + * Attempt to get strict C99 semantics for assignment with non-C99 compilers. + */ +#if FLT_EVAL_METHOD == 0 || __GNUC__ == 0 +#define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval)) +#else /* FLT_EVAL_METHOD == 0 || __GNUC__ == 0 */ +#define STRICT_ASSIGN(type, lval, rval) do { \ + volatile type __lval; \ + \ + if (sizeof(type) >= sizeof(long double)) \ + (lval) = (rval); \ + else { \ + __lval = (rval); \ + (lval) = __lval; \ + } \ +} while (0) +#endif /* FLT_EVAL_METHOD == 0 || __GNUC__ == 0 */ +#endif /* FLT_EVAL_METHOD */ + +/* fdlibm kernel function */ +extern int __ieee754_rem_pio2(double,double*); +extern double __kernel_sin(double,double,int); +extern double __kernel_cos(double,double); +extern double __kernel_tan(double,double,int); +extern int __kernel_rem_pio2(double*,double*,int,int,int); + +/* float versions of fdlibm kernel functions */ +extern int __ieee754_rem_pio2f(float,float*); +extern float __kernel_sinf(float,float,int); +extern float __kernel_cosf(float,float); +extern float __kernel_tanf(float,float,int); +extern int __kernel_rem_pio2f(float*,float*,int,int,int,const int*); + +/* long double precision kernel functions */ +long double __kernel_sinl(long double, long double, int); +long double __kernel_cosl(long double, long double); +long double __kernel_tanl(long double, long double, int); + +/* + * Common routine to process the arguments to nan(), nanf(), and nanl(). + */ +void _scan_nan(uint32_t *__words, int __num_words, const char *__s); + +/* + * TRUNC() is a macro that sets the trailing 27 bits in the mantissa + * of an IEEE double variable to zero. It must be expression-like + * for syntactic reasons, and we implement this expression using + * an inline function instead of a pure macro to avoid depending + * on the gcc feature of statement-expressions. + */ +#define TRUNC(d) (_b_trunc(&(d))) + +static __inline void +_b_trunc(volatile double *_dp) +{ + uint32_t _lw; + + GET_LOW_WORD(_lw, *_dp); + SET_LOW_WORD(*_dp, _lw & 0xf8000000); +} + +struct Double { + double a; + double b; +}; + +/* + * Functions internal to the math package, yet not static. + */ +double __exp__D(double, double); +struct Double __log__D(double); +long double __p1evll(long double, void *, int); +long double __polevll(long double, void *, int); + +#endif /* _MATH_PRIVATE_H_ */ diff --git a/liba/src/external/msun/s_copysignf.c b/liba/src/external/openbsd/s_copysignf.c similarity index 88% rename from liba/src/external/msun/s_copysignf.c rename to liba/src/external/openbsd/s_copysignf.c index 05ca1e368..267ed2da9 100644 --- a/liba/src/external/msun/s_copysignf.c +++ b/liba/src/external/openbsd/s_copysignf.c @@ -8,14 +8,11 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#include -__FBSDID("$FreeBSD$"); - /* * copysignf(float x, float y) * copysignf(x,y) returns a value with the magnitude of x and diff --git a/liba/src/external/msun/s_fabsf.c b/liba/src/external/openbsd/s_fabsf.c similarity index 86% rename from liba/src/external/msun/s_fabsf.c rename to liba/src/external/openbsd/s_fabsf.c index e9383d0db..71b9f755e 100644 --- a/liba/src/external/msun/s_fabsf.c +++ b/liba/src/external/openbsd/s_fabsf.c @@ -8,14 +8,11 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#include -__FBSDID("$FreeBSD$"); - /* * fabsf(x) returns the absolute value of x. */ diff --git a/liba/src/external/msun/s_scalbnf.c b/liba/src/external/openbsd/s_scalbnf.c similarity index 85% rename from liba/src/external/msun/s_scalbnf.c rename to liba/src/external/openbsd/s_scalbnf.c index 7666c74ab..5ab9208b5 100644 --- a/liba/src/external/msun/s_scalbnf.c +++ b/liba/src/external/openbsd/s_scalbnf.c @@ -8,17 +8,11 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ -#ifndef lint -static char rcsid[] = "$FreeBSD$"; -#endif - -#include - #include "math.h" #include "math_private.h" @@ -29,7 +23,7 @@ huge = 1.0e+30, tiny = 1.0e-30; float -scalbnf (float x, int n) +scalbnf(float x, int n) { int32_t k,ix; GET_FLOAT_WORD(ix,x); @@ -38,11 +32,11 @@ scalbnf (float x, int n) if ((ix&0x7fffffff)==0) return x; /* +-0 */ x *= two25; GET_FLOAT_WORD(ix,x); - k = ((ix&0x7f800000)>>23) - 25; + k = ((ix&0x7f800000)>>23) - 25; if (n< -50000) return tiny*x; /*underflow*/ } if (k==0xff) return x+x; /* NaN or Inf */ - k = k+n; + 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;} @@ -55,4 +49,8 @@ scalbnf (float x, int n) return x*twom25; } -__strong_reference(scalbnf, ldexpf); +float +ldexpf(float x, int n) +{ + return scalbnf(x, n); +} diff --git a/liba/src/external/openbsd/s_sinf.c b/liba/src/external/openbsd/s_sinf.c new file mode 100644 index 000000000..7bb90e67f --- /dev/null +++ b/liba/src/external/openbsd/s_sinf.c @@ -0,0 +1,45 @@ +/* s_sinf.c -- float version of s_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 "math.h" +#include "math_private.h" + +float +sinf(float x) +{ + 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_sinf(x,z,0); + + /* sin(Inf or NaN) is NaN */ + else if (ix>=0x7f800000) return x-x; + + /* argument reduction needed */ + else { + n = __ieee754_rem_pio2f(x,y); + switch(n&3) { + case 0: return __kernel_sinf(y[0],y[1],1); + case 1: return __kernel_cosf(y[0],y[1]); + case 2: return -__kernel_sinf(y[0],y[1],1); + default: + return -__kernel_cosf(y[0],y[1]); + } + } +}