sys/i386/gnu/fpemul/reg_u_div.s

   1         .file   "reg_u_div.S"
   2 /*
   3  *  reg_u_div.S
   4  *
   5  * Core division routines
   6  *
   7  *
   8  * Copyright (C) 1992,1993,1994
   9  *                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,
  10  *                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au
  11  * All rights reserved.
  12  *
  13  * This copyright notice covers the redistribution and use of the
  14  * FPU emulator developed by W. Metzenthen. It covers only its use
  15  * in the 386BSD, FreeBSD and NetBSD operating systems. Any other
  16  * use is not permitted under this copyright.
  17  *
  18  * Redistribution and use in source and binary forms, with or without
  19  * modification, are permitted provided that the following conditions
  20  * are met:
  21  * 1. Redistributions of source code must retain the above copyright
  22  *    notice, this list of conditions and the following disclaimer.
  23  * 2. Redistributions in binary form must include information specifying
  24  *    that source code for the emulator is freely available and include
  25  *    either:
  26  *      a) an offer to provide the source code for a nominal distribution
  27  *         fee, or
  28  *      b) list at least two alternative methods whereby the source
  29  *         can be obtained, e.g. a publically accessible bulletin board
  30  *         and an anonymous ftp site from which the software can be
  31  *         downloaded.
  32  * 3. All advertising materials specifically mentioning features or use of
  33  *    this emulator must acknowledge that it was developed by W. Metzenthen.
  34  * 4. The name of W. Metzenthen may not be used to endorse or promote
  35  *    products derived from this software without specific prior written
  36  *    permission.
  37  *
  38  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
  39  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
  40  * AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL
  41  * W. METZENTHEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  42  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  43  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  44  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  45  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  46  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  47  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  48  *
  49  *
  50  * The purpose of this copyright, based upon the Berkeley copyright, is to
  51  * ensure that the covered software remains freely available to everyone.
  52  *
  53  * The software (with necessary differences) is also available, but under
  54  * the terms of the GNU copyleft, for the Linux operating system and for
  55  * the djgpp ms-dos extender.
  56  *
  57  * W. Metzenthen   June 1994.
  58  *
  59  *
  60  * $FreeBSD: src/sys/gnu/i386/fpemul/reg_u_div.s,v 1.9 1999/08/28 00:42:58 peter Exp $
  61  *
  62  */
  63
  64 /*---------------------------------------------------------------------------+
  65  |  Kernel for the division routines.                                        |
  66  |                                                                           |
  67  |  void reg_u_div(FPU_REG *a, FPU_REG *a,                                   |
  68  |                 FPU_REG *dest, unsigned int control_word)                 |
  69  |                                                                           |
  70  |  Does not compute the destination exponent, but does adjust it.           |
  71  +---------------------------------------------------------------------------*/
  72
  73 #include <gnu/i386/fpemul/fpu_asm.h>
  74
  75
  76 /* #define      dSIGL(x)        (x) */
  77 /* #define      dSIGH(x)        4(x) */
  78
  79
  80 .data
  81 /*
  82         Local storage:
  83         Result:         accum_3:accum_2:accum_1:accum_0
  84         Overflow flag:  ovfl_flag
  85  */
  86         ALIGN_DATA
  87 accum_3:
  88         .long   0
  89 accum_2:
  90         .long   0
  91 accum_1:
  92         .long   0
  93 accum_0:
  94         .long   0
  95 result_1:
  96         .long   0
  97 result_2:
  98         .long   0
  99 ovfl_flag:
 100         .byte   0
 101
 102
 103 .text
 104
 105 .globl _divide_kernel
 106
 107 ENTRY(reg_u_div)
 108         pushl   %ebp
 109         movl    %esp,%ebp
 110
 111         pushl   %esi
 112         pushl   %edi
 113         pushl   %ebx
 114
 115         movl    PARAM1,%esi     /* pointer to num */
 116         movl    PARAM2,%ebx     /* pointer to denom */
 117         movl    PARAM3,%edi     /* pointer to answer */
 118
 119 #ifdef DENORM_OPERAND
 120         movl    EXP(%esi),%eax
 121         cmpl    EXP_UNDER,%eax
 122         jg      xOp1_not_denorm
 123
 124         call    _denormal_operand
 125         orl     %eax,%eax
 126         jnz     FPU_Arith_exit
 127
 128 xOp1_not_denorm:
 129         movl    EXP(%ebx),%eax
 130         cmpl    EXP_UNDER,%eax
 131         jg      xOp2_not_denorm
 132
 133         call    _denormal_operand
 134         orl     %eax,%eax
 135         jnz     FPU_Arith_exit
 136
 137 xOp2_not_denorm:
 138 #endif DENORM_OPERAND
 139
 140 _divide_kernel:
 141 #ifdef PARANOID
 142 /*      testl   $0x80000000, SIGH(%esi) *//* Dividend */
 143 /*      je      L_bugged */
 144         testl   $0x80000000, SIGH(%ebx) /* Divisor*/
 145         je      L_bugged
 146 #endif PARANOID
 147
 148 /* Check if the divisor can be treated as having just 32 bits */
 149         cmpl    $0,SIGL(%ebx)
 150         jnz     L_Full_Division /* Can't do a quick divide */
 151
 152 /* We should be able to zip through the division here */
 153         movl    SIGH(%ebx),%ecx /* The divisor */
 154         movl    SIGH(%esi),%edx /* Dividend */
 155         movl    SIGL(%esi),%eax /* Dividend */
 156
 157         cmpl    %ecx,%edx
 158         setaeb  ovfl_flag       /* Keep a record */
 159         jb      L_no_adjust
 160
 161         subl    %ecx,%edx       /* Prevent the overflow */
 162
 163 L_no_adjust:
 164         /* Divide the 64 bit number by the 32 bit denominator */
 165         divl    %ecx
 166         movl    %eax,result_2
 167
 168         /* Work on the remainder of the first division */
 169         xorl    %eax,%eax
 170         divl    %ecx
 171         movl    %eax,result_1
 172
 173         /* Work on the remainder of the 64 bit division */
 174         xorl    %eax,%eax
 175         divl    %ecx
 176
 177         testb   $255,ovfl_flag  /* was the num > denom ? */
 178         je      L_no_overflow
 179
 180         /* Do the shifting here */
 181         /* increase the exponent */
 182         incl    EXP(%edi)
 183
 184         /* shift the mantissa right one bit */
 185         stc                     /* To set the ms bit */
 186         rcrl    result_2
 187         rcrl    result_1
 188         rcrl    %eax
 189
 190 L_no_overflow:
 191         jmp     LRound_precision        /* Do the rounding as required*/
 192
 193
 194 /*---------------------------------------------------------------------------+
 195  |  Divide:   Return  arg1/arg2 to arg3.                                     |
 196  |                                                                           |
 197  |  This routine does not use the exponents of arg1 and arg2, but does       |
 198  |  adjust the exponent of arg3.                                             |
 199  |                                                                           |
 200  |  The maximum returned value is (ignoring exponents)                       |
 201  |               .ffffffff ffffffff                                          |
 202  |               ------------------  =  1.ffffffff fffffffe                  |
 203  |               .80000000 00000000                                          |
 204  | and the minimum is                                                        |
 205  |               .80000000 00000000                                          |
 206  |               ------------------  =  .80000000 00000001   (rounded)       |
 207  |               .ffffffff ffffffff                                          |
 208  |                                                                           |
 209  +---------------------------------------------------------------------------*/
 210
 211
 212 L_Full_Division:
 213         /* Save extended dividend in local register*/
 214         movl    SIGL(%esi),%eax
 215         movl    %eax,accum_2
 216         movl    SIGH(%esi),%eax
 217         movl    %eax,accum_3
 218         xorl    %eax,%eax
 219         movl    %eax,accum_1    /* zero the extension */
 220         movl    %eax,accum_0    /* zero the extension */
 221
 222         movl    SIGL(%esi),%eax /* Get the current num */
 223         movl    SIGH(%esi),%edx
 224
 225 /*----------------------------------------------------------------------*/
 226 /* Initialization done */
 227 /* Do the first 32 bits */
 228
 229         movb    $0,ovfl_flag
 230         cmpl    SIGH(%ebx),%edx /* Test for imminent overflow */
 231         jb      LLess_than_1
 232         ja      LGreater_than_1
 233
 234         cmpl    SIGL(%ebx),%eax
 235         jb      LLess_than_1
 236
 237 LGreater_than_1:
 238 /* The dividend is greater or equal, would cause overflow */
 239         setaeb  ovfl_flag               /* Keep a record */
 240
 241         subl    SIGL(%ebx),%eax
 242         sbbl    SIGH(%ebx),%edx /* Prevent the overflow */
 243         movl    %eax,accum_2
 244         movl    %edx,accum_3
 245
 246 LLess_than_1:
 247 /* At this point, we have a dividend < divisor, with a record of
 248    adjustment in ovfl_flag */
 249
 250         /* We will divide by a number which is too large */
 251         movl    SIGH(%ebx),%ecx
 252         addl    $1,%ecx
 253         jnc     LFirst_div_not_1
 254
 255         /* here we need to divide by 100000000h,
 256            i.e., no division at all.. */
 257         mov     %edx,%eax
 258         jmp     LFirst_div_done
 259
 260 LFirst_div_not_1:
 261         divl    %ecx            /* Divide the numerator by the augmented
 262                                    denom ms dw */
 263
 264 LFirst_div_done:
 265         movl    %eax,result_2   /* Put the result in the answer */
 266
 267         mull    SIGH(%ebx)      /* mul by the ms dw of the denom */
 268
 269         subl    %eax,accum_2    /* Subtract from the num local reg */
 270         sbbl    %edx,accum_3
 271
 272         movl    result_2,%eax   /* Get the result back */
 273         mull    SIGL(%ebx)      /* now mul the ls dw of the denom */
 274
 275         subl    %eax,accum_1    /* Subtract from the num local reg */
 276         sbbl    %edx,accum_2
 277         sbbl    $0,accum_3
 278         je      LDo_2nd_32_bits         /* Must check for non-zero result here */
 279
 280 #ifdef PARANOID
 281         jb      L_bugged_1
 282 #endif PARANOID
 283
 284         /* need to subtract another once of the denom */
 285         incl    result_2        /* Correct the answer */
 286
 287         movl    SIGL(%ebx),%eax
 288         movl    SIGH(%ebx),%edx
 289         subl    %eax,accum_1    /* Subtract from the num local reg */
 290         sbbl    %edx,accum_2
 291
 292 #ifdef PARANOID
 293         sbbl    $0,accum_3
 294         jne     L_bugged_1      /* Must check for non-zero result here */
 295 #endif PARANOID
 296
 297 /*----------------------------------------------------------------------*/
 298 /* Half of the main problem is done, there is just a reduced numerator
 299    to handle now */
 300 /* Work with the second 32 bits, accum_0 not used from now on */
 301 LDo_2nd_32_bits:
 302         movl    accum_2,%edx    /* get the reduced num */
 303         movl    accum_1,%eax
 304
 305         /* need to check for possible subsequent overflow */
 306         cmpl    SIGH(%ebx),%edx
 307         jb      LDo_2nd_div
 308         ja      LPrevent_2nd_overflow
 309
 310         cmpl    SIGL(%ebx),%eax
 311         jb      LDo_2nd_div
 312
 313 LPrevent_2nd_overflow:
 314 /* The numerator is greater or equal, would cause overflow */
 315         /* prevent overflow */
 316         subl    SIGL(%ebx),%eax
 317         sbbl    SIGH(%ebx),%edx
 318         movl    %edx,accum_2
 319         movl    %eax,accum_1
 320
 321         incl    result_2        /* Reflect the subtraction in the answer */
 322
 323 #ifdef PARANOID
 324         je      L_bugged_2      /* Can't bump the result to 1.0 */
 325 #endif PARANOID
 326
 327 LDo_2nd_div:
 328         cmpl    $0,%ecx         /* augmented denom msw*/
 329         jnz     LSecond_div_not_1
 330
 331         /* %ecx == 0, we are dividing by 1.0 */
 332         mov     %edx,%eax
 333         jmp     LSecond_div_done
 334
 335 LSecond_div_not_1:
 336         divl    %ecx            /* Divide the numerator by the denom ms dw */
 337
 338 LSecond_div_done:
 339         movl    %eax,result_1   /* Put the result in the answer */
 340
 341         mull    SIGH(%ebx)      /* mul by the ms dw of the denom */
 342
 343         subl    %eax,accum_1    /* Subtract from the num local reg */
 344         sbbl    %edx,accum_2
 345
 346 #ifdef PARANOID
 347         jc      L_bugged_2
 348 #endif PARANOID
 349
 350         movl    result_1,%eax   /* Get the result back */
 351         mull    SIGL(%ebx)      /* now mul the ls dw of the denom */
 352
 353         subl    %eax,accum_0    /* Subtract from the num local reg */
 354         sbbl    %edx,accum_1    /* Subtract from the num local reg */
 355         sbbl    $0,accum_2
 356
 357 #ifdef PARANOID
 358         jc      L_bugged_2
 359 #endif PARANOID
 360
 361         jz      LDo_3rd_32_bits
 362
 363 #ifdef PARANOID
 364         cmpl    $1,accum_2
 365         jne     L_bugged_2
 366 #endif PARANOID
 367
 368         /* need to subtract another once of the denom */
 369         movl    SIGL(%ebx),%eax
 370         movl    SIGH(%ebx),%edx
 371         subl    %eax,accum_0    /* Subtract from the num local reg */
 372         sbbl    %edx,accum_1
 373         sbbl    $0,accum_2
 374
 375 #ifdef PARANOID
 376         jc      L_bugged_2
 377         jne     L_bugged_2
 378 #endif PARANOID
 379
 380         addl    $1,result_1     /* Correct the answer */
 381         adcl    $0,result_2
 382
 383 #ifdef PARANOID
 384         jc      L_bugged_2      /* Must check for non-zero result here */
 385 #endif PARANOID
 386
 387 /*----------------------------------------------------------------------*/
 388 /* The division is essentially finished here, we just need to perform
 389    tidying operations. */
 390 /* deal with the 3rd 32 bits */
 391 LDo_3rd_32_bits:
 392         movl    accum_1,%edx            /* get the reduced num */
 393         movl    accum_0,%eax
 394
 395         /* need to check for possible subsequent overflow */
 396         cmpl    SIGH(%ebx),%edx /* denom*/
 397         jb      LRound_prep
 398         ja      LPrevent_3rd_overflow
 399
 400         cmpl    SIGL(%ebx),%eax /* denom */
 401         jb      LRound_prep
 402
 403 LPrevent_3rd_overflow:
 404         /* prevent overflow */
 405         subl    SIGL(%ebx),%eax
 406         sbbl    SIGH(%ebx),%edx
 407         movl    %edx,accum_1
 408         movl    %eax,accum_0
 409
 410         addl    $1,result_1     /* Reflect the subtraction in the answer */
 411         adcl    $0,result_2
 412         jne     LRound_prep
 413         jnc     LRound_prep
 414
 415         /* This is a tricky spot, there is an overflow of the answer */
 416         movb    $255,ovfl_flag          /* Overflow -> 1.000 */
 417
 418 LRound_prep:
 419 /* Prepare for rounding.
 420 // To test for rounding, we just need to compare 2*accum with the
 421 // denom. */
 422         movl    accum_0,%ecx
 423         movl    accum_1,%edx
 424         movl    %ecx,%eax
 425         orl     %edx,%eax
 426         jz      LRound_ovfl             /* The accumulator contains zero.*/
 427
 428         /* Multiply by 2 */
 429         clc
 430         rcll    $1,%ecx
 431         rcll    $1,%edx
 432         jc      LRound_large            /* No need to compare, denom smaller */
 433
 434         subl    SIGL(%ebx),%ecx
 435         sbbl    SIGH(%ebx),%edx
 436         jnc     LRound_not_small
 437
 438         movl    $0x70000000,%eax        /* Denom was larger */
 439         jmp     LRound_ovfl
 440
 441 LRound_not_small:
 442         jnz     LRound_large
 443
 444         movl    $0x80000000,%eax        /* Remainder was exactly 1/2 denom */
 445         jmp     LRound_ovfl
 446
 447 LRound_large:
 448         movl    $0xff000000,%eax        /* Denom was smaller */
 449
 450 LRound_ovfl:
 451 /* We are now ready to deal with rounding, but first we must get
 452    the bits properly aligned */
 453         testb   $255,ovfl_flag  /* was the num > denom ? */
 454         je      LRound_precision
 455
 456         incl    EXP(%edi)
 457
 458         /* shift the mantissa right one bit */
 459         stc                     /* Will set the ms bit */
 460         rcrl    result_2
 461         rcrl    result_1
 462         rcrl    %eax
 463
 464 /* Round the result as required */
 465 LRound_precision:
 466         decl    EXP(%edi)       /* binary point between 1st & 2nd bits */
 467
 468         movl    %eax,%edx
 469         movl    result_1,%ebx
 470         movl    result_2,%eax
 471         jmp     FPU_round
 472
 473
 474 #ifdef PARANOID
 475 /* The logic is wrong if we got here */
 476 L_bugged:
 477         pushl   EX_INTERNAL|0x202
 478         call    EXCEPTION
 479         pop     %ebx
 480         jmp     L_exit
 481
 482 L_bugged_1:
 483         pushl   EX_INTERNAL|0x203
 484         call    EXCEPTION
 485         pop     %ebx
 486         jmp     L_exit
 487
 488 L_bugged_2:
 489         pushl   EX_INTERNAL|0x204
 490         call    EXCEPTION
 491         pop     %ebx
 492         jmp     L_exit
 493
 494 L_exit:
 495         popl    %ebx
 496         popl    %edi
 497         popl    %esi
 498
 499         leave
 500         ret
 501 #endif PARANOID