lib/libc/sparc64/fpu/fpu_sqrt.c

   1 /*
   2  * Copyright (c) 1992, 1993
   3  *      The Regents of the University of California.  All rights reserved.
   4  *
   5  * This software was developed by the Computer Systems Engineering group
   6  * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
   7  * contributed to Berkeley.
   8  *
   9  * All advertising materials mentioning features or use of this software
  10  * must display the following acknowledgement:
  11  *      This product includes software developed by the University of
  12  *      California, Lawrence Berkeley Laboratory.
  13  *
  14  * Redistribution and use in source and binary forms, with or without
  15  * modification, are permitted provided that the following conditions
  16  * are met:
  17  * 1. Redistributions of source code must retain the above copyright
  18  *    notice, this list of conditions and the following disclaimer.
  19  * 2. Redistributions in binary form must reproduce the above copyright
  20  *    notice, this list of conditions and the following disclaimer in the
  21  *    documentation and/or other materials provided with the distribution.
  22  * 4. Neither the name of the University nor the names of its contributors
  23  *    may be used to endorse or promote products derived from this software
  24  *    without specific prior written permission.
  25  *
  26  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  27  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  28  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  29  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  30  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  31  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  32  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  33  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  34  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  35  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  36  * SUCH DAMAGE.
  37  *
  38  *      @(#)fpu_sqrt.c  8.1 (Berkeley) 6/11/93
  39  *      $NetBSD: fpu_sqrt.c,v 1.2 1994/11/20 20:52:46 deraadt Exp $
  40  */
  41
  42 #include <sys/cdefs.h>
  43 __FBSDID("$FreeBSD$");
  44
  45 /*
  46  * Perform an FPU square root (return sqrt(x)).
  47  */
  48
  49 #include <sys/types.h>
  50
  51 #include <machine/frame.h>
  52 #include <machine/fp.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  * equivalant 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 *
 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 `(1 << 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 = 1 << 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 = 1 << 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 = 1 << 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 }