.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 */