Subversion Repositories or1k

        .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 2,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