.file	"div_Xsig.S"
/*---------------------------------------------------------------------------+
 |  div_Xsig.S                                                               |
 |                                                                           |
 | Division subroutine for 96 bit quantities                                 |
 |                                                                           |
 | Copyright (C) 1994,1995                                                   |
 |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
 |                       Australia.  E-mail billm@jacobi.maths.monash.edu.au |
 |                                                                           |
 |                                                                           |
 +---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------+
 | Divide the 96 bit quantity pointed to by a, by that pointed to by b, and  |
 | put the 96 bit result at the location d.                                  |
 |                                                                           |
 | The result may not be accurate to 96 bits. It is intended for use where   |
 | a result better than 64 bits is required. The result should usually be    |
 | good to at least 94 bits.                                                 |
 | The returned result is actually divided by one half. This is done to      |
 | prevent overflow.                                                         |
 |                                                                           |
 |  .aaaaaaaaaaaaaa / .bbbbbbbbbbbbb  ->  .dddddddddddd                      |
 |                                                                           |
 |  void div_Xsig(Xsig *a, Xsig *b, Xsig *dest)                              |
 |                                                                           |
 +---------------------------------------------------------------------------*/

#include "exception.h"
#include "fpu_emu.h"


#define	XsigLL(x)	(x)
#define	XsigL(x)	4(x)
#define	XsigH(x)	8(x)


#ifndef NON_REENTRANT_FPU
/*
	Local storage on the stack:
	Accumulator:	FPU_accum_3:FPU_accum_2:FPU_accum_1:FPU_accum_0
 */
#define FPU_accum_3	-4(%ebp)
#define FPU_accum_2	-8(%ebp)
#define FPU_accum_1	-12(%ebp)
#define FPU_accum_0	-16(%ebp)
#define FPU_result_3	-20(%ebp)
#define FPU_result_2	-24(%ebp)
#define FPU_result_1	-28(%ebp)

#else
.data
/*
	Local storage in a static area:
	Accumulator:	FPU_accum_3:FPU_accum_2:FPU_accum_1:FPU_accum_0
 */
	.align 4,0
FPU_accum_3:
	.long	0
FPU_accum_2:
	.long	0
FPU_accum_1:
	.long	0
FPU_accum_0:
	.long	0
FPU_result_3:
	.long	0
FPU_result_2:
	.long	0
FPU_result_1:
	.long	0
#endif /* NON_REENTRANT_FPU */


.text
ENTRY(div_Xsig)
	pushl	%ebp
	movl	%esp,%ebp
#ifndef NON_REENTRANT_FPU
	subl	$28,%esp
#endif /* NON_REENTRANT_FPU */ 

	pushl	%esi
	pushl	%edi
	pushl	%ebx

	movl	PARAM1,%esi	/* pointer to num */
	movl	PARAM2,%ebx	/* pointer to denom */

#ifdef PARANOID
	testl	$0x80000000, XsigH(%ebx)	/* Divisor */
	je	L_bugged
#endif /* PARANOID */


/*---------------------------------------------------------------------------+
 |  Divide:   Return  arg1/arg2 to arg3.                                     |
 |                                                                           |
 |  The maximum returned value is (ignoring exponents)                       |
 |               .ffffffff ffffffff                                          |
 |               ------------------  =  1.ffffffff fffffffe                  |
 |               .80000000 00000000                                          |
 | and the minimum is                                                        |
 |               .80000000 00000000                                          |
 |               ------------------  =  .80000000 00000001   (rounded)       |
 |               .ffffffff ffffffff                                          |
 |                                                                           |
 +---------------------------------------------------------------------------*/

	/* Save extended dividend in local register */

	/* Divide by 2 to prevent overflow */
	clc
	movl	XsigH(%esi),%eax
	rcrl	%eax
	movl	%eax,FPU_accum_3
	movl	XsigL(%esi),%eax
	rcrl	%eax
	movl	%eax,FPU_accum_2
	movl	XsigLL(%esi),%eax
	rcrl	%eax
	movl	%eax,FPU_accum_1
	movl	$0,%eax
	rcrl	%eax
	movl	%eax,FPU_accum_0

	movl	FPU_accum_2,%eax	/* Get the current num */
	movl	FPU_accum_3,%edx

/*----------------------------------------------------------------------*/
/* Initialization done.
   Do the first 32 bits. */

	/* We will divide by a number which is too large */
	movl	XsigH(%ebx),%ecx
	addl	$1,%ecx
	jnc	LFirst_div_not_1

	/* here we need to divide by 100000000h,
	   i.e., no division at all.. */
	mov	%edx,%eax
	jmp	LFirst_div_done

LFirst_div_not_1:
	divl	%ecx		/* Divide the numerator by the augmented
				   denom ms dw */

LFirst_div_done:
	movl	%eax,FPU_result_3	/* Put the result in the answer */

	mull	XsigH(%ebx)	/* mul by the ms dw of the denom */

	subl	%eax,FPU_accum_2	/* Subtract from the num local reg */
	sbbl	%edx,FPU_accum_3

	movl	FPU_result_3,%eax	/* Get the result back */
	mull	XsigL(%ebx)	/* now mul the ls dw of the denom */

	subl	%eax,FPU_accum_1	/* Subtract from the num local reg */
	sbbl	%edx,FPU_accum_2
	sbbl	$0,FPU_accum_3
	je	LDo_2nd_32_bits		/* Must check for non-zero result here */

#ifdef PARANOID
	jb	L_bugged_1
#endif /* PARANOID */ 

	/* need to subtract another once of the denom */
	incl	FPU_result_3	/* Correct the answer */

	movl	XsigL(%ebx),%eax
	movl	XsigH(%ebx),%edx
	subl	%eax,FPU_accum_1	/* Subtract from the num local reg */
	sbbl	%edx,FPU_accum_2

#ifdef PARANOID
	sbbl	$0,FPU_accum_3
	jne	L_bugged_1	/* Must check for non-zero result here */
#endif /* PARANOID */ 

/*----------------------------------------------------------------------*/
/* Half of the main problem is done, there is just a reduced numerator
   to handle now.
   Work with the second 32 bits, FPU_accum_0 not used from now on */
LDo_2nd_32_bits:
	movl	FPU_accum_2,%edx	/* get the reduced num */
	movl	FPU_accum_1,%eax

	/* need to check for possible subsequent overflow */
	cmpl	XsigH(%ebx),%edx
	jb	LDo_2nd_div
	ja	LPrevent_2nd_overflow

	cmpl	XsigL(%ebx),%eax
	jb	LDo_2nd_div

LPrevent_2nd_overflow:
/* The numerator is greater or equal, would cause overflow */
	/* prevent overflow */
	subl	XsigL(%ebx),%eax
	sbbl	XsigH(%ebx),%edx
	movl	%edx,FPU_accum_2
	movl	%eax,FPU_accum_1

	incl	FPU_result_3	/* Reflect the subtraction in the answer */

#ifdef PARANOID
	je	L_bugged_2	/* Can't bump the result to 1.0 */
#endif /* PARANOID */ 

LDo_2nd_div:
	cmpl	$0,%ecx		/* augmented denom msw */
	jnz	LSecond_div_not_1

	/* %ecx == 0, we are dividing by 1.0 */
	mov	%edx,%eax
	jmp	LSecond_div_done

LSecond_div_not_1:
	divl	%ecx		/* Divide the numerator by the denom ms dw */

LSecond_div_done:
	movl	%eax,FPU_result_2	/* Put the result in the answer */

	mull	XsigH(%ebx)	/* mul by the ms dw of the denom */

	subl	%eax,FPU_accum_1	/* Subtract from the num local reg */
	sbbl	%edx,FPU_accum_2

#ifdef PARANOID
	jc	L_bugged_2
#endif /* PARANOID */

	movl	FPU_result_2,%eax	/* Get the result back */
	mull	XsigL(%ebx)	/* now mul the ls dw of the denom */

	subl	%eax,FPU_accum_0	/* Subtract from the num local reg */
	sbbl	%edx,FPU_accum_1	/* Subtract from the num local reg */
	sbbl	$0,FPU_accum_2

#ifdef PARANOID
	jc	L_bugged_2
#endif /* PARANOID */

	jz	LDo_3rd_32_bits

#ifdef PARANOID
	cmpl	$1,FPU_accum_2
	jne	L_bugged_2
#endif /* PARANOID */ 

	/* need to subtract another once of the denom */
	movl	XsigL(%ebx),%eax
	movl	XsigH(%ebx),%edx
	subl	%eax,FPU_accum_0	/* Subtract from the num local reg */
	sbbl	%edx,FPU_accum_1
	sbbl	$0,FPU_accum_2

#ifdef PARANOID
	jc	L_bugged_2
	jne	L_bugged_2
#endif /* PARANOID */ 

	addl	$1,FPU_result_2	/* Correct the answer */
	adcl	$0,FPU_result_3

#ifdef PARANOID
	jc	L_bugged_2	/* Must check for non-zero result here */
#endif /* PARANOID */ 

/*----------------------------------------------------------------------*/
/* The division is essentially finished here, we just need to perform
   tidying operations.
   Deal with the 3rd 32 bits */
LDo_3rd_32_bits:
	/* We use an approximation for the third 32 bits.
	To take account of the 3rd 32 bits of the divisor
	(call them del), we subtract  del * (a/b) */

	movl	FPU_result_3,%eax	/* a/b */
	mull	XsigLL(%ebx)		/* del */

	subl	%edx,FPU_accum_1

	/* A borrow indicates that the result is negative */
	jnb	LTest_over

	movl	XsigH(%ebx),%edx
	addl	%edx,FPU_accum_1

	subl	$1,FPU_result_2		/* Adjust the answer */
	sbbl	$0,FPU_result_3

	/* The above addition might not have been enough, check again. */
	movl	FPU_accum_1,%edx	/* get the reduced num */
	cmpl	XsigH(%ebx),%edx	/* denom */
	jb	LDo_3rd_div

	movl	XsigH(%ebx),%edx
	addl	%edx,FPU_accum_1

	subl	$1,FPU_result_2		/* Adjust the answer */
	sbbl	$0,FPU_result_3
	jmp	LDo_3rd_div

LTest_over:
	movl	FPU_accum_1,%edx	/* get the reduced num */

	/* need to check for possible subsequent overflow */
	cmpl	XsigH(%ebx),%edx	/* denom */
	jb	LDo_3rd_div

	/* prevent overflow */
	subl	XsigH(%ebx),%edx
	movl	%edx,FPU_accum_1

	addl	$1,FPU_result_2	/* Reflect the subtraction in the answer */
	adcl	$0,FPU_result_3

LDo_3rd_div:
	movl	FPU_accum_0,%eax
	movl	FPU_accum_1,%edx
	divl	XsigH(%ebx)

	movl    %eax,FPU_result_1       /* Rough estimate of third word */

	movl	PARAM3,%esi		/* pointer to answer */

	movl	FPU_result_1,%eax
	movl	%eax,XsigLL(%esi)
	movl	FPU_result_2,%eax
	movl	%eax,XsigL(%esi)
	movl	FPU_result_3,%eax
	movl	%eax,XsigH(%esi)

L_exit:
	popl	%ebx
	popl	%edi
	popl	%esi

	leave
	ret


#ifdef PARANOID
/* The logic is wrong if we got here */
L_bugged:
	pushl	EX_INTERNAL|0x240
	call	EXCEPTION
	pop	%ebx
	jmp	L_exit

L_bugged_1:
	pushl	EX_INTERNAL|0x241
	call	EXCEPTION
	pop	%ebx
	jmp	L_exit

L_bugged_2:
	pushl	EX_INTERNAL|0x242
	call	EXCEPTION
	pop	%ebx
	jmp	L_exit
#endif /* PARANOID */