rainlance/amath
1
0
Fork 0
mirror of https://gitlab.com/rnger/amath synced 2025-10-06 02:49:59 +00:00
Code Releases Wiki Activity
Files
7a63687d4b578a61d65e97df4cea881b6ae6a08d
amath/doc/man/man3/sqrt.c.3
Carsten Larsen 9e2c66af53 Update documentation
2017-01-24 22:03:15 +01:00

282 lines
6.8 KiB
Groff
Raw Blame History

.TH "lib/real/sqrt.c" 3 "Tue Jan 24 2017" "Version 1.6.2" "amath" \" -*- nroff -*-
.ad l
.nh
.SH NAME
lib/real/sqrt.c \-
.SH SYNOPSIS
.br
.PP
\fC#include 'prim\&.h'\fP
.br
\fC#include 'math\&.h'\fP
.br
.SS "Functions"
.in +1c
.ti -1c
.RI "double \fBsqrt\fP (double x)"
.br
.RI "\fISquare root function\&. \fP"
.in -1c
.SS "Variables"
.in +1c
.ti -1c
.RI "static const double \fBone\fP = 1\&.0"
.br
.ti -1c
.RI "static const double \fBtiny\fP = 1\&.0e\-300"
.br
.in -1c
.SH "Function Documentation"
.PP
.SS "double sqrt (double x)"
.PP
Square root function\&.
.PP
\fBVersion:\fP
.RS 4
1\&.3
.RE
.PP
\fBDate:\fP
.RS 4
95/01/18
.RE
.PP
.PP
.nf
Return correctly rounded sqrt\&.
------------------------------------------
| Use the hardware sqrt if you have one |
------------------------------------------
Method:
Bit by bit method using integer arithmetic\&. (Slow, but portable)
1\&. Normalization
Scale x to y in [1,4) with even powers of 2:
find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
sqrt(x) = 2^k * sqrt(y)
2\&. Bit by bit computation
Let q = sqrt(y) truncated to i bit after binary point (q = 1),
i 0
i+1 2
s = 2*q , and y = 2 * ( y - q )\&. (1)
i i i i
.fi
.PP
.PP
.PP
.nf
To compute q from q , one checks whether
i+1 i
.fi
.PP
.PP
.PP
.nf
-(i+1) 2
(q + 2 ) <= y\&. (2)
i
-(i+1)
If (2) is false, then q = q ; otherwise q = q + 2 \&.
i+1 i i+1 i
.fi
.PP
.PP
.PP
.nf
With some algebric manipulation, it is not difficult to see
that (2) is equivalent to
-(i+1)
s + 2 <= y (3)
i i
.fi
.PP
.PP
.PP
.nf
The advantage of (3) is that s and y can be computed by
i i
the following recurrence formula:
if (3) is false
.fi
.PP
.PP
.PP
.nf
s = s , y = y ; (4)
i+1 i i+1 i
.fi
.PP
.PP
.PP
.nf
otherwise,
-i -(i+1)
s = s + 2 , y = y - s - 2 (5)
i+1 i i+1 i i
.fi
.PP
.PP
.PP
.nf
One may easily use induction to prove (4) and (5)\&.
Note\&. Since the left hand side of (3) contain only i+2 bits,
it does not necessary to do a full (53-bit) comparison
in (3)\&.
3\&. Final rounding
After generating the 53 bits result, we compute one more bit\&.
Together with the remainder, we can decide whether the
result is exact, bigger than 1/2ulp, or less than 1/2ulp
(it will never equal to 1/2ulp)\&.
The rounding mode can be detected by checking whether
huge + tiny is equal to huge, and whether huge - tiny is
equal to huge for some floating point number 'huge' and 'tiny'\&.
.fi
.PP
.PP
.PP
.nf
Special cases:
sqrt(+-0) = +-0 \&.\&.\&. exact
sqrt(inf) = inf
sqrt(-ve) = NaN \&.\&.\&. with invalid signal
sqrt(NaN) = NaN \&.\&.\&. with invalid signal for signaling NaN
.fi
.PP
Other methods : see the \fBsquareroot\fP\&.
.PP
\fBCopyright:\fP
.RS 4
Copyright (C) 1993 by Sun Microsystems, Inc\&. All rights reserved\&. Developed at SunSoft, a Sun Microsystems, Inc\&. business\&. Permission to use, copy, modify, and distribute this software is freely granted, provided that this notice is preserved\&.
.RE
.PP
.PP
Definition at line 127 of file sqrt\&.c\&.
.PP
References one, and tiny\&.
.PP
Referenced by acos(), acosh(), asin(), asinh(), csqrt(), hypot(), pow(), and RealNumber::SquareRoot()\&.
.PP
.nf
128 {
129 double z;
130 sword sign = (int)0x80000000;
131 uword r,t1,s1,ix1,q1;
132 sword ix0,s0,q,m,t,i;
133
134 EXTRACT_WORDS(ix0,ix1,x);
135
136 /* take care of Inf and NaN */
137 if((ix0&0x7ff00000)==0x7ff00000) {
138 return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
139 sqrt(-inf)=sNaN */
140 }
141 /* take care of zero */
142 if(ix0<=0) {
143 if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
144 else if(ix0<0)
145 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
146 }
147 /* normalize x */
148 m = (ix0>>20);
149 if(m==0) { /* subnormal x */
150 while(ix0==0) {
151 m -= 21;
152 ix0 |= (ix1>>11);
153 ix1 <<= 21;
154 }
155 for(i=0; (ix0&0x00100000)==0; i++) ix0<<=1;
156 m -= i-1;
157 ix0 |= (ix1>>(32-i));
158 ix1 <<= i;
159 }
160 m -= 1023; /* unbias exponent */
161 ix0 = (ix0&0x000fffff)|0x00100000;
162 if(m&1) { /* odd m, double x to make it even */
163 ix0 += ix0 + ((ix1&sign)>>31);
164 ix1 += ix1;
165 }
166 m >>= 1; /* m = [m/2] */
167
168 /* generate sqrt(x) bit by bit */
169 ix0 += ix0 + ((ix1&sign)>>31);
170 ix1 += ix1;
171 q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
172 r = 0x00200000; /* r = moving bit from right to left */
173
174 while(r!=0) {
175 t = s0+r;
176 if(t<=ix0) {
177 s0 = t+r;
178 ix0 -= t;
179 q += r;
180 }
181 ix0 += ix0 + ((ix1&sign)>>31);
182 ix1 += ix1;
183 r>>=1;
184 }
185
186 r = sign;
187 while(r!=0) {
188 t1 = s1+r;
189 t = s0;
190 if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
191 s1 = t1+r;
192 if(((t1&sign)==(uword)sign)&&(s1&sign)==0) s0 += 1;
193 ix0 -= t;
194 if (ix1 < t1) ix0 -= 1;
195 ix1 -= t1;
196 q1 += r;
197 }
198 ix0 += ix0 + ((ix1&sign)>>31);
199 ix1 += ix1;
200 r>>=1;
201 }
202
203 /* use floating add to find out rounding direction */
204 if((ix0|ix1)!=0) {
205 z = one-tiny; /* trigger inexact flag */
206 if (z>=one) {
207 z = one+tiny;
208 if (q1==(uword)0xffffffff) {
209 q1=0;
210 q += 1;
211 }
212 else if (z>one) {
213 if (q1==(uword)0xfffffffe) q+=1;
214 q1+=2;
215 } else
216 q1 += (q1&1);
217 }
218 }
219 ix0 = (q>>1)+0x3fe00000;
220 ix1 = q1>>1;
221 if ((q&1)==1) ix1 |= sign;
222 ix0 += (m <<20);
223 INSERT_WORDS(z,ix0,ix1);
224 return z;
225 }
.fi
.SH "Variable Documentation"
.PP
.SS "const double one = 1\&.0\fC [static]\fP"
.PP
Definition at line 47 of file sqrt\&.c\&.
.PP
Referenced by sqrt()\&.
.SS "const double tiny = 1\&.0e\-300\fC [static]\fP"
.PP
Definition at line 47 of file sqrt\&.c\&.
.PP
Referenced by sqrt()\&.
.SH "Author"
.PP
Generated automatically by Doxygen for amath from the source code\&.
View Git Blame Copy Permalink