From 5d91a421b99fee96e1eda7bd383d7e685d803eeb Mon Sep 17 00:00:00 2001 From: Alexander von Gluck IV Date: Fri, 14 Aug 2015 23:37:49 -0500 Subject: [PATCH] glibc/arm: More generic math functions --- .../libroot/posix/glibc/arch/arm/Jamfile | 11 +- .../libroot/posix/glibc/arch/generic/e_logf.c | 38 ++----- .../libroot/posix/glibc/arch/generic/e_logl.c | 14 +++ .../posix/glibc/arch/generic/e_sqrtl.c | 104 ++++++++++++++++++ .../posix/glibc/arch/generic/s_atanl.c | 14 +++ 5 files changed, 152 insertions(+), 29 deletions(-) create mode 100644 src/system/libroot/posix/glibc/arch/generic/e_logl.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/e_sqrtl.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/s_atanl.c diff --git a/src/system/libroot/posix/glibc/arch/arm/Jamfile b/src/system/libroot/posix/glibc/arch/arm/Jamfile index 9bfb08399a..18924a7e2f 100644 --- a/src/system/libroot/posix/glibc/arch/arm/Jamfile +++ b/src/system/libroot/posix/glibc/arch/arm/Jamfile @@ -35,6 +35,7 @@ local genericSources = s_signbit.c s_signbitf.c s_signbitl.c s_nan.c s_nanf.c s_nanl.c + e_hypot.c e_hypotf.c e_hypotl.c w_hypot.c w_hypotf.c w_hypotl.c s_fpclassify.c s_fpclassifyf.c s_clog.c s_clogf.c s_clogl.c @@ -54,12 +55,16 @@ local genericSources = halfulp.c mpa.c mplog.c mpexp.c s_sin.c - s_atan.c s_atanf.c + s_atan.c s_atanf.c s_atanl.c s_tan.c - e_asin.c w_asin.c + e_asin.c e_asinl.c + w_asin.c w_asinl.c e_log10.c w_log10.c + e_logf.c e_logl.c e_acos.c w_acos.c - e_atan2.c w_atan2.c w_atan2l.c mpatan2.c mpatan.c mptan.c mpsqrt.c w_sqrt.c w_sqrtf.c + e_atan2.c e_atan2l.c + w_atan2.c w_atan2l.c mpatan2.c mpatan.c mptan.c mpsqrt.c w_sqrt.c w_sqrtf.c + e_sqrtl.c e_fmod.c w_fmod.c e_log.c w_log.c e_cosh.c w_cosh.c diff --git a/src/system/libroot/posix/glibc/arch/generic/e_logf.c b/src/system/libroot/posix/glibc/arch/generic/e_logf.c index de8f869df4..cf75e11781 100644 --- a/src/system/libroot/posix/glibc/arch/generic/e_logf.c +++ b/src/system/libroot/posix/glibc/arch/generic/e_logf.c @@ -13,18 +13,10 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_logf.c,v 1.4 1995/05/10 20:45:54 jtc Exp $"; -#endif +#include +#include -#include "math.h" -#include "math_private.h" - -#ifdef __STDC__ static const float -#else -static float -#endif ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ two25 = 3.355443200e+07, /* 0x4c000000 */ @@ -36,18 +28,10 @@ Lg5 = 1.8183572590e-01, /* 3E3A3325 */ Lg6 = 1.5313838422e-01, /* 3E1CD04F */ Lg7 = 1.4798198640e-01; /* 3E178897 */ -#ifdef __STDC__ static const float zero = 0.0; -#else -static float zero = 0.0; -#endif -#ifdef __STDC__ - float __ieee754_logf(float x) -#else - float __ieee754_logf(x) - float x; -#endif +float +__ieee754_logf(float x) { float hfsq,f,s,z,R,w,t1,t2,dk; int32_t k,ix,i,j; @@ -56,13 +40,14 @@ static float zero = 0.0; k=0; if (ix < 0x00800000) { /* x < 2**-126 */ - if ((ix&0x7fffffff)==0) - return -two25/(x-x); /* log(+-0)=-inf */ - if (ix<0) return (x-x)/(x-x); /* log(-#) = NaN */ + if (__builtin_expect((ix&0x7fffffff)==0, 0)) + return -two25/zero; /* log(+-0)=-inf */ + if (__builtin_expect(ix<0, 0)) + return (x-x)/(x-x); /* log(-#) = NaN */ k -= 25; x *= two25; /* subnormal number, scale up x */ GET_FLOAT_WORD(ix,x); } - if (ix >= 0x7f800000) return x+x; + if (__builtin_expect(ix >= 0x7f800000, 0)) return x+x; k += (ix>>23)-127; ix &= 0x007fffff; i = (ix+(0x95f64<<3))&0x800000; @@ -76,9 +61,9 @@ static float zero = 0.0; } 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);} + return dk*ln2_hi-((R-dk*ln2_lo)-f);} } - s = f/((float)2.0+f); + s = f/((float)2.0+f); dk = (float)k; z = s*s; i = ix-(0x6147a<<3); @@ -97,3 +82,4 @@ static float zero = 0.0; return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); } } +strong_alias (__ieee754_logf, __logf_finite) diff --git a/src/system/libroot/posix/glibc/arch/generic/e_logl.c b/src/system/libroot/posix/glibc/arch/generic/e_logl.c new file mode 100644 index 0000000000..7a4ea1b07f --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/e_logl.c @@ -0,0 +1,14 @@ +#include +#include +#include + +long double +__ieee754_logl (long double x) +{ + fputs ("__ieee754_logl not implemented\n", stderr); + __set_errno (ENOSYS); + return 0.0; +} +strong_alias (__ieee754_logl, __logl_finite) + +stub_warning (logl) diff --git a/src/system/libroot/posix/glibc/arch/generic/e_sqrtl.c b/src/system/libroot/posix/glibc/arch/generic/e_sqrtl.c new file mode 100644 index 0000000000..ac1b22175e --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/e_sqrtl.c @@ -0,0 +1,104 @@ +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; either version 2.1 of the License, or + * (at your option) any later version. + * + * 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, see . + */ +/*********************************************************************/ +/* MODULE_NAME: uroot.c */ +/* */ +/* FUNCTION: usqrt */ +/* */ +/* FILES NEEDED: dla.h endian.h mydefs.h uroot.h */ +/* uroot.tbl */ +/* */ +/* An ultimate sqrt routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of square */ +/* root of x. */ +/* Assumption: Machine arithmetic operations are performed in */ +/* round to nearest mode of IEEE 754 standard. */ +/* */ +/*********************************************************************/ + +#include + +typedef union {int64_t i[2]; long double x; double d[2]; } mynumber; + +static const double + t512 = 0x1p512, + tm256 = 0x1p-256, + two54 = 0x1p54, /* 0x4350000000000000 */ + twom54 = 0x1p-54; /* 0x3C90000000000000 */ + +/*********************************************************************/ +/* An ultimate sqrt routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of square */ +/* root of x. */ +/*********************************************************************/ +long double __ieee754_sqrtl(long double x) +{ + static const long double big = 134217728.0, big1 = 134217729.0; + long double t,s,i; + mynumber a,c; + uint64_t k, l; + int64_t m, n; + double d; + + a.x=x; + k=a.i[0] & INT64_C(0x7fffffffffffffff); + /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/ + if (k>INT64_C(0x000fffff00000000) && k> 53; + m = (a.i[1] >> 52) & 0x7ff; + if (m == 0) { + a.d[1] *= two54; + m = ((a.i[1] >> 52) & 0x7ff) - 54; + } + m += n; + if (m > 0) + a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52); + else if (m <= -54) { + a.i[1] &= INT64_C(0x8000000000000000); + } else { + m += 54; + a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52); + a.d[1] *= twom54; + } + } + a.i[0] = l; + s = a.x; + d = __ieee754_sqrt (a.d[0]); + c.i[0] = INT64_C(0x2000000000000000)+((k&INT64_C(0x7fe0000000000000))>>1); + c.i[1] = 0; + i = d; + t = 0.5L * (i + s / i); + i = 0.5L * (t + s / t); + return c.x * i; + } + else { + if (k>=INT64_C(0x7ff0000000000000)) { + if (a.i[0] == INT64_C(0xfff0000000000000)) + return (big1-big1)/(big-big); /* sqrt (-Inf) = NaN. */ + return x; /* sqrt (NaN) = NaN, sqrt (+Inf) = +Inf. */ + } + if (x == 0) return x; + if (x < 0) return (big1-big1)/(big-big); + return tm256*__ieee754_sqrtl(x*t512); + } +} +strong_alias (__ieee754_sqrtl, __sqrtl_finite) diff --git a/src/system/libroot/posix/glibc/arch/generic/s_atanl.c b/src/system/libroot/posix/glibc/arch/generic/s_atanl.c new file mode 100644 index 0000000000..2957d702d5 --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/s_atanl.c @@ -0,0 +1,14 @@ +#include +#include +#include + +long double +__atanl (long double x) +{ + fputs ("__atanl not implemented\n", stderr); + __set_errno (ENOSYS); + return 0.0; +} +weak_alias (__atanl, atanl) + +stub_warning (atanl)