mirror of
https://sourceware.org/git/glibc.git
synced 2024-12-13 14:50:17 +00:00
ca58f1dbeb
2001-03-12 Ulrich Drepper <drepper@redhat.com> * sysdeps/ieee754/dbl-64/e_remainder.c: Fix handling of boundary conditions. * sysdeps/ieee754/dbl-64/e_pow.c: Fix handling of boundary conditions. * sysdeps/ieee754/dbl-64/s_sin.c (__sin): Handle Inf and NaN correctly. (__cos): Likewise. * sysdeps/ieee754/dbl-64/e_asin.c (__ieee754_asin): Handle NaN correctly. (__ieee754_acos): Likewise. redefinition. * sysdeps/ieee754/dbl-64/endian.h: Define also one of BIG_ENDI and LITTLE_ENDI. * sysdeps/ieee754/dbl-64/MathLib.h (Init_Lib): Use void as parameter list.
103 lines
4.3 KiB
C
103 lines
4.3 KiB
C
|
|
/*
|
|
* IBM Accurate Mathematical Library
|
|
* Copyright (c) International Business Machines Corp., 2001
|
|
*
|
|
* 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 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 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);
|
|
}
|