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 #ifndefNON_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 2,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 #ifndefNON_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 #ifdefPARANOID 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 #ifdefPARANOID 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 #ifdefPARANOID 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 #ifdefPARANOID 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 #ifdefPARANOID 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 #ifdefPARANOID 243 jc L_bugged_2
244 #endif PARANOID
245
246 jz LDo_3rd_32_bits
247
248 #ifdefPARANOID 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 #ifdefPARANOID 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 #ifdefPARANOID 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 #ifdefPARANOID 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