diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/mpsqrt.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpsqrt.c | 111 |
1 files changed, 0 insertions, 111 deletions
diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c deleted file mode 100644 index be6d01eeef..0000000000 --- a/sysdeps/ieee754/dbl-64/mpsqrt.c +++ /dev/null @@ -1,111 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2017 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, see <http://www.gnu.org/licenses/>. - */ -/****************************************************************************/ -/* 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" - -#ifndef SECTION -# define SECTION -#endif - -#include "mpsqrt.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 */ -/****************************************************************************/ - -static double fastiroot (double); - -void -SECTION -__mpsqrt (mp_no *x, mp_no *y, int p) -{ - int i, m, ey; - double dx, dy; - static const mp_no mphalf = {0, {1.0, HALFRAD}}; - static const mp_no mp3halfs = {1, {1.0, 1.0, HALFRAD}}; - mp_no mpxn, mpz, mpu, mpt1, mpt2; - - 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 = __mpsqrt_mp[p]; - for (i = 0; i < m; i++) - { - __sqr (&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; -} - -/***********************************************************/ -/* Compute a double precision approximation for 1/sqrt(x) */ -/* with the relative error bounded by 2**-51. */ -/***********************************************************/ -static double -SECTION -fastiroot (double x) -{ - union - { - int i[2]; - double d; - } p, q; - double y, z, t; - int n; - static const double c0 = 0.99674, c1 = -0.53380; - static const double 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); -} |