]>
Commit | Line | Data |
---|---|---|
f7eac6eb | 1 | /* |
e4d82761 | 2 | * IBM Accurate Mathematical Library |
aeb25823 | 3 | * written by International Business Machines Corp. |
0ac5ae23 | 4 | * Copyright (C) 2001, 2011 Free Software Foundation |
f7eac6eb | 5 | * |
e4d82761 UD |
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 | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
6d52618b | 10 | * |
e4d82761 UD |
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 | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
6d52618b | 15 | * |
e4d82761 UD |
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. | |
f7eac6eb | 19 | */ |
e4d82761 UD |
20 | /*********************************************************************/ |
21 | /* MODULE_NAME: uroot.c */ | |
22 | /* */ | |
23 | /* FUNCTION: usqrt */ | |
24 | /* */ | |
25 | /* FILES NEEDED: dla.h endian.h mydefs.h uroot.h */ | |
26 | /* uroot.tbl */ | |
27 | /* */ | |
28 | /* An ultimate sqrt routine. Given an IEEE double machine number x */ | |
29 | /* it computes the correctly rounded (to nearest) value of square */ | |
30 | /* root of x. */ | |
31 | /* Assumption: Machine arithmetic operations are performed in */ | |
32 | /* round to nearest mode of IEEE 754 standard. */ | |
33 | /* */ | |
34 | /*********************************************************************/ | |
35 | ||
36 | #include "endian.h" | |
37 | #include "mydefs.h" | |
38 | #include "dla.h" | |
39 | #include "MathLib.h" | |
40 | #include "root.tbl" | |
e859d1d9 | 41 | #include "math_private.h" |
e4d82761 UD |
42 | |
43 | /*********************************************************************/ | |
7c370086 | 44 | /* An ultimate sqrt routine. Given an IEEE double machine number x */ |
e4d82761 UD |
45 | /* it computes the correctly rounded (to nearest) value of square */ |
46 | /* root of x. */ | |
47 | /*********************************************************************/ | |
48 | double __ieee754_sqrt(double x) { | |
49 | #include "uroot.h" | |
50 | static const double | |
51 | rt0 = 9.99999999859990725855365213134618E-01, | |
52 | rt1 = 4.99999999495955425917856814202739E-01, | |
53 | rt2 = 3.75017500867345182581453026130850E-01, | |
54 | rt3 = 3.12523626554518656309172508769531E-01; | |
7c370086 | 55 | static const double big = 134217728.0; |
e4d82761 | 56 | double y,t,del,res,res1,hy,z,zz,p,hx,tx,ty,s; |
7ea88fb4 AJ |
57 | mynumber a,c={{0,0}}; |
58 | int4 k; | |
e4d82761 UD |
59 | |
60 | a.x=x; | |
61 | k=a.i[HIGH_HALF]; | |
62 | a.i[HIGH_HALF]=(k&0x001fffff)|0x3fe00000; | |
63 | t=inroot[(k&0x001fffff)>>14]; | |
64 | s=a.x; | |
65 | /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/ | |
66 | if (k>0x000fffff && k<0x7ff00000) { | |
67 | y=1.0-t*(t*s); | |
68 | t=t*(rt0+y*(rt1+y*(rt2+y*rt3))); | |
69 | c.i[HIGH_HALF]=0x20000000+((k&0x7fe00000)>>1); | |
70 | y=t*s; | |
71 | hy=(y+big)-big; | |
72 | del=0.5*t*((s-hy*hy)-(y-hy)*(y+hy)); | |
73 | res=y+del; | |
74 | if (res == (res+1.002*((y-res)+del))) return res*c.x; | |
75 | else { | |
76 | res1=res+1.5*((y-res)+del); | |
77 | EMULV(res,res1,z,zz,p,hx,tx,hy,ty); /* (z+zz)=res*res1 */ | |
78 | return ((((z-s)+zz)<0)?max(res,res1):min(res,res1))*c.x; | |
79 | } | |
80 | } | |
81 | else { | |
7c370086 UD |
82 | if ((k & 0x7ff00000) == 0x7ff00000) |
83 | return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ | |
84 | if (x==0) return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */ | |
85 | if (k<0) return (x-x)/(x-x); /* sqrt(-ve)=sNaN */ | |
86 | return tm256.x*__ieee754_sqrt(x*t512.x); | |
e4d82761 | 87 | } |
f7eac6eb | 88 | } |
0ac5ae23 | 89 | strong_alias (__ieee754_sqrt, __sqrt_finite) |