Aren/haiku
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added more missing math functions for PPC from glibc 2.3.2.

Browse Source 
git-svn-id: file:///srv/svn/repos/haiku/haiku/trunk@15445 a95241bf-73f2-0310-859d-f6bbb57e9c96
This commit is contained in:
Ingo Weinhold 2005-12-09 14:16:39 +00:00
parent 7c360a1007
commit 379f146735
21 changed files with 2532 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

167
src/system/libroot/posix/glibc/arch/generic/atnat.h Normal file
View File

@ -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

66
src/system/libroot/posix/glibc/arch/generic/fraiseexcpt.c Normal file
View File

@ -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 <fenv_libc.h>
#include <bp-sym.h>
#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 <shlib-compat.h>
#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);

507
src/system/libroot/posix/glibc/arch/generic/mpa.c Normal file
View File

@ -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<m, the digits of x beyond the n'th are ignored. */
/* x=y is permissible. */
void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
int i,k;
EY = EX; k=MIN(m,n);
for (i=0; i <= k; i++) Y[i] = X[i];
for ( ; i <= n; i++) Y[i] = ZERO;
return;
}
/* Convert a multiple precision number *x into a double precision */
/* number *y, normalized case (|x| >= 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; i<EX; i++) c *= RADIX;
for (i=1; i>EX; 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]<TWO5))
{ *y=ZERO; return; }
if (p==1) {
if (EX==-42) {z[1]=X[1]+TWO10; z[2]=ZERO; z[3]=ZERO; k=3;}
else if (EX==-43) {z[1]= TWO10; z[2]=X[1]; z[3]=ZERO; k=2;}
else {z[1]= TWO10; z[2]=ZERO; z[3]=X[1]; k=1;}
}
else if (p==2) {
if (EX==-42) {z[1]=X[1]+TWO10; z[2]=X[2]; z[3]=ZERO; k=3;}
else if (EX==-43) {z[1]= TWO10; z[2]=X[1]; z[3]=X[2]; k=2;}
else {z[1]= TWO10; z[2]=ZERO; z[3]=X[1]; k=1;}
}
else {
if (EX==-42) {z[1]=X[1]+TWO10; z[2]=X[2]; k=3;}
else if (EX==-43) {z[1]= TWO10; z[2]=X[1]; k=2;}
else {z[1]= TWO10; z[2]=ZERO; k=1;}
z[3] = X[k];
}
u = (z[3] + TWO57) - TWO57;
if (u > 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<i2; i++,j--) Z[k] += X[i]*Y[j];
u = (Z[k] + CUTTER)-CUTTER;
if (u > 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; i<np1[p]; i++) {
__cpy(y,&w,p);
__mul(x,&w,y,p);
__sub(&mptwo,y,&z,p);
__mul(&w,&z,y,p);
}
return;
}
/* Divide one multiple precision number by another.Set *z = *x / *y. *x & *y */
/* are left unchanged. x&y may overlap but not x&z or y&z. */
/* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3 */
/* and 3.001*r**(1-p) for p>3. *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;
}

95
src/system/libroot/posix/glibc/arch/generic/mpa2.h Normal file
View File

@ -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

102
src/system/libroot/posix/glibc/arch/generic/mpatan.c Normal file
View File

@ -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; i<m; i++) {
__add(&mpone,&mpsm,&mpt1,p);
__mpsqrt(&mpt1,&mpt2,p);
__add(&mpt2,&mpt2,&mpt1,p);
__add(&mptwo,&mpsm,&mpt2,p);
__add(&mpt1,&mpt2,&mpt3,p);
__dvd(&mpsm,&mpt3,&mpt1,p);
__cpy(&mpt1,&mpsm,p);
}
__mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0];
}
/* Evaluate a truncated power series for Atan(s) */
n=np[p]; mptwoim1.d[1] = twonm1[p].d;
__dvd(&mpsm,&mptwoim1,&mpt,p);
for (i=n-1; i>1; 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;
}

174
src/system/libroot/posix/glibc/arch/generic/mpatan.h Normal file
View File

@ -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

70
src/system/libroot/posix/glibc/arch/generic/mpatan2.c Normal file
View File

@ -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;
}

103
src/system/libroot/posix/glibc/arch/generic/mpsqrt.c Normal file
View File

@ -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<m; i++) {
__mul(&mpu,&mpu,&mpt1,p);
__mul(&mpt1,&mpz,&mpt2,p);
__sub(&mp3halfs,&mpt2,&mpt1,p);
__mul(&mpu,&mpt1,&mpt2,p);
__cpy(&mpt2,&mpu,p);
}
__mul(&mpxn,&mpu,y,p); EY += ey;
return;
}
/***********************************************************/
/* Compute a double precision approximation for 1/sqrt(x) */
/* with the relative error bounded by 2**-51. */
/***********************************************************/
double fastiroot(double x) {
union {long i[2]; double d;} p,q;
double y,z, t;
long n;
static const double c0 = 0.99674, c1 = -0.53380, c2 = 0.45472, c3 = -0.21553;
p.d = x;
p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF ) | 0x3FE00000 ;
q.d = x;
y = p.d;
z = y -1.0;
n = (q.i[HIGH_HALF] - p.i[HIGH_HALF])>>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);
}

51
src/system/libroot/posix/glibc/arch/generic/mpsqrt.h Normal file
View File

@ -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

230
src/system/libroot/posix/glibc/arch/generic/s_atan.c Normal file
View File

@ -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<ZERO) ? -x : x;
if (u<C) {
if (u<B) {
if (u<A) { /* u < A */
return x; }
else { /* A <= u < B */
v=x*x; yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
if ((y=x+(yy-U1*x)) == x+(yy+U1*x)) return y;
EMULV(x,x,v,vv,t1,t2,t3,t4,t5) /* v+vv=x^2 */
s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1)) return y;
return atanMp(x,pr);
} }
else { /* B <= u < C */
i=(TWO52+TWO8*u)-TWO52; i-=16;
z=u-cij[i][0].d;
yy=z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
z*(cij[i][5].d+z* cij[i][6].d))));
t1=cij[i][1].d;
if (i<112) {
if (i<48) u2=U21; /* u < 1/4 */
else u2=U22; } /* 1/4 <= u < 1/2 */
else {
if (i<176) u2=U23; /* 1/2 <= u < 3/4 */
else u2=U24; } /* 3/4 <= u <= 1 */
if ((y=t1+(yy-u2*t1)) == t1+(yy+u2*t1)) return __signArctan(x,y);
z=u-hij[i][0].d;
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,ZERO,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,ZERO,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,ZERO,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,ZERO,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)
if ((y=s2+(ss2-U6*s2)) == s2+(ss2+U6*s2)) return __signArctan(x,y);
return atanMp(x,pr);
}
}
else {
if (u<D) { /* C <= u < D */
w=ONE/u;
EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
ww=w*((ONE-t1)-t2);
i=(TWO52+TWO8*w)-TWO52; i-=16;
z=(w-cij[i][0].d)+ww;
yy=HPI1-z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
z*(cij[i][5].d+z* cij[i][6].d))));
t1=HPI-cij[i][1].d;
if (i<112) u3=U31; /* w < 1/2 */
else u3=U32; /* w >= 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) { /* D <= u < E */
w=ONE/u; v=w*w;
EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
ww=w*((ONE-t1)-t2);
ESUB(HPI,w,t3,cor)
yy=((HPI1+cor)-ww)-yy;
if ((y=t3+(yy-U4)) == t3+(yy+U4)) return __signArctan(x,y);
DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
ADD2(w,ww,s2,ss2,s1,ss1,t1,t2)
SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2)
if ((y=s2+(ss2-U8)) == s2+(ss2+U8)) return __signArctan(x,y);
return atanMp(x,pr);
}
else {
/* 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<ZERO) return -y;
else return y;
}
/* Final stages. Compute atan(x) by multiple precision arithmetic */
static double atanMp(double x,const int pr[]){
mp_no mpx,mpy,mpy2,mperr,mpt1,mpy1;
double y1,y2;
int i,p;
for (i=0; i<M; i++) {
p = pr[i];
__dbl_mp(x,&mpx,p); __mpatan(&mpx,&mpy,p);
__dbl_mp(u9[i].d,&mpt1,p); __mul(&mpy,&mpt1,&mperr,p);
__add(&mpy,&mperr,&mpy1,p); __sub(&mpy,&mperr,&mpy2,p);
__mp_dbl(&mpy1,&y1,p); __mp_dbl(&mpy2,&y2,p);
if (y1==y2) return y1;
}
return y1; /*if unpossible to do exact computing */
}
#ifdef NO_LONG_DOUBLE
weak_alias (atan, atanl)
#endif

120
src/system/libroot/posix/glibc/arch/generic/s_atanf.c Normal file
View File

@ -0,0 +1,120 @@
/* s_atanf.c -- float version of s_atan.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_atanf.c,v 1.4 1995/05/10 20:46:47 jtc Exp $";
#endif
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
static const float atanhi[] = {
#else
static float atanhi[] = {
#endif
4.6364760399e-01, /* atan(0.5)hi 0x3eed6338 */
7.8539812565e-01, /* atan(1.0)hi 0x3f490fda */
9.8279368877e-01, /* atan(1.5)hi 0x3f7b985e */
1.5707962513e+00, /* atan(inf)hi 0x3fc90fda */
};
#ifdef __STDC__
static const float atanlo[] = {
#else
static float atanlo[] = {
#endif
5.0121582440e-09, /* atan(0.5)lo 0x31ac3769 */
3.7748947079e-08, /* atan(1.0)lo 0x33222168 */
3.4473217170e-08, /* atan(1.5)lo 0x33140fb4 */
7.5497894159e-08, /* atan(inf)lo 0x33a22168 */
};
#ifdef __STDC__
static const float aT[] = {
#else
static float aT[] = {
#endif
3.3333334327e-01, /* 0x3eaaaaaa */
-2.0000000298e-01, /* 0xbe4ccccd */
1.4285714924e-01, /* 0x3e124925 */
-1.1111110449e-01, /* 0xbde38e38 */
9.0908870101e-02, /* 0x3dba2e6e */
-7.6918758452e-02, /* 0xbd9d8795 */
6.6610731184e-02, /* 0x3d886b35 */
-5.8335702866e-02, /* 0xbd6ef16b */
4.9768779427e-02, /* 0x3d4bda59 */
-3.6531571299e-02, /* 0xbd15a221 */
1.6285819933e-02, /* 0x3c8569d7 */
};
#ifdef __STDC__
static const float
#else
static float
#endif
one = 1.0,
huge = 1.0e30;
#ifdef __STDC__
float __atanf(float x)
#else
float __atanf(x)
float x;
#endif
{
float w,s1,s2,z;
int32_t ix,hx,id;
GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff;
if(ix>=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)

14
src/system/libroot/posix/glibc/arch/ppc/Jamfile
View File

@ -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)
;

1
src/system/libroot/posix/glibc/arch/ppc/e_sqrtf.c Normal file
View File

@ -0,0 +1 @@
/* __ieee754_sqrtf is in w_sqrtf.c */

66
src/system/libroot/posix/glibc/arch/ppc/fraiseexcpt.c Normal file
View File

@ -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 <fenv_libc.h>
#include <bp-sym.h>
#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 <shlib-compat.h>
#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);

33
src/system/libroot/posix/glibc/arch/ppc/ftestexcept.c Normal file
View File

@ -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 <fenv_libc.h>
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;
}

144
src/system/libroot/posix/glibc/arch/ppc/t_sqrt.c Normal file
View File

@ -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 <math.h>
#include <stdio.h>
#include <assert.h>
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 */

146
src/system/libroot/posix/glibc/arch/ppc/w_sqrt.c Normal file
View File

@ -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 <math.h>
#include <math_private.h>
#include <fenv_libc.h>
#include <inttypes.h>
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

136
src/system/libroot/posix/glibc/arch/ppc/w_sqrtf.c Normal file
View File

@ -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 <math.h>
#include <math_private.h>
#include <fenv_libc.h>
#include <inttypes.h>
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)

145
src/system/libroot/posix/glibc/include/arch/ppc/bits/fenv.h Normal file
View File

@ -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 <bits/fenv.h> directly; include <fenv.h> 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

59
src/system/libroot/posix/glibc/include/arch/ppc/bits/fenvinline.h Normal file
View File

@ -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 */

106
src/system/libroot/posix/glibc/include/arch/ppc/fenv_libc.h Normal file
View File

@ -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 <fenv.h>
/* 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 */