sys/i386/gnu/fpemul/wm_sqrt.s

   1         .file   "wm_sqrt.S"
   2 /*
   3  *  wm_sqrt.S
   4  *
   5  * Fixed point arithmetic square root evaluation.
   6  *
   7  * Call from C as:
   8  *   void wm_sqrt(FPU_REG *n, unsigned int control_word)
   9  *
  10  *
  11  * Copyright (C) 1992,1993,1994
  12  *                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,
  13  *                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au
  14  * All rights reserved.
  15  *
  16  * This copyright notice covers the redistribution and use of the
  17  * FPU emulator developed by W. Metzenthen. It covers only its use
  18  * in the 386BSD, FreeBSD and NetBSD operating systems. Any other
  19  * use is not permitted under this copyright.
  20  *
  21  * Redistribution and use in source and binary forms, with or without
  22  * modification, are permitted provided that the following conditions
  23  * are met:
  24  * 1. Redistributions of source code must retain the above copyright
  25  *    notice, this list of conditions and the following disclaimer.
  26  * 2. Redistributions in binary form must include information specifying
  27  *    that source code for the emulator is freely available and include
  28  *    either:
  29  *      a) an offer to provide the source code for a nominal distribution
  30  *         fee, or
  31  *      b) list at least two alternative methods whereby the source
  32  *         can be obtained, e.g. a publically accessible bulletin board
  33  *         and an anonymous ftp site from which the software can be
  34  *         downloaded.
  35  * 3. All advertising materials specifically mentioning features or use of
  36  *    this emulator must acknowledge that it was developed by W. Metzenthen.
  37  * 4. The name of W. Metzenthen may not be used to endorse or promote
  38  *    products derived from this software without specific prior written
  39  *    permission.
  40  *
  41  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
  42  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
  43  * AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL
  44  * W. METZENTHEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  45  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  46  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  47  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  48  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  49  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  50  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  51  *
  52  *
  53  * The purpose of this copyright, based upon the Berkeley copyright, is to
  54  * ensure that the covered software remains freely available to everyone.
  55  *
  56  * The software (with necessary differences) is also available, but under
  57  * the terms of the GNU copyleft, for the Linux operating system and for
  58  * the djgpp ms-dos extender.
  59  *
  60  * W. Metzenthen   June 1994.
  61  *
  62  *
  63  * $FreeBSD: src/sys/gnu/i386/fpemul/wm_sqrt.s,v 1.9.2.1 2000/07/07 00:38:42 obrien Exp $
  64  *
  65  */
  66
  67
  68 /*---------------------------------------------------------------------------+
  69  |  wm_sqrt(FPU_REG *n, unsigned int control_word)                           |
  70  |    returns the square root of n in n.                                     |
  71  |                                                                           |
  72  |  Use Newton's method to compute the square root of a number, which must   |
  73  |  be in the range  [1.0 .. 4.0),  to 64 bits accuracy.                     |
  74  |  Does not check the sign or tag of the argument.                          |
  75  |  Sets the exponent, but not the sign or tag of the result.                |
  76  |                                                                           |
  77  |  The guess is kept in %esi:%edi                                           |
  78  +---------------------------------------------------------------------------*/
  79
  80 #include <gnu/i386/fpemul/fpu_asm.h>
  81
  82
  83 .data
  84 /*
  85         Local storage:
  86  */
  87         ALIGN_DATA
  88 accum_3:
  89         .long   0               /* ms word */
  90 accum_2:
  91         .long   0
  92 accum_1:
  93         .long   0
  94 accum_0:
  95         .long   0
  96
  97 /* The de-normalised argument:
  98 //                  sq_2                  sq_1              sq_0
  99 //        b b b b b b b ... b b b   b b b .... b b b   b 0 0 0 ... 0
 100 //           ^ binary point here */
 101 fsqrt_arg_2:
 102         .long   0               /* ms word */
 103 fsqrt_arg_1:
 104         .long   0
 105 fsqrt_arg_0:
 106         .long   0               /* ls word, at most the ms bit is set */
 107
 108 .text
 109
 110 ENTRY(wm_sqrt)
 111         pushl   %ebp
 112         movl    %esp,%ebp
 113         pushl   %esi
 114         pushl   %edi
 115         pushl   %ebx
 116
 117         movl    PARAM1,%esi
 118
 119         movl    SIGH(%esi),%eax
 120         movl    SIGL(%esi),%ecx
 121         xorl    %edx,%edx
 122
 123 /* We use a rough linear estimate for the first guess.. */
 124
 125         cmpl    EXP_BIAS,EXP(%esi)
 126         jnz     sqrt_arg_ge_2
 127
 128         shrl    $1,%eax                 /* arg is in the range  [1.0 .. 2.0) */
 129         rcrl    $1,%ecx
 130         rcrl    $1,%edx
 131
 132 sqrt_arg_ge_2:
 133 /* From here on, n is never accessed directly again until it is
 134 // replaced by the answer. */
 135
 136         movl    %eax,fsqrt_arg_2                /* ms word of n */
 137         movl    %ecx,fsqrt_arg_1
 138         movl    %edx,fsqrt_arg_0
 139
 140 /* Make a linear first estimate */
 141         shrl    $1,%eax
 142         addl    $0x40000000,%eax
 143         movl    $0xaaaaaaaa,%ecx
 144         mull    %ecx
 145         shll    %edx                    /* max result was 7fff... */
 146         testl   $0x80000000,%edx        /* but min was 3fff... */
 147         jnz     sqrt_prelim_no_adjust
 148
 149         movl    $0x80000000,%edx        /* round up */
 150
 151 sqrt_prelim_no_adjust:
 152         movl    %edx,%esi       /* Our first guess */
 153
 154 /* We have now computed (approx)   (2 + x) / 3, which forms the basis
 155    for a few iterations of Newton's method */
 156
 157         movl    fsqrt_arg_2,%ecx        /* ms word */
 158
 159 /* From our initial estimate, three iterations are enough to get us
 160 // to 30 bits or so. This will then allow two iterations at better
 161 // precision to complete the process.
 162
 163 // Compute  (g + n/g)/2  at each iteration (g is the guess). */
 164         shrl    %ecx            /* Doing this first will prevent a divide */
 165                                 /* overflow later. */
 166
 167         movl    %ecx,%edx       /* msw of the arg / 2 */
 168         divl    %esi            /* current estimate */
 169         shrl    %esi            /* divide by 2 */
 170         addl    %eax,%esi       /* the new estimate */
 171
 172         movl    %ecx,%edx
 173         divl    %esi
 174         shrl    %esi
 175         addl    %eax,%esi
 176
 177         movl    %ecx,%edx
 178         divl    %esi
 179         shrl    %esi
 180         addl    %eax,%esi
 181
 182 /* Now that an estimate accurate to about 30 bits has been obtained (in %esi),
 183 // we improve it to 60 bits or so.
 184
 185 // The strategy from now on is to compute new estimates from
 186 //      guess := guess + (n - guess^2) / (2 * guess) */
 187
 188 /* First, find the square of the guess */
 189         movl    %esi,%eax
 190         mull    %esi
 191 /* guess^2 now in %edx:%eax */
 192
 193         movl    fsqrt_arg_1,%ecx
 194         subl    %ecx,%eax
 195         movl    fsqrt_arg_2,%ecx        /* ms word of normalized n */
 196         sbbl    %ecx,%edx
 197         jnc     sqrt_stage_2_positive
 198 /* subtraction gives a negative result
 199 // negate the result before division */
 200         notl    %edx
 201         notl    %eax
 202         addl    $1,%eax
 203         adcl    $0,%edx
 204
 205         divl    %esi
 206         movl    %eax,%ecx
 207
 208         movl    %edx,%eax
 209         divl    %esi
 210         jmp     sqrt_stage_2_finish
 211
 212 sqrt_stage_2_positive:
 213         divl    %esi
 214         movl    %eax,%ecx
 215
 216         movl    %edx,%eax
 217         divl    %esi
 218
 219         notl    %ecx
 220         notl    %eax
 221         addl    $1,%eax
 222         adcl    $0,%ecx
 223
 224 sqrt_stage_2_finish:
 225         sarl    $1,%ecx         /* divide by 2 */
 226         rcrl    $1,%eax
 227
 228         /* Form the new estimate in %esi:%edi */
 229         movl    %eax,%edi
 230         addl    %ecx,%esi
 231
 232         jnz     sqrt_stage_2_done       /* result should be [1..2) */
 233
 234 #ifdef PARANOID
 235 /* It should be possible to get here only if the arg is ffff....ffff*/
 236         cmp     $0xffffffff,fsqrt_arg_1
 237         jnz     sqrt_stage_2_error
 238 #endif PARANOID
 239
 240 /* The best rounded result.*/
 241         xorl    %eax,%eax
 242         decl    %eax
 243         movl    %eax,%edi
 244         movl    %eax,%esi
 245         movl    $0x7fffffff,%eax
 246         jmp     sqrt_round_result
 247
 248 #ifdef PARANOID
 249 sqrt_stage_2_error:
 250         pushl   EX_INTERNAL|0x213
 251         call    EXCEPTION
 252 #endif PARANOID
 253
 254 sqrt_stage_2_done:
 255
 256 /* Now the square root has been computed to better than 60 bits */
 257
 258 /* Find the square of the guess*/
 259         movl    %edi,%eax               /* ls word of guess*/
 260         mull    %edi
 261         movl    %edx,accum_1
 262
 263         movl    %esi,%eax
 264         mull    %esi
 265         movl    %edx,accum_3
 266         movl    %eax,accum_2
 267
 268         movl    %edi,%eax
 269         mull    %esi
 270         addl    %eax,accum_1
 271         adcl    %edx,accum_2
 272         adcl    $0,accum_3
 273
 274 /*      movl    %esi,%eax*/
 275 /*      mull    %edi*/
 276         addl    %eax,accum_1
 277         adcl    %edx,accum_2
 278         adcl    $0,accum_3
 279
 280 /* guess^2 now in accum_3:accum_2:accum_1*/
 281
 282         movl    fsqrt_arg_0,%eax                /* get normalized n*/
 283         subl    %eax,accum_1
 284         movl    fsqrt_arg_1,%eax
 285         sbbl    %eax,accum_2
 286         movl    fsqrt_arg_2,%eax                /* ms word of normalized n*/
 287         sbbl    %eax,accum_3
 288         jnc     sqrt_stage_3_positive
 289
 290 /* subtraction gives a negative result*/
 291 /* negate the result before division */
 292         notl    accum_1
 293         notl    accum_2
 294         notl    accum_3
 295         addl    $1,accum_1
 296         adcl    $0,accum_2
 297
 298 #ifdef PARANOID
 299         adcl    $0,accum_3      /* This must be zero */
 300         jz      sqrt_stage_3_no_error
 301
 302 sqrt_stage_3_error:
 303         pushl   EX_INTERNAL|0x207
 304         call    EXCEPTION
 305
 306 sqrt_stage_3_no_error:
 307 #endif PARANOID
 308
 309         movl    accum_2,%edx
 310         movl    accum_1,%eax
 311         divl    %esi
 312         movl    %eax,%ecx
 313
 314         movl    %edx,%eax
 315         divl    %esi
 316
 317         sarl    $1,%ecx         /* divide by 2*/
 318         rcrl    $1,%eax
 319
 320         /* prepare to round the result*/
 321
 322         addl    %ecx,%edi
 323         adcl    $0,%esi
 324
 325         jmp     sqrt_stage_3_finished
 326
 327 sqrt_stage_3_positive:
 328         movl    accum_2,%edx
 329         movl    accum_1,%eax
 330         divl    %esi
 331         movl    %eax,%ecx
 332
 333         movl    %edx,%eax
 334         divl    %esi
 335
 336         sarl    $1,%ecx         /* divide by 2*/
 337         rcrl    $1,%eax
 338
 339         /* prepare to round the result*/
 340
 341         notl    %eax            /* Negate the correction term*/
 342         notl    %ecx
 343         addl    $1,%eax
 344         adcl    $0,%ecx         /* carry here ==> correction == 0*/
 345         adcl    $0xffffffff,%esi
 346
 347         addl    %ecx,%edi
 348         adcl    $0,%esi
 349
 350 sqrt_stage_3_finished:
 351
 352 /* The result in %esi:%edi:%esi should be good to about 90 bits here,
 353 // and the rounding information here does not have sufficient accuracy
 354 // in a few rare cases. */
 355         cmpl    $0xffffffe0,%eax
 356         ja      sqrt_near_exact_x
 357
 358         cmpl    $0x00000020,%eax
 359         jb      sqrt_near_exact
 360
 361         cmpl    $0x7fffffe0,%eax
 362         jb      sqrt_round_result
 363
 364         cmpl    $0x80000020,%eax
 365         jb      sqrt_get_more_precision
 366
 367 sqrt_round_result:
 368 /* Set up for rounding operations*/
 369         movl    %eax,%edx
 370         movl    %esi,%eax
 371         movl    %edi,%ebx
 372         movl    PARAM1,%edi
 373         movl    EXP_BIAS,EXP(%edi)      /* Result is in  [1.0 .. 2.0)*/
 374         movl    PARAM2,%ecx
 375         jmp     FPU_round_sqrt
 376
 377
 378 sqrt_near_exact_x:
 379 /* First, the estimate must be rounded up.*/
 380         addl    $1,%edi
 381         adcl    $0,%esi
 382
 383 sqrt_near_exact:
 384 /* This is an easy case because x^1/2 is monotonic.
 385 // We need just find the square of our estimate, compare it
 386 // with the argument, and deduce whether our estimate is
 387 // above, below, or exact. We use the fact that the estimate
 388 // is known to be accurate to about 90 bits. */
 389         movl    %edi,%eax               /* ls word of guess*/
 390         mull    %edi
 391         movl    %edx,%ebx               /* 2nd ls word of square*/
 392         movl    %eax,%ecx               /* ls word of square*/
 393
 394         movl    %edi,%eax
 395         mull    %esi
 396         addl    %eax,%ebx
 397         addl    %eax,%ebx
 398
 399 #ifdef PARANOID
 400         cmp     $0xffffffb0,%ebx
 401         jb      sqrt_near_exact_ok
 402
 403         cmp     $0x00000050,%ebx
 404         ja      sqrt_near_exact_ok
 405
 406         pushl   EX_INTERNAL|0x214
 407         call    EXCEPTION
 408
 409 sqrt_near_exact_ok:
 410 #endif PARANOID
 411
 412         or      %ebx,%ebx
 413         js      sqrt_near_exact_small
 414
 415         jnz     sqrt_near_exact_large
 416
 417         or      %ebx,%edx
 418         jnz     sqrt_near_exact_large
 419
 420 /* Our estimate is exactly the right answer*/
 421         xorl    %eax,%eax
 422         jmp     sqrt_round_result
 423
 424 sqrt_near_exact_small:
 425 /* Our estimate is too small*/
 426         movl    $0x000000ff,%eax
 427         jmp     sqrt_round_result
 428
 429 sqrt_near_exact_large:
 430 /* Our estimate is too large, we need to decrement it*/
 431         subl    $1,%edi
 432         sbbl    $0,%esi
 433         movl    $0xffffff00,%eax
 434         jmp     sqrt_round_result
 435
 436
 437 sqrt_get_more_precision:
 438 /* This case is almost the same as the above, except we start*/
 439 /* with an extra bit of precision in the estimate.*/
 440         stc                     /* The extra bit.*/
 441         rcll    $1,%edi         /* Shift the estimate left one bit*/
 442         rcll    $1,%esi
 443
 444         movl    %edi,%eax               /* ls word of guess*/
 445         mull    %edi
 446         movl    %edx,%ebx               /* 2nd ls word of square*/
 447         movl    %eax,%ecx               /* ls word of square*/
 448
 449         movl    %edi,%eax
 450         mull    %esi
 451         addl    %eax,%ebx
 452         addl    %eax,%ebx
 453
 454 /* Put our estimate back to its original value*/
 455         stc                     /* The ms bit.*/
 456         rcrl    $1,%esi         /* Shift the estimate left one bit*/
 457         rcrl    $1,%edi
 458
 459 #ifdef PARANOID
 460         cmp     $0xffffff60,%ebx
 461         jb      sqrt_more_prec_ok
 462
 463         cmp     $0x000000a0,%ebx
 464         ja      sqrt_more_prec_ok
 465
 466         pushl   EX_INTERNAL|0x215
 467         call    EXCEPTION
 468
 469 sqrt_more_prec_ok:
 470 #endif PARANOID
 471
 472         or      %ebx,%ebx
 473         js      sqrt_more_prec_small
 474
 475         jnz     sqrt_more_prec_large
 476
 477         or      %ebx,%ecx
 478         jnz     sqrt_more_prec_large
 479
 480 /* Our estimate is exactly the right answer*/
 481         movl    $0x80000000,%eax
 482         jmp     sqrt_round_result
 483
 484 sqrt_more_prec_small:
 485 /* Our estimate is too small*/
 486         movl    $0x800000ff,%eax
 487         jmp     sqrt_round_result
 488
 489 sqrt_more_prec_large:
 490 /* Our estimate is too large*/
 491         movl    $0x7fffff00,%eax
 492         jmp     sqrt_round_result