arch/i386/math-emu/div_Xsig.S

   1         .file   "div_Xsig.S"
   2 /*---------------------------------------------------------------------------+
   3  |  div_Xsig.S                                                               |
   4  |                                                                           |
   5  | Division subroutine for 96 bit quantities                                 |
   6  |                                                                           |
   7  | Copyright (C) 1994,1995                                                   |
   8  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
   9  |                       Australia.  E-mail billm@jacobi.maths.monash.edu.au |
  10  |                                                                           |
  11  |                                                                           |
  12  +---------------------------------------------------------------------------*/
  13
  14 /*---------------------------------------------------------------------------+
  15  | Divide the 96 bit quantity pointed to by a, by that pointed to by b, and  |
  16  | put the 96 bit result at the location d.                                  |
  17  |                                                                           |
  18  | The result may not be accurate to 96 bits. It is intended for use where   |
  19  | a result better than 64 bits is required. The result should usually be    |
  20  | good to at least 94 bits.                                                 |
  21  | The returned result is actually divided by one half. This is done to      |
  22  | prevent overflow.                                                         |
  23  |                                                                           |
  24  |  .aaaaaaaaaaaaaa / .bbbbbbbbbbbbb  ->  .dddddddddddd                      |
  25  |                                                                           |
  26  |  void div_Xsig(Xsig *a, Xsig *b, Xsig *dest)                              |
  27  |                                                                           |
  28  +---------------------------------------------------------------------------*/
  29
  30 #include "exception.h"
  31 #include "fpu_emu.h"
  32
  33
  34 #define XsigLL(x)       (x)
  35 #define XsigL(x)        4(x)
  36 #define XsigH(x)        8(x)
  37
  38
  39 #ifndef NON_REENTRANT_FPU
  40 /*
  41         Local storage on the stack:
  42         Accumulator:    FPU_accum_3:FPU_accum_2:FPU_accum_1:FPU_accum_0
  43  */
  44 #define FPU_accum_3     -4(%ebp)
  45 #define FPU_accum_2     -8(%ebp)
  46 #define FPU_accum_1     -12(%ebp)
  47 #define FPU_accum_0     -16(%ebp)
  48 #define FPU_result_3    -20(%ebp)
  49 #define FPU_result_2    -24(%ebp)
  50 #define FPU_result_1    -28(%ebp)
  51
  52 #else
  53 .data
  54 /*
  55         Local storage in a static area:
  56         Accumulator:    FPU_accum_3:FPU_accum_2:FPU_accum_1:FPU_accum_0
  57  */
  58         .align 4,0
  59 FPU_accum_3:
  60         .long   0
  61 FPU_accum_2:
  62         .long   0
  63 FPU_accum_1:
  64         .long   0
  65 FPU_accum_0:
  66         .long   0
  67 FPU_result_3:
  68         .long   0
  69 FPU_result_2:
  70         .long   0
  71 FPU_result_1:
  72         .long   0
  73 #endif /* NON_REENTRANT_FPU */
  74
  75
  76 .text
  77 ENTRY(div_Xsig)
  78         pushl   %ebp
  79         movl    %esp,%ebp
  80 #ifndef NON_REENTRANT_FPU
  81         subl    $28,%esp
  82 #endif /* NON_REENTRANT_FPU */
  83
  84         pushl   %esi
  85         pushl   %edi
  86         pushl   %ebx
  87
  88         movl    PARAM1,%esi     /* pointer to num */
  89         movl    PARAM2,%ebx     /* pointer to denom */
  90
  91 #ifdef PARANOID
  92         testl   $0x80000000, XsigH(%ebx)        /* Divisor */
  93         je      L_bugged
  94 #endif /* PARANOID */
  95
  96
  97 /*---------------------------------------------------------------------------+
  98  |  Divide:   Return  arg1/arg2 to arg3.                                     |
  99  |                                                                           |
 100  |  The maximum returned value is (ignoring exponents)                       |
 101  |               .ffffffff ffffffff                                          |
 102  |               ------------------  =  1.ffffffff fffffffe                  |
 103  |               .80000000 00000000                                          |
 104  | and the minimum is                                                        |
 105  |               .80000000 00000000                                          |
 106  |               ------------------  =  .80000000 00000001   (rounded)       |
 107  |               .ffffffff ffffffff                                          |
 108  |                                                                           |
 109  +---------------------------------------------------------------------------*/
 110
 111         /* Save extended dividend in local register */
 112
 113         /* Divide by 2 to prevent overflow */
 114         clc
 115         movl    XsigH(%esi),%eax
 116         rcrl    %eax
 117         movl    %eax,FPU_accum_3
 118         movl    XsigL(%esi),%eax
 119         rcrl    %eax
 120         movl    %eax,FPU_accum_2
 121         movl    XsigLL(%esi),%eax
 122         rcrl    %eax
 123         movl    %eax,FPU_accum_1
 124         movl    $0,%eax
 125         rcrl    %eax
 126         movl    %eax,FPU_accum_0
 127
 128         movl    FPU_accum_2,%eax        /* Get the current num */
 129         movl    FPU_accum_3,%edx
 130
 131 /*----------------------------------------------------------------------*/
 132 /* Initialization done.
 133    Do the first 32 bits. */
 134
 135         /* We will divide by a number which is too large */
 136         movl    XsigH(%ebx),%ecx
 137         addl    $1,%ecx
 138         jnc     LFirst_div_not_1
 139
 140         /* here we need to divide by 100000000h,
 141            i.e., no division at all.. */
 142         mov     %edx,%eax
 143         jmp     LFirst_div_done
 144
 145 LFirst_div_not_1:
 146         divl    %ecx            /* Divide the numerator by the augmented
 147                                    denom ms dw */
 148
 149 LFirst_div_done:
 150         movl    %eax,FPU_result_3       /* Put the result in the answer */
 151
 152         mull    XsigH(%ebx)     /* mul by the ms dw of the denom */
 153
 154         subl    %eax,FPU_accum_2        /* Subtract from the num local reg */
 155         sbbl    %edx,FPU_accum_3
 156
 157         movl    FPU_result_3,%eax       /* Get the result back */
 158         mull    XsigL(%ebx)     /* now mul the ls dw of the denom */
 159
 160         subl    %eax,FPU_accum_1        /* Subtract from the num local reg */
 161         sbbl    %edx,FPU_accum_2
 162         sbbl    $0,FPU_accum_3
 163         je      LDo_2nd_32_bits         /* Must check for non-zero result here */
 164
 165 #ifdef PARANOID
 166         jb      L_bugged_1
 167 #endif /* PARANOID */
 168
 169         /* need to subtract another once of the denom */
 170         incl    FPU_result_3    /* Correct the answer */
 171
 172         movl    XsigL(%ebx),%eax
 173         movl    XsigH(%ebx),%edx
 174         subl    %eax,FPU_accum_1        /* Subtract from the num local reg */
 175         sbbl    %edx,FPU_accum_2
 176
 177 #ifdef PARANOID
 178         sbbl    $0,FPU_accum_3
 179         jne     L_bugged_1      /* Must check for non-zero result here */
 180 #endif /* PARANOID */
 181
 182 /*----------------------------------------------------------------------*/
 183 /* Half of the main problem is done, there is just a reduced numerator
 184    to handle now.
 185    Work with the second 32 bits, FPU_accum_0 not used from now on */
 186 LDo_2nd_32_bits:
 187         movl    FPU_accum_2,%edx        /* get the reduced num */
 188         movl    FPU_accum_1,%eax
 189
 190         /* need to check for possible subsequent overflow */
 191         cmpl    XsigH(%ebx),%edx
 192         jb      LDo_2nd_div
 193         ja      LPrevent_2nd_overflow
 194
 195         cmpl    XsigL(%ebx),%eax
 196         jb      LDo_2nd_div
 197
 198 LPrevent_2nd_overflow:
 199 /* The numerator is greater or equal, would cause overflow */
 200         /* prevent overflow */
 201         subl    XsigL(%ebx),%eax
 202         sbbl    XsigH(%ebx),%edx
 203         movl    %edx,FPU_accum_2
 204         movl    %eax,FPU_accum_1
 205
 206         incl    FPU_result_3    /* Reflect the subtraction in the answer */
 207
 208 #ifdef PARANOID
 209         je      L_bugged_2      /* Can't bump the result to 1.0 */
 210 #endif /* PARANOID */
 211
 212 LDo_2nd_div:
 213         cmpl    $0,%ecx         /* augmented denom msw */
 214         jnz     LSecond_div_not_1
 215
 216         /* %ecx == 0, we are dividing by 1.0 */
 217         mov     %edx,%eax
 218         jmp     LSecond_div_done
 219
 220 LSecond_div_not_1:
 221         divl    %ecx            /* Divide the numerator by the denom ms dw */
 222
 223 LSecond_div_done:
 224         movl    %eax,FPU_result_2       /* Put the result in the answer */
 225
 226         mull    XsigH(%ebx)     /* mul by the ms dw of the denom */
 227
 228         subl    %eax,FPU_accum_1        /* Subtract from the num local reg */
 229         sbbl    %edx,FPU_accum_2
 230
 231 #ifdef PARANOID
 232         jc      L_bugged_2
 233 #endif /* PARANOID */
 234
 235         movl    FPU_result_2,%eax       /* Get the result back */
 236         mull    XsigL(%ebx)     /* now mul the ls dw of the denom */
 237
 238         subl    %eax,FPU_accum_0        /* Subtract from the num local reg */
 239         sbbl    %edx,FPU_accum_1        /* Subtract from the num local reg */
 240         sbbl    $0,FPU_accum_2
 241
 242 #ifdef PARANOID
 243         jc      L_bugged_2
 244 #endif /* PARANOID */
 245
 246         jz      LDo_3rd_32_bits
 247
 248 #ifdef PARANOID
 249         cmpl    $1,FPU_accum_2
 250         jne     L_bugged_2
 251 #endif /* PARANOID */
 252
 253         /* need to subtract another once of the denom */
 254         movl    XsigL(%ebx),%eax
 255         movl    XsigH(%ebx),%edx
 256         subl    %eax,FPU_accum_0        /* Subtract from the num local reg */
 257         sbbl    %edx,FPU_accum_1
 258         sbbl    $0,FPU_accum_2
 259
 260 #ifdef PARANOID
 261         jc      L_bugged_2
 262         jne     L_bugged_2
 263 #endif /* PARANOID */
 264
 265         addl    $1,FPU_result_2 /* Correct the answer */
 266         adcl    $0,FPU_result_3
 267
 268 #ifdef PARANOID
 269         jc      L_bugged_2      /* Must check for non-zero result here */
 270 #endif /* PARANOID */
 271
 272 /*----------------------------------------------------------------------*/
 273 /* The division is essentially finished here, we just need to perform
 274    tidying operations.
 275    Deal with the 3rd 32 bits */
 276 LDo_3rd_32_bits:
 277         /* We use an approximation for the third 32 bits.
 278         To take account of the 3rd 32 bits of the divisor
 279         (call them del), we subtract  del * (a/b) */
 280
 281         movl    FPU_result_3,%eax       /* a/b */
 282         mull    XsigLL(%ebx)            /* del */
 283
 284         subl    %edx,FPU_accum_1
 285
 286         /* A borrow indicates that the result is negative */
 287         jnb     LTest_over
 288
 289         movl    XsigH(%ebx),%edx
 290         addl    %edx,FPU_accum_1
 291
 292         subl    $1,FPU_result_2         /* Adjust the answer */
 293         sbbl    $0,FPU_result_3
 294
 295         /* The above addition might not have been enough, check again. */
 296         movl    FPU_accum_1,%edx        /* get the reduced num */
 297         cmpl    XsigH(%ebx),%edx        /* denom */
 298         jb      LDo_3rd_div
 299
 300         movl    XsigH(%ebx),%edx
 301         addl    %edx,FPU_accum_1
 302
 303         subl    $1,FPU_result_2         /* Adjust the answer */
 304         sbbl    $0,FPU_result_3
 305         jmp     LDo_3rd_div
 306
 307 LTest_over:
 308         movl    FPU_accum_1,%edx        /* get the reduced num */
 309
 310         /* need to check for possible subsequent overflow */
 311         cmpl    XsigH(%ebx),%edx        /* denom */
 312         jb      LDo_3rd_div
 313
 314         /* prevent overflow */
 315         subl    XsigH(%ebx),%edx
 316         movl    %edx,FPU_accum_1
 317
 318         addl    $1,FPU_result_2 /* Reflect the subtraction in the answer */
 319         adcl    $0,FPU_result_3
 320
 321 LDo_3rd_div:
 322         movl    FPU_accum_0,%eax
 323         movl    FPU_accum_1,%edx
 324         divl    XsigH(%ebx)
 325
 326         movl    %eax,FPU_result_1       /* Rough estimate of third word */
 327
 328         movl    PARAM3,%esi             /* pointer to answer */
 329
 330         movl    FPU_result_1,%eax
 331         movl    %eax,XsigLL(%esi)
 332         movl    FPU_result_2,%eax
 333         movl    %eax,XsigL(%esi)
 334         movl    FPU_result_3,%eax
 335         movl    %eax,XsigH(%esi)
 336
 337 L_exit:
 338         popl    %ebx
 339         popl    %edi
 340         popl    %esi
 341
 342         leave
 343         ret
 344
 345
 346 #ifdef PARANOID
 347 /* The logic is wrong if we got here */
 348 L_bugged:
 349         pushl   EX_INTERNAL|0x240
 350         call    EXCEPTION
 351         pop     %ebx
 352         jmp     L_exit
 353
 354 L_bugged_1:
 355         pushl   EX_INTERNAL|0x241
 356         call    EXCEPTION
 357         pop     %ebx
 358         jmp     L_exit
 359
 360 L_bugged_2:
 361         pushl   EX_INTERNAL|0x242
 362         call    EXCEPTION
 363         pop     %ebx
 364         jmp     L_exit
 365 #endif /* PARANOID */