]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/mpsqrt.c
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2019 Free Software Foundation, Inc.
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <https://www.gnu.org/licenses/>.
19 /****************************************************************************/
20 /* MODULE_NAME:mpsqrt.c */
25 /* FILES NEEDED:endian.h mpa.h mpsqrt.h */
27 /* Multi-Precision square root function subroutine for precision p >= 4. */
28 /* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */
30 /****************************************************************************/
40 /****************************************************************************/
41 /* Multi-Precision square root function subroutine for precision p >= 4. */
42 /* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */
43 /* Routine receives two pointers to Multi Precision numbers: */
44 /* x (left argument) and y (next argument). Routine also receives precision */
45 /* p as integer. Routine computes sqrt(*x) and stores result in *y */
46 /****************************************************************************/
48 static double fastiroot (double);
52 __mpsqrt (mp_no
*x
, mp_no
*y
, int p
)
56 static const mp_no mphalf
= {0, {1.0, HALFRAD
}};
57 static const mp_no mp3halfs
= {1, {1.0, 1.0, HALFRAD
}};
58 mp_no mpxn
, mpz
, mpu
, mpt1
, mpt2
;
63 __mp_dbl (&mpxn
, &dx
, p
);
65 __dbl_mp (dy
, &mpu
, p
);
66 __mul (&mpxn
, &mphalf
, &mpz
, p
);
69 for (i
= 0; i
< m
; i
++)
71 __sqr (&mpu
, &mpt1
, p
);
72 __mul (&mpt1
, &mpz
, &mpt2
, p
);
73 __sub (&mp3halfs
, &mpt2
, &mpt1
, p
);
74 __mul (&mpu
, &mpt1
, &mpt2
, p
);
75 __cpy (&mpt2
, &mpu
, p
);
77 __mul (&mpxn
, &mpu
, y
, p
);
81 /***********************************************************/
82 /* Compute a double precision approximation for 1/sqrt(x) */
83 /* with the relative error bounded by 2**-51. */
84 /***********************************************************/
96 static const double c0
= 0.99674, c1
= -0.53380;
97 static const double c2
= 0.45472, c3
= -0.21553;
100 p
.i
[HIGH_HALF
] = (p
.i
[HIGH_HALF
] & 0x3FFFFFFF) | 0x3FE00000;
104 n
= (q
.i
[HIGH_HALF
] - p
.i
[HIGH_HALF
]) >> 1;
105 z
= ((c3
* z
+ c2
) * z
+ c1
) * z
+ c0
; /* 2**-7 */
106 z
= z
* (1.5 - 0.5 * y
* z
* z
); /* 2**-14 */
107 p
.d
= z
* (1.5 - 0.5 * y
* z
* z
); /* 2**-28 */
110 return p
.d
* (1.5 - 0.5 * p
.d
* t
);