1 /* $OpenBSD: fpu_sqrt.c,v 1.7 2024/03/29 21:02:11 miod Exp $ */
2
3 /*
4 * Copyright (c) 1992, 1993
5 * The Regents of the University of California. All rights reserved.
6 *
7 * This software was developed by the Computer Systems Engineering group
8 * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
9 * contributed to Berkeley.
10 *
11 * All advertising materials mentioning features or use of this software
12 * must display the following acknowledgement:
13 * This product includes software developed by the University of
14 * California, Lawrence Berkeley Laboratory.
15 *
16 * Redistribution and use in source and binary forms, with or without
17 * modification, are permitted provided that the following conditions
18 * are met:
19 * 1. Redistributions of source code must retain the above copyright
20 * notice, this list of conditions and the following disclaimer.
21 * 2. Redistributions in binary form must reproduce the above copyright
22 * notice, this list of conditions and the following disclaimer in the
23 * documentation and/or other materials provided with the distribution.
24 * 3. All advertising materials mentioning features or use of this software
25 * must display the following acknowledgement:
26 * This product includes software developed by the University of
27 * California, Berkeley and its contributors.
28 * 4. Neither the name of the University nor the names of its contributors
29 * may be used to endorse or promote products derived from this software
30 * without specific prior written permission.
31 *
32 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
33 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
34 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
35 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
36 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
37 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
38 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
39 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
40 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
41 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
42 * SUCH DAMAGE.
43 *
44 * @(#)fpu_sqrt.c 8.1 (Berkeley) 6/11/93
45 * $NetBSD: fpu_sqrt.c,v 1.2 1994/11/20 20:52:46 deraadt Exp $
46 */
47
48 /*
49 * Perform an FPU square root (return sqrt(x)).
50 */
51
52 #include <sys/types.h>
53
54 #include "fpu_arith.h"
55 #include "fpu_emu.h"
56 #include "fpu_extern.h"
57
58 /*
59 * Our task is to calculate the square root of a floating point number x0.
60 * This number x normally has the form:
61 *
62 * exp
63 * x = mant * 2 (where 1 <= mant < 2 and exp is an integer)
64 *
65 * This can be left as it stands, or the mantissa can be doubled and the
66 * exponent decremented:
67 *
68 * exp-1
69 * x = (2 * mant) * 2 (where 2 <= 2 * mant < 4)
70 *
71 * If the exponent `exp' is even, the square root of the number is best
72 * handled using the first form, and is by definition equal to:
73 *
74 * exp/2
75 * sqrt(x) = sqrt(mant) * 2
76 *
77 * If exp is odd, on the other hand, it is convenient to use the second
78 * form, giving:
79 *
80 * (exp-1)/2
81 * sqrt(x) = sqrt(2 * mant) * 2
82 *
83 * In the first case, we have
84 *
85 * 1 <= mant < 2
86 *
87 * and therefore
88 *
89 * sqrt(1) <= sqrt(mant) < sqrt(2)
90 *
91 * while in the second case we have
92 *
93 * 2 <= 2*mant < 4
94 *
95 * and therefore
96 *
97 * sqrt(2) <= sqrt(2*mant) < sqrt(4)
98 *
99 * so that in any case, we are sure that
100 *
101 * sqrt(1) <= sqrt(n * mant) < sqrt(4), n = 1 or 2
102 *
103 * or
104 *
105 * 1 <= sqrt(n * mant) < 2, n = 1 or 2.
106 *
107 * This root is therefore a properly formed mantissa for a floating
108 * point number. The exponent of sqrt(x) is either exp/2 or (exp-1)/2
109 * as above. This leaves us with the problem of finding the square root
110 * of a fixed-point number in the range [1..4).
111 *
112 * Though it may not be instantly obvious, the following square root
113 * algorithm works for any integer x of an even number of bits, provided
114 * that no overflows occur:
115 *
116 * let q = 0
117 * for k = NBITS-1 to 0 step -1 do -- for each digit in the answer...
118 * x *= 2 -- multiply by radix, for next digit
119 * if x >= 2q + 2^k then -- if adding 2^k does not
120 * x -= 2q + 2^k -- exceed the correct root,
121 * q += 2^k -- add 2^k and adjust x
122 * fi
123 * done
124 * sqrt = q / 2^(NBITS/2) -- (and any remainder is in x)
125 *
126 * If NBITS is odd (so that k is initially even), we can just add another
127 * zero bit at the top of x. Doing so means that q is not going to acquire
128 * a 1 bit in the first trip around the loop (since x0 < 2^NBITS). If the
129 * final value in x is not needed, or can be off by a factor of 2, this is
130 * equivalent to moving the `x *= 2' step to the bottom of the loop:
131 *
132 * for k = NBITS-1 to 0 step -1 do if ... fi; x *= 2; done
133 *
134 * and the result q will then be sqrt(x0) * 2^floor(NBITS / 2).
135 * (Since the algorithm is destructive on x, we will call x's initial
136 * value, for which q is some power of two times its square root, x0.)
137 *
138 * If we insert a loop invariant y = 2q, we can then rewrite this using
139 * C notation as:
140 *
141 * q = y = 0; x = x0;
142 * for (k = NBITS; --k >= 0;) {
143 * #if (NBITS is even)
144 * x *= 2;
145 * #endif
146 * t = y + (1 << k);
147 * if (x >= t) {
148 * x -= t;
149 * q += 1 << k;
150 * y += 1 << (k + 1);
151 * }
152 * #if (NBITS is odd)
153 * x *= 2;
154 * #endif
155 * }
156 *
157 * If x0 is fixed point, rather than an integer, we can simply alter the
158 * scale factor between q and sqrt(x0). As it happens, we can easily arrange
159 * for the scale factor to be 2**0 or 1, so that sqrt(x0) == q.
160 *
161 * In our case, however, x0 (and therefore x, y, q, and t) are multiword
162 * integers, which adds some complication. But note that q is built one
163 * bit at a time, from the top down, and is not used itself in the loop
164 * (we use 2q as held in y instead). This means we can build our answer
165 * in an integer, one word at a time, which saves a bit of work. Also,
166 * since 1 << k is always a `new' bit in q, 1 << k and 1 << (k+1) are
167 * `new' bits in y and we can set them with an `or' operation rather than
168 * a full-blown multiword add.
169 *
170 * We are almost done, except for one snag. We must prove that none of our
171 * intermediate calculations can overflow. We know that x0 is in [1..4)
172 * and therefore the square root in q will be in [1..2), but what about x,
173 * y, and t?
174 *
175 * We know that y = 2q at the beginning of each loop. (The relation only
176 * fails temporarily while y and q are being updated.) Since q < 2, y < 4.
177 * The sum in t can, in our case, be as much as y+(1<<1) = y+2 < 6, and.
178 * Furthermore, we can prove with a bit of work that x never exceeds y by
179 * more than 2, so that even after doubling, 0 <= x < 8. (This is left as
180 * an exercise to the reader, mostly because I have become tired of working
181 * on this comment.)
182 *
183 * If our floating point mantissas (which are of the form 1.frac) occupy
184 * B+1 bits, our largest intermediary needs at most B+3 bits, or two extra.
185 * In fact, we want even one more bit (for a carry, to avoid compares), or
186 * three extra. There is a comment in fpu_emu.h reminding maintainers of
187 * this, so we have some justification in assuming it.
188 */
189 struct fpn *
__fpu_sqrt(fe)190 __fpu_sqrt(fe)
191 struct fpemu *fe;
192 {
193 struct fpn *x = &fe->fe_f1;
194 u_int bit, q, tt;
195 u_int x0, x1, x2, x3;
196 u_int y0, y1, y2, y3;
197 u_int d0, d1, d2, d3;
198 int e;
199
200 /*
201 * Take care of special cases first. In order:
202 *
203 * sqrt(NaN) = NaN
204 * sqrt(+0) = +0
205 * sqrt(-0) = -0
206 * sqrt(x < 0) = NaN (including sqrt(-Inf))
207 * sqrt(+Inf) = +Inf
208 *
209 * Then all that remains are numbers with mantissas in [1..2).
210 */
211 if (ISNAN(x) || ISZERO(x))
212 return (x);
213 if (x->fp_sign)
214 return (__fpu_newnan(fe));
215 if (ISINF(x))
216 return (x);
217
218 /*
219 * Calculate result exponent. As noted above, this may involve
220 * doubling the mantissa. We will also need to double x each
221 * time around the loop, so we define a macro for this here, and
222 * we break out the multiword mantissa.
223 */
224 #ifdef FPU_SHL1_BY_ADD
225 #define DOUBLE_X { \
226 FPU_ADDS(x3, x3, x3); FPU_ADDCS(x2, x2, x2); \
227 FPU_ADDCS(x1, x1, x1); FPU_ADDC(x0, x0, x0); \
228 }
229 #else
230 #define DOUBLE_X { \
231 x0 = (x0 << 1) | (x1 >> 31); x1 = (x1 << 1) | (x2 >> 31); \
232 x2 = (x2 << 1) | (x3 >> 31); x3 <<= 1; \
233 }
234 #endif
235 #if (FP_NMANT & 1) != 0
236 # define ODD_DOUBLE DOUBLE_X
237 # define EVEN_DOUBLE /* nothing */
238 #else
239 # define ODD_DOUBLE /* nothing */
240 # define EVEN_DOUBLE DOUBLE_X
241 #endif
242 x0 = x->fp_mant[0];
243 x1 = x->fp_mant[1];
244 x2 = x->fp_mant[2];
245 x3 = x->fp_mant[3];
246 e = x->fp_exp;
247 if (e & 1) /* exponent is odd; use sqrt(2mant) */
248 DOUBLE_X;
249 /* THE FOLLOWING ASSUMES THAT RIGHT SHIFT DOES SIGN EXTENSION */
250 x->fp_exp = e >> 1; /* calculates (e&1 ? (e-1)/2 : e/2 */
251
252 /*
253 * Now calculate the mantissa root. Since x is now in [1..4),
254 * we know that the first trip around the loop will definitely
255 * set the top bit in q, so we can do that manually and start
256 * the loop at the next bit down instead. We must be sure to
257 * double x correctly while doing the `known q=1.0'.
258 *
259 * We do this one mantissa-word at a time, as noted above, to
260 * save work. To avoid `(1U << 31) << 1', we also do the top bit
261 * outside of each per-word loop.
262 *
263 * The calculation `t = y + bit' breaks down into `t0 = y0, ...,
264 * t3 = y3, t? |= bit' for the appropriate word. Since the bit
265 * is always a `new' one, this means that three of the `t?'s are
266 * just the corresponding `y?'; we use `#define's here for this.
267 * The variable `tt' holds the actual `t?' variable.
268 */
269
270 /* calculate q0 */
271 #define t0 tt
272 bit = FP_1;
273 EVEN_DOUBLE;
274 /* if (x >= (t0 = y0 | bit)) { */ /* always true */
275 q = bit;
276 x0 -= bit;
277 y0 = bit << 1;
278 /* } */
279 ODD_DOUBLE;
280 while ((bit >>= 1) != 0) { /* for remaining bits in q0 */
281 EVEN_DOUBLE;
282 t0 = y0 | bit; /* t = y + bit */
283 if (x0 >= t0) { /* if x >= t then */
284 x0 -= t0; /* x -= t */
285 q |= bit; /* q += bit */
286 y0 |= bit << 1; /* y += bit << 1 */
287 }
288 ODD_DOUBLE;
289 }
290 x->fp_mant[0] = q;
291 #undef t0
292
293 /* calculate q1. note (y0&1)==0. */
294 #define t0 y0
295 #define t1 tt
296 q = 0;
297 y1 = 0;
298 bit = 1U << 31;
299 EVEN_DOUBLE;
300 t1 = bit;
301 FPU_SUBS(d1, x1, t1);
302 FPU_SUBC(d0, x0, t0); /* d = x - t */
303 if ((int)d0 >= 0) { /* if d >= 0 (i.e., x >= t) then */
304 x0 = d0, x1 = d1; /* x -= t */
305 q = bit; /* q += bit */
306 y0 |= 1; /* y += bit << 1 */
307 }
308 ODD_DOUBLE;
309 while ((bit >>= 1) != 0) { /* for remaining bits in q1 */
310 EVEN_DOUBLE; /* as before */
311 t1 = y1 | bit;
312 FPU_SUBS(d1, x1, t1);
313 FPU_SUBC(d0, x0, t0);
314 if ((int)d0 >= 0) {
315 x0 = d0, x1 = d1;
316 q |= bit;
317 y1 |= bit << 1;
318 }
319 ODD_DOUBLE;
320 }
321 x->fp_mant[1] = q;
322 #undef t1
323
324 /* calculate q2. note (y1&1)==0; y0 (aka t0) is fixed. */
325 #define t1 y1
326 #define t2 tt
327 q = 0;
328 y2 = 0;
329 bit = 1U << 31;
330 EVEN_DOUBLE;
331 t2 = bit;
332 FPU_SUBS(d2, x2, t2);
333 FPU_SUBCS(d1, x1, t1);
334 FPU_SUBC(d0, x0, t0);
335 if ((int)d0 >= 0) {
336 x0 = d0, x1 = d1, x2 = d2;
337 q |= bit;
338 y1 |= 1; /* now t1, y1 are set in concrete */
339 }
340 ODD_DOUBLE;
341 while ((bit >>= 1) != 0) {
342 EVEN_DOUBLE;
343 t2 = y2 | bit;
344 FPU_SUBS(d2, x2, t2);
345 FPU_SUBCS(d1, x1, t1);
346 FPU_SUBC(d0, x0, t0);
347 if ((int)d0 >= 0) {
348 x0 = d0, x1 = d1, x2 = d2;
349 q |= bit;
350 y2 |= bit << 1;
351 }
352 ODD_DOUBLE;
353 }
354 x->fp_mant[2] = q;
355 #undef t2
356
357 /* calculate q3. y0, t0, y1, t1 all fixed; y2, t2, almost done. */
358 #define t2 y2
359 #define t3 tt
360 q = 0;
361 y3 = 0;
362 bit = 1U << 31;
363 EVEN_DOUBLE;
364 t3 = bit;
365 FPU_SUBS(d3, x3, t3);
366 FPU_SUBCS(d2, x2, t2);
367 FPU_SUBCS(d1, x1, t1);
368 FPU_SUBC(d0, x0, t0);
369 if ((int)d0 >= 0) {
370 x0 = d0, x1 = d1, x2 = d2, x3 = d3;
371 q |= bit;
372 y2 |= 1;
373 }
374 ODD_DOUBLE;
375 while ((bit >>= 1) != 0) {
376 EVEN_DOUBLE;
377 t3 = y3 | bit;
378 FPU_SUBS(d3, x3, t3);
379 FPU_SUBCS(d2, x2, t2);
380 FPU_SUBCS(d1, x1, t1);
381 FPU_SUBC(d0, x0, t0);
382 if ((int)d0 >= 0) {
383 x0 = d0, x1 = d1, x2 = d2, x3 = d3;
384 q |= bit;
385 y3 |= bit << 1;
386 }
387 ODD_DOUBLE;
388 }
389 x->fp_mant[3] = q;
390
391 /*
392 * The result, which includes guard and round bits, is exact iff
393 * x is now zero; any nonzero bits in x represent sticky bits.
394 */
395 x->fp_sticky = x0 | x1 | x2 | x3;
396 return (x);
397 }
398