|
|
1.1 root 1: .file "reg_u_div.S"
2: /*---------------------------------------------------------------------------+
3: | reg_u_div.S |
4: | |
5: | Core division routines |
6: | |
7: | Copyright (C) 1992 W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
8: | Australia. E-mail [email protected] |
9: | |
10: | |
11: +---------------------------------------------------------------------------*/
12:
13: /*---------------------------------------------------------------------------+
14: | Kernel for the division routines. |
15: | |
16: | void reg_u_div(unsigned long long *a, unsigned long long *a, |
17: | FPU_REG *dest) |
18: | |
19: | Does not compute the destination exponent, but does adjust it. |
20: +---------------------------------------------------------------------------*/
21:
22: #include "exception.h"
23: #include "fpu_asm.h"
24:
25:
26: #define dSIGL(x) (x)
27: #define dSIGH(x) 4(x)
28:
29:
30: .data
31: /*
32: Local storage:
33: Result: accum_3:accum_2:accum_1:accum_0
34: Overflow flag: ovfl_flag
35: */
36: .align 2,0
37: accum_3:
38: .long 0
39: accum_2:
40: .long 0
41: accum_1:
42: .long 0
43: accum_0:
44: .long 0
45: result_1:
46: .long 0
47: result_2:
48: .long 0
49: ovfl_flag:
50: .byte 0
51:
52:
53: .text
54: .align 2,144
55:
56: .globl reg_u_div
57:
58: .globl divide_kernel
59:
60: reg_u_div:
61: pushl %ebp
62: movl %esp,%ebp
63:
64: pushl %esi
65: pushl %edi
66: pushl %ebx
67:
68: movl PARAM1,%esi /* pointer to num */
69: movl PARAM2,%ebx /* pointer to denom */
70: movl PARAM3,%edi /* pointer to answer */
71:
72: divide_kernel:
73: #ifdef PARANOID
74: // testl $0x80000000, dSIGH(%esi)
75: // je xL_bugged
76: testl $0x80000000, dSIGH(%ebx)
77: je xL_bugged
78: #endif PARANOID
79:
80: /* Check if the denominator can be treated as having just 32 bits */
81: cmpl $0,dSIGL(%ebx)
82: jnz L_Full_Division /* Can't do a quick divide */
83:
84: /* We should be able to zip through the division here */
85: movl dSIGH(%ebx),%ecx /* The denominator */
86: movl dSIGH(%esi),%edx /* Get the current num */
87: movl dSIGL(%esi),%eax /* Get the current num */
88:
89: cmpl %ecx,%edx
90: setae ovfl_flag /* Keep a record */
91: jb xL_no_adjust
92:
93: subl %ecx,%edx /* Prevent the overflow */
94:
95: xL_no_adjust:
96: /* Divide the 64 bit number by the 32 bit denominator */
97: divl %ecx
98: movl %eax,SIGH(%edi) /* Put the result in the answer */
99:
100: /* Work on the remainder of the first division */
101: xorl %eax,%eax
102: divl %ecx
103: movl %eax,SIGL(%edi) /* Put the result in the answer */
104:
105: /* Work on the remainder of the 64 bit division */
106: xorl %eax,%eax
107: divl %ecx
108:
109: testb $255,ovfl_flag /* was the num > denom ? */
110: je xL_no_overflow
111:
112: /* Do the shifting here */
113: /* increase the exponent */
114: incl EXP(%edi)
115:
116: /* shift the mantissa right one bit */
117: stc /* To set the ms bit */
118: rcrl $1,SIGH(%edi)
119: rcrl $1,SIGL(%edi)
120: rcrl $1,%eax
121:
122: xL_no_overflow:
123: cmpl $0x80000000,%eax
124: jc xL_no_round
125: jnz xL_round_up
126:
127: /* "round to even" used */
128: testb $1,SIGL(%edi)
129: jz xL_no_round
130:
131: xL_round_up:
132: addl $1,SIGL(%edi)
133: adcl $0,SIGH(%edi)
134:
135: #ifdef PARANOID
136: jc xL_bugged2
137: #endif PARANOID
138:
139: xL_no_round:
140: jmp xL_done
141:
142:
143: #ifdef PARANOID
144: xL_bugged:
145: pushl EX_INTERNAL|0x203
146: call EXCEPTION
147: pop %ebx
148: jmp xL_exit
149:
150: xL_bugged2:
151: pushl EX_INTERNAL|0x204
152: call EXCEPTION
153: pop %ebx
154: jmp xL_exit
155: #endif PARANOID
156:
157:
158: /*---------------------------------------------------------------------------+
159: | Divide: Return arg1/arg2 to arg3. |
160: | |
161: | This routine does not use the exponents of arg1 and arg2, but does |
162: | adjust the exponent of arg3. |
163: | |
164: | The maximum returned value is (ignoring exponents) |
165: | .ffffffff ffffffff |
166: | ------------------ = 1.ffffffff fffffffe |
167: | .80000000 00000000 |
168: | and the minimum is |
169: | .80000000 00000000 |
170: | ------------------ = .80000000 00000001 (rounded) |
171: | .ffffffff ffffffff |
172: | |
173: +---------------------------------------------------------------------------*/
174:
175:
176: L_Full_Division:
177: movl dSIGL(%esi),%eax /* Save extended num in local register */
178: movl %eax,accum_2
179: movl dSIGH(%esi),%eax
180: movl %eax,accum_3
181: xorl %eax,%eax
182: movl %eax,accum_1 /* zero the extension */
183: movl %eax,accum_0 /* zero the extension */
184:
185: movl dSIGL(%esi),%eax /* Get the current num */
186: movl dSIGH(%esi),%edx
187:
188: /*----------------------------------------------------------------------*/
189: /* Initialization done */
190: /* Do the first 32 bits */
191:
192: movb $0,ovfl_flag
193: cmpl dSIGH(%ebx),%edx /* Test for imminent overflow */
194: jb L02
195: ja L01
196:
197: cmpl dSIGL(%ebx),%eax
198: jb L02
199:
200: L01:
201: /* The numerator is greater or equal, would cause overflow */
202: setae ovfl_flag /* Keep a record */
203:
204: subl dSIGL(%ebx),%eax
205: sbbl dSIGH(%ebx),%edx /* Prevent the overflow */
206: movl %eax,accum_2
207: movl %edx,accum_3
208:
209: L02:
210: /* At this point, we have a num < denom, with a record of
211: adjustment in ovfl_flag */
212:
213: /* We will divide by a number which is too large */
214: movl dSIGH(%ebx),%ecx
215: addl $1,%ecx
216: jnc L04
217:
218: /* here we need to divide by 100000000h,
219: i.e., no division at all.. */
220:
221: mov %edx,%eax
222: jmp L05
223:
224: L04:
225: divl %ecx /* Divide the numerator by the augmented
226: denom ms dw */
227:
228: L05:
229: movl %eax,result_2 /* Put the result in the answer */
230:
231: mull dSIGH(%ebx) /* mul by the ms dw of the denom */
232:
233: subl %eax,accum_2 /* Subtract from the num local reg */
234: sbbl %edx,accum_3
235:
236: movl result_2,%eax /* Get the result back */
237: mull dSIGL(%ebx) /* now mul the ls dw of the denom */
238:
239: subl %eax,accum_1 /* Subtract from the num local reg */
240: sbbl %edx,accum_2
241: sbbl $0,accum_3
242: je L10 /* Must check for non-zero result here */
243:
244: #ifdef PARANOID
245: jb L_bugged
246: #endif PARANOID
247:
248: /* need to subtract another once of the denom */
249: incl result_2 /* Correct the answer */
250:
251: movl dSIGL(%ebx),%eax
252: movl dSIGH(%ebx),%edx
253: subl %eax,accum_1 /* Subtract from the num local reg */
254: sbbl %edx,accum_2
255:
256: #ifdef PARANOID
257: sbbl $0,accum_3
258: jne L_bugged /* Must check for non-zero result here */
259: #endif PARANOID
260:
261: /*----------------------------------------------------------------------*/
262: /* Half of the main problem is done, there is just a reduced numerator
263: to handle now */
264: /* Work with the second 32 bits, accum_0 not used from now on */
265: L10:
266: movl accum_2,%edx /* get the reduced num */
267: movl accum_1,%eax
268:
269: /* need to check for possible subsequent overflow */
270: cmpl dSIGH(%ebx),%edx
271: jb L22
272: ja L21
273:
274: cmpl dSIGL(%ebx),%eax
275: jb L22
276:
277: L21:
278: /* The numerator is greater or equal, would cause overflow */
279: /* prevent overflow */
280: subl dSIGL(%ebx),%eax
281: sbbl dSIGH(%ebx),%edx
282: movl %edx,accum_2
283: movl %eax,accum_1
284:
285: incl result_2 /* Reflect the subtraction in the answer */
286:
287: #ifdef PARANOID
288: je L_bugged
289: #endif PARANOID
290:
291: L22:
292: cmpl $0,%ecx
293: jnz L24
294:
295: /* %ecx == 0, we are dividing by 1.0 */
296: mov %edx,%eax
297: jmp L25
298:
299: L24:
300: divl %ecx /* Divide the numerator by the denom ms dw */
301:
302: L25:
303: movl %eax,result_1 /* Put the result in the answer */
304:
305: mull dSIGH(%ebx) /* mul by the ms dw of the denom */
306:
307: subl %eax,accum_1 /* Subtract from the num local reg */
308: sbbl %edx,accum_2
309:
310: #ifdef PARANOID
311: jc L_bugged
312: #endif PARANOID
313:
314: movl result_1,%eax /* Get the result back */
315: mull dSIGL(%ebx) /* now mul the ls dw of the denom */
316:
317: /* Here we are throwing away some ls bits */
318: subl %eax,accum_0 /* Subtract from the num local reg */
319: sbbl %edx,accum_1 /* Subtract from the num local reg */
320: sbbl $0,accum_2
321:
322: #ifdef PARANOID
323: jc L_bugged
324: #endif PARANOID
325:
326: jz L35 /* Just deal with rounding now */
327:
328: #ifdef PARANOID
329: cmpl $1,accum_2
330: jne L_bugged
331: #endif PARANOID
332:
333: L32:
334: /* need to subtract another once of the denom */
335: movl dSIGL(%ebx),%eax
336: movl dSIGH(%ebx),%edx
337: subl %eax,accum_0 /* Subtract from the num local reg */
338: sbbl %edx,accum_1
339: sbbl $0,accum_2
340:
341: #ifdef PARANOID
342: jc L_bugged
343: jne L_bugged
344: #endif PARANOID
345:
346: addl $1,result_1 /* Correct the answer */
347: adcl $0,result_2
348:
349: #ifdef PARANOID
350: /* Do we ever really need this check? ***** */
351: jc L_bugged /* Must check for non-zero result here */
352: #endif PARANOID
353:
354: /*----------------------------------------------------------------------*/
355: /* The division is essentially finished here, we just need to perform
356: tidying operations. */
357: /* deal with the 3rd 32 bits */
358: L35:
359: movl accum_1,%edx /* get the reduced num */
360: movl accum_0,%eax
361:
362: /* need to check for possible subsequent overflow */
363: cmpl dSIGH(%ebx),%edx
364: jb L42
365: ja L41
366:
367: cmpl dSIGL(%ebx),%eax
368: jb L42
369:
370: L41:
371: /* prevent overflow */
372: subl dSIGL(%ebx),%eax
373: sbbl dSIGH(%ebx),%edx
374: movl %edx,accum_1
375: movl %eax,accum_0 /* ***** not needed unless extended rounding */
376:
377: addl $1,result_1 /* Reflect the subtraction in the answer */
378: adcl $0,result_2
379: jne L42
380: jnc L42
381:
382: /* This is a tricky spot, there is an overflow of the answer */
383: movb $255,ovfl_flag /* Overflow -> 1.000 */
384:
385: L42:
386: /* Now test for rounding */
387: movl accum_1,%edx /* ms byte of accum extension */
388:
389: L44:
390: cmpl $0,%ecx
391: jnz L241
392:
393: /* %ecx == 0, we are dividing by 1.0 */
394: mov %edx,%eax
395: jmp L251
396:
397: L241:
398: divl %ecx /* Divide the numerator by the denom ms dw */
399:
400: L251:
401: /* We are now ready to deal with rounding, but first we must get
402: the bits properly aligned */
403: testb $255,ovfl_flag /* was the num > denom ? */
404: je L45
405:
406: incl EXP(%edi)
407:
408: /* shift the mantissa right one bit */
409: stc // Will set the ms bit
410: rcrl $1,result_2
411: rcrl $1,result_1
412: rcrl $1,%eax
413:
414: L45:
415: cmpl $0x80000000,%eax
416: jc xL_no_round_2 // No round up
417: jnz xL_round_up_2
418:
419: /* "round to even" used here for now... */
420: testb $1,result_1
421: jz xL_no_round_2 // No round up
422:
423: xL_round_up_2:
424: addl $1,result_1
425: adcl $0,result_2
426:
427: #ifdef PARANOID
428: jc L_bugged
429: #endif PARANOID
430:
431: xL_no_round_2:
432: movl result_1,%eax
433: movl %eax,SIGL(%edi)
434: movl result_2,%eax
435: movl %eax,SIGH(%edi)
436:
437: xL_done:
438: decl EXP(%edi) /* binary point between 1st & 2nd bits */
439:
440: xL_check_exponent:
441: cmpl EXP_OVER,EXP(%edi)
442: jge xL_overflow
443:
444: cmpl EXP_UNDER,EXP(%edi)
445: jle xL_underflow
446:
447: xL_exit:
448: popl %ebx
449: popl %edi
450: popl %esi
451:
452: leave
453: ret
454:
455:
456: xL_overflow:
457: pushl %edi
458: call arith_overflow
459: popl %ebx
460: jmp xL_exit
461:
462: xL_underflow:
463: pushl %edi
464: call arith_underflow
465: popl %ebx
466: jmp xL_exit
467:
468:
469: #ifdef PARANOID
470: /* The logic is wrong if we got here */
471: L_bugged:
472: pushl EX_INTERNAL|0x202
473: call EXCEPTION
474: pop %ebx
475: jmp xL_exit
476: #endif PARANOID
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.