x86/math-emu/reg_u_div.S

0001 /* SPDX-License-Identifier: GPL-2.0 */
0002     .file   "reg_u_div.S"
0003 /*---------------------------------------------------------------------------+
0004  |  reg_u_div.S                                                              |
0005  |                                                                           |
0006  | Divide one FPU_REG by another and put the result in a destination FPU_REG.|
0007  |                                                                           |
0008  | Copyright (C) 1992,1993,1995,1997                                         |
0009  |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
0010  |                  E-mail   billm@suburbia.net                              |
0011  |                                                                           |
0012  |                                                                           |
0013  +---------------------------------------------------------------------------*/
0014
0015 /*---------------------------------------------------------------------------+
0016  | Call from C as:                                                           |
0017  |    int FPU_u_div(FPU_REG *a, FPU_REG *b, FPU_REG *dest,                   |
0018  |                unsigned int control_word, char *sign)                     |
0019  |                                                                           |
0020  |  Does not compute the destination exponent, but does adjust it.           |
0021  |                                                                           |
0022  |    Return value is the tag of the answer, or-ed with FPU_Exception if     |
0023  |    one was raised, or -1 on internal error.                               |
0024  +---------------------------------------------------------------------------*/
0025
0026 #include "exception.h"
0027 #include "fpu_emu.h"
0028 #include "control_w.h"
0029
0030
0031 /* #define  dSIGL(x)    (x) */
0032 /* #define  dSIGH(x)    4(x) */
0033
0034
0035 #ifndef NON_REENTRANT_FPU
0036 /*
0037     Local storage on the stack:
0038     Result:     FPU_accum_3:FPU_accum_2:FPU_accum_1:FPU_accum_0
0039     Overflow flag:  ovfl_flag
0040  */
0041 #define FPU_accum_3 -4(%ebp)
0042 #define FPU_accum_2 -8(%ebp)
0043 #define FPU_accum_1 -12(%ebp)
0044 #define FPU_accum_0 -16(%ebp)
0045 #define FPU_result_1    -20(%ebp)
0046 #define FPU_result_2    -24(%ebp)
0047 #define FPU_ovfl_flag   -28(%ebp)
0048
0049 #else
0050 .data
0051 /*
0052     Local storage in a static area:
0053     Result:     FPU_accum_3:FPU_accum_2:FPU_accum_1:FPU_accum_0
0054     Overflow flag:  ovfl_flag
0055  */
0056     .align 4,0
0057 FPU_accum_3:
0058     .long   0
0059 FPU_accum_2:
0060     .long   0
0061 FPU_accum_1:
0062     .long   0
0063 FPU_accum_0:
0064     .long   0
0065 FPU_result_1:
0066     .long   0
0067 FPU_result_2:
0068     .long   0
0069 FPU_ovfl_flag:
0070     .byte   0
0071 #endif /* NON_REENTRANT_FPU */
0072
0073 #define REGA    PARAM1
0074 #define REGB    PARAM2
0075 #define DEST    PARAM3
0076
0077 .text
0078 SYM_FUNC_START(FPU_u_div)
0079     pushl   %ebp
0080     movl    %esp,%ebp
0081 #ifndef NON_REENTRANT_FPU
0082     subl    $28,%esp
0083 #endif /* NON_REENTRANT_FPU */
0084
0085     pushl   %esi
0086     pushl   %edi
0087     pushl   %ebx
0088
0089     movl    REGA,%esi
0090     movl    REGB,%ebx
0091     movl    DEST,%edi
0092
0093     movswl  EXP(%esi),%edx
0094     movswl  EXP(%ebx),%eax
0095     subl    %eax,%edx
0096     addl    EXP_BIAS,%edx
0097
0098     /* A denormal and a large number can cause an exponent underflow */
0099     cmpl    EXP_WAY_UNDER,%edx
0100     jg  xExp_not_underflow
0101
0102     /* Set to a really low value allow correct handling */
0103     movl    EXP_WAY_UNDER,%edx
0104
0105 xExp_not_underflow:
0106
0107     movw    %dx,EXP(%edi)
0108
0109 #ifdef PARANOID
0110 /*  testl   $0x80000000, SIGH(%esi) // Dividend */
0111 /*  je  L_bugged */
0112     testl   $0x80000000, SIGH(%ebx) /* Divisor */
0113     je  L_bugged
0114 #endif /* PARANOID */
0115
0116 /* Check if the divisor can be treated as having just 32 bits */
0117     cmpl    $0,SIGL(%ebx)
0118     jnz L_Full_Division /* Can't do a quick divide */
0119
0120 /* We should be able to zip through the division here */
0121     movl    SIGH(%ebx),%ecx /* The divisor */
0122     movl    SIGH(%esi),%edx /* Dividend */
0123     movl    SIGL(%esi),%eax /* Dividend */
0124
0125     cmpl    %ecx,%edx
0126     setaeb  FPU_ovfl_flag   /* Keep a record */
0127     jb  L_no_adjust
0128
0129     subl    %ecx,%edx   /* Prevent the overflow */
0130
0131 L_no_adjust:
0132     /* Divide the 64 bit number by the 32 bit denominator */
0133     divl    %ecx
0134     movl    %eax,FPU_result_2
0135
0136     /* Work on the remainder of the first division */
0137     xorl    %eax,%eax
0138     divl    %ecx
0139     movl    %eax,FPU_result_1
0140
0141     /* Work on the remainder of the 64 bit division */
0142     xorl    %eax,%eax
0143     divl    %ecx
0144
0145     testb   $255,FPU_ovfl_flag  /* was the num > denom ? */
0146     je  L_no_overflow
0147
0148     /* Do the shifting here */
0149     /* increase the exponent */
0150     incw    EXP(%edi)
0151
0152     /* shift the mantissa right one bit */
0153     stc         /* To set the ms bit */
0154     rcrl    FPU_result_2
0155     rcrl    FPU_result_1
0156     rcrl    %eax
0157
0158 L_no_overflow:
0159     jmp LRound_precision    /* Do the rounding as required */
0160
0161
0162 /*---------------------------------------------------------------------------+
0163  |  Divide:   Return  arg1/arg2 to arg3.                                     |
0164  |                                                                           |
0165  |  This routine does not use the exponents of arg1 and arg2, but does       |
0166  |  adjust the exponent of arg3.                                             |
0167  |                                                                           |
0168  |  The maximum returned value is (ignoring exponents)                       |
0169  |               .ffffffff ffffffff                                          |
0170  |               ------------------  =  1.ffffffff fffffffe                  |
0171  |               .80000000 00000000                                          |
0172  | and the minimum is                                                        |
0173  |               .80000000 00000000                                          |
0174  |               ------------------  =  .80000000 00000001   (rounded)       |
0175  |               .ffffffff ffffffff                                          |
0176  |                                                                           |
0177  +---------------------------------------------------------------------------*/
0178
0179
0180 L_Full_Division:
0181     /* Save extended dividend in local register */
0182     movl    SIGL(%esi),%eax
0183     movl    %eax,FPU_accum_2
0184     movl    SIGH(%esi),%eax
0185     movl    %eax,FPU_accum_3
0186     xorl    %eax,%eax
0187     movl    %eax,FPU_accum_1    /* zero the extension */
0188     movl    %eax,FPU_accum_0    /* zero the extension */
0189
0190     movl    SIGL(%esi),%eax /* Get the current num */
0191     movl    SIGH(%esi),%edx
0192
0193 /*----------------------------------------------------------------------*/
0194 /* Initialization done.
0195    Do the first 32 bits. */
0196
0197     movb    $0,FPU_ovfl_flag
0198     cmpl    SIGH(%ebx),%edx /* Test for imminent overflow */
0199     jb  LLess_than_1
0200     ja  LGreater_than_1
0201
0202     cmpl    SIGL(%ebx),%eax
0203     jb  LLess_than_1
0204
0205 LGreater_than_1:
0206 /* The dividend is greater or equal, would cause overflow */
0207     setaeb  FPU_ovfl_flag       /* Keep a record */
0208
0209     subl    SIGL(%ebx),%eax
0210     sbbl    SIGH(%ebx),%edx /* Prevent the overflow */
0211     movl    %eax,FPU_accum_2
0212     movl    %edx,FPU_accum_3
0213
0214 LLess_than_1:
0215 /* At this point, we have a dividend < divisor, with a record of
0216    adjustment in FPU_ovfl_flag */
0217
0218     /* We will divide by a number which is too large */
0219     movl    SIGH(%ebx),%ecx
0220     addl    $1,%ecx
0221     jnc LFirst_div_not_1
0222
0223     /* here we need to divide by 100000000h,
0224        i.e., no division at all.. */
0225     mov %edx,%eax
0226     jmp LFirst_div_done
0227
0228 LFirst_div_not_1:
0229     divl    %ecx        /* Divide the numerator by the augmented
0230                    denom ms dw */
0231
0232 LFirst_div_done:
0233     movl    %eax,FPU_result_2   /* Put the result in the answer */
0234
0235     mull    SIGH(%ebx)  /* mul by the ms dw of the denom */
0236
0237     subl    %eax,FPU_accum_2    /* Subtract from the num local reg */
0238     sbbl    %edx,FPU_accum_3
0239
0240     movl    FPU_result_2,%eax   /* Get the result back */
0241     mull    SIGL(%ebx)  /* now mul the ls dw of the denom */
0242
0243     subl    %eax,FPU_accum_1    /* Subtract from the num local reg */
0244     sbbl    %edx,FPU_accum_2
0245     sbbl    $0,FPU_accum_3
0246     je  LDo_2nd_32_bits     /* Must check for non-zero result here */
0247
0248 #ifdef PARANOID
0249     jb  L_bugged_1
0250 #endif /* PARANOID */
0251
0252     /* need to subtract another once of the denom */
0253     incl    FPU_result_2    /* Correct the answer */
0254
0255     movl    SIGL(%ebx),%eax
0256     movl    SIGH(%ebx),%edx
0257     subl    %eax,FPU_accum_1    /* Subtract from the num local reg */
0258     sbbl    %edx,FPU_accum_2
0259
0260 #ifdef PARANOID
0261     sbbl    $0,FPU_accum_3
0262     jne L_bugged_1  /* Must check for non-zero result here */
0263 #endif /* PARANOID */
0264
0265 /*----------------------------------------------------------------------*/
0266 /* Half of the main problem is done, there is just a reduced numerator
0267    to handle now.
0268    Work with the second 32 bits, FPU_accum_0 not used from now on */
0269 LDo_2nd_32_bits:
0270     movl    FPU_accum_2,%edx    /* get the reduced num */
0271     movl    FPU_accum_1,%eax
0272
0273     /* need to check for possible subsequent overflow */
0274     cmpl    SIGH(%ebx),%edx
0275     jb  LDo_2nd_div
0276     ja  LPrevent_2nd_overflow
0277
0278     cmpl    SIGL(%ebx),%eax
0279     jb  LDo_2nd_div
0280
0281 LPrevent_2nd_overflow:
0282 /* The numerator is greater or equal, would cause overflow */
0283     /* prevent overflow */
0284     subl    SIGL(%ebx),%eax
0285     sbbl    SIGH(%ebx),%edx
0286     movl    %edx,FPU_accum_2
0287     movl    %eax,FPU_accum_1
0288
0289     incl    FPU_result_2    /* Reflect the subtraction in the answer */
0290
0291 #ifdef PARANOID
0292     je  L_bugged_2  /* Can't bump the result to 1.0 */
0293 #endif /* PARANOID */
0294
0295 LDo_2nd_div:
0296     cmpl    $0,%ecx     /* augmented denom msw */
0297     jnz LSecond_div_not_1
0298
0299     /* %ecx == 0, we are dividing by 1.0 */
0300     mov %edx,%eax
0301     jmp LSecond_div_done
0302
0303 LSecond_div_not_1:
0304     divl    %ecx        /* Divide the numerator by the denom ms dw */
0305
0306 LSecond_div_done:
0307     movl    %eax,FPU_result_1   /* Put the result in the answer */
0308
0309     mull    SIGH(%ebx)  /* mul by the ms dw of the denom */
0310
0311     subl    %eax,FPU_accum_1    /* Subtract from the num local reg */
0312     sbbl    %edx,FPU_accum_2
0313
0314 #ifdef PARANOID
0315     jc  L_bugged_2
0316 #endif /* PARANOID */
0317
0318     movl    FPU_result_1,%eax   /* Get the result back */
0319     mull    SIGL(%ebx)  /* now mul the ls dw of the denom */
0320
0321     subl    %eax,FPU_accum_0    /* Subtract from the num local reg */
0322     sbbl    %edx,FPU_accum_1    /* Subtract from the num local reg */
0323     sbbl    $0,FPU_accum_2
0324
0325 #ifdef PARANOID
0326     jc  L_bugged_2
0327 #endif /* PARANOID */
0328
0329     jz  LDo_3rd_32_bits
0330
0331 #ifdef PARANOID
0332     cmpl    $1,FPU_accum_2
0333     jne L_bugged_2
0334 #endif /* PARANOID */
0335
0336     /* need to subtract another once of the denom */
0337     movl    SIGL(%ebx),%eax
0338     movl    SIGH(%ebx),%edx
0339     subl    %eax,FPU_accum_0    /* Subtract from the num local reg */
0340     sbbl    %edx,FPU_accum_1
0341     sbbl    $0,FPU_accum_2
0342
0343 #ifdef PARANOID
0344     jc  L_bugged_2
0345     jne L_bugged_2
0346 #endif /* PARANOID */
0347
0348     addl    $1,FPU_result_1 /* Correct the answer */
0349     adcl    $0,FPU_result_2
0350
0351 #ifdef PARANOID
0352     jc  L_bugged_2  /* Must check for non-zero result here */
0353 #endif /* PARANOID */
0354
0355 /*----------------------------------------------------------------------*/
0356 /* The division is essentially finished here, we just need to perform
0357    tidying operations.
0358    Deal with the 3rd 32 bits */
0359 LDo_3rd_32_bits:
0360     movl    FPU_accum_1,%edx        /* get the reduced num */
0361     movl    FPU_accum_0,%eax
0362
0363     /* need to check for possible subsequent overflow */
0364     cmpl    SIGH(%ebx),%edx /* denom */
0365     jb  LRound_prep
0366     ja  LPrevent_3rd_overflow
0367
0368     cmpl    SIGL(%ebx),%eax /* denom */
0369     jb  LRound_prep
0370
0371 LPrevent_3rd_overflow:
0372     /* prevent overflow */
0373     subl    SIGL(%ebx),%eax
0374     sbbl    SIGH(%ebx),%edx
0375     movl    %edx,FPU_accum_1
0376     movl    %eax,FPU_accum_0
0377
0378     addl    $1,FPU_result_1 /* Reflect the subtraction in the answer */
0379     adcl    $0,FPU_result_2
0380     jne LRound_prep
0381     jnc LRound_prep
0382
0383     /* This is a tricky spot, there is an overflow of the answer */
0384     movb    $255,FPU_ovfl_flag      /* Overflow -> 1.000 */
0385
0386 LRound_prep:
0387 /*
0388  * Prepare for rounding.
0389  * To test for rounding, we just need to compare 2*accum with the
0390  * denom.
0391  */
0392     movl    FPU_accum_0,%ecx
0393     movl    FPU_accum_1,%edx
0394     movl    %ecx,%eax
0395     orl %edx,%eax
0396     jz  LRound_ovfl     /* The accumulator contains zero. */
0397
0398     /* Multiply by 2 */
0399     clc
0400     rcll    $1,%ecx
0401     rcll    $1,%edx
0402     jc  LRound_large        /* No need to compare, denom smaller */
0403
0404     subl    SIGL(%ebx),%ecx
0405     sbbl    SIGH(%ebx),%edx
0406     jnc LRound_not_small
0407
0408     movl    $0x70000000,%eax    /* Denom was larger */
0409     jmp LRound_ovfl
0410
0411 LRound_not_small:
0412     jnz LRound_large
0413
0414     movl    $0x80000000,%eax    /* Remainder was exactly 1/2 denom */
0415     jmp LRound_ovfl
0416
0417 LRound_large:
0418     movl    $0xff000000,%eax    /* Denom was smaller */
0419
0420 LRound_ovfl:
0421 /* We are now ready to deal with rounding, but first we must get
0422    the bits properly aligned */
0423     testb   $255,FPU_ovfl_flag  /* was the num > denom ? */
0424     je  LRound_precision
0425
0426     incw    EXP(%edi)
0427
0428     /* shift the mantissa right one bit */
0429     stc         /* Will set the ms bit */
0430     rcrl    FPU_result_2
0431     rcrl    FPU_result_1
0432     rcrl    %eax
0433
0434 /* Round the result as required */
0435 LRound_precision:
0436     decw    EXP(%edi)   /* binary point between 1st & 2nd bits */
0437
0438     movl    %eax,%edx
0439     movl    FPU_result_1,%ebx
0440     movl    FPU_result_2,%eax
0441     jmp fpu_reg_round
0442
0443
0444 #ifdef PARANOID
0445 /* The logic is wrong if we got here */
0446 L_bugged:
0447     pushl   EX_INTERNAL|0x202
0448     call    EXCEPTION
0449     pop %ebx
0450     jmp L_exit
0451
0452 L_bugged_1:
0453     pushl   EX_INTERNAL|0x203
0454     call    EXCEPTION
0455     pop %ebx
0456     jmp L_exit
0457
0458 L_bugged_2:
0459     pushl   EX_INTERNAL|0x204
0460     call    EXCEPTION
0461     pop %ebx
0462     jmp L_exit
0463
0464 L_exit:
0465     movl    $-1,%eax
0466     popl    %ebx
0467     popl    %edi
0468     popl    %esi
0469
0470     leave
0471     RET
0472 #endif /* PARANOID */
0473
0474 SYM_FUNC_END(FPU_u_div)