From 379f1467353183335fc90b76bd09141cef93d9e3 Mon Sep 17 00:00:00 2001 From: Ingo Weinhold Date: Fri, 9 Dec 2005 14:16:39 +0000 Subject: [PATCH] Added more missing math functions for PPC from glibc 2.3.2. git-svn-id: file:///srv/svn/repos/haiku/haiku/trunk@15445 a95241bf-73f2-0310-859d-f6bbb57e9c96 --- .../libroot/posix/glibc/arch/generic/atnat.h | 167 ++++++ .../posix/glibc/arch/generic/fraiseexcpt.c | 66 +++ .../libroot/posix/glibc/arch/generic/mpa.c | 507 ++++++++++++++++++ .../libroot/posix/glibc/arch/generic/mpa2.h | 95 ++++ .../libroot/posix/glibc/arch/generic/mpatan.c | 102 ++++ .../libroot/posix/glibc/arch/generic/mpatan.h | 174 ++++++ .../posix/glibc/arch/generic/mpatan2.c | 70 +++ .../libroot/posix/glibc/arch/generic/mpsqrt.c | 103 ++++ .../libroot/posix/glibc/arch/generic/mpsqrt.h | 51 ++ .../libroot/posix/glibc/arch/generic/s_atan.c | 230 ++++++++ .../posix/glibc/arch/generic/s_atanf.c | 120 +++++ .../libroot/posix/glibc/arch/ppc/Jamfile | 14 +- .../libroot/posix/glibc/arch/ppc/e_sqrtf.c | 1 + .../posix/glibc/arch/ppc/fraiseexcpt.c | 66 +++ .../posix/glibc/arch/ppc/ftestexcept.c | 33 ++ .../libroot/posix/glibc/arch/ppc/t_sqrt.c | 144 +++++ .../libroot/posix/glibc/arch/ppc/w_sqrt.c | 146 +++++ .../libroot/posix/glibc/arch/ppc/w_sqrtf.c | 136 +++++ .../posix/glibc/include/arch/ppc/bits/fenv.h | 145 +++++ .../glibc/include/arch/ppc/bits/fenvinline.h | 59 ++ .../posix/glibc/include/arch/ppc/fenv_libc.h | 106 ++++ 21 files changed, 2532 insertions(+), 3 deletions(-) create mode 100644 src/system/libroot/posix/glibc/arch/generic/atnat.h create mode 100644 src/system/libroot/posix/glibc/arch/generic/fraiseexcpt.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/mpa.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/mpa2.h create mode 100644 src/system/libroot/posix/glibc/arch/generic/mpatan.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/mpatan.h create mode 100644 src/system/libroot/posix/glibc/arch/generic/mpatan2.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/mpsqrt.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/mpsqrt.h create mode 100644 src/system/libroot/posix/glibc/arch/generic/s_atan.c create mode 100644 src/system/libroot/posix/glibc/arch/generic/s_atanf.c create mode 100644 src/system/libroot/posix/glibc/arch/ppc/e_sqrtf.c create mode 100644 src/system/libroot/posix/glibc/arch/ppc/fraiseexcpt.c create mode 100644 src/system/libroot/posix/glibc/arch/ppc/ftestexcept.c create mode 100644 src/system/libroot/posix/glibc/arch/ppc/t_sqrt.c create mode 100644 src/system/libroot/posix/glibc/arch/ppc/w_sqrt.c create mode 100644 src/system/libroot/posix/glibc/arch/ppc/w_sqrtf.c create mode 100644 src/system/libroot/posix/glibc/include/arch/ppc/bits/fenv.h create mode 100644 src/system/libroot/posix/glibc/include/arch/ppc/bits/fenvinline.h create mode 100644 src/system/libroot/posix/glibc/include/arch/ppc/fenv_libc.h diff --git a/src/system/libroot/posix/glibc/arch/generic/atnat.h b/src/system/libroot/posix/glibc/arch/generic/atnat.h new file mode 100644 index 0000000000..72f67f01fc --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/atnat.h @@ -0,0 +1,167 @@ +/* + * IBM Accurate Mathematical Library + * Written by International Business Machines Corp. + * Copyright (C) 2001 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ + +/************************************************************************/ +/* MODULE_NAME: atnat.h */ +/* */ +/* */ +/* common data and variables definition for BIG or LITTLE ENDIAN */ +/************************************************************************/ +#ifndef ATNAT_H +#define ATNAT_H + +#define M 4 + +#ifdef BIG_ENDI + static const number + /* polynomial I */ +/**/ d3 = {{0xbfd55555, 0x55555555} }, /* -0.333... */ +/**/ d5 = {{0x3fc99999, 0x999997fd} }, /* 0.199... */ +/**/ d7 = {{0xbfc24924, 0x923f7603} }, /* -0.142... */ +/**/ d9 = {{0x3fbc71c6, 0xe5129a3b} }, /* 0.111... */ +/**/ d11 = {{0xbfb74580, 0x22b13c25} }, /* -0.090... */ +/**/ d13 = {{0x3fb375f0, 0x8b31cbce} }, /* 0.076... */ + /* polynomial II */ +/**/ f3 = {{0xbfd55555, 0x55555555} }, /* -1/3 */ +/**/ ff3 = {{0xbc755555, 0x55555555} }, /* -1/3-f3 */ +/**/ f5 = {{0x3fc99999, 0x9999999a} }, /* 1/5 */ +/**/ ff5 = {{0xbc699999, 0x9999999a} }, /* 1/5-f5 */ +/**/ f7 = {{0xbfc24924, 0x92492492} }, /* -1/7 */ +/**/ ff7 = {{0xbc624924, 0x92492492} }, /* -1/7-f7 */ +/**/ f9 = {{0x3fbc71c7, 0x1c71c71c} }, /* 1/9 */ +/**/ ff9 = {{0x3c5c71c7, 0x1c71c71c} }, /* 1/9-f9 */ +/**/ f11 = {{0xbfb745d1, 0x745d1746} }, /* -1/11 */ +/**/ f13 = {{0x3fb3b13b, 0x13b13b14} }, /* 1/13 */ +/**/ f15 = {{0xbfb11111, 0x11111111} }, /* -1/15 */ +/**/ f17 = {{0x3fae1e1e, 0x1e1e1e1e} }, /* 1/17 */ +/**/ f19 = {{0xbfaaf286, 0xbca1af28} }, /* -1/19 */ + /* constants */ +/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */ +/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ a = {{0x3e4bb67a, 0x00000000} }, /* 1.290e-8 */ +/**/ b = {{0x3fb00000, 0x00000000} }, /* 1/16 */ +/**/ c = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ d = {{0x40300000, 0x00000000} }, /* 16 */ +/**/ e = {{0x43349ff2, 0x00000000} }, /* 5.805e15 */ +/**/ hpi = {{0x3ff921fb, 0x54442d18} }, /* pi/2 */ +/**/ mhpi = {{0xbff921fb, 0x54442d18} }, /* -pi/2 */ +/**/ hpi1 = {{0x3c91a626, 0x33145c07} }, /* pi/2-hpi */ +/**/ u1 = {{0x3c2d3382, 0x00000000} }, /* 7.915e-19 */ +/**/ u21 = {{0x3c6dffc0, 0x00000000} }, /* 1.301e-17 */ +/**/ u22 = {{0x3c527bd0, 0x00000000} }, /* 4.008e-18 */ +/**/ u23 = {{0x3c3cd057, 0x00000000} }, /* 1.562e-18 */ +/**/ u24 = {{0x3c329cdf, 0x00000000} }, /* 1.009e-18 */ +/**/ u31 = {{0x3c3a1edf, 0x00000000} }, /* 1.416e-18 */ +/**/ u32 = {{0x3c33f0e1, 0x00000000} }, /* 1.081e-18 */ +/**/ u4 = {{0x3bf955e4, 0x00000000} }, /* 8.584e-20 */ +/**/ u5 = {{0x3aaef2d1, 0x00000000} }, /* 5e-26 */ +/**/ u6 = {{0x3a98c56d, 0x00000000} }, /* 2.001e-26 */ +/**/ u7 = {{0x3a9375de, 0x00000000} }, /* 1.572e-26 */ +/**/ u8 = {{0x3a6eeb36, 0x00000000} }, /* 3.122e-27 */ +/**/ u9[M] ={{{0x38c1aa5b, 0x00000000} }, /* 2.658e-35 */ +/**/ {{0x35c1aa4d, 0x00000000} }, /* 9.443e-50 */ +/**/ {{0x32c1aa88, 0x00000000} }, /* 3.355e-64 */ +/**/ {{0x11c1aa56, 0x00000000} }},/* 3.818e-223 */ +/**/ two8 = {{0x40700000, 0x00000000} }, /* 2**8=256 */ +/**/ two52 = {{0x43300000, 0x00000000} }; /* 2**52 */ + +#else +#ifdef LITTLE_ENDI + static const number + /* polynomial I */ +/**/ d3 = {{0x55555555, 0xbfd55555} }, /* -0.333... */ +/**/ d5 = {{0x999997fd, 0x3fc99999} }, /* 0.199... */ +/**/ d7 = {{0x923f7603, 0xbfc24924} }, /* -0.142... */ +/**/ d9 = {{0xe5129a3b, 0x3fbc71c6} }, /* 0.111... */ +/**/ d11 = {{0x22b13c25, 0xbfb74580} }, /* -0.090... */ +/**/ d13 = {{0x8b31cbce, 0x3fb375f0} }, /* 0.076... */ + /* polynomial II */ +/**/ f3 = {{0x55555555, 0xbfd55555} }, /* -1/3 */ +/**/ ff3 = {{0x55555555, 0xbc755555} }, /* -1/3-f3 */ +/**/ f5 = {{0x9999999a, 0x3fc99999} }, /* 1/5 */ +/**/ ff5 = {{0x9999999a, 0xbc699999} }, /* 1/5-f5 */ +/**/ f7 = {{0x92492492, 0xbfc24924} }, /* -1/7 */ +/**/ ff7 = {{0x92492492, 0xbc624924} }, /* -1/7-f7 */ +/**/ f9 = {{0x1c71c71c, 0x3fbc71c7} }, /* 1/9 */ +/**/ ff9 = {{0x1c71c71c, 0x3c5c71c7} }, /* 1/9-f9 */ +/**/ f11 = {{0x745d1746, 0xbfb745d1} }, /* -1/11 */ +/**/ f13 = {{0x13b13b14, 0x3fb3b13b} }, /* 1/13 */ +/**/ f15 = {{0x11111111, 0xbfb11111} }, /* -1/15 */ +/**/ f17 = {{0x1e1e1e1e, 0x3fae1e1e} }, /* 1/17 */ +/**/ f19 = {{0xbca1af28, 0xbfaaf286} }, /* -1/19 */ + /* constants */ +/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */ +/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ a = {{0x00000000, 0x3e4bb67a} }, /* 1.290e-8 */ +/**/ b = {{0x00000000, 0x3fb00000} }, /* 1/16 */ +/**/ c = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ d = {{0x00000000, 0x40300000} }, /* 16 */ +/**/ e = {{0x00000000, 0x43349ff2} }, /* 5.805e15 */ +/**/ hpi = {{0x54442d18, 0x3ff921fb} }, /* pi/2 */ +/**/ mhpi = {{0x54442d18, 0xbff921fb} }, /* -pi/2 */ +/**/ hpi1 = {{0x33145c07, 0x3c91a626} }, /* pi/2-hpi */ +/**/ u1 = {{0x00000000, 0x3c2d3382} }, /* 7.915e-19 */ +/**/ u21 = {{0x00000000, 0x3c6dffc0} }, /* 1.301e-17 */ +/**/ u22 = {{0x00000000, 0x3c527bd0} }, /* 4.008e-18 */ +/**/ u23 = {{0x00000000, 0x3c3cd057} }, /* 1.562e-18 */ +/**/ u24 = {{0x00000000, 0x3c329cdf} }, /* 1.009e-18 */ +/**/ u31 = {{0x00000000, 0x3c3a1edf} }, /* 1.416e-18 */ +/**/ u32 = {{0x00000000, 0x3c33f0e1} }, /* 1.081e-18 */ +/**/ u4 = {{0x00000000, 0x3bf955e4} }, /* 8.584e-20 */ +/**/ u5 = {{0x00000000, 0x3aaef2d1} }, /* 5e-26 */ +/**/ u6 = {{0x00000000, 0x3a98c56d} }, /* 2.001e-26 */ +/**/ u7 = {{0x00000000, 0x3a9375de} }, /* 1.572e-26 */ +/**/ u8 = {{0x00000000, 0x3a6eeb36} }, /* 3.122e-27 */ +/**/ u9[M] ={{{0x00000000, 0x38c1aa5b} }, /* 2.658e-35 */ +/**/ {{0x00000000, 0x35c1aa4d} }, /* 9.443e-50 */ +/**/ {{0x00000000, 0x32c1aa88} }, /* 3.355e-64 */ +/**/ {{0x00000000, 0x11c1aa56} }},/* 3.818e-223 */ +/**/ two8 = {{0x00000000, 0x40700000} }, /* 2**8=256 */ +/**/ two52 = {{0x00000000, 0x43300000} }; /* 2**52 */ + +#endif +#endif + +#define ZERO zero.d +#define ONE one.d +#define A a.d +#define B b.d +#define C c.d +#define D d.d +#define E e.d +#define HPI hpi.d +#define MHPI mhpi.d +#define HPI1 hpi1.d +#define U1 u1.d +#define U21 u21.d +#define U22 u22.d +#define U23 u23.d +#define U24 u24.d +#define U31 u31.d +#define U32 u32.d +#define U4 u4.d +#define U5 u5.d +#define U6 u6.d +#define U7 u7.d +#define U8 u8.d +#define TWO8 two8.d +#define TWO52 two52.d + +#endif diff --git a/src/system/libroot/posix/glibc/arch/generic/fraiseexcpt.c b/src/system/libroot/posix/glibc/arch/generic/fraiseexcpt.c new file mode 100644 index 0000000000..dbe36c3d5a --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/fraiseexcpt.c @@ -0,0 +1,66 @@ +/* Raise given exceptions. + Copyright (C) 1997,99,2000,01,02 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library 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. + + The GNU C Library 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 the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include +#include + +#undef feraiseexcept +int +__feraiseexcept (int excepts) +{ + fenv_union_t u; + + /* Raise exceptions represented by EXCEPTS. It is the responsibility of + the OS to ensure that if multiple exceptions occur they are fed back + to this process in the proper way; this can happen in hardware, + anyway (in particular, inexact with overflow or underflow). */ + + /* Get the current state. */ + u.fenv = fegetenv_register (); + + /* Add the exceptions */ + u.l[1] = (u.l[1] + | (excepts & FPSCR_STICKY_BITS) + /* Turn FE_INVALID into FE_INVALID_SOFTWARE. */ + | (excepts >> ((31 - FPSCR_VX) - (31 - FPSCR_VXSOFT)) + & FE_INVALID_SOFTWARE)); + + /* Store the new status word (along with the rest of the environment), + triggering any appropriate exceptions. */ + fesetenv_register (u.fenv); + + if ((excepts & FE_INVALID) + /* For some reason, some PowerPC chips (the 601, in particular) + don't have FE_INVALID_SOFTWARE implemented. Detect this + case and raise FE_INVALID_SNAN instead. */ + && !fetestexcept (FE_INVALID)) + set_fpscr_bit (FPSCR_VXSNAN); + + /* Success. */ + return 0; +} + +#include +#if SHLIB_COMPAT (libm, GLIBC_2_1, GLIBC_2_2) +strong_alias (__feraiseexcept, __old_feraiseexcept) +compat_symbol (libm, BP_SYM (__old_feraiseexcept), BP_SYM (feraiseexcept), GLIBC_2_1); +#endif + +libm_hidden_ver (__feraiseexcept, feraiseexcept) +versioned_symbol (libm, BP_SYM (__feraiseexcept), BP_SYM (feraiseexcept), GLIBC_2_2); diff --git a/src/system/libroot/posix/glibc/arch/generic/mpa.c b/src/system/libroot/posix/glibc/arch/generic/mpa.c new file mode 100644 index 0000000000..acd7c4102f --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/mpa.c @@ -0,0 +1,507 @@ + +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001 Free Software Foundation + * + * 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ +/************************************************************************/ +/* MODULE_NAME: mpa.c */ +/* */ +/* FUNCTIONS: */ +/* mcr */ +/* acr */ +/* cr */ +/* cpy */ +/* cpymn */ +/* norm */ +/* denorm */ +/* mp_dbl */ +/* dbl_mp */ +/* add_magnitudes */ +/* sub_magnitudes */ +/* add */ +/* sub */ +/* mul */ +/* inv */ +/* dvd */ +/* */ +/* Arithmetic functions for multiple precision numbers. */ +/* Relative errors are bounded */ +/************************************************************************/ + + +#include "endian.h" +#include "mpa.h" +#include "mpa2.h" +/* mcr() compares the sizes of the mantissas of two multiple precision */ +/* numbers. Mantissas are compared regardless of the signs of the */ +/* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also */ +/* disregarded. */ +static int mcr(const mp_no *x, const mp_no *y, int p) { + int i; + for (i=1; i<=p; i++) { + if (X[i] == Y[i]) continue; + else if (X[i] > Y[i]) return 1; + else return -1; } + return 0; +} + + + +/* acr() compares the absolute values of two multiple precision numbers */ +int __acr(const mp_no *x, const mp_no *y, int p) { + int i; + + if (X[0] == ZERO) { + if (Y[0] == ZERO) i= 0; + else i=-1; + } + else if (Y[0] == ZERO) i= 1; + else { + if (EX > EY) i= 1; + else if (EX < EY) i=-1; + else i= mcr(x,y,p); + } + + return i; +} + + +/* cr90 compares the values of two multiple precision numbers */ +int __cr(const mp_no *x, const mp_no *y, int p) { + int i; + + if (X[0] > Y[0]) i= 1; + else if (X[0] < Y[0]) i=-1; + else if (X[0] < ZERO ) i= __acr(y,x,p); + else i= __acr(x,y,p); + + return i; +} + + +/* Copy a multiple precision number. Set *y=*x. x=y is permissible. */ +void __cpy(const mp_no *x, mp_no *y, int p) { + int i; + + EY = EX; + for (i=0; i <= p; i++) Y[i] = X[i]; + + return; +} + + +/* Copy a multiple precision number x of precision m into a */ +/* multiple precision number y of precision n. In case n>m, */ +/* the digits of y beyond the m'th are set to zero. In case */ +/* n= 2**(-1022))) */ +static void norm(const mp_no *x, double *y, int p) +{ + #define R radixi.d + int i; +#if 0 + int k; +#endif + double a,c,u,v,z[5]; + if (p<5) { + if (p==1) c = X[1]; + else if (p==2) c = X[1] + R* X[2]; + else if (p==3) c = X[1] + R*(X[2] + R* X[3]); + else if (p==4) c =(X[1] + R* X[2]) + R*R*(X[3] + R*X[4]); + } + else { + for (a=ONE, z[1]=X[1]; z[1] < TWO23; ) + {a *= TWO; z[1] *= TWO; } + + for (i=2; i<5; i++) { + z[i] = X[i]*a; + u = (z[i] + CUTTER)-CUTTER; + if (u > z[i]) u -= RADIX; + z[i] -= u; + z[i-1] += u*RADIXI; + } + + u = (z[3] + TWO71) - TWO71; + if (u > z[3]) u -= TWO19; + v = z[3]-u; + + if (v == TWO18) { + if (z[4] == ZERO) { + for (i=5; i <= p; i++) { + if (X[i] == ZERO) continue; + else {z[3] += ONE; break; } + } + } + else z[3] += ONE; + } + + c = (z[1] + R *(z[2] + R * z[3]))/a; + } + + c *= X[0]; + + for (i=1; iEX; i--) c *= RADIXI; + + *y = c; + return; +#undef R +} + +/* Convert a multiple precision number *x into a double precision */ +/* number *y, denormalized case (|x| < 2**(-1022))) */ +static void denorm(const mp_no *x, double *y, int p) +{ + int i,k; + double c,u,z[5]; +#if 0 + double a,v; +#endif + +#define R radixi.d + if (EX<-44 || (EX==-44 && X[1] z[3]) u -= TWO5; + + if (u==z[3]) { + for (i=k+1; i <= p; i++) { + if (X[i] == ZERO) continue; + else {z[3] += ONE; break; } + } + } + + c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10); + + *y = c*TWOM1032; + return; + +#undef R +} + +/* Convert a multiple precision number *x into a double precision number *y. */ +/* The result is correctly rounded to the nearest/even. *x is left unchanged */ + +void __mp_dbl(const mp_no *x, double *y, int p) { +#if 0 + int i,k; + double a,c,u,v,z[5]; +#endif + + if (X[0] == ZERO) {*y = ZERO; return; } + + if (EX> -42) norm(x,y,p); + else if (EX==-42 && X[1]>=TWO10) norm(x,y,p); + else denorm(x,y,p); +} + + +/* dbl_mp() converts a double precision number x into a multiple precision */ +/* number *y. If the precision p is too small the result is truncated. x is */ +/* left unchanged. */ + +void __dbl_mp(double x, mp_no *y, int p) { + + int i,n; + double u; + + /* Sign */ + if (x == ZERO) {Y[0] = ZERO; return; } + else if (x > ZERO) Y[0] = ONE; + else {Y[0] = MONE; x=-x; } + + /* Exponent */ + for (EY=ONE; x >= RADIX; EY += ONE) x *= RADIXI; + for ( ; x < ONE; EY -= ONE) x *= RADIX; + + /* Digits */ + n=MIN(p,4); + for (i=1; i<=n; i++) { + u = (x + TWO52) - TWO52; + if (u>x) u -= ONE; + Y[i] = u; x -= u; x *= RADIX; } + for ( ; i<=p; i++) Y[i] = ZERO; + return; +} + + +/* add_magnitudes() adds the magnitudes of *x & *y assuming that */ +/* abs(*x) >= abs(*y) > 0. */ +/* The sign of the sum *z is undefined. x&y may overlap but not x&z or y&z. */ +/* No guard digit is used. The result equals the exact sum, truncated. */ +/* *x & *y are left unchanged. */ + +static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { + + int i,j,k; + + EZ = EX; + + i=p; j=p+ EY - EX; k=p+1; + + if (j<1) + {__cpy(x,z,p); return; } + else Z[k] = ZERO; + + for (; j>0; i--,j--) { + Z[k] += X[i] + Y[j]; + if (Z[k] >= RADIX) { + Z[k] -= RADIX; + Z[--k] = ONE; } + else + Z[--k] = ZERO; + } + + for (; i>0; i--) { + Z[k] += X[i]; + if (Z[k] >= RADIX) { + Z[k] -= RADIX; + Z[--k] = ONE; } + else + Z[--k] = ZERO; + } + + if (Z[1] == ZERO) { + for (i=1; i<=p; i++) Z[i] = Z[i+1]; } + else EZ += ONE; +} + + +/* sub_magnitudes() subtracts the magnitudes of *x & *y assuming that */ +/* abs(*x) > abs(*y) > 0. */ +/* The sign of the difference *z is undefined. x&y may overlap but not x&z */ +/* or y&z. One guard digit is used. The error is less than one ulp. */ +/* *x & *y are left unchanged. */ + +static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { + + int i,j,k; + + EZ = EX; + + if (EX == EY) { + i=j=k=p; + Z[k] = Z[k+1] = ZERO; } + else { + j= EX - EY; + if (j > p) {__cpy(x,z,p); return; } + else { + i=p; j=p+1-j; k=p; + if (Y[j] > ZERO) { + Z[k+1] = RADIX - Y[j--]; + Z[k] = MONE; } + else { + Z[k+1] = ZERO; + Z[k] = ZERO; j--;} + } + } + + for (; j>0; i--,j--) { + Z[k] += (X[i] - Y[j]); + if (Z[k] < ZERO) { + Z[k] += RADIX; + Z[--k] = MONE; } + else + Z[--k] = ZERO; + } + + for (; i>0; i--) { + Z[k] += X[i]; + if (Z[k] < ZERO) { + Z[k] += RADIX; + Z[--k] = MONE; } + else + Z[--k] = ZERO; + } + + for (i=1; Z[i] == ZERO; i++) ; + EZ = EZ - i + 1; + for (k=1; i <= p+1; ) + Z[k++] = Z[i++]; + for (; k <= p; ) + Z[k++] = ZERO; + + return; +} + + +/* Add two multiple precision numbers. Set *z = *x + *y. x&y may overlap */ +/* but not x&z or y&z. One guard digit is used. The error is less than */ +/* one ulp. *x & *y are left unchanged. */ + +void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { + + int n; + + if (X[0] == ZERO) {__cpy(y,z,p); return; } + else if (Y[0] == ZERO) {__cpy(x,z,p); return; } + + if (X[0] == Y[0]) { + if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; } + else {add_magnitudes(y,x,z,p); Z[0] = Y[0]; } + } + else { + if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; } + else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; } + else Z[0] = ZERO; + } + return; +} + + +/* Subtract two multiple precision numbers. *z is set to *x - *y. x&y may */ +/* overlap but not x&z or y&z. One guard digit is used. The error is */ +/* less than one ulp. *x & *y are left unchanged. */ + +void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { + + int n; + + if (X[0] == ZERO) {__cpy(y,z,p); Z[0] = -Z[0]; return; } + else if (Y[0] == ZERO) {__cpy(x,z,p); return; } + + if (X[0] != Y[0]) { + if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; } + else {add_magnitudes(y,x,z,p); Z[0] = -Y[0]; } + } + else { + if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; } + else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = -Y[0]; } + else Z[0] = ZERO; + } + return; +} + + +/* Multiply two multiple precision numbers. *z is set to *x * *y. x&y */ +/* may overlap but not x&z or y&z. In case p=1,2,3 the exact result is */ +/* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp. */ +/* *x & *y are left unchanged. */ + +void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { + + int i, i1, i2, j, k, k2; + double u; + + /* Is z=0? */ + if (X[0]*Y[0]==ZERO) + { Z[0]=ZERO; return; } + + /* Multiply, add and carry */ + k2 = (p<3) ? p+p : p+3; + Z[k2]=ZERO; + for (k=k2; k>1; ) { + if (k > p) {i1=k-p; i2=p+1; } + else {i1=1; i2=k; } + for (i=i1,j=i2-1; i Z[k]) u -= RADIX; + Z[k] -= u; + Z[--k] = u*RADIXI; + } + + /* Is there a carry beyond the most significant digit? */ + if (Z[1] == ZERO) { + for (i=1; i<=p; i++) Z[i]=Z[i+1]; + EZ = EX + EY - 1; } + else + EZ = EX + EY; + + Z[0] = X[0] * Y[0]; + return; +} + + +/* Invert a multiple precision number. Set *y = 1 / *x. */ +/* Relative error bound = 1.001*r**(1-p) for p=2, 1.063*r**(1-p) for p=3, */ +/* 2.001*r**(1-p) for p>3. */ +/* *x=0 is not permissible. *x is left unchanged. */ + +void __inv(const mp_no *x, mp_no *y, int p) { + int i; +#if 0 + int l; +#endif + double t; + mp_no z,w; + static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3, + 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4}; + const mp_no mptwo = {1,{1.0,2.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + + __cpy(x,&z,p); z.e=0; __mp_dbl(&z,&t,p); + t=ONE/t; __dbl_mp(t,y,p); EY -= EX; + + for (i=0; i3. *y=0 is not permissible. */ + +void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) { + + mp_no w; + + if (X[0] == ZERO) Z[0] = ZERO; + else {__inv(y,&w,p); __mul(x,&w,z,p);} + return; +} diff --git a/src/system/libroot/posix/glibc/arch/generic/mpa2.h b/src/system/libroot/posix/glibc/arch/generic/mpa2.h new file mode 100644 index 0000000000..ba93e8fc13 --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/mpa2.h @@ -0,0 +1,95 @@ + +/* + * IBM Accurate Mathematical Library + * Written by International Business Machines Corp. + * Copyright (C) 2001 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ + +/**************************************************************************/ +/* */ +/* MODULE_NAME:mpa2.h */ +/* */ +/* */ +/* variables prototype and definition according to type of processor */ +/* types definition */ +/**************************************************************************/ + +#ifndef MPA2_H +#define MPA2_H + + +#ifdef BIG_ENDI +static const number +/**/ radix = {{0x41700000, 0x00000000} }, /* 2**24 */ +/**/ radixi = {{0x3e700000, 0x00000000} }, /* 2**-24 */ +/**/ cutter = {{0x44b00000, 0x00000000} }, /* 2**76 */ +/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */ +/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ mone = {{0xbff00000, 0x00000000} }, /* -1 */ +/**/ two = {{0x40000000, 0x00000000} }, /* 2 */ +/**/ two5 = {{0x40400000, 0x00000000} }, /* 2**5 */ +/**/ two10 = {{0x40900000, 0x00000000} }, /* 2**10 */ +/**/ two18 = {{0x41100000, 0x00000000} }, /* 2**18 */ +/**/ two19 = {{0x41200000, 0x00000000} }, /* 2**19 */ +/**/ two23 = {{0x41600000, 0x00000000} }, /* 2**23 */ +/**/ two52 = {{0x43300000, 0x00000000} }, /* 2**52 */ +/**/ two57 = {{0x43800000, 0x00000000} }, /* 2**57 */ +/**/ two71 = {{0x44600000, 0x00000000} }, /* 2**71 */ +/**/ twom1032 = {{0x00000400, 0x00000000} }; /* 2**-1032 */ + +#else +#ifdef LITTLE_ENDI +static const number +/**/ radix = {{0x00000000, 0x41700000} }, /* 2**24 */ +/**/ radixi = {{0x00000000, 0x3e700000} }, /* 2**-24 */ +/**/ cutter = {{0x00000000, 0x44b00000} }, /* 2**76 */ +/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */ +/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ mone = {{0x00000000, 0xbff00000} }, /* -1 */ +/**/ two = {{0x00000000, 0x40000000} }, /* 2 */ +/**/ two5 = {{0x00000000, 0x40400000} }, /* 2**5 */ +/**/ two10 = {{0x00000000, 0x40900000} }, /* 2**10 */ +/**/ two18 = {{0x00000000, 0x41100000} }, /* 2**18 */ +/**/ two19 = {{0x00000000, 0x41200000} }, /* 2**19 */ +/**/ two23 = {{0x00000000, 0x41600000} }, /* 2**23 */ +/**/ two52 = {{0x00000000, 0x43300000} }, /* 2**52 */ +/**/ two57 = {{0x00000000, 0x43800000} }, /* 2**57 */ +/**/ two71 = {{0x00000000, 0x44600000} }, /* 2**71 */ +/**/ twom1032 = {{0x00000000, 0x00000400} }; /* 2**-1032 */ + +#endif +#endif + +#define RADIX radix.d +#define RADIXI radixi.d +#define CUTTER cutter.d +#define ZERO zero.d +#define ONE one.d +#define MONE mone.d +#define TWO two.d +#define TWO5 two5.d +#define TWO10 two10.d +#define TWO18 two18.d +#define TWO19 two19.d +#define TWO23 two23.d +#define TWO52 two52.d +#define TWO57 two57.d +#define TWO71 two71.d +#define TWOM1032 twom1032.d + + +#endif diff --git a/src/system/libroot/posix/glibc/arch/generic/mpatan.c b/src/system/libroot/posix/glibc/arch/generic/mpatan.c new file mode 100644 index 0000000000..ee21c25138 --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/mpatan.c @@ -0,0 +1,102 @@ + +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001 Free Software Foundation + * + * 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ +/******************************************************************/ +/* */ +/* MODULE_NAME:mpatan.c */ +/* */ +/* FUNCTIONS:mpatan */ +/* */ +/* FILES NEEDED: mpa.h endian.h mpatan.h */ +/* mpa.c */ +/* */ +/* Multi-Precision Atan function subroutine, for precision p >= 4.*/ +/* The relative error of the result is bounded by 34.32*r**(1-p), */ +/* where r=2**24. */ +/******************************************************************/ + +#include "endian.h" +#include "mpa.h" +void __mpsqrt(mp_no *, mp_no *, int); + +void __mpatan(mp_no *x, mp_no *y, int p) { +#include "mpatan.h" + + int i,m,n; + double dx; + mp_no + mpone = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, + mptwo = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, + mptwoim1 = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + + mp_no mps,mpsm,mpt,mpt1,mpt2,mpt3; + + /* Choose m and initiate mpone, mptwo & mptwoim1 */ + if (EX>0) m=7; + else if (EX<0) m=0; + else { + __mp_dbl(x,&dx,p); dx=ABS(dx); + for (m=6; m>0; m--) + {if (dx>xm[m].d) break;} + } + mpone.e = mptwo.e = mptwoim1.e = 1; + mpone.d[0] = mpone.d[1] = mptwo.d[0] = mptwoim1.d[0] = ONE; + mptwo.d[1] = TWO; + + /* Reduce x m times */ + __mul(x,x,&mpsm,p); + if (m==0) __cpy(x,&mps,p); + else { + for (i=0; i1; i--) { + mptwoim1.d[1] -= TWO; + __dvd(&mpsm,&mptwoim1,&mpt1,p); + __mul(&mpsm,&mpt,&mpt2,p); + __sub(&mpt1,&mpt2,&mpt,p); + } + __mul(&mps,&mpt,&mpt1,p); + __sub(&mps,&mpt1,&mpt,p); + + /* Compute Atan(x) */ + mptwoim1.d[1] = twom[m].d; + __mul(&mptwoim1,&mpt,y,p); + + return; +} diff --git a/src/system/libroot/posix/glibc/arch/generic/mpatan.h b/src/system/libroot/posix/glibc/arch/generic/mpatan.h new file mode 100644 index 0000000000..d420ff3408 --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/mpatan.h @@ -0,0 +1,174 @@ + +/* + * IBM Accurate Mathematical Library + * Written by International Business Machines Corp. + * Copyright (C) 2001 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ + +/******************************************************************/ +/* */ +/* MODULE_NAME:mpatan.h */ +/* */ +/* common data and variables prototype and definition */ +/******************************************************************/ + +#ifndef MPATAN_H +#define MPATAN_H + +#ifdef BIG_ENDI + static const number + xm[8] = { /* x[m] */ +/**/ {{0x00000000, 0x00000000} }, /* 0.0 */ +/**/ {{0x3f8930be, 0x00000000} }, /* 0.0123 */ +/**/ {{0x3f991687, 0x00000000} }, /* 0.0245 */ +/**/ {{0x3fa923a2, 0x00000000} }, /* 0.0491 */ +/**/ {{0x3fb930be, 0x00000000} }, /* 0.0984 */ +/**/ {{0x3fc95810, 0x00000000} }, /* 0.198 */ +/**/ {{0x3fda7ef9, 0x00000000} }, /* 0.414 */ +/**/ {{0x3ff00000, 0x00000000} }, /* 1.0 */ + }; + static const number + twonm1[33] = { /* 2n-1 */ +/**/ {{0x00000000, 0x00000000} }, /* 0 */ +/**/ {{0x00000000, 0x00000000} }, /* 0 */ +/**/ {{0x00000000, 0x00000000} }, /* 0 */ +/**/ {{0x00000000, 0x00000000} }, /* 0 */ +/**/ {{0x40260000, 0x00000000} }, /* 11 */ +/**/ {{0x402e0000, 0x00000000} }, /* 15 */ +/**/ {{0x40330000, 0x00000000} }, /* 19 */ +/**/ {{0x40350000, 0x00000000} }, /* 21 */ +/**/ {{0x40390000, 0x00000000} }, /* 25 */ +/**/ {{0x403d0000, 0x00000000} }, /* 29 */ +/**/ {{0x40408000, 0x00000000} }, /* 33 */ +/**/ {{0x40428000, 0x00000000} }, /* 37 */ +/**/ {{0x40448000, 0x00000000} }, /* 41 */ +/**/ {{0x40468000, 0x00000000} }, /* 45 */ +/**/ {{0x40488000, 0x00000000} }, /* 49 */ +/**/ {{0x404a8000, 0x00000000} }, /* 53 */ +/**/ {{0x404b8000, 0x00000000} }, /* 55 */ +/**/ {{0x404d8000, 0x00000000} }, /* 59 */ +/**/ {{0x404f8000, 0x00000000} }, /* 63 */ +/**/ {{0x4050c000, 0x00000000} }, /* 67 */ +/**/ {{0x4051c000, 0x00000000} }, /* 71 */ +/**/ {{0x4052c000, 0x00000000} }, /* 75 */ +/**/ {{0x4053c000, 0x00000000} }, /* 79 */ +/**/ {{0x4054c000, 0x00000000} }, /* 83 */ +/**/ {{0x40554000, 0x00000000} }, /* 85 */ +/**/ {{0x40564000, 0x00000000} }, /* 89 */ +/**/ {{0x40574000, 0x00000000} }, /* 93 */ +/**/ {{0x40584000, 0x00000000} }, /* 97 */ +/**/ {{0x40594000, 0x00000000} }, /* 101 */ +/**/ {{0x405a4000, 0x00000000} }, /* 105 */ +/**/ {{0x405b4000, 0x00000000} }, /* 109 */ +/**/ {{0x405c4000, 0x00000000} }, /* 113 */ +/**/ {{0x405d4000, 0x00000000} }, /* 117 */ + }; + + static const number + twom[8] = { /* 2**m */ +/**/ {{0x3ff00000, 0x00000000} }, /* 1.0 */ +/**/ {{0x40000000, 0x00000000} }, /* 2.0 */ +/**/ {{0x40100000, 0x00000000} }, /* 4.0 */ +/**/ {{0x40200000, 0x00000000} }, /* 8.0 */ +/**/ {{0x40300000, 0x00000000} }, /* 16.0 */ +/**/ {{0x40400000, 0x00000000} }, /* 32.0 */ +/**/ {{0x40500000, 0x00000000} }, /* 64.0 */ +/**/ {{0x40600000, 0x00000000} }, /* 128.0 */ + }; + + static const number +/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ two = {{0x40000000, 0x00000000} }; /* 2 */ + +#else +#ifdef LITTLE_ENDI + + static const number + xm[8] = { /* x[m] */ +/**/ {{0x00000000, 0x00000000} }, /* 0.0 */ +/**/ {{0x00000000, 0x3f8930be} }, /* 0.0123 */ +/**/ {{0x00000000, 0x3f991687} }, /* 0.0245 */ +/**/ {{0x00000000, 0x3fa923a2} }, /* 0.0491 */ +/**/ {{0x00000000, 0x3fb930be} }, /* 0.0984 */ +/**/ {{0x00000000, 0x3fc95810} }, /* 0.198 */ +/**/ {{0x00000000, 0x3fda7ef9} }, /* 0.414 */ +/**/ {{0x00000000, 0x3ff00000} }, /* 1.0 */ + }; + static const number + twonm1[33] = { /* 2n-1 */ +/**/ {{0x00000000, 0x00000000} }, /* 0 */ +/**/ {{0x00000000, 0x00000000} }, /* 0 */ +/**/ {{0x00000000, 0x00000000} }, /* 0 */ +/**/ {{0x00000000, 0x00000000} }, /* 0 */ +/**/ {{0x00000000, 0x40260000} }, /* 11 */ +/**/ {{0x00000000, 0x402e0000} }, /* 15 */ +/**/ {{0x00000000, 0x40330000} }, /* 19 */ +/**/ {{0x00000000, 0x40350000} }, /* 21 */ +/**/ {{0x00000000, 0x40390000} }, /* 25 */ +/**/ {{0x00000000, 0x403d0000} }, /* 29 */ +/**/ {{0x00000000, 0x40408000} }, /* 33 */ +/**/ {{0x00000000, 0x40428000} }, /* 37 */ +/**/ {{0x00000000, 0x40448000} }, /* 41 */ +/**/ {{0x00000000, 0x40468000} }, /* 45 */ +/**/ {{0x00000000, 0x40488000} }, /* 49 */ +/**/ {{0x00000000, 0x404a8000} }, /* 53 */ +/**/ {{0x00000000, 0x404b8000} }, /* 55 */ +/**/ {{0x00000000, 0x404d8000} }, /* 59 */ +/**/ {{0x00000000, 0x404f8000} }, /* 63 */ +/**/ {{0x00000000, 0x4050c000} }, /* 67 */ +/**/ {{0x00000000, 0x4051c000} }, /* 71 */ +/**/ {{0x00000000, 0x4052c000} }, /* 75 */ +/**/ {{0x00000000, 0x4053c000} }, /* 79 */ +/**/ {{0x00000000, 0x4054c000} }, /* 83 */ +/**/ {{0x00000000, 0x40554000} }, /* 85 */ +/**/ {{0x00000000, 0x40564000} }, /* 89 */ +/**/ {{0x00000000, 0x40574000} }, /* 93 */ +/**/ {{0x00000000, 0x40584000} }, /* 97 */ +/**/ {{0x00000000, 0x40594000} }, /* 101 */ +/**/ {{0x00000000, 0x405a4000} }, /* 105 */ +/**/ {{0x00000000, 0x405b4000} }, /* 109 */ +/**/ {{0x00000000, 0x405c4000} }, /* 113 */ +/**/ {{0x00000000, 0x405d4000} }, /* 117 */ + }; + + static const number + twom[8] = { /* 2**m */ +/**/ {{0x00000000, 0x3ff00000} }, /* 1.0 */ +/**/ {{0x00000000, 0x40000000} }, /* 2.0 */ +/**/ {{0x00000000, 0x40100000} }, /* 4.0 */ +/**/ {{0x00000000, 0x40200000} }, /* 8.0 */ +/**/ {{0x00000000, 0x40300000} }, /* 16.0 */ +/**/ {{0x00000000, 0x40400000} }, /* 32.0 */ +/**/ {{0x00000000, 0x40500000} }, /* 64.0 */ +/**/ {{0x00000000, 0x40600000} }, /* 128.0 */ + }; + + static const number +/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ two = {{0x00000000, 0x40000000} }; /* 2 */ + +#endif +#endif + +#define ONE one.d +#define TWO two.d + + static const int + np[33] = { 0, 0, 0, 0, 6, 8,10,11,13,15,17,19,21,23,25,27,28, + 30,32,34,36,38,40,42,43,45,47,49,51,53,55,57,59}; + +#endif diff --git a/src/system/libroot/posix/glibc/arch/generic/mpatan2.c b/src/system/libroot/posix/glibc/arch/generic/mpatan2.c new file mode 100644 index 0000000000..8977ec9042 --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/mpatan2.c @@ -0,0 +1,70 @@ + +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001 Free Software Foundation + * + * 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ +/******************************************************************/ +/* MODULE_NAME: mpatan2.c */ +/* */ +/* FUNCTIONS:mpatan2 */ +/* */ +/* FILES NEEDED: mpa.h */ +/* mpa.c mpatan.c mpsqrt.c */ +/* */ +/* Multi-Precision Atan2(y,x) function subroutine, */ +/* for precision p >= 4. */ +/* y=0 is not permitted if x<=0. No error messages are given. */ +/* The relative error of the result is bounded by 44.84*r**(1-p) */ +/* if x <= 0, y != 0 and by 37.33*r**(1-p) if x>0. here r=2**24. */ +/* */ +/******************************************************************/ + + + +#include "mpa.h" + +void __mpsqrt(mp_no *, mp_no *, int); +void __mpatan(mp_no *, mp_no *, int); + +/* Multi-Precision Atan2(y,x) function subroutine, for p >= 4. */ +/* y=0 is not permitted if x<=0. No error messages are given. */ +void __mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) { + + static const double ZERO = 0.0, ONE = 1.0; + + mp_no mpone = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + mp_no mpt1,mpt2,mpt3; + + + if (X[0] <= ZERO) { + mpone.e = 1; mpone.d[0] = mpone.d[1] = ONE; + __dvd(x,y,&mpt1,p); __mul(&mpt1,&mpt1,&mpt2,p); + if (mpt1.d[0] != ZERO) mpt1.d[0] = ONE; + __add(&mpt2,&mpone,&mpt3,p); __mpsqrt(&mpt3,&mpt2,p); + __add(&mpt1,&mpt2,&mpt3,p); mpt3.d[0]=Y[0]; + __mpatan(&mpt3,&mpt1,p); __add(&mpt1,&mpt1,z,p); + } + else + { __dvd(y,x,&mpt1,p); + __mpatan(&mpt1,z,p); + } + + return; +} diff --git a/src/system/libroot/posix/glibc/arch/generic/mpsqrt.c b/src/system/libroot/posix/glibc/arch/generic/mpsqrt.c new file mode 100644 index 0000000000..2ad060be20 --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/mpsqrt.c @@ -0,0 +1,103 @@ + +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001 Free Software Foundation + * + * 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ +/****************************************************************************/ +/* MODULE_NAME:mpsqrt.c */ +/* */ +/* FUNCTION:mpsqrt */ +/* fastiroot */ +/* */ +/* FILES NEEDED:endian.h mpa.h mpsqrt.h */ +/* mpa.c */ +/* Multi-Precision square root function subroutine for precision p >= 4. */ +/* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */ +/* */ +/****************************************************************************/ +#include "endian.h" +#include "mpa.h" + +/****************************************************************************/ +/* Multi-Precision square root function subroutine for precision p >= 4. */ +/* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */ +/* Routine receives two pointers to Multi Precision numbers: */ +/* x (left argument) and y (next argument). Routine also receives precision */ +/* p as integer. Routine computes sqrt(*x) and stores result in *y */ +/****************************************************************************/ + +double fastiroot(double); + +void __mpsqrt(mp_no *x, mp_no *y, int p) { +#include "mpsqrt.h" + + int i,m,ex,ey; + double dx,dy; + mp_no + mphalf = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, + mp3halfs = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + mp_no mpxn,mpz,mpu,mpt1,mpt2; + + /* Prepare multi-precision 1/2 and 3/2 */ + mphalf.e =0; mphalf.d[0] =ONE; mphalf.d[1] =HALFRAD; + mp3halfs.e=1; mp3halfs.d[0]=ONE; mp3halfs.d[1]=ONE; mp3halfs.d[2]=HALFRAD; + + ex=EX; ey=EX/2; __cpy(x,&mpxn,p); mpxn.e -= (ey+ey); + __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p); + __mul(&mpxn,&mphalf,&mpz,p); + + m=mp[p]; + for (i=0; i>1; + z = ((c3*z + c2)*z + c1)*z + c0; /* 2**-7 */ + z = z*(1.5 - 0.5*y*z*z); /* 2**-14 */ + p.d = z*(1.5 - 0.5*y*z*z); /* 2**-28 */ + p.i[HIGH_HALF] -= n; + t = x*p.d; + return p.d*(1.5 - 0.5*p.d*t); +} diff --git a/src/system/libroot/posix/glibc/arch/generic/mpsqrt.h b/src/system/libroot/posix/glibc/arch/generic/mpsqrt.h new file mode 100644 index 0000000000..729d57af2c --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/mpsqrt.h @@ -0,0 +1,51 @@ +/* + * IBM Accurate Mathematical Library + * Written by International Business Machines Corp. + * Copyright (C) 2001 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ + +/******************************************************************/ +/* */ +/* MODULE_NAME:mpatan.h */ +/* */ +/* common data and variables prototype and definition */ +/******************************************************************/ + +#ifndef MPSQRT_H +#define MPSQRT_H + +#ifdef BIG_ENDI + static const number +/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ halfrad = {{0x41600000, 0x00000000} }; /* 2**23 */ + +#else +#ifdef LITTLE_ENDI + static const number +/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ halfrad = {{0x00000000, 0x41600000} }; /* 2**23 */ + +#endif +#endif + +#define ONE one.d +#define HALFRAD halfrad.d + + static const int mp[33] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4, + 4,4,4,4,4,4,4,4,4}; + +#endif diff --git a/src/system/libroot/posix/glibc/arch/generic/s_atan.c b/src/system/libroot/posix/glibc/arch/generic/s_atan.c new file mode 100644 index 0000000000..556f5b216d --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/generic/s_atan.c @@ -0,0 +1,230 @@ +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001 Free Software Foundation + * + * 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ +/************************************************************************/ +/* MODULE_NAME: atnat.c */ +/* */ +/* FUNCTIONS: uatan */ +/* atanMp */ +/* signArctan */ +/* */ +/* */ +/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat.h */ +/* mpatan.c mpatan2.c mpsqrt.c */ +/* uatan.tbl */ +/* */ +/* An ultimate atan() routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of atan(x). */ +/* */ +/* Assumption: Machine arithmetic operations are performed in */ +/* round to nearest mode of IEEE 754 standard. */ +/* */ +/************************************************************************/ + +#include "dla.h" +#include "mpa.h" +#include "MathLib.h" +#include "uatan.tbl" +#include "atnat.h" +#include "math.h" + +void __mpatan(mp_no *,mp_no *,int); /* see definition in mpatan.c */ +static double atanMp(double,const int[]); +double __signArctan(double,double); +/* An ultimate atan() routine. Given an IEEE double machine number x, */ +/* routine computes the correctly rounded (to nearest) value of atan(x). */ +double atan(double x) { + + + double cor,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,u,u2,u3, + v,vv,w,ww,y,yy,z,zz; +#if 0 + double y1,y2; +#endif + int i,ux,dx; +#if 0 + int p; +#endif + static const int pr[M]={6,8,10,32}; + number num; +#if 0 + mp_no mpt1,mpx,mpy,mpy1,mpy2,mperr; +#endif + + num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF]; + + /* x=NaN */ + if (((ux&0x7ff00000)==0x7ff00000) && (((ux&0x000fffff)|dx)!=0x00000000)) + return x+x; + + /* Regular values of x, including denormals +-0 and +-INF */ + u = (x= 1/2 */ + if ((y=t1+(yy-u3)) == t1+(yy+u3)) return __signArctan(x,y); + + DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) + t1=w-hij[i][0].d; + EADD(t1,ww,z,zz) + s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+ + z*(hij[i][14].d+z* hij[i][15].d)))); + ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) + MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) + MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2) + MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2) + MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2) + SUB2(HPI,HPI1,s2,ss2,s1,ss1,t1,t2) + if ((y=s1+(ss1-U7)) == s1+(ss1+U7)) return __signArctan(x,y); + + return atanMp(x,pr); + } + else { + if (u= E */ + if (x>0) return HPI; + else return MHPI; } + } + } + +} + + + /* Fix the sign of y and return */ +double __signArctan(double x,double y){ + + if (x=0x50800000) { /* if |x| >= 2^34 */ + if(ix>0x7f800000) + return x+x; /* NaN */ + if(hx>0) return atanhi[3]+atanlo[3]; + else return -atanhi[3]-atanlo[3]; + } if (ix < 0x3ee00000) { /* |x| < 0.4375 */ + if (ix < 0x31000000) { /* |x| < 2^-29 */ + if(huge+x>one) return x; /* raise inexact */ + } + id = -1; + } else { + x = fabsf(x); + if (ix < 0x3f980000) { /* |x| < 1.1875 */ + if (ix < 0x3f300000) { /* 7/16 <=|x|<11/16 */ + id = 0; x = ((float)2.0*x-one)/((float)2.0+x); + } else { /* 11/16<=|x|< 19/16 */ + id = 1; x = (x-one)/(x+one); + } + } else { + if (ix < 0x401c0000) { /* |x| < 2.4375 */ + id = 2; x = (x-(float)1.5)/(one+(float)1.5*x); + } else { /* 2.4375 <= |x| < 2^66 */ + id = 3; x = -(float)1.0/x; + } + }} + /* end of argument reduction */ + z = x*x; + w = z*z; + /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ + s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10]))))); + s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9])))); + if (id<0) return x - x*(s1+s2); + else { + z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); + return (hx<0)? -z:z; + } +} +weak_alias (__atanf, atanf) diff --git a/src/system/libroot/posix/glibc/arch/ppc/Jamfile b/src/system/libroot/posix/glibc/arch/ppc/Jamfile index 33b0e194e0..a20d50bf5d 100644 --- a/src/system/libroot/posix/glibc/arch/ppc/Jamfile +++ b/src/system/libroot/posix/glibc/arch/ppc/Jamfile @@ -19,10 +19,12 @@ SubDirCcFlags -D_GNU_SOURCE -D_IEEE_LIBM ; # Note: There is no *l() support yet. Our compiler says sizeof(long double) = 8, # while there are only 96 and 128 bit implementation in glibc. local genericSources = - cmp.c dbl2mpn.c divrem.c mpn2dbl.c mul.c mul_n.c + cmp.c dbl2mpn.c divrem.c mpa.c mpatan.c mpatan2.c mpn2dbl.c mpsqrt.c + mul.c mul_n.c e_acos.c e_acosf.c # e_acosl.c e_atan2.c e_atan2f.c # e_atan2l.c e_fmod.c e_fmodf.c # e_fmodl.c + s_atan.c s_atanf.c # s_atanl.c s_ceil.c s_ceilf.c # s_ceill.c s_finite.c s_finitef.c # s_finitel.c s_floor.c s_floorf.c # s_floorl.c @@ -30,12 +32,15 @@ local genericSources = s_logb.c s_logbf.c # s_logbl.c s_modf.c s_modff.c s_modfl.c s_scalbn.c s_scalbnf.c # s_scalbnl.c - w_atan2.c w_atan2f.c w_atan2l.c - w_fmod.c w_fmodl.c w_fmodf.c + w_atan2.c w_atan2f.c # w_atan2l.c + w_fmod.c w_fmodf.c # w_fmodl.c ; MergeObject posix_gnu_arch_$(TARGET_ARCH).o : add_n.S addmul_1.S + e_sqrtf.c + fraiseexcpt.c + ftestexcept.c ldbl2mpn.c mul_1.S lshift.S @@ -54,6 +59,9 @@ MergeObject posix_gnu_arch_$(TARGET_ARCH).o : # s_rintl.c sub_n.S submul_1.S + t_sqrt.c + w_sqrt.c + w_sqrtf.c $(genericSources) ; diff --git a/src/system/libroot/posix/glibc/arch/ppc/e_sqrtf.c b/src/system/libroot/posix/glibc/arch/ppc/e_sqrtf.c new file mode 100644 index 0000000000..01c76d6757 --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/ppc/e_sqrtf.c @@ -0,0 +1 @@ +/* __ieee754_sqrtf is in w_sqrtf.c */ diff --git a/src/system/libroot/posix/glibc/arch/ppc/fraiseexcpt.c b/src/system/libroot/posix/glibc/arch/ppc/fraiseexcpt.c new file mode 100644 index 0000000000..dbe36c3d5a --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/ppc/fraiseexcpt.c @@ -0,0 +1,66 @@ +/* Raise given exceptions. + Copyright (C) 1997,99,2000,01,02 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library 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. + + The GNU C Library 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 the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include +#include + +#undef feraiseexcept +int +__feraiseexcept (int excepts) +{ + fenv_union_t u; + + /* Raise exceptions represented by EXCEPTS. It is the responsibility of + the OS to ensure that if multiple exceptions occur they are fed back + to this process in the proper way; this can happen in hardware, + anyway (in particular, inexact with overflow or underflow). */ + + /* Get the current state. */ + u.fenv = fegetenv_register (); + + /* Add the exceptions */ + u.l[1] = (u.l[1] + | (excepts & FPSCR_STICKY_BITS) + /* Turn FE_INVALID into FE_INVALID_SOFTWARE. */ + | (excepts >> ((31 - FPSCR_VX) - (31 - FPSCR_VXSOFT)) + & FE_INVALID_SOFTWARE)); + + /* Store the new status word (along with the rest of the environment), + triggering any appropriate exceptions. */ + fesetenv_register (u.fenv); + + if ((excepts & FE_INVALID) + /* For some reason, some PowerPC chips (the 601, in particular) + don't have FE_INVALID_SOFTWARE implemented. Detect this + case and raise FE_INVALID_SNAN instead. */ + && !fetestexcept (FE_INVALID)) + set_fpscr_bit (FPSCR_VXSNAN); + + /* Success. */ + return 0; +} + +#include +#if SHLIB_COMPAT (libm, GLIBC_2_1, GLIBC_2_2) +strong_alias (__feraiseexcept, __old_feraiseexcept) +compat_symbol (libm, BP_SYM (__old_feraiseexcept), BP_SYM (feraiseexcept), GLIBC_2_1); +#endif + +libm_hidden_ver (__feraiseexcept, feraiseexcept) +versioned_symbol (libm, BP_SYM (__feraiseexcept), BP_SYM (feraiseexcept), GLIBC_2_2); diff --git a/src/system/libroot/posix/glibc/arch/ppc/ftestexcept.c b/src/system/libroot/posix/glibc/arch/ppc/ftestexcept.c new file mode 100644 index 0000000000..64406c41dc --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/ppc/ftestexcept.c @@ -0,0 +1,33 @@ +/* Test exception in current environment. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library 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. + + The GNU C Library 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 the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include + +int +fetestexcept (int excepts) +{ + fenv_union_t u; + + /* Get the current state. */ + u.fenv = fegetenv_register (); + + /* The FE_INVALID bit is dealt with correctly by the hardware, so we can + just: */ + return u.l[1] & excepts; +} diff --git a/src/system/libroot/posix/glibc/arch/ppc/t_sqrt.c b/src/system/libroot/posix/glibc/arch/ppc/t_sqrt.c new file mode 100644 index 0000000000..c49380c0fd --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/ppc/t_sqrt.c @@ -0,0 +1,144 @@ +const float __t_sqrt[1024] = { +0.7078,0.7064, 0.7092,0.7050, 0.7106,0.7037, 0.7119,0.7023, 0.7133,0.7010, +0.7147,0.6996, 0.7160,0.6983, 0.7174,0.6970, 0.7187,0.6957, 0.7201,0.6943, +0.7215,0.6930, 0.7228,0.6917, 0.7242,0.6905, 0.7255,0.6892, 0.7269,0.6879, +0.7282,0.6866, 0.7295,0.6854, 0.7309,0.6841, 0.7322,0.6829, 0.7335,0.6816, +0.7349,0.6804, 0.7362,0.6792, 0.7375,0.6779, 0.7388,0.6767, 0.7402,0.6755, +0.7415,0.6743, 0.7428,0.6731, 0.7441,0.6719, 0.7454,0.6708, 0.7467,0.6696, +0.7480,0.6684, 0.7493,0.6672, 0.7507,0.6661, 0.7520,0.6649, 0.7532,0.6638, +0.7545,0.6627, 0.7558,0.6615, 0.7571,0.6604, 0.7584,0.6593, 0.7597,0.6582, +0.7610,0.6570, 0.7623,0.6559, 0.7635,0.6548, 0.7648,0.6537, 0.7661,0.6527, +0.7674,0.6516, 0.7686,0.6505, 0.7699,0.6494, 0.7712,0.6484, 0.7725,0.6473, +0.7737,0.6462, 0.7750,0.6452, 0.7762,0.6441, 0.7775,0.6431, 0.7787,0.6421, +0.7800,0.6410, 0.7812,0.6400, 0.7825,0.6390, 0.7837,0.6380, 0.7850,0.6370, +0.7862,0.6359, 0.7875,0.6349, 0.7887,0.6339, 0.7900,0.6330, 0.7912,0.6320, +0.7924,0.6310, 0.7937,0.6300, 0.7949,0.6290, 0.7961,0.6281, 0.7973,0.6271, +0.7986,0.6261, 0.7998,0.6252, 0.8010,0.6242, 0.8022,0.6233, 0.8034,0.6223, +0.8046,0.6214, 0.8059,0.6205, 0.8071,0.6195, 0.8083,0.6186, 0.8095,0.6177, +0.8107,0.6168, 0.8119,0.6158, 0.8131,0.6149, 0.8143,0.6140, 0.8155,0.6131, +0.8167,0.6122, 0.8179,0.6113, 0.8191,0.6104, 0.8203,0.6096, 0.8215,0.6087, +0.8227,0.6078, 0.8238,0.6069, 0.8250,0.6060, 0.8262,0.6052, 0.8274,0.6043, +0.8286,0.6035, 0.8297,0.6026, 0.8309,0.6017, 0.8321,0.6009, 0.8333,0.6000, +0.8344,0.5992, 0.8356,0.5984, 0.8368,0.5975, 0.8379,0.5967, 0.8391,0.5959, +0.8403,0.5950, 0.8414,0.5942, 0.8426,0.5934, 0.8437,0.5926, 0.8449,0.5918, +0.8461,0.5910, 0.8472,0.5902, 0.8484,0.5894, 0.8495,0.5886, 0.8507,0.5878, +0.8518,0.5870, 0.8530,0.5862, 0.8541,0.5854, 0.8552,0.5846, 0.8564,0.5838, +0.8575,0.5831, 0.8587,0.5823, 0.8598,0.5815, 0.8609,0.5808, 0.8621,0.5800, +0.8632,0.5792, 0.8643,0.5785, 0.8655,0.5777, 0.8666,0.5770, 0.8677,0.5762, +0.8688,0.5755, 0.8700,0.5747, 0.8711,0.5740, 0.8722,0.5733, 0.8733,0.5725, +0.8744,0.5718, 0.8756,0.5711, 0.8767,0.5703, 0.8778,0.5696, 0.8789,0.5689, +0.8800,0.5682, 0.8811,0.5675, 0.8822,0.5667, 0.8833,0.5660, 0.8844,0.5653, +0.8855,0.5646, 0.8866,0.5639, 0.8877,0.5632, 0.8888,0.5625, 0.8899,0.5618, +0.8910,0.5611, 0.8921,0.5605, 0.8932,0.5598, 0.8943,0.5591, 0.8954,0.5584, +0.8965,0.5577, 0.8976,0.5570, 0.8987,0.5564, 0.8998,0.5557, 0.9008,0.5550, +0.9019,0.5544, 0.9030,0.5537, 0.9041,0.5530, 0.9052,0.5524, 0.9062,0.5517, +0.9073,0.5511, 0.9084,0.5504, 0.9095,0.5498, 0.9105,0.5491, 0.9116,0.5485, +0.9127,0.5478, 0.9138,0.5472, 0.9148,0.5465, 0.9159,0.5459, 0.9170,0.5453, +0.9180,0.5446, 0.9191,0.5440, 0.9202,0.5434, 0.9212,0.5428, 0.9223,0.5421, +0.9233,0.5415, 0.9244,0.5409, 0.9254,0.5403, 0.9265,0.5397, 0.9276,0.5391, +0.9286,0.5384, 0.9297,0.5378, 0.9307,0.5372, 0.9318,0.5366, 0.9328,0.5360, +0.9338,0.5354, 0.9349,0.5348, 0.9359,0.5342, 0.9370,0.5336, 0.9380,0.5330, +0.9391,0.5324, 0.9401,0.5319, 0.9411,0.5313, 0.9422,0.5307, 0.9432,0.5301, +0.9442,0.5295, 0.9453,0.5289, 0.9463,0.5284, 0.9473,0.5278, 0.9484,0.5272, +0.9494,0.5266, 0.9504,0.5261, 0.9515,0.5255, 0.9525,0.5249, 0.9535,0.5244, +0.9545,0.5238, 0.9556,0.5233, 0.9566,0.5227, 0.9576,0.5221, 0.9586,0.5216, +0.9596,0.5210, 0.9607,0.5205, 0.9617,0.5199, 0.9627,0.5194, 0.9637,0.5188, +0.9647,0.5183, 0.9657,0.5177, 0.9667,0.5172, 0.9677,0.5167, 0.9687,0.5161, +0.9698,0.5156, 0.9708,0.5151, 0.9718,0.5145, 0.9728,0.5140, 0.9738,0.5135, +0.9748,0.5129, 0.9758,0.5124, 0.9768,0.5119, 0.9778,0.5114, 0.9788,0.5108, +0.9798,0.5103, 0.9808,0.5098, 0.9818,0.5093, 0.9828,0.5088, 0.9838,0.5083, +0.9847,0.5077, 0.9857,0.5072, 0.9867,0.5067, 0.9877,0.5062, 0.9887,0.5057, +0.9897,0.5052, 0.9907,0.5047, 0.9917,0.5042, 0.9926,0.5037, 0.9936,0.5032, +0.9946,0.5027, 0.9956,0.5022, 0.9966,0.5017, 0.9976,0.5012, 0.9985,0.5007, +0.9995,0.5002, +1.0010,0.4995, 1.0029,0.4985, 1.0049,0.4976, 1.0068,0.4966, 1.0088,0.4957, +1.0107,0.4947, 1.0126,0.4938, 1.0145,0.4928, 1.0165,0.4919, 1.0184,0.4910, +1.0203,0.4901, 1.0222,0.4891, 1.0241,0.4882, 1.0260,0.4873, 1.0279,0.4864, +1.0298,0.4855, 1.0317,0.4846, 1.0336,0.4837, 1.0355,0.4829, 1.0374,0.4820, +1.0393,0.4811, 1.0411,0.4802, 1.0430,0.4794, 1.0449,0.4785, 1.0468,0.4777, +1.0486,0.4768, 1.0505,0.4760, 1.0523,0.4751, 1.0542,0.4743, 1.0560,0.4735, +1.0579,0.4726, 1.0597,0.4718, 1.0616,0.4710, 1.0634,0.4702, 1.0653,0.4694, +1.0671,0.4686, 1.0689,0.4678, 1.0707,0.4670, 1.0726,0.4662, 1.0744,0.4654, +1.0762,0.4646, 1.0780,0.4638, 1.0798,0.4630, 1.0816,0.4623, 1.0834,0.4615, +1.0852,0.4607, 1.0870,0.4600, 1.0888,0.4592, 1.0906,0.4585, 1.0924,0.4577, +1.0942,0.4570, 1.0960,0.4562, 1.0978,0.4555, 1.0995,0.4547, 1.1013,0.4540, +1.1031,0.4533, 1.1049,0.4525, 1.1066,0.4518, 1.1084,0.4511, 1.1101,0.4504, +1.1119,0.4497, 1.1137,0.4490, 1.1154,0.4483, 1.1172,0.4476, 1.1189,0.4469, +1.1207,0.4462, 1.1224,0.4455, 1.1241,0.4448, 1.1259,0.4441, 1.1276,0.4434, +1.1293,0.4427, 1.1311,0.4421, 1.1328,0.4414, 1.1345,0.4407, 1.1362,0.4401, +1.1379,0.4394, 1.1397,0.4387, 1.1414,0.4381, 1.1431,0.4374, 1.1448,0.4368, +1.1465,0.4361, 1.1482,0.4355, 1.1499,0.4348, 1.1516,0.4342, 1.1533,0.4335, +1.1550,0.4329, 1.1567,0.4323, 1.1584,0.4316, 1.1600,0.4310, 1.1617,0.4304, +1.1634,0.4298, 1.1651,0.4292, 1.1668,0.4285, 1.1684,0.4279, 1.1701,0.4273, +1.1718,0.4267, 1.1734,0.4261, 1.1751,0.4255, 1.1768,0.4249, 1.1784,0.4243, +1.1801,0.4237, 1.1817,0.4231, 1.1834,0.4225, 1.1850,0.4219, 1.1867,0.4213, +1.1883,0.4208, 1.1900,0.4202, 1.1916,0.4196, 1.1932,0.4190, 1.1949,0.4185, +1.1965,0.4179, 1.1981,0.4173, 1.1998,0.4167, 1.2014,0.4162, 1.2030,0.4156, +1.2046,0.4151, 1.2063,0.4145, 1.2079,0.4139, 1.2095,0.4134, 1.2111,0.4128, +1.2127,0.4123, 1.2143,0.4117, 1.2159,0.4112, 1.2175,0.4107, 1.2192,0.4101, +1.2208,0.4096, 1.2224,0.4090, 1.2239,0.4085, 1.2255,0.4080, 1.2271,0.4075, +1.2287,0.4069, 1.2303,0.4064, 1.2319,0.4059, 1.2335,0.4054, 1.2351,0.4048, +1.2366,0.4043, 1.2382,0.4038, 1.2398,0.4033, 1.2414,0.4028, 1.2429,0.4023, +1.2445,0.4018, 1.2461,0.4013, 1.2477,0.4008, 1.2492,0.4003, 1.2508,0.3998, +1.2523,0.3993, 1.2539,0.3988, 1.2555,0.3983, 1.2570,0.3978, 1.2586,0.3973, +1.2601,0.3968, 1.2617,0.3963, 1.2632,0.3958, 1.2648,0.3953, 1.2663,0.3949, +1.2678,0.3944, 1.2694,0.3939, 1.2709,0.3934, 1.2725,0.3929, 1.2740,0.3925, +1.2755,0.3920, 1.2771,0.3915, 1.2786,0.3911, 1.2801,0.3906, 1.2816,0.3901, +1.2832,0.3897, 1.2847,0.3892, 1.2862,0.3887, 1.2877,0.3883, 1.2892,0.3878, +1.2907,0.3874, 1.2923,0.3869, 1.2938,0.3865, 1.2953,0.3860, 1.2968,0.3856, +1.2983,0.3851, 1.2998,0.3847, 1.3013,0.3842, 1.3028,0.3838, 1.3043,0.3834, +1.3058,0.3829, 1.3073,0.3825, 1.3088,0.3820, 1.3103,0.3816, 1.3118,0.3812, +1.3132,0.3807, 1.3147,0.3803, 1.3162,0.3799, 1.3177,0.3794, 1.3192,0.3790, +1.3207,0.3786, 1.3221,0.3782, 1.3236,0.3778, 1.3251,0.3773, 1.3266,0.3769, +1.3280,0.3765, 1.3295,0.3761, 1.3310,0.3757, 1.3324,0.3753, 1.3339,0.3748, +1.3354,0.3744, 1.3368,0.3740, 1.3383,0.3736, 1.3397,0.3732, 1.3412,0.3728, +1.3427,0.3724, 1.3441,0.3720, 1.3456,0.3716, 1.3470,0.3712, 1.3485,0.3708, +1.3499,0.3704, 1.3514,0.3700, 1.3528,0.3696, 1.3542,0.3692, 1.3557,0.3688, +1.3571,0.3684, 1.3586,0.3680, 1.3600,0.3676, 1.3614,0.3673, 1.3629,0.3669, +1.3643,0.3665, 1.3657,0.3661, 1.3672,0.3657, 1.3686,0.3653, 1.3700,0.3650, +1.3714,0.3646, 1.3729,0.3642, 1.3743,0.3638, 1.3757,0.3634, 1.3771,0.3631, +1.3785,0.3627, 1.3800,0.3623, 1.3814,0.3620, 1.3828,0.3616, 1.3842,0.3612, +1.3856,0.3609, 1.3870,0.3605, 1.3884,0.3601, 1.3898,0.3598, 1.3912,0.3594, +1.3926,0.3590, 1.3940,0.3587, 1.3954,0.3583, 1.3968,0.3580, 1.3982,0.3576, +1.3996,0.3572, 1.4010,0.3569, 1.4024,0.3565, 1.4038,0.3562, 1.4052,0.3558, +1.4066,0.3555, 1.4080,0.3551, 1.4094,0.3548, 1.4108,0.3544, 1.4121,0.3541, +1.4135,0.3537 +}; + + +/* Generated by: */ +#if 0 +#include +#include +#include + +int +main(int argc, char **argv) +{ + int i, j; + + printf ("const float __t_sqrt[1024] = {"); + for (i = 0; i < 2; i++) + { + putchar('\n'); + for (j = 0; j < 256; j++) + { + double mval = j/512.0 + 0.5; + double eval = i==0 ? 1.0 : 2.0; + double ls = sqrt(mval*eval); + double hs = sqrt((mval+1/512.0)*eval); + double as = (ls+hs)*0.5; + double lx = 1/(2.0*ls); + double hx = 1/(2.0*hs); + double ax = (lx+hx)*0.5; + + printf("%.4f,%.4f%s",as,ax, + i*j==255 ? "\n" : j % 5 == 4 ? ",\n" : ", "); + assert((hs-ls)/as < 1/256.0); + assert((hx-lx)/ax < 1/256.0); + } + } + printf ("};\n"); + return 0; +} +#endif /* 0 */ diff --git a/src/system/libroot/posix/glibc/arch/ppc/w_sqrt.c b/src/system/libroot/posix/glibc/arch/ppc/w_sqrt.c new file mode 100644 index 0000000000..8141f361a7 --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/ppc/w_sqrt.c @@ -0,0 +1,146 @@ +/* Single-precision floating point square root. + Copyright (C) 1997, 2002 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library 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. + + The GNU C Library 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 the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include +#include +#include +#include + +static const double almost_half = 0.5000000000000001; /* 0.5 + 2^-53 */ +static const uint32_t a_nan = 0x7fc00000; +static const uint32_t a_inf = 0x7f800000; +static const float two108 = 3.245185536584267269e+32; +static const float twom54 = 5.551115123125782702e-17; +extern const float __t_sqrt[1024]; + +/* The method is based on a description in + Computation of elementary functions on the IBM RISC System/6000 processor, + P. W. Markstein, IBM J. Res. Develop, 34(1) 1990. + Basically, it consists of two interleaved Newton-Rhapson approximations, + one to find the actual square root, and one to find its reciprocal + without the expense of a division operation. The tricky bit here + is the use of the POWER/PowerPC multiply-add operation to get the + required accuracy with high speed. + + The argument reduction works by a combination of table lookup to + obtain the initial guesses, and some careful modification of the + generated guesses (which mostly runs on the integer unit, while the + Newton-Rhapson is running on the FPU). */ +double +__sqrt(double x) +{ + const float inf = *(const float *)&a_inf; + /* x = f_wash(x); *//* This ensures only one exception for SNaN. */ + if (x > 0) + { + if (x != inf) + { + /* Variables named starting with 's' exist in the + argument-reduced space, so that 2 > sx >= 0.5, + 1.41... > sg >= 0.70.., 0.70.. >= sy > 0.35... . + Variables named ending with 'i' are integer versions of + floating-point values. */ + double sx; /* The value of which we're trying to find the + square root. */ + double sg,g; /* Guess of the square root of x. */ + double sd,d; /* Difference between the square of the guess and x. */ + double sy; /* Estimate of 1/2g (overestimated by 1ulp). */ + double sy2; /* 2*sy */ + double e; /* Difference between y*g and 1/2 (se = e * fsy). */ + double shx; /* == sx * fsg */ + double fsg; /* sg*fsg == g. */ + fenv_t fe; /* Saved floating-point environment (stores rounding + mode and whether the inexact exception is + enabled). */ + uint32_t xi0, xi1, sxi, fsgi; + const float *t_sqrt; + + fe = fegetenv_register(); + EXTRACT_WORDS (xi0,xi1,x); + relax_fenv_state(); + sxi = (xi0 & 0x3fffffff) | 0x3fe00000; + INSERT_WORDS (sx, sxi, xi1); + t_sqrt = __t_sqrt + (xi0 >> (52-32-8-1) & 0x3fe); + sg = t_sqrt[0]; + sy = t_sqrt[1]; + + /* Here we have three Newton-Rhapson iterations each of a + division and a square root and the remainder of the + argument reduction, all interleaved. */ + sd = -(sg*sg - sx); + fsgi = (xi0 + 0x40000000) >> 1 & 0x7ff00000; + sy2 = sy + sy; + sg = sy*sd + sg; /* 16-bit approximation to sqrt(sx). */ + INSERT_WORDS (fsg, fsgi, 0); + e = -(sy*sg - almost_half); + sd = -(sg*sg - sx); + if ((xi0 & 0x7ff00000) == 0) + goto denorm; + sy = sy + e*sy2; + sg = sg + sy*sd; /* 32-bit approximation to sqrt(sx). */ + sy2 = sy + sy; + e = -(sy*sg - almost_half); + sd = -(sg*sg - sx); + sy = sy + e*sy2; + shx = sx * fsg; + sg = sg + sy*sd; /* 64-bit approximation to sqrt(sx), + but perhaps rounded incorrectly. */ + sy2 = sy + sy; + g = sg * fsg; + e = -(sy*sg - almost_half); + d = -(g*sg - shx); + sy = sy + e*sy2; + fesetenv_register (fe); + return g + sy*d; + denorm: + /* For denormalised numbers, we normalise, calculate the + square root, and return an adjusted result. */ + fesetenv_register (fe); + return __sqrt(x * two108) * twom54; + } + } + else if (x < 0) + { +#ifdef FE_INVALID_SQRT + feraiseexcept (FE_INVALID_SQRT); + /* For some reason, some PowerPC processors don't implement + FE_INVALID_SQRT. I guess no-one ever thought they'd be + used for square roots... :-) */ + if (!fetestexcept (FE_INVALID)) +#endif + feraiseexcept (FE_INVALID); +#ifndef _IEEE_LIBM + if (_LIB_VERSION != _IEEE_) + x = __kernel_standard(x,x,26); + else +#endif + x = *(const float*)&a_nan; + } + return f_wash(x); +} + +weak_alias (__sqrt, sqrt) +/* Strictly, this is wrong, but the only places where _ieee754_sqrt is + used will not pass in a negative result. */ +strong_alias(__sqrt,__ieee754_sqrt) + +#ifdef NO_LONG_DOUBLE +weak_alias (__sqrt, __sqrtl) +weak_alias (__sqrt, sqrtl) +#endif diff --git a/src/system/libroot/posix/glibc/arch/ppc/w_sqrtf.c b/src/system/libroot/posix/glibc/arch/ppc/w_sqrtf.c new file mode 100644 index 0000000000..20fbebc06e --- /dev/null +++ b/src/system/libroot/posix/glibc/arch/ppc/w_sqrtf.c @@ -0,0 +1,136 @@ +/* Single-precision floating point square root. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library 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. + + The GNU C Library 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 the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include +#include +#include +#include + +static const float almost_half = 0.50000006; /* 0.5 + 2^-24 */ +static const uint32_t a_nan = 0x7fc00000; +static const uint32_t a_inf = 0x7f800000; +static const float two48 = 281474976710656.0; +static const float twom24 = 5.9604644775390625e-8; +extern const float __t_sqrt[1024]; + +/* The method is based on a description in + Computation of elementary functions on the IBM RISC System/6000 processor, + P. W. Markstein, IBM J. Res. Develop, 34(1) 1990. + Basically, it consists of two interleaved Newton-Rhapson approximations, + one to find the actual square root, and one to find its reciprocal + without the expense of a division operation. The tricky bit here + is the use of the POWER/PowerPC multiply-add operation to get the + required accuracy with high speed. + + The argument reduction works by a combination of table lookup to + obtain the initial guesses, and some careful modification of the + generated guesses (which mostly runs on the integer unit, while the + Newton-Rhapson is running on the FPU). */ +float +__sqrtf(float x) +{ + const float inf = *(const float *)&a_inf; + /* x = f_washf(x); *//* This ensures only one exception for SNaN. */ + if (x > 0) + { + if (x != inf) + { + /* Variables named starting with 's' exist in the + argument-reduced space, so that 2 > sx >= 0.5, + 1.41... > sg >= 0.70.., 0.70.. >= sy > 0.35... . + Variables named ending with 'i' are integer versions of + floating-point values. */ + float sx; /* The value of which we're trying to find the square + root. */ + float sg,g; /* Guess of the square root of x. */ + float sd,d; /* Difference between the square of the guess and x. */ + float sy; /* Estimate of 1/2g (overestimated by 1ulp). */ + float sy2; /* 2*sy */ + float e; /* Difference between y*g and 1/2 (note that e==se). */ + float shx; /* == sx * fsg */ + float fsg; /* sg*fsg == g. */ + fenv_t fe; /* Saved floating-point environment (stores rounding + mode and whether the inexact exception is + enabled). */ + uint32_t xi, sxi, fsgi; + const float *t_sqrt; + + GET_FLOAT_WORD (xi, x); + fe = fegetenv_register (); + relax_fenv_state (); + sxi = (xi & 0x3fffffff) | 0x3f000000; + SET_FLOAT_WORD (sx, sxi); + t_sqrt = __t_sqrt + (xi >> (23-8-1) & 0x3fe); + sg = t_sqrt[0]; + sy = t_sqrt[1]; + + /* Here we have three Newton-Rhapson iterations each of a + division and a square root and the remainder of the + argument reduction, all interleaved. */ + sd = -(sg*sg - sx); + fsgi = (xi + 0x40000000) >> 1 & 0x7f800000; + sy2 = sy + sy; + sg = sy*sd + sg; /* 16-bit approximation to sqrt(sx). */ + e = -(sy*sg - almost_half); + SET_FLOAT_WORD (fsg, fsgi); + sd = -(sg*sg - sx); + sy = sy + e*sy2; + if ((xi & 0x7f800000) == 0) + goto denorm; + shx = sx * fsg; + sg = sg + sy*sd; /* 32-bit approximation to sqrt(sx), + but perhaps rounded incorrectly. */ + sy2 = sy + sy; + g = sg * fsg; + e = -(sy*sg - almost_half); + d = -(g*sg - shx); + sy = sy + e*sy2; + fesetenv_register (fe); + return g + sy*d; + denorm: + /* For denormalised numbers, we normalise, calculate the + square root, and return an adjusted result. */ + fesetenv_register (fe); + return __sqrtf(x * two48) * twom24; + } + } + else if (x < 0) + { +#ifdef FE_INVALID_SQRT + feraiseexcept (FE_INVALID_SQRT); + /* For some reason, some PowerPC processors don't implement + FE_INVALID_SQRT. I guess no-one ever thought they'd be + used for square roots... :-) */ + if (!fetestexcept (FE_INVALID)) +#endif + feraiseexcept (FE_INVALID); +#ifndef _IEEE_LIBM + if (_LIB_VERSION != _IEEE_) + x = __kernel_standard(x,x,126); + else +#endif + x = *(const float*)&a_nan; + } + return f_washf(x); +} + +weak_alias (__sqrtf, sqrtf) +/* Strictly, this is wrong, but the only places where _ieee754_sqrt is + used will not pass in a negative result. */ +strong_alias(__sqrtf,__ieee754_sqrtf) diff --git a/src/system/libroot/posix/glibc/include/arch/ppc/bits/fenv.h b/src/system/libroot/posix/glibc/include/arch/ppc/bits/fenv.h new file mode 100644 index 0000000000..8509b4b0c3 --- /dev/null +++ b/src/system/libroot/posix/glibc/include/arch/ppc/bits/fenv.h @@ -0,0 +1,145 @@ +/* Copyright (C) 1997, 1998, 1999 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library 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. + + The GNU C Library 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 the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#ifndef _FENV_H +# error "Never use directly; include instead." +#endif + + +/* Define bits representing the exception. We use the bit positions of + the appropriate bits in the FPSCR... */ +enum + { + FE_INEXACT = 1 << (31 - 6), +#define FE_INEXACT FE_INEXACT + FE_DIVBYZERO = 1 << (31 - 5), +#define FE_DIVBYZERO FE_DIVBYZERO + FE_UNDERFLOW = 1 << (31 - 4), +#define FE_UNDERFLOW FE_UNDERFLOW + FE_OVERFLOW = 1 << (31 - 3), +#define FE_OVERFLOW FE_OVERFLOW + + /* ... except for FE_INVALID, for which we use bit 31. FE_INVALID + actually corresponds to bits 7 through 12 and 21 through 23 + in the FPSCR, but we can't use that because the current draft + says that it must be a power of 2. Instead we use bit 2 which + is the summary bit for all the FE_INVALID exceptions, which + kind of makes sense. */ + FE_INVALID = 1 << (31 - 2), +#define FE_INVALID FE_INVALID + +#ifdef __USE_GNU + /* Breakdown of the FE_INVALID bits. Setting FE_INVALID on an + input to a routine is equivalent to setting all of these bits; + FE_INVALID will be set on output from a routine iff one of + these bits is set. Note, though, that you can't disable or + enable these exceptions individually. */ + + /* Operation with SNaN. */ + FE_INVALID_SNAN = 1 << (31 - 7), +# define FE_INVALID_SNAN FE_INVALID_SNAN + + /* Inf - Inf */ + FE_INVALID_ISI = 1 << (31 - 8), +# define FE_INVALID_ISI FE_INVALID_ISI + + /* Inf / Inf */ + FE_INVALID_IDI = 1 << (31 - 9), +# define FE_INVALID_IDI FE_INVALID_IDI + + /* 0 / 0 */ + FE_INVALID_ZDZ = 1 << (31 - 10), +# define FE_INVALID_ZDZ FE_INVALID_ZDZ + + /* Inf * 0 */ + FE_INVALID_IMZ = 1 << (31 - 11), +# define FE_INVALID_IMZ FE_INVALID_IMZ + + /* Comparison with NaN or SNaN. */ + FE_INVALID_COMPARE = 1 << (31 - 12), +# define FE_INVALID_COMPARE FE_INVALID_COMPARE + + /* Invalid operation flag for software (not set by hardware). */ + /* Note that some chips don't have this implemented, presumably + because no-one expected anyone to write software for them %-). */ + FE_INVALID_SOFTWARE = 1 << (31 - 21), +# define FE_INVALID_SOFTWARE FE_INVALID_SOFTWARE + + /* Square root of negative number (including -Inf). */ + /* Note that some chips don't have this implemented. */ + FE_INVALID_SQRT = 1 << (31 - 22), +# define FE_INVALID_SQRT FE_INVALID_SQRT + + /* Conversion-to-integer of a NaN or a number too large or too small. */ + FE_INVALID_INTEGER_CONVERSION = 1 << (31 - 23) +# define FE_INVALID_INTEGER_CONVERSION FE_INVALID_INTEGER_CONVERSION + +# define FE_ALL_INVALID \ + (FE_INVALID_SNAN | FE_INVALID_ISI | FE_INVALID_IDI | FE_INVALID_ZDZ \ + | FE_INVALID_IMZ | FE_INVALID_COMPARE | FE_INVALID_SOFTWARE \ + | FE_INVALID_SQRT | FE_INVALID_INTEGER_CONVERSION) +#endif + }; + +#define FE_ALL_EXCEPT \ + (FE_INEXACT | FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID) + +/* PowerPC chips support all of the four defined rounding modes. We + use the bit pattern in the FPSCR as the values for the + appropriate macros. */ +enum + { + FE_TONEAREST = 0, +#define FE_TONEAREST FE_TONEAREST + FE_TOWARDZERO = 1, +#define FE_TOWARDZERO FE_TOWARDZERO + FE_UPWARD = 2, +#define FE_UPWARD FE_UPWARD + FE_DOWNWARD = 3 +#define FE_DOWNWARD FE_DOWNWARD + }; + +/* Type representing exception flags. */ +typedef unsigned int fexcept_t; + +/* Type representing floating-point environment. We leave it as 'double' + for efficiency reasons (rather than writing it to a 32-bit integer). */ +typedef double fenv_t; + +/* If the default argument is used we use this value. */ +extern const fenv_t __fe_dfl_env; +#define FE_DFL_ENV (&__fe_dfl_env) + +#ifdef __USE_GNU +/* Floating-point environment where all exceptions are enabled. Note that + this is not sufficient to give you SIGFPE. */ +extern const fenv_t __fe_enabled_env; +# define FE_ENABLED_ENV (&__fe_enabled_env) + +/* Floating-point environment with (processor-dependent) non-IEEE floating + point. */ +extern const fenv_t __fe_nonieee_env; +# define FE_NONIEEE_ENV (&__fe_nonieee_env) + +/* Floating-point environment with all exceptions enabled. Note that + just evaluating this value will set the processor into 'FPU + exceptions imprecise recoverable' mode, which may cause a significant + performance penalty (but have no other visible effect). */ +extern const fenv_t *__fe_nomask_env (void); +# define FE_NOMASK_ENV (__fe_nomask_env ()) +#endif diff --git a/src/system/libroot/posix/glibc/include/arch/ppc/bits/fenvinline.h b/src/system/libroot/posix/glibc/include/arch/ppc/bits/fenvinline.h new file mode 100644 index 0000000000..552c8c9db7 --- /dev/null +++ b/src/system/libroot/posix/glibc/include/arch/ppc/bits/fenvinline.h @@ -0,0 +1,59 @@ +/* Inline floating-point environment handling functions for powerpc. + Copyright (C) 1995, 1996, 1997, 1998, 1999 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library 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. + + The GNU C Library 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 the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#if defined __GNUC__ && !defined _SOFT_FLOAT && !defined __NO_MATH_INLINES + +/* Inline definition for fegetround. */ +# define fegetround() \ + (__extension__ ({ int __fegetround_result; \ + __asm__ ("mcrfs 7,7 ; mfcr %0" \ + : "=r"(__fegetround_result) : : "cr7"); \ + __fegetround_result & 3; })) + +/* The weird 'i#*X' constraints on the following suppress a gcc + warning when __excepts is not a constant. Otherwise, they mean the + same as just plain 'i'. */ + +/* Inline definition for feraiseexcept. */ +# define feraiseexcept(__excepts) \ + ((__builtin_constant_p (__excepts) \ + && ((__excepts) & ((__excepts)-1)) == 0 \ + && (__excepts) != FE_INVALID) \ + ? ((__excepts) != 0 \ + ? (__extension__ ({ __asm__ __volatile__ \ + ("mtfsb1 %s0" \ + : : "i#*X"(__builtin_ffs (__excepts))); \ + 0; })) \ + : 0) \ + : (feraiseexcept) (__excepts)) + +/* Inline definition for feclearexcept. */ +# define feclearexcept(__excepts) \ + ((__builtin_constant_p (__excepts) \ + && ((__excepts) & ((__excepts)-1)) == 0 \ + && (__excepts) != FE_INVALID) \ + ? ((__excepts) != 0 \ + ? (__extension__ ({ __asm__ __volatile__ \ + ("mtfsb0 %s0" \ + : : "i#*X"(__builtin_ffs (__excepts))); \ + 0; })) \ + : 0) \ + : (feclearexcept) (__excepts)) + +#endif /* __GNUC__ && !_SOFT_FLOAT */ diff --git a/src/system/libroot/posix/glibc/include/arch/ppc/fenv_libc.h b/src/system/libroot/posix/glibc/include/arch/ppc/fenv_libc.h new file mode 100644 index 0000000000..7ae12a7d2b --- /dev/null +++ b/src/system/libroot/posix/glibc/include/arch/ppc/fenv_libc.h @@ -0,0 +1,106 @@ +/* Internal libc stuff for floating point environment routines. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library 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. + + The GNU C Library 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 the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#ifndef _FENV_LIBC_H +#define _FENV_LIBC_H 1 + +#include + +/* The sticky bits in the FPSCR indicating exceptions have occurred. */ +#define FPSCR_STICKY_BITS ((FE_ALL_EXCEPT | FE_ALL_INVALID) & ~FE_INVALID) + +/* Equivalent to fegetenv, but returns a fenv_t instead of taking a + pointer. */ +#define fegetenv_register() \ + ({ fenv_t env; asm volatile ("mffs %0" : "=f" (env)); env; }) + +/* Equivalent to fesetenv, but takes a fenv_t instead of a pointer. */ +#define fesetenv_register(env) \ + ({ double d = (env); asm volatile ("mtfsf 0xff,%0" : : "f" (d)); }) + +/* This very handy macro: + - Sets the rounding mode to 'round to nearest'; + - Sets the processor into IEEE mode; and + - Prevents exceptions from being raised for inexact results. + These things happen to be exactly what you need for typical elementary + functions. */ +#define relax_fenv_state() asm ("mtfsfi 7,0") + +/* Set/clear a particular FPSCR bit (for instance, + reset_fpscr_bit(FPSCR_VE); + prevents INVALID exceptions from being raised). */ +#define set_fpscr_bit(x) asm volatile ("mtfsb1 %0" : : "i"(x)) +#define reset_fpscr_bit(x) asm volatile ("mtfsb0 %0" : : "i"(x)) + +typedef union +{ + fenv_t fenv; + unsigned int l[2]; +} fenv_union_t; + +/* Definitions of all the FPSCR bit numbers */ +enum { + FPSCR_FX = 0, /* exception summary */ + FPSCR_FEX, /* enabled exception summary */ + FPSCR_VX, /* invalid operation summary */ + FPSCR_OX, /* overflow */ + FPSCR_UX, /* underflow */ + FPSCR_ZX, /* zero divide */ + FPSCR_XX, /* inexact */ + FPSCR_VXSNAN, /* invalid operation for SNaN */ + FPSCR_VXISI, /* invalid operation for Inf-Inf */ + FPSCR_VXIDI, /* invalid operation for Inf/Inf */ + FPSCR_VXZDZ, /* invalid operation for 0/0 */ + FPSCR_VXIMZ, /* invalid operation for Inf*0 */ + FPSCR_VXVC, /* invalid operation for invalid compare */ + FPSCR_FR, /* fraction rounded [fraction was incremented by round] */ + FPSCR_FI, /* fraction inexact */ + FPSCR_FPRF_C, /* result class descriptor */ + FPSCR_FPRF_FL, /* result less than (usually, less than 0) */ + FPSCR_FPRF_FG, /* result greater than */ + FPSCR_FPRF_FE, /* result equal to */ + FPSCR_FPRF_FU, /* result unordered */ + FPSCR_20, /* reserved */ + FPSCR_VXSOFT, /* invalid operation set by software */ + FPSCR_VXSQRT, /* invalid operation for square root */ + FPSCR_VXCVI, /* invalid operation for invalid integer convert */ + FPSCR_VE, /* invalid operation exception enable */ + FPSCR_OE, /* overflow exception enable */ + FPSCR_UE, /* underflow exception enable */ + FPSCR_ZE, /* zero divide exception enable */ + FPSCR_XE, /* inexact exception enable */ + FPSCR_NI /* non-IEEE mode (typically, no denormalised numbers) */ + /* the remaining two least-significant bits keep the rounding mode */ +}; + +/* This operation (i) sets the appropriate FPSCR bits for its + parameter, (ii) converts SNaN to the corresponding NaN, and (iii) + otherwise passes its parameter through unchanged (in particular, -0 + and +0 stay as they were). The `obvious' way to do this is optimised + out by gcc. */ +#define f_wash(x) \ + ({ double d; asm volatile ("fmul %0,%1,%2" \ + : "=f"(d) \ + : "f" (x), "f"((float)1.0)); d; }) +#define f_washf(x) \ + ({ float f; asm volatile ("fmuls %0,%1,%2" \ + : "=f"(f) \ + : "f" (x), "f"((float)1.0)); f; }) + +#endif /* fenv_libc.h */