3 * IBM Accurate Mathematical Library
4 * Copyright (c) International Business Machines Corp., 2001
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 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 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, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 /****************************************************************************/
21 /* MODULE_NAME:mpsqrt.c */
26 /* FILES NEEDED:endian.h mpa.h mpsqrt.h */
28 /* Multi-Precision square root function subroutine for precision p >= 4. */
29 /* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */
31 /****************************************************************************/
35 /****************************************************************************/
36 /* Multi-Precision square root function subroutine for precision p >= 4. */
37 /* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */
38 /* Routine receives two pointers to Multi Precision numbers: */
39 /* x (left argument) and y (next argument). Routine also receives precision */
40 /* p as integer. Routine computes sqrt(*x) and stores result in *y */
41 /****************************************************************************/
43 double fastiroot(double);
45 void mpsqrt(mp_no *x, mp_no *y, int p) {
51 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,
52 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,
53 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}},
54 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,
55 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,
56 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
57 mp_no mpxn,mpz,mpu,mpt1,mpt2;
59 /* Prepare multi-precision 1/2 and 3/2 */
60 mphalf.e =0; mphalf.d[0] =ONE; mphalf.d[1] =HALFRAD;
61 mp3halfs.e=1; mp3halfs.d[0]=ONE; mp3halfs.d[1]=ONE; mp3halfs.d[2]=HALFRAD;
63 ex=EX; ey=EX/2; cpy(x,&mpxn,p); mpxn.e -= (ey+ey);
64 __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); dbl_mp(dy,&mpu,p);
65 mul(&mpxn,&mphalf,&mpz,p);
69 mul(&mpu,&mpu,&mpt1,p);
70 mul(&mpt1,&mpz,&mpt2,p);
71 sub(&mp3halfs,&mpt2,&mpt1,p);
72 mul(&mpu,&mpt1,&mpt2,p);
75 mul(&mpxn,&mpu,y,p); EY += ey;
80 /***********************************************************/
81 /* Compute a double precision approximation for 1/sqrt(x) */
82 /* with the relative error bounded by 2**-51. */
83 /***********************************************************/
84 double fastiroot(double x) {
85 union {long i[2]; double d;} p,q;
88 static const double c0 = 0.99674, c1 = -0.53380, c2 = 0.45472, c3 = -0.21553;
91 p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF ) | 0x3FE00000 ;
95 n = (q.i[HIGH_HALF] - p.i[HIGH_HALF])>>1;
96 z = ((c3*z + c2)*z + c1)*z + c0; /* 2**-7 */
97 z = z*(1.5 - 0.5*y*z*z); /* 2**-14 */
98 p.d = z*(1.5 - 0.5*y*z*z); /* 2**-28 */
101 return p.d*(1.5 - 0.5*p.d*t);