summaryrefslogtreecommitdiffstats
path: root/arch/x86/math-emu/div_Xsig.S
diff options
context:
space:
mode:
Diffstat (limited to 'arch/x86/math-emu/div_Xsig.S')
-rw-r--r--arch/x86/math-emu/div_Xsig.S365
1 files changed, 365 insertions, 0 deletions
diff --git a/arch/x86/math-emu/div_Xsig.S b/arch/x86/math-emu/div_Xsig.S
new file mode 100644
index 00000000000..f77ba3058b3
--- /dev/null
+++ b/arch/x86/math-emu/div_Xsig.S
@@ -0,0 +1,365 @@
+ .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 */