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  * $DragonFly: src/sys/i386/gnu/fpemul/Attic/reg_u_div.s,v 1.4 2003/08/07 21:17:20 dillon Exp $
  62  *
  63  */
  64
  65 /*---------------------------------------------------------------------------+
  66  |  Kernel for the division routines.                                        |
  67  |                                                                           |
  68  |  void reg_u_div(FPU_REG *a, FPU_REG *a,                                   |
  69  |                 FPU_REG *dest, unsigned int control_word)                 |
  70  |                                                                           |
  71  |  Does not compute the destination exponent, but does adjust it.           |
  72  +---------------------------------------------------------------------------*/
  73
  74 #include "fpu_asm.h"
  75
  76
  77 /* #define      dSIGL(x)        (x) */
  78 /* #define      dSIGH(x)        4(x) */
  79
  80
  81 .data
  82 /*
  83         Local storage:
  84         Result:         accum_3:accum_2:accum_1:accum_0
  85         Overflow flag:  ovfl_flag
  86  */
  87         ALIGN_DATA
  88 accum_3:
  89         .long   0
  90 accum_2:
  91         .long   0
  92 accum_1:
  93         .long   0
  94 accum_0:
  95         .long   0
  96 result_1:
  97         .long   0
  98 result_2:
  99         .long   0
 100 ovfl_flag:
 101         .byte   0
 102
 103
 104 .text
 105
 106 .globl _divide_kernel
 107
 108 ENTRY(reg_u_div)
 109         pushl   %ebp
 110         movl    %esp,%ebp
 111
 112         pushl   %esi
 113         pushl   %edi
 114         pushl   %ebx
 115
 116         movl    PARAM1,%esi     /* pointer to num */
 117         movl    PARAM2,%ebx     /* pointer to denom */
 118         movl    PARAM3,%edi     /* pointer to answer */
 119
 120 #ifdef DENORM_OPERAND
 121         movl    EXP(%esi),%eax
 122         cmpl    EXP_UNDER,%eax
 123         jg      xOp1_not_denorm
 124
 125         call    denormal_operand
 126         orl     %eax,%eax
 127         jnz     FPU_Arith_exit
 128
 129 xOp1_not_denorm:
 130         movl    EXP(%ebx),%eax
 131         cmpl    EXP_UNDER,%eax
 132         jg      xOp2_not_denorm
 133
 134         call    denormal_operand
 135         orl     %eax,%eax
 136         jnz     FPU_Arith_exit
 137
 138 xOp2_not_denorm:
 139 #endif DENORM_OPERAND
 140
 141 _divide_kernel:
 142 #ifdef PARANOID
 143 /*      testl   $0x80000000, SIGH(%esi) *//* Dividend */
 144 /*      je      L_bugged */
 145         testl   $0x80000000, SIGH(%ebx) /* Divisor*/
 146         je      L_bugged
 147 #endif PARANOID
 148
 149 /* Check if the divisor can be treated as having just 32 bits */
 150         cmpl    $0,SIGL(%ebx)
 151         jnz     L_Full_Division /* Can't do a quick divide */
 152
 153 /* We should be able to zip through the division here */
 154         movl    SIGH(%ebx),%ecx /* The divisor */
 155         movl    SIGH(%esi),%edx /* Dividend */
 156         movl    SIGL(%esi),%eax /* Dividend */
 157
 158         cmpl    %ecx,%edx
 159         setaeb  ovfl_flag       /* Keep a record */
 160         jb      L_no_adjust
 161
 162         subl    %ecx,%edx       /* Prevent the overflow */
 163
 164 L_no_adjust:
 165         /* Divide the 64 bit number by the 32 bit denominator */
 166         divl    %ecx
 167         movl    %eax,result_2
 168
 169         /* Work on the remainder of the first division */
 170         xorl    %eax,%eax
 171         divl    %ecx
 172         movl    %eax,result_1
 173
 174         /* Work on the remainder of the 64 bit division */
 175         xorl    %eax,%eax
 176         divl    %ecx
 177
 178         testb   $255,ovfl_flag  /* was the num > denom ? */
 179         je      L_no_overflow
 180
 181         /* Do the shifting here */
 182         /* increase the exponent */
 183         incl    EXP(%edi)
 184
 185         /* shift the mantissa right one bit */
 186         stc                     /* To set the ms bit */
 187         rcrl    result_2
 188         rcrl    result_1
 189         rcrl    %eax
 190
 191 L_no_overflow:
 192         jmp     LRound_precision        /* Do the rounding as required*/
 193
 194
 195 /*---------------------------------------------------------------------------+
 196  |  Divide:   Return  arg1/arg2 to arg3.                                     |
 197  |                                                                           |
 198  |  This routine does not use the exponents of arg1 and arg2, but does       |
 199  |  adjust the exponent of arg3.                                             |
 200  |                                                                           |
 201  |  The maximum returned value is (ignoring exponents)                       |
 202  |               .ffffffff ffffffff                                          |
 203  |               ------------------  =  1.ffffffff fffffffe                  |
 204  |               .80000000 00000000                                          |
 205  | and the minimum is                                                        |
 206  |               .80000000 00000000                                          |
 207  |               ------------------  =  .80000000 00000001   (rounded)       |
 208  |               .ffffffff ffffffff                                          |
 209  |                                                                           |
 210  +---------------------------------------------------------------------------*/
 211
 212
 213 L_Full_Division:
 214         /* Save extended dividend in local register*/
 215         movl    SIGL(%esi),%eax
 216         movl    %eax,accum_2
 217         movl    SIGH(%esi),%eax
 218         movl    %eax,accum_3
 219         xorl    %eax,%eax
 220         movl    %eax,accum_1    /* zero the extension */
 221         movl    %eax,accum_0    /* zero the extension */
 222
 223         movl    SIGL(%esi),%eax /* Get the current num */
 224         movl    SIGH(%esi),%edx
 225
 226 /*----------------------------------------------------------------------*/
 227 /* Initialization done */
 228 /* Do the first 32 bits */
 229
 230         movb    $0,ovfl_flag
 231         cmpl    SIGH(%ebx),%edx /* Test for imminent overflow */
 232         jb      LLess_than_1
 233         ja      LGreater_than_1
 234
 235         cmpl    SIGL(%ebx),%eax
 236         jb      LLess_than_1
 237
 238 LGreater_than_1:
 239 /* The dividend is greater or equal, would cause overflow */
 240         setaeb  ovfl_flag               /* Keep a record */
 241
 242         subl    SIGL(%ebx),%eax
 243         sbbl    SIGH(%ebx),%edx /* Prevent the overflow */
 244         movl    %eax,accum_2
 245         movl    %edx,accum_3
 246
 247 LLess_than_1:
 248 /* At this point, we have a dividend < divisor, with a record of
 249    adjustment in ovfl_flag */
 250
 251         /* We will divide by a number which is too large */
 252         movl    SIGH(%ebx),%ecx
 253         addl    $1,%ecx
 254         jnc     LFirst_div_not_1
 255
 256         /* here we need to divide by 100000000h,
 257            i.e., no division at all.. */
 258         mov     %edx,%eax
 259         jmp     LFirst_div_done
 260
 261 LFirst_div_not_1:
 262         divl    %ecx            /* Divide the numerator by the augmented
 263                                    denom ms dw */
 264
 265 LFirst_div_done:
 266         movl    %eax,result_2   /* Put the result in the answer */
 267
 268         mull    SIGH(%ebx)      /* mul by the ms dw of the denom */
 269
 270         subl    %eax,accum_2    /* Subtract from the num local reg */
 271         sbbl    %edx,accum_3
 272
 273         movl    result_2,%eax   /* Get the result back */
 274         mull    SIGL(%ebx)      /* now mul the ls dw of the denom */
 275
 276         subl    %eax,accum_1    /* Subtract from the num local reg */
 277         sbbl    %edx,accum_2
 278         sbbl    $0,accum_3
 279         je      LDo_2nd_32_bits         /* Must check for non-zero result here */
 280
 281 #ifdef PARANOID
 282         jb      L_bugged_1
 283 #endif PARANOID
 284
 285         /* need to subtract another once of the denom */
 286         incl    result_2        /* Correct the answer */
 287
 288         movl    SIGL(%ebx),%eax
 289         movl    SIGH(%ebx),%edx
 290         subl    %eax,accum_1    /* Subtract from the num local reg */
 291         sbbl    %edx,accum_2
 292
 293 #ifdef PARANOID
 294         sbbl    $0,accum_3
 295         jne     L_bugged_1      /* Must check for non-zero result here */
 296 #endif PARANOID
 297
 298 /*----------------------------------------------------------------------*/
 299 /* Half of the main problem is done, there is just a reduced numerator
 300    to handle now */
 301 /* Work with the second 32 bits, accum_0 not used from now on */
 302 LDo_2nd_32_bits:
 303         movl    accum_2,%edx    /* get the reduced num */
 304         movl    accum_1,%eax
 305
 306         /* need to check for possible subsequent overflow */
 307         cmpl    SIGH(%ebx),%edx
 308         jb      LDo_2nd_div
 309         ja      LPrevent_2nd_overflow
 310
 311         cmpl    SIGL(%ebx),%eax
 312         jb      LDo_2nd_div
 313
 314 LPrevent_2nd_overflow:
 315 /* The numerator is greater or equal, would cause overflow */
 316         /* prevent overflow */
 317         subl    SIGL(%ebx),%eax
 318         sbbl    SIGH(%ebx),%edx
 319         movl    %edx,accum_2
 320         movl    %eax,accum_1
 321
 322         incl    result_2        /* Reflect the subtraction in the answer */
 323
 324 #ifdef PARANOID
 325         je      L_bugged_2      /* Can't bump the result to 1.0 */
 326 #endif PARANOID
 327
 328 LDo_2nd_div:
 329         cmpl    $0,%ecx         /* augmented denom msw*/
 330         jnz     LSecond_div_not_1
 331
 332         /* %ecx == 0, we are dividing by 1.0 */
 333         mov     %edx,%eax
 334         jmp     LSecond_div_done
 335
 336 LSecond_div_not_1:
 337         divl    %ecx            /* Divide the numerator by the denom ms dw */
 338
 339 LSecond_div_done:
 340         movl    %eax,result_1   /* Put the result in the answer */
 341
 342         mull    SIGH(%ebx)      /* mul by the ms dw of the denom */
 343
 344         subl    %eax,accum_1    /* Subtract from the num local reg */
 345         sbbl    %edx,accum_2
 346
 347 #ifdef PARANOID
 348         jc      L_bugged_2
 349 #endif PARANOID
 350
 351         movl    result_1,%eax   /* Get the result back */
 352         mull    SIGL(%ebx)      /* now mul the ls dw of the denom */
 353
 354         subl    %eax,accum_0    /* Subtract from the num local reg */
 355         sbbl    %edx,accum_1    /* Subtract from the num local reg */
 356         sbbl    $0,accum_2
 357
 358 #ifdef PARANOID
 359         jc      L_bugged_2
 360 #endif PARANOID
 361
 362         jz      LDo_3rd_32_bits
 363
 364 #ifdef PARANOID
 365         cmpl    $1,accum_2
 366         jne     L_bugged_2
 367 #endif PARANOID
 368
 369         /* need to subtract another once of the denom */
 370         movl    SIGL(%ebx),%eax
 371         movl    SIGH(%ebx),%edx
 372         subl    %eax,accum_0    /* Subtract from the num local reg */
 373         sbbl    %edx,accum_1
 374         sbbl    $0,accum_2
 375
 376 #ifdef PARANOID
 377         jc      L_bugged_2
 378         jne     L_bugged_2
 379 #endif PARANOID
 380
 381         addl    $1,result_1     /* Correct the answer */
 382         adcl    $0,result_2
 383
 384 #ifdef PARANOID
 385         jc      L_bugged_2      /* Must check for non-zero result here */
 386 #endif PARANOID
 387
 388 /*----------------------------------------------------------------------*/
 389 /* The division is essentially finished here, we just need to perform
 390    tidying operations. */
 391 /* deal with the 3rd 32 bits */
 392 LDo_3rd_32_bits:
 393         movl    accum_1,%edx            /* get the reduced num */
 394         movl    accum_0,%eax
 395
 396         /* need to check for possible subsequent overflow */
 397         cmpl    SIGH(%ebx),%edx /* denom*/
 398         jb      LRound_prep
 399         ja      LPrevent_3rd_overflow
 400
 401         cmpl    SIGL(%ebx),%eax /* denom */
 402         jb      LRound_prep
 403
 404 LPrevent_3rd_overflow:
 405         /* prevent overflow */
 406         subl    SIGL(%ebx),%eax
 407         sbbl    SIGH(%ebx),%edx
 408         movl    %edx,accum_1
 409         movl    %eax,accum_0
 410
 411         addl    $1,result_1     /* Reflect the subtraction in the answer */
 412         adcl    $0,result_2
 413         jne     LRound_prep
 414         jnc     LRound_prep
 415
 416         /* This is a tricky spot, there is an overflow of the answer */
 417         movb    $255,ovfl_flag          /* Overflow -> 1.000 */
 418
 419 LRound_prep:
 420 /* Prepare for rounding.
 421 // To test for rounding, we just need to compare 2*accum with the
 422 // denom. */
 423         movl    accum_0,%ecx
 424         movl    accum_1,%edx
 425         movl    %ecx,%eax
 426         orl     %edx,%eax
 427         jz      LRound_ovfl             /* The accumulator contains zero.*/
 428
 429         /* Multiply by 2 */
 430         clc
 431         rcll    $1,%ecx
 432         rcll    $1,%edx
 433         jc      LRound_large            /* No need to compare, denom smaller */
 434
 435         subl    SIGL(%ebx),%ecx
 436         sbbl    SIGH(%ebx),%edx
 437         jnc     LRound_not_small
 438
 439         movl    $0x70000000,%eax        /* Denom was larger */
 440         jmp     LRound_ovfl
 441
 442 LRound_not_small:
 443         jnz     LRound_large
 444
 445         movl    $0x80000000,%eax        /* Remainder was exactly 1/2 denom */
 446         jmp     LRound_ovfl
 447
 448 LRound_large:
 449         movl    $0xff000000,%eax        /* Denom was smaller */
 450
 451 LRound_ovfl:
 452 /* We are now ready to deal with rounding, but first we must get
 453    the bits properly aligned */
 454         testb   $255,ovfl_flag  /* was the num > denom ? */
 455         je      LRound_precision
 456
 457         incl    EXP(%edi)
 458
 459         /* shift the mantissa right one bit */
 460         stc                     /* Will set the ms bit */
 461         rcrl    result_2
 462         rcrl    result_1
 463         rcrl    %eax
 464
 465 /* Round the result as required */
 466 LRound_precision:
 467         decl    EXP(%edi)       /* binary point between 1st & 2nd bits */
 468
 469         movl    %eax,%edx
 470         movl    result_1,%ebx
 471         movl    result_2,%eax
 472         jmp     FPU_round
 473
 474
 475 #ifdef PARANOID
 476 /* The logic is wrong if we got here */
 477 L_bugged:
 478         pushl   EX_INTERNAL|0x202
 479         call    EXCEPTION
 480         pop     %ebx
 481         jmp     L_exit
 482
 483 L_bugged_1:
 484         pushl   EX_INTERNAL|0x203
 485         call    EXCEPTION
 486         pop     %ebx
 487         jmp     L_exit
 488
 489 L_bugged_2:
 490         pushl   EX_INTERNAL|0x204
 491         call    EXCEPTION
 492         pop     %ebx
 493         jmp     L_exit
 494
 495 L_exit:
 496         popl    %ebx
 497         popl    %edi
 498         popl    %esi
 499
 500         leave
 501         ret
 502 #endif PARANOID