diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/mpsqrt.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpsqrt.c | 102 |
1 files changed, 102 insertions, 0 deletions
diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c new file mode 100644 index 0000000000..6f73fae137 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/mpsqrt.c @@ -0,0 +1,102 @@ + +/* + * 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); +} |