Add some more PPC math functions from glibc 2.3.2.
git-svn-id: file:///srv/svn/repos/haiku/haiku/trunk@15406 a95241bf-73f2-0310-859d-f6bbb57e9c96
This commit is contained in:
Ingo Weinhold
@@ -0,0 +1,103 @@
|
||||
/*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
/********************************************************************/
|
||||
/* Ultimate math functions. Each function computes the exact */
|
||||
/* theoretical value of its argument rounded to nearest or even. */
|
||||
/* */
|
||||
/* Assumption: Machine arithmetic operations are performed in */
|
||||
/* round nearest mode of IEEE 754 standard. */
|
||||
/********************************************************************/
|
||||
|
||||
#ifndef UMATH_LIB
|
||||
#define UMATH_LIB
|
||||
/********************************************************************/
|
||||
/* Function changes the precision mode to IEEE 754 double precision */
|
||||
/* and the rounding mode to nearest or even. */
|
||||
/* It returns the original status of these modes. */
|
||||
/* See further explanations of usage in DPChange.h */
|
||||
/********************************************************************/
|
||||
unsigned short Init_Lib(void);
|
||||
|
||||
/********************************************************************/
|
||||
/* Function that changes the precision and rounding modes to the */
|
||||
/* specified by the argument received. See further explanations in */
|
||||
/* DPChange.h */
|
||||
/********************************************************************/
|
||||
void Exit_Lib(unsigned short);
|
||||
|
||||
|
||||
/* The asin() function calculates the arc sine of its argument. */
|
||||
/* The function returns the arc sine in radians */
|
||||
/* (between -PI/2 and PI/2). */
|
||||
/* If the argument is greater than 1 or less than -1 it returns */
|
||||
/* a NaN. */
|
||||
double uasin(double );
|
||||
|
||||
|
||||
/* The acos() function calculates the arc cosine of its argument. */
|
||||
/* The function returns the arc cosine in radians */
|
||||
/* (between -PI/2 and PI/2). */
|
||||
/* If the argument is greater than 1 or less than -1 it returns */
|
||||
/* a NaN. */
|
||||
double uacos(double );
|
||||
|
||||
/* The atan() function calculates the arctanget of its argument. */
|
||||
/* The function returns the arc tangent in radians */
|
||||
/* (between -PI/2 and PI/2). */
|
||||
double uatan(double );
|
||||
|
||||
|
||||
/* The uatan2() function calculates the arc tangent of the two arguments x */
|
||||
/* and y (x is the right argument and y is the left one).The signs of both */
|
||||
/* arguments are used to determine the quadrant of the result. */
|
||||
/* The function returns the result in radians, which is between -PI and PI */
|
||||
double uatan2(double ,double );
|
||||
|
||||
/* Compute log(x). The base of log is e (natural logarithm) */
|
||||
double ulog(double );
|
||||
|
||||
/* Compute e raised to the power of argument x. */
|
||||
double uexp(double );
|
||||
|
||||
/* Compute sin(x). The argument x is assumed to be given in radians.*/
|
||||
double usin(double );
|
||||
|
||||
/* Compute cos(x). The argument x is assumed to be given in radians.*/
|
||||
double ucos(double );
|
||||
|
||||
/* Compute tan(x). The argument x is assumed to be given in radians.*/
|
||||
double utan(double );
|
||||
|
||||
/* Compute the square root of non-negative argument x. */
|
||||
/* If x is negative the returned value is NaN. */
|
||||
double usqrt(double );
|
||||
|
||||
/* Compute x raised to the power of y, where x is the left argument */
|
||||
/* and y is the right argument. The function returns a NaN if x<0. */
|
||||
/* If x equals zero it returns -inf */
|
||||
double upow(double , double );
|
||||
|
||||
/* Computing x mod y, where x is the left argument and y is the */
|
||||
/* right one. */
|
||||
double uremainder(double , double );
|
||||
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,185 @@
|
||||
|
||||
/*
|
||||
* 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: atnat2.h */
|
||||
/* */
|
||||
/* */
|
||||
/* common data and variables definition for BIG or LITTLE ENDIAN */
|
||||
/************************************************************************/
|
||||
|
||||
|
||||
|
||||
#ifndef ATNAT2_H
|
||||
#define ATNAT2_H
|
||||
|
||||
|
||||
#define MM 5
|
||||
#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 */
|
||||
/**/ mzero = {{0x80000000, 0x00000000} }, /* -0 */
|
||||
/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */
|
||||
/**/ inv16 = {{0x3fb00000, 0x00000000} }, /* 1/16 */
|
||||
/**/ opi = {{0x400921fb, 0x54442d18} }, /* pi */
|
||||
/**/ opi1 = {{0x3ca1a626, 0x33145c07} }, /* pi-opi */
|
||||
/**/ mopi = {{0xc00921fb, 0x54442d18} }, /* -pi */
|
||||
/**/ hpi = {{0x3ff921fb, 0x54442d18} }, /* pi/2 */
|
||||
/**/ hpi1 = {{0x3c91a626, 0x33145c07} }, /* pi/2-hpi */
|
||||
/**/ mhpi = {{0xbff921fb, 0x54442d18} }, /* -pi/2 */
|
||||
/**/ qpi = {{0x3fe921fb, 0x54442d18} }, /* pi/4 */
|
||||
/**/ qpi1 = {{0x3c81a626, 0x33145c07} }, /* pi/4-qpi */
|
||||
/**/ mqpi = {{0xbfe921fb, 0x54442d18} }, /* -pi/4 */
|
||||
/**/ tqpi = {{0x4002d97c, 0x7f3321d2} }, /* 3pi/4 */
|
||||
/**/ tqpi1 = {{0x3c9a7939, 0x4c9e8a0a} }, /* 3pi/4-tqpi */
|
||||
/**/ mtqpi = {{0xc002d97c, 0x7f3321d2} }, /* -3pi/4 */
|
||||
/**/ u1 = {{0x3c314c2a, 0x00000000} }, /* 9.377e-19 */
|
||||
/**/ u2 = {{0x3bf955e4, 0x00000000} }, /* 8.584e-20 */
|
||||
/**/ u3 = {{0x3bf955e4, 0x00000000} }, /* 8.584e-20 */
|
||||
/**/ u4 = {{0x3bf955e4, 0x00000000} }, /* 8.584e-20 */
|
||||
/**/ u5 = {{0x3aaef2d1, 0x00000000} }, /* 5e-26 */
|
||||
/**/ u6 = {{0x3a6eeb36, 0x00000000} }, /* 3.122e-27 */
|
||||
/**/ u7 = {{0x3a6eeb36, 0x00000000} }, /* 3.122e-27 */
|
||||
/**/ u8 = {{0x3a6eeb36, 0x00000000} }, /* 3.122e-27 */
|
||||
/**/ u91 = {{0x3c6dffc0, 0x00000000} }, /* 1.301e-17 */
|
||||
/**/ u92 = {{0x3c527bd0, 0x00000000} }, /* 4.008e-18 */
|
||||
/**/ u93 = {{0x3c3cd057, 0x00000000} }, /* 1.562e-18 */
|
||||
/**/ u94 = {{0x3c329cdf, 0x00000000} }, /* 1.009e-18 */
|
||||
/**/ ua1 = {{0x3c3a1edf, 0x00000000} }, /* 1.416e-18 */
|
||||
/**/ ua2 = {{0x3c33f0e1, 0x00000000} }, /* 1.081e-18 */
|
||||
/**/ ub = {{0x3a98c56d, 0x00000000} }, /* 2.001e-26 */
|
||||
/**/ uc = {{0x3a9375de, 0x00000000} }, /* 1.572e-26 */
|
||||
/**/ ud[MM] ={{{0x38c6eddf, 0x00000000} }, /* 3.450e-35 */
|
||||
/**/ {{0x35c6ef60, 0x00000000} }, /* 1.226e-49 */
|
||||
/**/ {{0x32c6ed2f, 0x00000000} }, /* 4.354e-64 */
|
||||
/**/ {{0x23c6eee8, 0x00000000} }, /* 2.465e-136 */
|
||||
/**/ {{0x11c6ed16, 0x00000000} }},/* 4.955e-223 */
|
||||
/**/ ue = {{0x38900e9d, 0x00000000} }, /* 3.02e-36 */
|
||||
/**/ two8 = {{0x40700000, 0x00000000} }, /* 2**8=256 */
|
||||
/**/ two52 = {{0x43300000, 0x00000000} }, /* 2**52 */
|
||||
/**/ two500 = {{0x5f300000, 0x00000000} }, /* 2**500 */
|
||||
/**/ twom500 = {{0x20b00000, 0x00000000} }, /* 2**(-500) */
|
||||
/**/ twom1022 = {{0x00100000, 0x00000000} }; /* 2**(-1022) */
|
||||
|
||||
#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 */
|
||||
/**/ mzero = {{0x00000000, 0x80000000} }, /* -0 */
|
||||
/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */
|
||||
/**/ inv16 = {{0x00000000, 0x3fb00000} }, /* 1/16 */
|
||||
/**/ opi = {{0x54442d18, 0x400921fb} }, /* pi */
|
||||
/**/ opi1 = {{0x33145c07, 0x3ca1a626} }, /* pi-opi */
|
||||
/**/ mopi = {{0x54442d18, 0xc00921fb} }, /* -pi */
|
||||
/**/ hpi = {{0x54442d18, 0x3ff921fb} }, /* pi/2 */
|
||||
/**/ hpi1 = {{0x33145c07, 0x3c91a626} }, /* pi/2-hpi */
|
||||
/**/ mhpi = {{0x54442d18, 0xbff921fb} }, /* -pi/2 */
|
||||
/**/ qpi = {{0x54442d18, 0x3fe921fb} }, /* pi/4 */
|
||||
/**/ qpi1 = {{0x33145c07, 0x3c81a626} }, /* pi/4-qpi */
|
||||
/**/ mqpi = {{0x54442d18, 0xbfe921fb} }, /* -pi/4 */
|
||||
/**/ tqpi = {{0x7f3321d2, 0x4002d97c} }, /* 3pi/4 */
|
||||
/**/ tqpi1 = {{0x4c9e8a0a, 0x3c9a7939} }, /* 3pi/4-tqpi */
|
||||
/**/ mtqpi = {{0x7f3321d2, 0xc002d97c} }, /* -3pi/4 */
|
||||
/**/ u1 = {{0x00000000, 0x3c314c2a} }, /* 9.377e-19 */
|
||||
/**/ u2 = {{0x00000000, 0x3bf955e4} }, /* 8.584e-20 */
|
||||
/**/ u3 = {{0x00000000, 0x3bf955e4} }, /* 8.584e-20 */
|
||||
/**/ u4 = {{0x00000000, 0x3bf955e4} }, /* 8.584e-20 */
|
||||
/**/ u5 = {{0x00000000, 0x3aaef2d1} }, /* 5e-26 */
|
||||
/**/ u6 = {{0x00000000, 0x3a6eeb36} }, /* 3.122e-27 */
|
||||
/**/ u7 = {{0x00000000, 0x3a6eeb36} }, /* 3.122e-27 */
|
||||
/**/ u8 = {{0x00000000, 0x3a6eeb36} }, /* 3.122e-27 */
|
||||
/**/ u91 = {{0x00000000, 0x3c6dffc0} }, /* 1.301e-17 */
|
||||
/**/ u92 = {{0x00000000, 0x3c527bd0} }, /* 4.008e-18 */
|
||||
/**/ u93 = {{0x00000000, 0x3c3cd057} }, /* 1.562e-18 */
|
||||
/**/ u94 = {{0x00000000, 0x3c329cdf} }, /* 1.009e-18 */
|
||||
/**/ ua1 = {{0x00000000, 0x3c3a1edf} }, /* 1.416e-18 */
|
||||
/**/ ua2 = {{0x00000000, 0x3c33f0e1} }, /* 1.081e-18 */
|
||||
/**/ ub = {{0x00000000, 0x3a98c56d} }, /* 2.001e-26 */
|
||||
/**/ uc = {{0x00000000, 0x3a9375de} }, /* 1.572e-26 */
|
||||
/**/ ud[MM] ={{{0x00000000, 0x38c6eddf} }, /* 3.450e-35 */
|
||||
/**/ {{0x00000000, 0x35c6ef60} }, /* 1.226e-49 */
|
||||
/**/ {{0x00000000, 0x32c6ed2f} }, /* 4.354e-64 */
|
||||
/**/ {{0x00000000, 0x23c6eee8} }, /* 2.465e-136 */
|
||||
/**/ {{0x00000000, 0x11c6ed16} }},/* 4.955e-223 */
|
||||
/**/ ue = {{0x00000000, 0x38900e9d} }, /* 3.02e-36 */
|
||||
/**/ two8 = {{0x00000000, 0x40700000} }, /* 2**8=256 */
|
||||
/**/ two52 = {{0x00000000, 0x43300000} }, /* 2**52 */
|
||||
/**/ two500 = {{0x00000000, 0x5f300000} }, /* 2**500 */
|
||||
/**/ twom500 = {{0x00000000, 0x20b00000} }, /* 2**(-500) */
|
||||
/**/ twom1022 = {{0x00000000, 0x00100000} }; /* 2**(-1022) */
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#define ZERO zero.d
|
||||
#define MZERO mzero.d
|
||||
#define ONE one.d
|
||||
#define TWO8 two8.d
|
||||
#define TWO52 two52.d
|
||||
#define TWOM1022 twom1022.d
|
||||
|
||||
#endif
|
||||
@@ -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: dla.h */
|
||||
/* */
|
||||
/* This file holds C language macros for 'Double Length Floating Point */
|
||||
/* Arithmetic'. The macros are based on the paper: */
|
||||
/* T.J.Dekker, "A floating-point Technique for extending the */
|
||||
/* Available Precision", Number. Math. 18, 224-242 (1971). */
|
||||
/* A Double-Length number is defined by a pair (r,s), of IEEE double */
|
||||
/* precision floating point numbers that satisfy, */
|
||||
/* */
|
||||
/* abs(s) <= abs(r+s)*2**(-53)/(1+2**(-53)). */
|
||||
/* */
|
||||
/* The computer arithmetic assumed is IEEE double precision in */
|
||||
/* round to nearest mode. All variables in the macros must be of type */
|
||||
/* IEEE double. */
|
||||
/***********************************************************************/
|
||||
|
||||
/* CN = 1+2**27 = '41a0000002000000' IEEE double format */
|
||||
#define CN 134217729.0
|
||||
|
||||
|
||||
/* Exact addition of two single-length floating point numbers, Dekker. */
|
||||
/* The macro produces a double-length number (z,zz) that satisfies */
|
||||
/* z+zz = x+y exactly. */
|
||||
|
||||
#define EADD(x,y,z,zz) \
|
||||
z=(x)+(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x));
|
||||
|
||||
|
||||
/* Exact subtraction of two single-length floating point numbers, Dekker. */
|
||||
/* The macro produces a double-length number (z,zz) that satisfies */
|
||||
/* z+zz = x-y exactly. */
|
||||
|
||||
#define ESUB(x,y,z,zz) \
|
||||
z=(x)-(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z)));
|
||||
|
||||
|
||||
/* Exact multiplication of two single-length floating point numbers, */
|
||||
/* Veltkamp. The macro produces a double-length number (z,zz) that */
|
||||
/* satisfies z+zz = x*y exactly. p,hx,tx,hy,ty are temporary */
|
||||
/* storage variables of type double. */
|
||||
|
||||
#define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \
|
||||
p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \
|
||||
p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \
|
||||
z=(x)*(y); zz=(((hx*hy-z)+hx*ty)+tx*hy)+tx*ty;
|
||||
|
||||
|
||||
/* Exact multiplication of two single-length floating point numbers, Dekker. */
|
||||
/* The macro produces a nearly double-length number (z,zz) (see Dekker) */
|
||||
/* that satisfies z+zz = x*y exactly. p,hx,tx,hy,ty,q are temporary */
|
||||
/* storage variables of type double. */
|
||||
|
||||
#define MUL12(x,y,z,zz,p,hx,tx,hy,ty,q) \
|
||||
p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \
|
||||
p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \
|
||||
p=hx*hy; q=hx*ty+tx*hy; z=p+q; zz=((p-z)+q)+tx*ty;
|
||||
|
||||
|
||||
/* Double-length addition, Dekker. The macro produces a double-length */
|
||||
/* number (z,zz) which satisfies approximately z+zz = x+xx + y+yy. */
|
||||
/* An error bound: (abs(x+xx)+abs(y+yy))*4.94e-32. (x,xx), (y,yy) */
|
||||
/* are assumed to be double-length numbers. r,s are temporary */
|
||||
/* storage variables of type double. */
|
||||
|
||||
#define ADD2(x,xx,y,yy,z,zz,r,s) \
|
||||
r=(x)+(y); s=(ABS(x)>ABS(y)) ? \
|
||||
(((((x)-r)+(y))+(yy))+(xx)) : \
|
||||
(((((y)-r)+(x))+(xx))+(yy)); \
|
||||
z=r+s; zz=(r-z)+s;
|
||||
|
||||
|
||||
/* Double-length subtraction, Dekker. The macro produces a double-length */
|
||||
/* number (z,zz) which satisfies approximately z+zz = x+xx - (y+yy). */
|
||||
/* An error bound: (abs(x+xx)+abs(y+yy))*4.94e-32. (x,xx), (y,yy) */
|
||||
/* are assumed to be double-length numbers. r,s are temporary */
|
||||
/* storage variables of type double. */
|
||||
|
||||
#define SUB2(x,xx,y,yy,z,zz,r,s) \
|
||||
r=(x)-(y); s=(ABS(x)>ABS(y)) ? \
|
||||
(((((x)-r)-(y))-(yy))+(xx)) : \
|
||||
((((x)-((y)+r))+(xx))-(yy)); \
|
||||
z=r+s; zz=(r-z)+s;
|
||||
|
||||
|
||||
/* Double-length multiplication, Dekker. The macro produces a double-length */
|
||||
/* number (z,zz) which satisfies approximately z+zz = (x+xx)*(y+yy). */
|
||||
/* An error bound: abs((x+xx)*(y+yy))*1.24e-31. (x,xx), (y,yy) */
|
||||
/* are assumed to be double-length numbers. p,hx,tx,hy,ty,q,c,cc are */
|
||||
/* temporary storage variables of type double. */
|
||||
|
||||
#define MUL2(x,xx,y,yy,z,zz,p,hx,tx,hy,ty,q,c,cc) \
|
||||
MUL12(x,y,c,cc,p,hx,tx,hy,ty,q) \
|
||||
cc=((x)*(yy)+(xx)*(y))+cc; z=c+cc; zz=(c-z)+cc;
|
||||
|
||||
|
||||
/* Double-length division, Dekker. The macro produces a double-length */
|
||||
/* number (z,zz) which satisfies approximately z+zz = (x+xx)/(y+yy). */
|
||||
/* An error bound: abs((x+xx)/(y+yy))*1.50e-31. (x,xx), (y,yy) */
|
||||
/* are assumed to be double-length numbers. p,hx,tx,hy,ty,q,c,cc,u,uu */
|
||||
/* are temporary storage variables of type double. */
|
||||
|
||||
#define DIV2(x,xx,y,yy,z,zz,p,hx,tx,hy,ty,q,c,cc,u,uu) \
|
||||
c=(x)/(y); MUL12(c,y,u,uu,p,hx,tx,hy,ty,q) \
|
||||
cc=(((((x)-u)-uu)+(xx))-c*(yy))/(y); z=c+cc; zz=(c-z)+cc;
|
||||
|
||||
|
||||
/* Double-length addition, slower but more accurate than ADD2. */
|
||||
/* The macro produces a double-length */
|
||||
/* number (z,zz) which satisfies approximately z+zz = (x+xx)+(y+yy). */
|
||||
/* An error bound: abs(x+xx + y+yy)*1.50e-31. (x,xx), (y,yy) */
|
||||
/* are assumed to be double-length numbers. r,rr,s,ss,u,uu,w */
|
||||
/* are temporary storage variables of type double. */
|
||||
|
||||
#define ADD2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w) \
|
||||
r=(x)+(y); \
|
||||
if (ABS(x)>ABS(y)) { rr=((x)-r)+(y); s=(rr+(yy))+(xx); } \
|
||||
else { rr=((y)-r)+(x); s=(rr+(xx))+(yy); } \
|
||||
if (rr!=0.0) { \
|
||||
z=r+s; zz=(r-z)+s; } \
|
||||
else { \
|
||||
ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)+(yy)) : (((yy)-s)+(xx)); \
|
||||
u=r+s; \
|
||||
uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \
|
||||
w=uu+ss; z=u+w; \
|
||||
zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; }
|
||||
|
||||
|
||||
/* Double-length subtraction, slower but more accurate than SUB2. */
|
||||
/* The macro produces a double-length */
|
||||
/* number (z,zz) which satisfies approximately z+zz = (x+xx)-(y+yy). */
|
||||
/* An error bound: abs(x+xx - (y+yy))*1.50e-31. (x,xx), (y,yy) */
|
||||
/* are assumed to be double-length numbers. r,rr,s,ss,u,uu,w */
|
||||
/* are temporary storage variables of type double. */
|
||||
|
||||
#define SUB2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w) \
|
||||
r=(x)-(y); \
|
||||
if (ABS(x)>ABS(y)) { rr=((x)-r)-(y); s=(rr-(yy))+(xx); } \
|
||||
else { rr=(x)-((y)+r); s=(rr+(xx))-(yy); } \
|
||||
if (rr!=0.0) { \
|
||||
z=r+s; zz=(r-z)+s; } \
|
||||
else { \
|
||||
ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)-(yy)) : ((xx)-((yy)+s)); \
|
||||
u=r+s; \
|
||||
uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \
|
||||
w=uu+ss; z=u+w; \
|
||||
zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; }
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1 @@
|
||||
/* In e_asin.c */
|
||||
@@ -0,0 +1,89 @@
|
||||
/* e_acosf.c -- float version of e_acos.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* 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: e_acosf.c,v 1.5 1995/05/12 04:57:16 jtc Exp $";
|
||||
#endif
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const float
|
||||
#else
|
||||
static float
|
||||
#endif
|
||||
one = 1.0000000000e+00, /* 0x3F800000 */
|
||||
pi = 3.1415925026e+00, /* 0x40490fda */
|
||||
pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */
|
||||
pio2_lo = 7.5497894159e-08, /* 0x33a22168 */
|
||||
pS0 = 1.6666667163e-01, /* 0x3e2aaaab */
|
||||
pS1 = -3.2556581497e-01, /* 0xbea6b090 */
|
||||
pS2 = 2.0121252537e-01, /* 0x3e4e0aa8 */
|
||||
pS3 = -4.0055535734e-02, /* 0xbd241146 */
|
||||
pS4 = 7.9153501429e-04, /* 0x3a4f7f04 */
|
||||
pS5 = 3.4793309169e-05, /* 0x3811ef08 */
|
||||
qS1 = -2.4033949375e+00, /* 0xc019d139 */
|
||||
qS2 = 2.0209457874e+00, /* 0x4001572d */
|
||||
qS3 = -6.8828397989e-01, /* 0xbf303361 */
|
||||
qS4 = 7.7038154006e-02; /* 0x3d9dc62e */
|
||||
|
||||
#ifdef __STDC__
|
||||
float __ieee754_acosf(float x)
|
||||
#else
|
||||
float __ieee754_acosf(x)
|
||||
float x;
|
||||
#endif
|
||||
{
|
||||
float z,p,q,r,w,s,c,df;
|
||||
int32_t hx,ix;
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
ix = hx&0x7fffffff;
|
||||
if(ix==0x3f800000) { /* |x|==1 */
|
||||
if(hx>0) return 0.0; /* acos(1) = 0 */
|
||||
else return pi+(float)2.0*pio2_lo; /* acos(-1)= pi */
|
||||
} else if(ix>0x3f800000) { /* |x| >= 1 */
|
||||
return (x-x)/(x-x); /* acos(|x|>1) is NaN */
|
||||
}
|
||||
if(ix<0x3f000000) { /* |x| < 0.5 */
|
||||
if(ix<=0x23000000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
|
||||
z = x*x;
|
||||
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
|
||||
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
|
||||
r = p/q;
|
||||
return pio2_hi - (x - (pio2_lo-x*r));
|
||||
} else if (hx<0) { /* x < -0.5 */
|
||||
z = (one+x)*(float)0.5;
|
||||
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
|
||||
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
|
||||
s = __ieee754_sqrtf(z);
|
||||
r = p/q;
|
||||
w = r*s-pio2_lo;
|
||||
return pi - (float)2.0*(s+w);
|
||||
} else { /* x > 0.5 */
|
||||
int32_t idf;
|
||||
z = (one-x)*(float)0.5;
|
||||
s = __ieee754_sqrtf(z);
|
||||
df = s;
|
||||
GET_FLOAT_WORD(idf,df);
|
||||
SET_FLOAT_WORD(df,idf&0xfffff000);
|
||||
c = (z-df*df)/(s+df);
|
||||
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
|
||||
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
|
||||
r = p/q;
|
||||
w = r*s+c;
|
||||
return (float)2.0*(df+w);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,406 @@
|
||||
/*
|
||||
* 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: atnat2.c */
|
||||
/* */
|
||||
/* FUNCTIONS: uatan2 */
|
||||
/* atan2Mp */
|
||||
/* signArctan2 */
|
||||
/* normalized */
|
||||
/* */
|
||||
/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat2.h */
|
||||
/* mpatan.c mpatan2.c mpsqrt.c */
|
||||
/* uatan.tbl */
|
||||
/* */
|
||||
/* An ultimate atan2() routine. Given two IEEE double machine numbers y,*/
|
||||
/* x it computes the correctly rounded (to nearest) value of atan2(y,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 "atnat2.h"
|
||||
#include "math_private.h"
|
||||
|
||||
/************************************************************************/
|
||||
/* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */
|
||||
/* it computes the correctly rounded (to nearest) value of atan2(y,x). */
|
||||
/* Assumption: Machine arithmetic operations are performed in */
|
||||
/* round to nearest mode of IEEE 754 standard. */
|
||||
/************************************************************************/
|
||||
static double atan2Mp(double ,double ,const int[]);
|
||||
static double signArctan2(double ,double);
|
||||
static double normalized(double ,double,double ,double);
|
||||
void __mpatan2(mp_no *,mp_no *,mp_no *,int);
|
||||
|
||||
double __ieee754_atan2(double y,double x) {
|
||||
|
||||
int i,de,ux,dx,uy,dy;
|
||||
#if 0
|
||||
int p;
|
||||
#endif
|
||||
static const int pr[MM]={6,8,10,20,32};
|
||||
double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t4,t5,t6,t7,t8,
|
||||
z,zz,cor,s1,ss1,s2,ss2;
|
||||
#if 0
|
||||
double z1,z2;
|
||||
#endif
|
||||
number num;
|
||||
#if 0
|
||||
mp_no mperr,mpt1,mpx,mpy,mpz,mpz1,mpz2;
|
||||
#endif
|
||||
|
||||
static const int ep= 59768832, /* 57*16**5 */
|
||||
em=-59768832; /* -57*16**5 */
|
||||
|
||||
/* x=NaN or y=NaN */
|
||||
num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF];
|
||||
if ((ux&0x7ff00000) ==0x7ff00000) {
|
||||
if (((ux&0x000fffff)|dx)!=0x00000000) return x+x; }
|
||||
num.d = y; uy = num.i[HIGH_HALF]; dy = num.i[LOW_HALF];
|
||||
if ((uy&0x7ff00000) ==0x7ff00000) {
|
||||
if (((uy&0x000fffff)|dy)!=0x00000000) return y+y; }
|
||||
|
||||
/* y=+-0 */
|
||||
if (uy==0x00000000) {
|
||||
if (dy==0x00000000) {
|
||||
if ((ux&0x80000000)==0x00000000) return ZERO;
|
||||
else return opi.d; } }
|
||||
else if (uy==0x80000000) {
|
||||
if (dy==0x00000000) {
|
||||
if ((ux&0x80000000)==0x00000000) return MZERO;
|
||||
else return mopi.d;} }
|
||||
|
||||
/* x=+-0 */
|
||||
if (x==ZERO) {
|
||||
if ((uy&0x80000000)==0x00000000) return hpi.d;
|
||||
else return mhpi.d; }
|
||||
|
||||
/* x=+-INF */
|
||||
if (ux==0x7ff00000) {
|
||||
if (dx==0x00000000) {
|
||||
if (uy==0x7ff00000) {
|
||||
if (dy==0x00000000) return qpi.d; }
|
||||
else if (uy==0xfff00000) {
|
||||
if (dy==0x00000000) return mqpi.d; }
|
||||
else {
|
||||
if ((uy&0x80000000)==0x00000000) return ZERO;
|
||||
else return MZERO; }
|
||||
}
|
||||
}
|
||||
else if (ux==0xfff00000) {
|
||||
if (dx==0x00000000) {
|
||||
if (uy==0x7ff00000) {
|
||||
if (dy==0x00000000) return tqpi.d; }
|
||||
else if (uy==0xfff00000) {
|
||||
if (dy==0x00000000) return mtqpi.d; }
|
||||
else {
|
||||
if ((uy&0x80000000)==0x00000000) return opi.d;
|
||||
else return mopi.d; }
|
||||
}
|
||||
}
|
||||
|
||||
/* y=+-INF */
|
||||
if (uy==0x7ff00000) {
|
||||
if (dy==0x00000000) return hpi.d; }
|
||||
else if (uy==0xfff00000) {
|
||||
if (dy==0x00000000) return mhpi.d; }
|
||||
|
||||
/* either x/y or y/x is very close to zero */
|
||||
ax = (x<ZERO) ? -x : x; ay = (y<ZERO) ? -y : y;
|
||||
de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
|
||||
if (de>=ep) { return ((y>ZERO) ? hpi.d : mhpi.d); }
|
||||
else if (de<=em) {
|
||||
if (x>ZERO) {
|
||||
if ((z=ay/ax)<TWOM1022) return normalized(ax,ay,y,z);
|
||||
else return signArctan2(y,z); }
|
||||
else { return ((y>ZERO) ? opi.d : mopi.d); } }
|
||||
|
||||
/* if either x or y is extremely close to zero, scale abs(x), abs(y). */
|
||||
if (ax<twom500.d || ay<twom500.d) { ax*=two500.d; ay*=two500.d; }
|
||||
|
||||
/* x,y which are neither special nor extreme */
|
||||
if (ay<ax) {
|
||||
u=ay/ax;
|
||||
EMULV(ax,u,v,vv,t1,t2,t3,t4,t5)
|
||||
du=((ay-v)-vv)/ax; }
|
||||
else {
|
||||
u=ax/ay;
|
||||
EMULV(ay,u,v,vv,t1,t2,t3,t4,t5)
|
||||
du=((ax-v)-vv)/ay; }
|
||||
|
||||
if (x>ZERO) {
|
||||
|
||||
/* (i) x>0, abs(y)< abs(x): atan(ay/ax) */
|
||||
if (ay<ax) {
|
||||
if (u<inv16.d) {
|
||||
v=u*u; zz=du+u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
|
||||
if ((z=u+(zz-u1.d*u)) == u+(zz+u1.d*u)) return signArctan2(y,z);
|
||||
|
||||
MUL2(u,du,u,du,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(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
|
||||
ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
|
||||
if ((z=s1+(ss1-u5.d*s1)) == s1+(ss1+u5.d*s1)) return signArctan2(y,z);
|
||||
return atan2Mp(x,y,pr);
|
||||
}
|
||||
else {
|
||||
i=(TWO52+TWO8*u)-TWO52; i-=16;
|
||||
t3=u-cij[i][0].d;
|
||||
EADD(t3,du,v,dv)
|
||||
t1=cij[i][1].d; t2=cij[i][2].d;
|
||||
zz=v*t2+(dv*t2+v*v*(cij[i][3].d+v*(cij[i][4].d+
|
||||
v*(cij[i][5].d+v* cij[i][6].d))));
|
||||
if (i<112) {
|
||||
if (i<48) u9=u91.d; /* u < 1/4 */
|
||||
else u9=u92.d; } /* 1/4 <= u < 1/2 */
|
||||
else {
|
||||
if (i<176) u9=u93.d; /* 1/2 <= u < 3/4 */
|
||||
else u9=u94.d; } /* 3/4 <= u <= 1 */
|
||||
if ((z=t1+(zz-u9*t1)) == t1+(zz+u9*t1)) return signArctan2(y,z);
|
||||
|
||||
t1=u-hij[i][0].d;
|
||||
EADD(t1,du,v,vv)
|
||||
s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
|
||||
v*(hij[i][14].d+v* hij[i][15].d))));
|
||||
ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
|
||||
MUL2(v,vv,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(v,vv,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(v,vv,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(v,vv,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 ((z=s2+(ss2-ub.d*s2)) == s2+(ss2+ub.d*s2)) return signArctan2(y,z);
|
||||
return atan2Mp(x,y,pr);
|
||||
}
|
||||
}
|
||||
|
||||
/* (ii) x>0, abs(x)<=abs(y): pi/2-atan(ax/ay) */
|
||||
else {
|
||||
if (u<inv16.d) {
|
||||
v=u*u;
|
||||
zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
|
||||
ESUB(hpi.d,u,t2,cor)
|
||||
t3=((hpi1.d+cor)-du)-zz;
|
||||
if ((z=t2+(t3-u2.d)) == t2+(t3+u2.d)) return signArctan2(y,z);
|
||||
|
||||
MUL2(u,du,u,du,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(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
|
||||
ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
|
||||
SUB2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
|
||||
if ((z=s2+(ss2-u6.d)) == s2+(ss2+u6.d)) return signArctan2(y,z);
|
||||
return atan2Mp(x,y,pr);
|
||||
}
|
||||
else {
|
||||
i=(TWO52+TWO8*u)-TWO52; i-=16;
|
||||
v=(u-cij[i][0].d)+du;
|
||||
zz=hpi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
|
||||
v*(cij[i][5].d+v* cij[i][6].d))));
|
||||
t1=hpi.d-cij[i][1].d;
|
||||
if (i<112) ua=ua1.d; /* w < 1/2 */
|
||||
else ua=ua2.d; /* w >= 1/2 */
|
||||
if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z);
|
||||
|
||||
t1=u-hij[i][0].d;
|
||||
EADD(t1,du,v,vv)
|
||||
s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
|
||||
v*(hij[i][14].d+v* hij[i][15].d))));
|
||||
ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
|
||||
MUL2(v,vv,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(v,vv,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(v,vv,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(v,vv,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.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
|
||||
if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z);
|
||||
return atan2Mp(x,y,pr);
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
|
||||
/* (iii) x<0, abs(x)< abs(y): pi/2+atan(ax/ay) */
|
||||
if (ax<ay) {
|
||||
if (u<inv16.d) {
|
||||
v=u*u;
|
||||
zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
|
||||
EADD(hpi.d,u,t2,cor)
|
||||
t3=((hpi1.d+cor)+du)+zz;
|
||||
if ((z=t2+(t3-u3.d)) == t2+(t3+u3.d)) return signArctan2(y,z);
|
||||
|
||||
MUL2(u,du,u,du,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(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
|
||||
ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
|
||||
ADD2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
|
||||
if ((z=s2+(ss2-u7.d)) == s2+(ss2+u7.d)) return signArctan2(y,z);
|
||||
return atan2Mp(x,y,pr);
|
||||
}
|
||||
else {
|
||||
i=(TWO52+TWO8*u)-TWO52; i-=16;
|
||||
v=(u-cij[i][0].d)+du;
|
||||
zz=hpi1.d+v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
|
||||
v*(cij[i][5].d+v* cij[i][6].d))));
|
||||
t1=hpi.d+cij[i][1].d;
|
||||
if (i<112) ua=ua1.d; /* w < 1/2 */
|
||||
else ua=ua2.d; /* w >= 1/2 */
|
||||
if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z);
|
||||
|
||||
t1=u-hij[i][0].d;
|
||||
EADD(t1,du,v,vv)
|
||||
s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
|
||||
v*(hij[i][14].d+v* hij[i][15].d))));
|
||||
ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
|
||||
MUL2(v,vv,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(v,vv,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(v,vv,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(v,vv,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)
|
||||
ADD2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
|
||||
if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z);
|
||||
return atan2Mp(x,y,pr);
|
||||
}
|
||||
}
|
||||
|
||||
/* (iv) x<0, abs(y)<=abs(x): pi-atan(ax/ay) */
|
||||
else {
|
||||
if (u<inv16.d) {
|
||||
v=u*u;
|
||||
zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
|
||||
ESUB(opi.d,u,t2,cor)
|
||||
t3=((opi1.d+cor)-du)-zz;
|
||||
if ((z=t2+(t3-u4.d)) == t2+(t3+u4.d)) return signArctan2(y,z);
|
||||
|
||||
MUL2(u,du,u,du,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(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
|
||||
ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
|
||||
SUB2(opi.d,opi1.d,s1,ss1,s2,ss2,t1,t2)
|
||||
if ((z=s2+(ss2-u8.d)) == s2+(ss2+u8.d)) return signArctan2(y,z);
|
||||
return atan2Mp(x,y,pr);
|
||||
}
|
||||
else {
|
||||
i=(TWO52+TWO8*u)-TWO52; i-=16;
|
||||
v=(u-cij[i][0].d)+du;
|
||||
zz=opi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
|
||||
v*(cij[i][5].d+v* cij[i][6].d))));
|
||||
t1=opi.d-cij[i][1].d;
|
||||
if (i<112) ua=ua1.d; /* w < 1/2 */
|
||||
else ua=ua2.d; /* w >= 1/2 */
|
||||
if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z);
|
||||
|
||||
t1=u-hij[i][0].d;
|
||||
EADD(t1,du,v,vv)
|
||||
s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
|
||||
v*(hij[i][14].d+v* hij[i][15].d))));
|
||||
ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
|
||||
MUL2(v,vv,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(v,vv,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(v,vv,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(v,vv,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(opi.d,opi1.d,s2,ss2,s1,ss1,t1,t2)
|
||||
if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z);
|
||||
return atan2Mp(x,y,pr);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
/* Treat the Denormalized case */
|
||||
static double normalized(double ax,double ay,double y, double z)
|
||||
{ int p;
|
||||
mp_no mpx,mpy,mpz,mperr,mpz2,mpt1;
|
||||
p=6;
|
||||
__dbl_mp(ax,&mpx,p); __dbl_mp(ay,&mpy,p); __dvd(&mpy,&mpx,&mpz,p);
|
||||
__dbl_mp(ue.d,&mpt1,p); __mul(&mpz,&mpt1,&mperr,p);
|
||||
__sub(&mpz,&mperr,&mpz2,p); __mp_dbl(&mpz2,&z,p);
|
||||
return signArctan2(y,z);
|
||||
}
|
||||
/* Fix the sign and return after stage 1 or stage 2 */
|
||||
static double signArctan2(double y,double z)
|
||||
{
|
||||
return ((y<ZERO) ? -z : z);
|
||||
}
|
||||
/* Stage 3: Perform a multi-Precision computation */
|
||||
static double atan2Mp(double x,double y,const int pr[])
|
||||
{
|
||||
double z1,z2;
|
||||
int i,p;
|
||||
mp_no mpx,mpy,mpz,mpz1,mpz2,mperr,mpt1;
|
||||
for (i=0; i<MM; i++) {
|
||||
p = pr[i];
|
||||
__dbl_mp(x,&mpx,p); __dbl_mp(y,&mpy,p);
|
||||
__mpatan2(&mpy,&mpx,&mpz,p);
|
||||
__dbl_mp(ud[i].d,&mpt1,p); __mul(&mpz,&mpt1,&mperr,p);
|
||||
__add(&mpz,&mperr,&mpz1,p); __sub(&mpz,&mperr,&mpz2,p);
|
||||
__mp_dbl(&mpz1,&z1,p); __mp_dbl(&mpz2,&z2,p);
|
||||
if (z1==z2) return z1;
|
||||
}
|
||||
return z1; /*if unpossible to do exact computing */
|
||||
}
|
||||
@@ -0,0 +1,105 @@
|
||||
/* e_atan2f.c -- float version of e_atan2.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* 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: e_atan2f.c,v 1.4 1995/05/10 20:44:53 jtc Exp $";
|
||||
#endif
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const float
|
||||
#else
|
||||
static float
|
||||
#endif
|
||||
tiny = 1.0e-30,
|
||||
zero = 0.0,
|
||||
pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
|
||||
pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
|
||||
pi = 3.1415927410e+00, /* 0x40490fdb */
|
||||
pi_lo = -8.7422776573e-07; /* 0xb3bbbd2e */
|
||||
|
||||
#ifdef __STDC__
|
||||
float __ieee754_atan2f(float y, float x)
|
||||
#else
|
||||
float __ieee754_atan2f(y,x)
|
||||
float y,x;
|
||||
#endif
|
||||
{
|
||||
float z;
|
||||
int32_t k,m,hx,hy,ix,iy;
|
||||
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
ix = hx&0x7fffffff;
|
||||
GET_FLOAT_WORD(hy,y);
|
||||
iy = hy&0x7fffffff;
|
||||
if((ix>0x7f800000)||
|
||||
(iy>0x7f800000)) /* x or y is NaN */
|
||||
return x+y;
|
||||
if(hx==0x3f800000) return __atanf(y); /* x=1.0 */
|
||||
m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
|
||||
|
||||
/* when y = 0 */
|
||||
if(iy==0) {
|
||||
switch(m) {
|
||||
case 0:
|
||||
case 1: return y; /* atan(+-0,+anything)=+-0 */
|
||||
case 2: return pi+tiny;/* atan(+0,-anything) = pi */
|
||||
case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
|
||||
}
|
||||
}
|
||||
/* when x = 0 */
|
||||
if(ix==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
|
||||
|
||||
/* when x is INF */
|
||||
if(ix==0x7f800000) {
|
||||
if(iy==0x7f800000) {
|
||||
switch(m) {
|
||||
case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
|
||||
case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
|
||||
case 2: return (float)3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
|
||||
case 3: return (float)-3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
|
||||
}
|
||||
} else {
|
||||
switch(m) {
|
||||
case 0: return zero ; /* atan(+...,+INF) */
|
||||
case 1: return -zero ; /* atan(-...,+INF) */
|
||||
case 2: return pi+tiny ; /* atan(+...,-INF) */
|
||||
case 3: return -pi-tiny ; /* atan(-...,-INF) */
|
||||
}
|
||||
}
|
||||
}
|
||||
/* when y is INF */
|
||||
if(iy==0x7f800000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
|
||||
|
||||
/* compute y/x */
|
||||
k = (iy-ix)>>23;
|
||||
if(k > 60) z=pi_o_2+(float)0.5*pi_lo; /* |y/x| > 2**60 */
|
||||
else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */
|
||||
else z=__atanf(fabsf(y/x)); /* safe to do y/x */
|
||||
switch (m) {
|
||||
case 0: return z ; /* atan(+,+) */
|
||||
case 1: {
|
||||
u_int32_t zh;
|
||||
GET_FLOAT_WORD(zh,z);
|
||||
SET_FLOAT_WORD(z,zh ^ 0x80000000);
|
||||
}
|
||||
return z ; /* atan(-,+) */
|
||||
case 2: return pi-(z-pi_lo);/* atan(+,-) */
|
||||
default: /* case 3 */
|
||||
return (z-pi_lo)-pi;/* atan(-,-) */
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,140 @@
|
||||
/* @(#)e_fmod.c 5.1 93/09/24 */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: e_fmod.c,v 1.8 1995/05/10 20:45:07 jtc Exp $";
|
||||
#endif
|
||||
|
||||
/*
|
||||
* __ieee754_fmod(x,y)
|
||||
* Return x mod y in exact arithmetic
|
||||
* Method: shift and subtract
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const double one = 1.0, Zero[] = {0.0, -0.0,};
|
||||
#else
|
||||
static double one = 1.0, Zero[] = {0.0, -0.0,};
|
||||
#endif
|
||||
|
||||
#ifdef __STDC__
|
||||
double __ieee754_fmod(double x, double y)
|
||||
#else
|
||||
double __ieee754_fmod(x,y)
|
||||
double x,y ;
|
||||
#endif
|
||||
{
|
||||
int32_t n,hx,hy,hz,ix,iy,sx,i;
|
||||
u_int32_t lx,ly,lz;
|
||||
|
||||
EXTRACT_WORDS(hx,lx,x);
|
||||
EXTRACT_WORDS(hy,ly,y);
|
||||
sx = hx&0x80000000; /* sign of x */
|
||||
hx ^=sx; /* |x| */
|
||||
hy &= 0x7fffffff; /* |y| */
|
||||
|
||||
/* purge off exception values */
|
||||
if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
|
||||
((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
|
||||
return (x*y)/(x*y);
|
||||
if(hx<=hy) {
|
||||
if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
|
||||
if(lx==ly)
|
||||
return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
|
||||
}
|
||||
|
||||
/* determine ix = ilogb(x) */
|
||||
if(hx<0x00100000) { /* subnormal x */
|
||||
if(hx==0) {
|
||||
for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
|
||||
} else {
|
||||
for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
|
||||
}
|
||||
} else ix = (hx>>20)-1023;
|
||||
|
||||
/* determine iy = ilogb(y) */
|
||||
if(hy<0x00100000) { /* subnormal y */
|
||||
if(hy==0) {
|
||||
for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
|
||||
} else {
|
||||
for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
|
||||
}
|
||||
} else iy = (hy>>20)-1023;
|
||||
|
||||
/* set up {hx,lx}, {hy,ly} and align y to x */
|
||||
if(ix >= -1022)
|
||||
hx = 0x00100000|(0x000fffff&hx);
|
||||
else { /* subnormal x, shift x to normal */
|
||||
n = -1022-ix;
|
||||
if(n<=31) {
|
||||
hx = (hx<<n)|(lx>>(32-n));
|
||||
lx <<= n;
|
||||
} else {
|
||||
hx = lx<<(n-32);
|
||||
lx = 0;
|
||||
}
|
||||
}
|
||||
if(iy >= -1022)
|
||||
hy = 0x00100000|(0x000fffff&hy);
|
||||
else { /* subnormal y, shift y to normal */
|
||||
n = -1022-iy;
|
||||
if(n<=31) {
|
||||
hy = (hy<<n)|(ly>>(32-n));
|
||||
ly <<= n;
|
||||
} else {
|
||||
hy = ly<<(n-32);
|
||||
ly = 0;
|
||||
}
|
||||
}
|
||||
|
||||
/* fix point fmod */
|
||||
n = ix - iy;
|
||||
while(n--) {
|
||||
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
|
||||
if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
|
||||
else {
|
||||
if((hz|lz)==0) /* return sign(x)*0 */
|
||||
return Zero[(u_int32_t)sx>>31];
|
||||
hx = hz+hz+(lz>>31); lx = lz+lz;
|
||||
}
|
||||
}
|
||||
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
|
||||
if(hz>=0) {hx=hz;lx=lz;}
|
||||
|
||||
/* convert back to floating value and restore the sign */
|
||||
if((hx|lx)==0) /* return sign(x)*0 */
|
||||
return Zero[(u_int32_t)sx>>31];
|
||||
while(hx<0x00100000) { /* normalize x */
|
||||
hx = hx+hx+(lx>>31); lx = lx+lx;
|
||||
iy -= 1;
|
||||
}
|
||||
if(iy>= -1022) { /* normalize output */
|
||||
hx = ((hx-0x00100000)|((iy+1023)<<20));
|
||||
INSERT_WORDS(x,hx|sx,lx);
|
||||
} else { /* subnormal output */
|
||||
n = -1022 - iy;
|
||||
if(n<=20) {
|
||||
lx = (lx>>n)|((u_int32_t)hx<<(32-n));
|
||||
hx >>= n;
|
||||
} else if (n<=31) {
|
||||
lx = (hx<<(32-n))|(lx>>n); hx = sx;
|
||||
} else {
|
||||
lx = hx>>(n-32); hx = sx;
|
||||
}
|
||||
INSERT_WORDS(x,hx|sx,lx);
|
||||
x *= one; /* create necessary signal */
|
||||
}
|
||||
return x; /* exact output */
|
||||
}
|
||||
@@ -0,0 +1,113 @@
|
||||
/* e_fmodf.c -- float version of e_fmod.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* 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: e_fmodf.c,v 1.4 1995/05/10 20:45:10 jtc Exp $";
|
||||
#endif
|
||||
|
||||
/*
|
||||
* __ieee754_fmodf(x,y)
|
||||
* Return x mod y in exact arithmetic
|
||||
* Method: shift and subtract
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const float one = 1.0, Zero[] = {0.0, -0.0,};
|
||||
#else
|
||||
static float one = 1.0, Zero[] = {0.0, -0.0,};
|
||||
#endif
|
||||
|
||||
#ifdef __STDC__
|
||||
float __ieee754_fmodf(float x, float y)
|
||||
#else
|
||||
float __ieee754_fmodf(x,y)
|
||||
float x,y ;
|
||||
#endif
|
||||
{
|
||||
int32_t n,hx,hy,hz,ix,iy,sx,i;
|
||||
|
||||
GET_FLOAT_WORD(hx,x);
|
||||
GET_FLOAT_WORD(hy,y);
|
||||
sx = hx&0x80000000; /* sign of x */
|
||||
hx ^=sx; /* |x| */
|
||||
hy &= 0x7fffffff; /* |y| */
|
||||
|
||||
/* purge off exception values */
|
||||
if(hy==0||(hx>=0x7f800000)|| /* y=0,or x not finite */
|
||||
(hy>0x7f800000)) /* or y is NaN */
|
||||
return (x*y)/(x*y);
|
||||
if(hx<hy) return x; /* |x|<|y| return x */
|
||||
if(hx==hy)
|
||||
return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
|
||||
|
||||
/* determine ix = ilogb(x) */
|
||||
if(hx<0x00800000) { /* subnormal x */
|
||||
for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
|
||||
} else ix = (hx>>23)-127;
|
||||
|
||||
/* determine iy = ilogb(y) */
|
||||
if(hy<0x00800000) { /* subnormal y */
|
||||
for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1;
|
||||
} else iy = (hy>>23)-127;
|
||||
|
||||
/* set up {hx,lx}, {hy,ly} and align y to x */
|
||||
if(ix >= -126)
|
||||
hx = 0x00800000|(0x007fffff&hx);
|
||||
else { /* subnormal x, shift x to normal */
|
||||
n = -126-ix;
|
||||
hx = hx<<n;
|
||||
}
|
||||
if(iy >= -126)
|
||||
hy = 0x00800000|(0x007fffff&hy);
|
||||
else { /* subnormal y, shift y to normal */
|
||||
n = -126-iy;
|
||||
hy = hy<<n;
|
||||
}
|
||||
|
||||
/* fix point fmod */
|
||||
n = ix - iy;
|
||||
while(n--) {
|
||||
hz=hx-hy;
|
||||
if(hz<0){hx = hx+hx;}
|
||||
else {
|
||||
if(hz==0) /* return sign(x)*0 */
|
||||
return Zero[(u_int32_t)sx>>31];
|
||||
hx = hz+hz;
|
||||
}
|
||||
}
|
||||
hz=hx-hy;
|
||||
if(hz>=0) {hx=hz;}
|
||||
|
||||
/* convert back to floating value and restore the sign */
|
||||
if(hx==0) /* return sign(x)*0 */
|
||||
return Zero[(u_int32_t)sx>>31];
|
||||
while(hx<0x00800000) { /* normalize x */
|
||||
hx = hx+hx;
|
||||
iy -= 1;
|
||||
}
|
||||
if(iy>= -126) { /* normalize output */
|
||||
hx = ((hx-0x00800000)|((iy+127)<<23));
|
||||
SET_FLOAT_WORD(x,hx|sx);
|
||||
} else { /* subnormal output */
|
||||
n = -126 - iy;
|
||||
hx >>= n;
|
||||
SET_FLOAT_WORD(x,hx|sx);
|
||||
x *= one; /* create necessary signal */
|
||||
}
|
||||
return x; /* exact output */
|
||||
}
|
||||
@@ -0,0 +1,87 @@
|
||||
|
||||
/*
|
||||
* 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: mpa.h */
|
||||
/* */
|
||||
/* FUNCTIONS: */
|
||||
/* mcr */
|
||||
/* acr */
|
||||
/* cr */
|
||||
/* cpy */
|
||||
/* cpymn */
|
||||
/* mp_dbl */
|
||||
/* dbl_mp */
|
||||
/* add */
|
||||
/* sub */
|
||||
/* mul */
|
||||
/* inv */
|
||||
/* dvd */
|
||||
/* */
|
||||
/* Arithmetic functions for multiple precision numbers. */
|
||||
/* Common types and definition */
|
||||
/************************************************************************/
|
||||
|
||||
|
||||
typedef struct {/* This structure holds the details of a multi-precision */
|
||||
int e; /* floating point number, x: d[0] holds its sign (-1,0 or 1) */
|
||||
double d[40]; /* e holds its exponent (...,-2,-1,0,1,2,...) and */
|
||||
} mp_no; /* d[1]...d[p] hold its mantissa digits. The value of x is, */
|
||||
/* x = d[1]*r**(e-1) + d[2]*r**(e-2) + ... + d[p]*r**(e-p). */
|
||||
/* Here r = 2**24, 0 <= d[i] < r and 1 <= p <= 32. */
|
||||
/* p is a global variable. A multi-precision number is */
|
||||
/* always normalized. Namely, d[1] > 0. An exception is */
|
||||
/* a zero which is characterized by d[0] = 0. The terms */
|
||||
/* d[p+1], d[p+2], ... of a none zero number have no */
|
||||
/* significance and so are the terms e, d[1],d[2],... */
|
||||
/* of a zero. */
|
||||
|
||||
typedef union { int i[2]; double d; } number;
|
||||
|
||||
#define X x->d
|
||||
#define Y y->d
|
||||
#define Z z->d
|
||||
#define EX x->e
|
||||
#define EY y->e
|
||||
#define EZ z->e
|
||||
|
||||
#define MAX(x,y) ((x) < (y) ? (y) : (x))
|
||||
#define MIN(x,y) ((x) < (y) ? (x) : (y))
|
||||
#define ABS(x) ((x) < 0 ? -(x) : (x))
|
||||
|
||||
int __acr(const mp_no *, const mp_no *, int);
|
||||
int __cr(const mp_no *, const mp_no *, int);
|
||||
void __cpy(const mp_no *, mp_no *, int);
|
||||
void __cpymn(const mp_no *, int, mp_no *, int);
|
||||
void __mp_dbl(const mp_no *, double *, int);
|
||||
void __dbl_mp(double, mp_no *, int);
|
||||
void __add(const mp_no *, const mp_no *, mp_no *, int);
|
||||
void __sub(const mp_no *, const mp_no *, mp_no *, int);
|
||||
void __mul(const mp_no *, const mp_no *, mp_no *, int);
|
||||
void __inv(const mp_no *, mp_no *, int);
|
||||
void __dvd(const mp_no *, const mp_no *, mp_no *, int);
|
||||
|
||||
extern void __mpatan (mp_no *, mp_no *, int);
|
||||
extern void __mpatan2 (mp_no *, mp_no *, mp_no *, int);
|
||||
extern void __mpsqrt (mp_no *, mp_no *, int);
|
||||
extern void __mpexp (mp_no *, mp_no *__y, int);
|
||||
extern void __c32 (mp_no *, mp_no *, mp_no *, int);
|
||||
extern int __mpranred (double, mp_no *, int);
|
||||
@@ -0,0 +1,85 @@
|
||||
/* @(#)s_ceil.c 5.1 93/09/24 */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: s_ceil.c,v 1.8 1995/05/10 20:46:53 jtc Exp $";
|
||||
#endif
|
||||
|
||||
/*
|
||||
* ceil(x)
|
||||
* Return x rounded toward -inf to integral value
|
||||
* Method:
|
||||
* Bit twiddling.
|
||||
* Exception:
|
||||
* Inexact flag raised if x not equal to ceil(x).
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const double huge = 1.0e300;
|
||||
#else
|
||||
static double huge = 1.0e300;
|
||||
#endif
|
||||
|
||||
#ifdef __STDC__
|
||||
double __ceil(double x)
|
||||
#else
|
||||
double __ceil(x)
|
||||
double x;
|
||||
#endif
|
||||
{
|
||||
int32_t i0,i1,j0;
|
||||
u_int32_t i,j;
|
||||
EXTRACT_WORDS(i0,i1,x);
|
||||
j0 = ((i0>>20)&0x7ff)-0x3ff;
|
||||
if(j0<20) {
|
||||
if(j0<0) { /* raise inexact if x != 0 */
|
||||
if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
|
||||
if(i0<0) {i0=0x80000000;i1=0;}
|
||||
else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;}
|
||||
}
|
||||
} else {
|
||||
i = (0x000fffff)>>j0;
|
||||
if(((i0&i)|i1)==0) return x; /* x is integral */
|
||||
if(huge+x>0.0) { /* raise inexact flag */
|
||||
if(i0>0) i0 += (0x00100000)>>j0;
|
||||
i0 &= (~i); i1=0;
|
||||
}
|
||||
}
|
||||
} else if (j0>51) {
|
||||
if(j0==0x400) return x+x; /* inf or NaN */
|
||||
else return x; /* x is integral */
|
||||
} else {
|
||||
i = ((u_int32_t)(0xffffffff))>>(j0-20);
|
||||
if((i1&i)==0) return x; /* x is integral */
|
||||
if(huge+x>0.0) { /* raise inexact flag */
|
||||
if(i0>0) {
|
||||
if(j0==20) i0+=1;
|
||||
else {
|
||||
j = i1 + (1<<(52-j0));
|
||||
if(j<i1) i0+=1; /* got a carry */
|
||||
i1 = j;
|
||||
}
|
||||
}
|
||||
i1 &= (~i);
|
||||
}
|
||||
}
|
||||
INSERT_WORDS(x,i0,i1);
|
||||
return x;
|
||||
}
|
||||
weak_alias (__ceil, ceil)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__ceil, __ceill)
|
||||
weak_alias (__ceil, ceill)
|
||||
#endif
|
||||
@@ -0,0 +1,41 @@
|
||||
/* @(#)s_finite.c 5.1 93/09/24 */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: s_finite.c,v 1.8 1995/05/10 20:47:17 jtc Exp $";
|
||||
#endif
|
||||
|
||||
/*
|
||||
* finite(x) returns 1 is x is finite, else 0;
|
||||
* no branching!
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
int __finite(double x)
|
||||
#else
|
||||
int __finite(x)
|
||||
double x;
|
||||
#endif
|
||||
{
|
||||
int32_t hx;
|
||||
GET_HIGH_WORD(hx,x);
|
||||
return (int)((u_int32_t)((hx&0x7fffffff)-0x7ff00000)>>31);
|
||||
}
|
||||
hidden_def (__finite)
|
||||
weak_alias (__finite, finite)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__finite, __finitel)
|
||||
weak_alias (__finite, finitel)
|
||||
#endif
|
||||
@@ -0,0 +1,40 @@
|
||||
/* s_finitef.c -- float version of s_finite.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* 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_finitef.c,v 1.4 1995/05/10 20:47:18 jtc Exp $";
|
||||
#endif
|
||||
|
||||
/*
|
||||
* finitef(x) returns 1 is x is finite, else 0;
|
||||
* no branching!
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
int __finitef(float x)
|
||||
#else
|
||||
int __finitef(x)
|
||||
float x;
|
||||
#endif
|
||||
{
|
||||
int32_t ix;
|
||||
GET_FLOAT_WORD(ix,x);
|
||||
return (int)((u_int32_t)((ix&0x7fffffff)-0x7f800000)>>31);
|
||||
}
|
||||
hidden_def (__finitef)
|
||||
weak_alias (__finitef, finitef)
|
||||
@@ -0,0 +1,86 @@
|
||||
/* @(#)s_floor.c 5.1 93/09/24 */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: s_floor.c,v 1.8 1995/05/10 20:47:20 jtc Exp $";
|
||||
#endif
|
||||
|
||||
/*
|
||||
* floor(x)
|
||||
* Return x rounded toward -inf to integral value
|
||||
* Method:
|
||||
* Bit twiddling.
|
||||
* Exception:
|
||||
* Inexact flag raised if x not equal to floor(x).
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const double huge = 1.0e300;
|
||||
#else
|
||||
static double huge = 1.0e300;
|
||||
#endif
|
||||
|
||||
#ifdef __STDC__
|
||||
double __floor(double x)
|
||||
#else
|
||||
double __floor(x)
|
||||
double x;
|
||||
#endif
|
||||
{
|
||||
int32_t i0,i1,j0;
|
||||
u_int32_t i,j;
|
||||
EXTRACT_WORDS(i0,i1,x);
|
||||
j0 = ((i0>>20)&0x7ff)-0x3ff;
|
||||
if(j0<20) {
|
||||
if(j0<0) { /* raise inexact if x != 0 */
|
||||
if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
|
||||
if(i0>=0) {i0=i1=0;}
|
||||
else if(((i0&0x7fffffff)|i1)!=0)
|
||||
{ i0=0xbff00000;i1=0;}
|
||||
}
|
||||
} else {
|
||||
i = (0x000fffff)>>j0;
|
||||
if(((i0&i)|i1)==0) return x; /* x is integral */
|
||||
if(huge+x>0.0) { /* raise inexact flag */
|
||||
if(i0<0) i0 += (0x00100000)>>j0;
|
||||
i0 &= (~i); i1=0;
|
||||
}
|
||||
}
|
||||
} else if (j0>51) {
|
||||
if(j0==0x400) return x+x; /* inf or NaN */
|
||||
else return x; /* x is integral */
|
||||
} else {
|
||||
i = ((u_int32_t)(0xffffffff))>>(j0-20);
|
||||
if((i1&i)==0) return x; /* x is integral */
|
||||
if(huge+x>0.0) { /* raise inexact flag */
|
||||
if(i0<0) {
|
||||
if(j0==20) i0+=1;
|
||||
else {
|
||||
j = i1+(1<<(52-j0));
|
||||
if(j<i1) i0 +=1 ; /* got a carry */
|
||||
i1=j;
|
||||
}
|
||||
}
|
||||
i1 &= (~i);
|
||||
}
|
||||
}
|
||||
INSERT_WORDS(x,i0,i1);
|
||||
return x;
|
||||
}
|
||||
weak_alias (__floor, floor)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__floor, __floorl)
|
||||
weak_alias (__floor, floorl)
|
||||
#endif
|
||||
@@ -0,0 +1,47 @@
|
||||
/* @(#)s_logb.c 5.1 93/09/24 */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: s_logb.c,v 1.8 1995/05/10 20:47:50 jtc Exp $";
|
||||
#endif
|
||||
|
||||
/*
|
||||
* double logb(x)
|
||||
* IEEE 754 logb. Included to pass IEEE test suite. Not recommend.
|
||||
* Use ilogb instead.
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
double __logb(double x)
|
||||
#else
|
||||
double __logb(x)
|
||||
double x;
|
||||
#endif
|
||||
{
|
||||
int32_t lx,ix;
|
||||
EXTRACT_WORDS(ix,lx,x);
|
||||
ix &= 0x7fffffff; /* high |x| */
|
||||
if((ix|lx)==0) return -1.0/fabs(x);
|
||||
if(ix>=0x7ff00000) return x*x;
|
||||
if((ix>>=20)==0) /* IEEE 754 logb */
|
||||
return -1022.0;
|
||||
else
|
||||
return (double) (ix-1023);
|
||||
}
|
||||
weak_alias (__logb, logb)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__logb, __logbl)
|
||||
weak_alias (__logb, logbl)
|
||||
#endif
|
||||
@@ -0,0 +1,40 @@
|
||||
/* s_logbf.c -- float version of s_logb.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* 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_logbf.c,v 1.4 1995/05/10 20:47:51 jtc Exp $";
|
||||
#endif
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
float __logbf(float x)
|
||||
#else
|
||||
float __logbf(x)
|
||||
float x;
|
||||
#endif
|
||||
{
|
||||
int32_t ix;
|
||||
GET_FLOAT_WORD(ix,x);
|
||||
ix &= 0x7fffffff; /* high |x| */
|
||||
if(ix==0) return (float)-1.0/fabsf(x);
|
||||
if(ix>=0x7f800000) return x*x;
|
||||
if((ix>>=23)==0) /* IEEE 754 logb */
|
||||
return -126.0;
|
||||
else
|
||||
return (float) (ix-127);
|
||||
}
|
||||
weak_alias (__logbf, logbf)
|
||||
@@ -0,0 +1,70 @@
|
||||
/* @(#)s_scalbn.c 5.1 93/09/24 */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: s_scalbn.c,v 1.8 1995/05/10 20:48:08 jtc Exp $";
|
||||
#endif
|
||||
|
||||
/*
|
||||
* scalbn (double x, int n)
|
||||
* scalbn(x,n) returns x* 2**n computed by exponent
|
||||
* manipulation rather than by actually performing an
|
||||
* exponentiation or a multiplication.
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const double
|
||||
#else
|
||||
static double
|
||||
#endif
|
||||
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
|
||||
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
|
||||
huge = 1.0e+300,
|
||||
tiny = 1.0e-300;
|
||||
|
||||
#ifdef __STDC__
|
||||
double __scalbn (double x, int n)
|
||||
#else
|
||||
double __scalbn (x,n)
|
||||
double x; int n;
|
||||
#endif
|
||||
{
|
||||
int32_t k,hx,lx;
|
||||
EXTRACT_WORDS(hx,lx,x);
|
||||
k = (hx&0x7ff00000)>>20; /* extract exponent */
|
||||
if (k==0) { /* 0 or subnormal x */
|
||||
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
|
||||
x *= two54;
|
||||
GET_HIGH_WORD(hx,x);
|
||||
k = ((hx&0x7ff00000)>>20) - 54;
|
||||
}
|
||||
if (k==0x7ff) return x+x; /* NaN or Inf */
|
||||
k = k+n;
|
||||
if (n> 50000 || k > 0x7fe)
|
||||
return huge*__copysign(huge,x); /* overflow */
|
||||
if (n< -50000) return tiny*__copysign(tiny,x); /*underflow*/
|
||||
if (k > 0) /* normal result */
|
||||
{SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
|
||||
if (k <= -54)
|
||||
return tiny*__copysign(tiny,x); /*underflow*/
|
||||
k += 54; /* subnormal result */
|
||||
SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
|
||||
return x*twom54;
|
||||
}
|
||||
weak_alias (__scalbn, scalbn)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__scalbn, __scalbnl)
|
||||
weak_alias (__scalbn, scalbnl)
|
||||
#endif
|
||||
@@ -0,0 +1,63 @@
|
||||
/* s_scalbnf.c -- float version of s_scalbn.c.
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
|
||||
*/
|
||||
|
||||
/*
|
||||
* ====================================================
|
||||
* 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_scalbnf.c,v 1.4 1995/05/10 20:48:10 jtc Exp $";
|
||||
#endif
|
||||
|
||||
#include "math.h"
|
||||
#include "math_private.h"
|
||||
|
||||
#ifdef __STDC__
|
||||
static const float
|
||||
#else
|
||||
static float
|
||||
#endif
|
||||
two25 = 3.355443200e+07, /* 0x4c000000 */
|
||||
twom25 = 2.9802322388e-08, /* 0x33000000 */
|
||||
huge = 1.0e+30,
|
||||
tiny = 1.0e-30;
|
||||
|
||||
#ifdef __STDC__
|
||||
float __scalbnf (float x, int n)
|
||||
#else
|
||||
float __scalbnf (x,n)
|
||||
float x; int n;
|
||||
#endif
|
||||
{
|
||||
int32_t k,ix;
|
||||
GET_FLOAT_WORD(ix,x);
|
||||
k = (ix&0x7f800000)>>23; /* extract exponent */
|
||||
if (k==0) { /* 0 or subnormal x */
|
||||
if ((ix&0x7fffffff)==0) return x; /* +-0 */
|
||||
x *= two25;
|
||||
GET_FLOAT_WORD(ix,x);
|
||||
k = ((ix&0x7f800000)>>23) - 25;
|
||||
}
|
||||
if (k==0xff) return x+x; /* NaN or Inf */
|
||||
k = k+n;
|
||||
if (n> 50000 || k > 0xfe)
|
||||
return huge*copysignf(huge,x); /* overflow */
|
||||
if (n< -50000)
|
||||
return tiny*copysignf(tiny,x); /*underflow*/
|
||||
if (k > 0) /* normal result */
|
||||
{SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
|
||||
if (k <= -25)
|
||||
return tiny*copysignf(tiny,x); /*underflow*/
|
||||
k += 25; /* subnormal result */
|
||||
SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
|
||||
return x*twom25;
|
||||
}
|
||||
weak_alias (__scalbnf, scalbnf)
|
||||
File diff suppressed because it is too large
Load Diff
@@ -16,53 +16,42 @@ DEBUG = 0 ;
|
||||
|
||||
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
|
||||
s_fpclassify.c s_fpclassifyf.c s_fpclassifyl.c s_modf.c s_modff.c s_modfl.c
|
||||
w_atan2.c w_atan2f.c w_atan2l.c w_fmod.c w_fmodl.c w_fmodf.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_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
|
||||
s_fpclassify.c s_fpclassifyf.c s_fpclassifyl.c
|
||||
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
|
||||
;
|
||||
|
||||
MergeObject posix_gnu_arch_$(TARGET_ARCH).o :
|
||||
add_n.S
|
||||
addmul_1.S
|
||||
# e_acos.S
|
||||
# e_atan2.S
|
||||
# e_atan2f.S
|
||||
# e_atan2l.c
|
||||
# e_fmod.S
|
||||
# e_fmodf.S
|
||||
# e_fmodl.c
|
||||
ldbl2mpn.c
|
||||
mul_1.S
|
||||
lshift.S
|
||||
rshift.S
|
||||
# s_ceil.S
|
||||
s_ceilf.c
|
||||
# s_ceill.S
|
||||
# s_copysign.S
|
||||
# s_copysignf.S
|
||||
s_copysign.S
|
||||
s_copysignf.S
|
||||
# s_copysignl.S
|
||||
# s_fdim.S
|
||||
# s_fdimf.S
|
||||
s_fdim.c
|
||||
s_fdimf.c
|
||||
# s_fdiml.S
|
||||
# s_fabs.S
|
||||
# s_fabsf.S
|
||||
s_fabs.S
|
||||
s_fabsf.S
|
||||
# s_fabsl.S
|
||||
# s_finite.S
|
||||
# s_finitef.S
|
||||
# s_finitel.S
|
||||
# s_floor.S
|
||||
s_floorf.c
|
||||
# s_floorl.S
|
||||
# s_logb.S
|
||||
# s_logbf.S
|
||||
# s_logbl.c
|
||||
# s_rint.S
|
||||
# s_rintf.S
|
||||
s_rint.c
|
||||
s_rintf.c
|
||||
# s_rintl.c
|
||||
# s_scalbn.S
|
||||
# s_scalbnf.S
|
||||
# s_scalbnl.S
|
||||
sub_n.S
|
||||
submul_1.S
|
||||
|
||||
|
||||
@@ -0,0 +1,50 @@
|
||||
/* Copy a sign bit between floating-point values.
|
||||
Copyright (C) 1997, 1999, 2000 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. */
|
||||
|
||||
/* This has been coded in assembler because GCC makes such a mess of it
|
||||
when it's coded in C. */
|
||||
|
||||
#include <sysdep.h>
|
||||
|
||||
ENTRY(__copysign)
|
||||
/* double [f1] copysign (double [f1] x, double [f2] y);
|
||||
copysign(x,y) returns a value with the magnitude of x and
|
||||
with the sign bit of y. */
|
||||
stwu r1,-16(r1)
|
||||
stfd fp2,8(r1)
|
||||
lwz r3,8(r1)
|
||||
cmpwi r3,0
|
||||
addi r1,r1,16
|
||||
blt L(0)
|
||||
fabs fp1,fp1
|
||||
blr
|
||||
L(0): fnabs fp1,fp1
|
||||
blr
|
||||
END (__copysign)
|
||||
|
||||
weak_alias(__copysign,copysign)
|
||||
|
||||
/* It turns out that it's safe to use this code even for single-precision. */
|
||||
weak_alias(__copysign,copysignf)
|
||||
strong_alias(__copysign,__copysignf)
|
||||
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
weak_alias(__copysign,copysignl)
|
||||
strong_alias(__copysign,__copysignl)
|
||||
#endif
|
||||
@@ -0,0 +1 @@
|
||||
/* __copysignf is in s_copysign.S */
|
||||
@@ -0,0 +1,37 @@
|
||||
/* Floating-point absolute value. PowerPC version.
|
||||
Copyright (C) 1997, 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. */
|
||||
|
||||
#include <sysdep.h>
|
||||
|
||||
ENTRY(__fabs)
|
||||
/* double [f1] fabs (double [f1] x); */
|
||||
fabs fp1,fp1
|
||||
blr
|
||||
END(__fabs)
|
||||
|
||||
weak_alias(__fabs,fabs)
|
||||
|
||||
/* It turns out that it's safe to use this code even for single-precision. */
|
||||
strong_alias(__fabs,__fabsf)
|
||||
weak_alias(__fabs,fabsf)
|
||||
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
weak_alias(__fabs,__fabsl)
|
||||
weak_alias(__fabs,fabsl)
|
||||
#endif
|
||||
@@ -0,0 +1 @@
|
||||
/* __fabsf is in s_fabs.S */
|
||||
@@ -0,0 +1,31 @@
|
||||
/* Return positive difference between arguments.
|
||||
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>
|
||||
|
||||
double
|
||||
__fdim (double x, double y)
|
||||
{
|
||||
return x < y ? 0 : x - y;
|
||||
}
|
||||
weak_alias (__fdim, fdim)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__fdim, __fdiml)
|
||||
weak_alias (__fdim, fdiml)
|
||||
#endif
|
||||
@@ -0,0 +1,27 @@
|
||||
/* Return positive difference between arguments.
|
||||
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>
|
||||
|
||||
float
|
||||
__fdimf (float x, float y)
|
||||
{
|
||||
return x < y ? 0 : x - y;
|
||||
}
|
||||
weak_alias (__fdimf, fdimf)
|
||||
@@ -0,0 +1,47 @@
|
||||
/* Round a 64-bit floating point value to the nearest integer.
|
||||
Copyright (C) 1997, 2000 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"
|
||||
|
||||
double
|
||||
__rint (double x)
|
||||
{
|
||||
static const float TWO52 = 4503599627370496.0;
|
||||
|
||||
if (fabs (x) < TWO52)
|
||||
{
|
||||
if (x > 0.0)
|
||||
{
|
||||
x += TWO52;
|
||||
x -= TWO52;
|
||||
}
|
||||
else if (x < 0.0)
|
||||
{
|
||||
x = TWO52 - x;
|
||||
x = -(x - TWO52);
|
||||
}
|
||||
}
|
||||
|
||||
return x;
|
||||
}
|
||||
weak_alias (__rint, rint)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__rint, __rintl)
|
||||
weak_alias (__rint, rintl)
|
||||
#endif
|
||||
@@ -0,0 +1,43 @@
|
||||
/* Round a 32-bit floating point value to the nearest integer.
|
||||
Copyright (C) 1997, 2000 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"
|
||||
|
||||
float
|
||||
__rintf (float x)
|
||||
{
|
||||
static const float TWO23 = 8388608.0;
|
||||
|
||||
if (fabsf (x) < TWO23)
|
||||
{
|
||||
if (x > 0.0)
|
||||
{
|
||||
x += TWO23;
|
||||
x -= TWO23;
|
||||
}
|
||||
else if (x < 0.0)
|
||||
{
|
||||
x = TWO23 - x;
|
||||
x = -(x - TWO23);
|
||||
}
|
||||
}
|
||||
|
||||
return x;
|
||||
}
|
||||
weak_alias (__rintf, rintf)
|
||||
Reference in New Issue
Block a user