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
59FPU_accum_3:
60	.long	0
61FPU_accum_2:
62	.long	0
63FPU_accum_1:
64	.long	0
65FPU_accum_0:
66	.long	0
67FPU_result_3:
68	.long	0
69FPU_result_2:
70	.long	0
71FPU_result_1:
72	.long	0
73#endif /* NON_REENTRANT_FPU */
74
75
76.text
77ENTRY(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
145LFirst_div_not_1:
146	divl	%ecx		/* Divide the numerator by the augmented
147				   denom ms dw */
148
149LFirst_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 */
186LDo_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
198LPrevent_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
212LDo_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
220LSecond_div_not_1:
221	divl	%ecx		/* Divide the numerator by the denom ms dw */
222
223LSecond_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 */
276LDo_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
307LTest_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
321LDo_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
337L_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 */
348L_bugged:
349	pushl	EX_INTERNAL|0x240
350	call	EXCEPTION
351	pop	%ebx
352	jmp	L_exit
353
354L_bugged_1:
355	pushl	EX_INTERNAL|0x241
356	call	EXCEPTION
357	pop	%ebx
358	jmp	L_exit
359
360L_bugged_2:
361	pushl	EX_INTERNAL|0x242
362	call	EXCEPTION
363	pop	%ebx
364	jmp	L_exit
365#endif /* PARANOID */
366