]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/mpsqrt.c
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
);