Subversion Repositories openrisc

/* libgcc routines for the Texas Instruments TMS320C[34]x
   Copyright (C) 1997,98, 1999 Free Software Foundation, Inc.

 Contributed by Michael Hayes (m.hayes@elec.canterbury.ac.nz)
            and Herman Ten Brugge (Haj.Ten.Brugge@net.HCC.nl).

        
This file is part of GCC.

GCC is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.

In addition to the permissions in the GNU General Public License, the
Free Software Foundation gives you unlimited permission to link the
compiled version of this file into combinations with other programs,
and to distribute those combinations without any restriction coming
from the use of this file.  (The General Public License restrictions
do apply in other respects; for example, they cover modification of
the file, and distribution when not linked into a combine
executable.)

This file is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, write to
the Free Software Foundation, 51 Franklin Street, Fifth Floor,
Boston, MA 02110-1301, USA.  */

; These routines are called using the standard TI register argument
; passing model.
; The following registers do not have to be saved:
; r0, r1, r2, r3, ar0, ar1, ar2, ir0, ir1, bk, rs, rc, re, (r9, r10, r11)
;
; Perform floating point divqf3
;
; This routine performs a reciprocal of the divisor using the method
; described in the C30/C40 user manuals.  It then multiplies that
; result by the dividend.
; 
; Let r be the reciprocal of the divisor v and let the ith estimate
; of r be denoted by r[i].  An iterative approach can be used to
; improve the estimate of r, given an initial estimate r[0], where
;
; r[i + 1] = r[i] * (2.0 - v * r[i])
;
; The normalized error e[i] at the ith iteration is
;
; e[i] = (r - r[i]) / r = (1 / v - r[i]) * v = (1 - v * r[i])
;
; Note that 
;
; e[i + 1]  = (1 - v * r[i + 1]) = 1 - 2 * v * r[i] + v^2 + (r[i])^2
;           = (1 - v * r[i])^2 = (e[i])^2

; r2 dividend, r3 divisor, r0 quotient
; clobbers r1, ar1
#ifdef L_divsf3
        .text
        .global ___divqf3
___divqf3:

#ifdef _TMS320C4x
        .if .REGPARM == 0
        lda     sp,ar0
        ldf     *-ar0(2), r3
        .endif

        pop     ar1             ; Pop return address

; r0 = estimate of r, r1 = tmp, r2 = dividend, r3 = divisor
        rcpf    r3, r0          ; Compute initial estimate r[0]

        mpyf3   r0, r3, r1      ; r1 = r[0] * v
        subrf   2.0, r1         ; r1 = 2.0 - r[0] * v
        mpyf    r1, r0          ; r0 = r[0] * (2.0 - r[0] * v) = r[1]
; End of 1st iteration (16 bits accuracy)

        mpyf3   r0, r3, r1      ; r1 = r[1] * v
        subrf   2.0, r1         ; r1 = 2.0 - r[1] * v

        bud     ar1             ; Delayed branch
        mpyf    r1, r0          ; r0 = r[1] * (2.0 - r[1] * v) = r[2]
; End of 2nd iteration (32 bits accuracy)
        .if .REGPARM == 0
        mpyf    *-ar0(1), r0    ; Multiply by the dividend
        .else
        mpyf    r2, r0          ; Multiply by the dividend
        .endif
        rnd     r0
        ; Branch occurs here
#else
        .if .REGPARM == 0
        ldiu    sp,ar0
        ldf     *-ar0(2), r3
        .endif

        pop     ar1             ; Pop return address

; Initial estimate       r[0] = 1.0 * 2^(-e - 1)
; where                  v = m * 2^e

; r0 = estimate of r, r1 = tmp, r2 = dividend, r3 = divisor

; Calculate initial estimate r[0]
        pushf   r3
        pop     r0
        not     r0              ; r0 = -e
                                ; complement exponent = -e -1
                                ; complement sign (side effect)
                                ; complement mantissa (almost 3 bit accurate)
        push    r0
        popf    r0              ; r0 = 1.0 * e^(-e - 1) + inverted mantissa
        ldf     -1.0, r1        ; undo complement sign bit
        xor     r1, r0

        mpyf3   r0, r3, r1      ; r1 = r[0] * v
        subrf   2.0, r1         ; r1 = 2.0 - r[0] * v
        mpyf    r1, r0          ; r0 = r[0] * (2.0 - r[0] * v) = r[1]
; End of 1st iteration

        mpyf3   r0, r3, r1      ; r1 = r[1] * v
        subrf   2.0, r1         ; r1 = 2.0 - r[1] * v
        mpyf    r1, r0          ; r0 = r[1] * (2.0 - r[1] * v) = r[2]
; End of 2nd iteration

        mpyf3   r0, r3, r1      ; r1 = r[2] * v
        subrf   2.0, r1         ; r1 = 2.0 - r[2] * v
        mpyf    r1, r0          ; r0 = r[2] * (2.0 - r[2] * v) = r[3]
; End of 3rd iteration

        rnd     r0              ; Minimize error in x[3]'s LSBs

; Use modified last iteration
; r[4] = (r[3] * (1.0 - (v * r[3]))) + r[3]
        mpyf3   r0, r3, r1      ; r1 = r[3] * v
        subrf   1.0, r1         ; r1 = 1.0 - r[3] * v
        mpyf    r0, r1          ; r1 = r[3] * (1.0 - r[3] * v)
        addf    r1, r0          ; r0 = r[3] * (1.0 - r[3] * v) + r[3] = r[4]

        rnd     r0              ; Minimize error in x[4]'s LSBs

        bud     ar1             ; Delayed branch

        .if .REGPARM == 0
        ldfu    *-ar0(1), r2    ; Dividend in mem has only 24 bits significance
        .else
        rnd     r2              ; Minimize error in reg dividend's LSBs
                                ; since this may have 32 bit significance
        .endif

        mpyf    r2, r0          ; Multiply by the dividend
        rnd     r0              ; Round result to 32 bits

        ; Branch occurs here
#endif

#endif
;
; Integer signed division
;
; ar2 dividend, r2 divisor, r0 quotient
; clobbers r1, r3, ar0, ar1, ir0, ir1, rc, rs, re
#ifdef L_divsi3
        .text
        .global ___divqi3
        .ref    udivqi3n
___divqi3:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldi     *-ar0(1), ar2
        ldi     *-ar0(2), r2
        .endif

        xor3    ar2, r2, r3     ; Get the sign
        absi    ar2, r0
        bvd     divq32
        ldi     r0, ar2
        absi    r2, r2
        cmpi    ar2, r2         ; Divisor > dividend?

        pop     ir1
        bhid    zero            ; If so, return 0

;
; Normalize oeprands.  Use difference exponents as shift count
; for divisor, and as repeat count for "subc"
;
        float   ar2, r1         ; Normalize dividend
        pushf   r1              ; Get as integer
        pop     ar0
        lsh     -24, ar0        ; Get exponent
        
        float   r2, r1          ; Normalize divisor
        pushf   r1              ; Get as integer
        pop     ir0
        lsh     -24, ir0        ; Get exponent

        subi    ir0, ar0        ; Get difference of exponents
        lsh     ar0, r2         ; Align divisor with dividend

;
; Do count + 1 subtracts and shifts
;
        rpts    ar0
                subc    r2, ar2

;
; Mask off the lower count+1 bits of ar2
;
        subri   31, ar0         ; Shift count is (32 - (ar0 + 1))
        lsh     ar0, ar2        ; Shift left
        negi    ar0, ar0
        lsh3    ar0, ar2, r0    ; Shift right and put result in r0

;
; Check sign and negate result if necessary
;
        bud     ir1             ; Delayed return
        negi    r0, r1          ; Negate result
        ash     -31, r3         ; Check sign
        ldinz   r1, r0          ; If set, use negative result
        ; Branch occurs here

zero:   bud     ir1             ; Delayed branch
        ldi     0, r0
        nop
        nop
        ; Branch occurs here
;
; special case where ar2 = abs(ar2) = 0x80000000.  We handle this by
; calling unsigned divide and negating the result if necessary.
;
divq32:
        push    r3              ; Save sign
        call    udivqi3n
        pop     r3
        pop     ir1
        bd      ir1
        negi    r0, r1          ; Negate result
        ash     -31, r3         ; Check sign
        ldinz   r1, r0          ; If set, use negative result
        ; Branch occurs here
#endif
;
;
; ar2 dividend, r2 divisor, r0 quotient, 
; clobbers r1, r3, ar0, ar1, ir0, ir1, rc, rs, re
#ifdef L_udivsi3
        .text
        .global ___udivqi3
        .global udivqi3n
___udivqi3:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldi     *-ar0(1), ar2
        ldi     *-ar0(2), r2
        .endif

udivqi3n:
        pop     ir1

        cmpi    ar2, r2         ; If divisor > dividend
        bhi     qzero           ; return zero
        ldi     r2, ar1         ; Store divisor in ar1

        tstb    ar2, ar2        ; Check top bit, jump if set to special handler
        bld     div_32          ; Delayed branch

;
; Get divisor exponent
;
        float   ar1, r1         ; Normalize the divisor
        pushf   r1              ; Get into int register
        pop     rc
        ; branch occurs here

        bzd     qzero           ; if (float) divisor zero, return zero

        float   ar2, r1         ; Normalize the dividend
        pushf   r1              ; Get into int register
        pop     ar0
        lsh     -24, ar0        ; Get both the exponents
        lsh     -24, rc

        subi    rc, ar0         ; Get the difference between the exponents
        lsh     ar0, ar1        ; Normalize the divisor with the dividend

;
; Do count_1 subtracts and shifts
;
        rpts    ar0
                subc    ar1, ar2

;
; mask off the lower count+1 bits
;
        subri   31, ar0         ; Shift count (31 - (ar0+1))
        bud     ir1             ; Delayed return
        lsh3    ar0, ar2, r0
        negi    ar0, ar0
        lsh     ar0, r0
        ; Branch occurs here

;
; Handle a full 32-bit dividend
;
div_32: tstb    ar1, ar1
        bld     qone            ; if divisor high bit is one, the result is one
        lsh     -24, rc
        subri   31, rc
        lsh     rc, ar1         ; Line up the divisor

;
; Now divisor and dividend are aligned.  Do first SUBC by hand, save
; of the forst quotient digit.  Then, shift divisor right rather
; than shifting dividend left.  This leaves a zero in the top bit of
; the divident
;
        ldi     1, ar0          ; Initizialize MSB of quotient
        lsh     rc, ar0         ; create a mask for MSBs
        subi    1, ar0          ; mask is (2 << count) - 1

        subi3   ar1, ar2, r1
        ldihs   r1, ar2
        ldihs   1, r1
        ldilo   0, r1
        lsh     rc, r1

        lsh     -1, ar1
        subi    1, rc
;
; do the rest of the shifts and subtracts
;
        rpts    rc
                subc    ar1, ar2

        bud     ir1
        and     ar0, ar2
        or3     r1, ar2, r0
        nop

qone:
        bud     ir1
        ldi     1, r0
        nop
        nop

qzero:
        bud     ir1
        ldi     0, r0
        nop
        nop
#endif

#ifdef L_umodsi3
        .text
        .global ___umodqi3
        .global umodqi3n
___umodqi3:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldi     *-ar0(1), ar2
        ldi     *-ar0(2), r2
        .endif

umodqi3n:
        pop     ir1             ; return address
        cmpi    ar2, r2         ; divisor > dividend ? 
        bhi     uzero           ;    if so, return dividend
        ldi     r2, ar1         ; load divisor
;
; If top bit of dividend is set, handle specially.
;
        tstb    ar2, ar2        ; check top bit
        bld     umod_32         ; get divisor exponent, then jump.
;
; Get divisor exponent by converting to float.
;
        float   ar1, r1         ; normalize divisor
        pushf   r1              ; push as float
        pop     rc              ; pop as int to get exponent
        bzd     uzero           ; if (float)divisor was zero, return
;
; 31 or less bits in dividend.  Get dividend exponent.
;
        float   ar2, r1         ; normalize dividend
        pushf   r1              ; push as float
        pop     ar0             ; pop as int to get exponent
;
; Use difference in exponents as shift count to line up MSBs.
;
        lsh     -24, rc         ; divisor exponent
        lsh     -24, ar0        ; dividend exponent
        subi    rc, ar0         ; difference
        lsh     ar0, ar1        ; shift divisor up
; 
; Do COUNT+1 subtract & shifts.
;
        rpts    ar0
                subc    ar1, ar2  
;
;  Remainder is in upper 31-COUNT bits.
;
        bud     ir1             ; delayed branch to return
        addi    1, ar0          ; shift count is COUNT+1
        negi    ar0, ar0        ; negate for right shift
        lsh3    ar0, ar2, r0    ; shift to get result
        ; Return occurs here

;
; The following code handles cases of a full 32-bit dividend.  Before
; SUBC can be used, the top bit must be cleared (otherwise SUBC can
; possibly shift a significant 1 out the top of the dividend).  This
; is accomplished by first doing a normal subtraction, then proceeding
; with SUBCs. 
;
umod_32:
;
; If the top bit of the divisor is set too, the remainder is simply
; the difference between the dividend and divisor.  Otherwise, shift 
; the divisor up to line up the MSBs.
;
        tstb    ar1, ar1        ; check divisor
        bld     uone            ; if negative, remainder is diff

        lsh     -24, rc         ; divisor exponent
        subri   31, rc          ; shift count = 31 - exp
        negi    rc, ar0         ; used later as shift count
        lsh     rc, ar1         ; shift up to line up MSBs
;
; Now MSBs are aligned.  Do first SUBC by hand using a plain subtraction.
; Then, shift divisor right rather than shifting dividend left.  This leaves
; a 0 in the top bit of the dividend.
;
        subi3   ar1, ar2, r1    ; subtract 
        ldihs   r1, ar2         ; if positive, replace dividend
        subi    1, rc           ; first iteration is done
        lsh     -1, ar1         ; shift divisor down
; 
; Do EXP subtract & shifts.
;
        rpts    rc  
                subc    ar1, ar2   
;
;  Quotient is in EXP+1 LSBs; shift remainder (in MSBs) down.
;
        bud     ir1
        lsh3    ar0, ar2, r0    ; COUNT contains -(EXP+1)
        nop
        nop
;
;  Return (dividend - divisor).
;
uone:   bud     ir1
        subi3   r2, ar2, r0  
        nop
        nop
;
;  Return dividend.
;
uzero:  bud     ir1
        ldi     ar2, r0         ; set status from result
        nop
        nop
#endif

#ifdef L_modsi3
        .text
        .global ___modqi3
        .ref umodqi3n
___modqi3:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldi     *-ar0(1), ar2
        ldi     *-ar0(2), r2
        .endif

; 
; Determine sign of result.  Get absolute value of operands.
; 
        ldi     ar2, ar0        ; sign of result same as dividend
        absi    ar2, r0         ; make dividend positive
        bvd     mod_32          ; if still negative, escape
        absi    r2, r1          ; make divisor positive
        ldi     r1, ar1         ; save in ar1       
        cmpi    r0, ar1         ; divisor > dividend ? 

        pop     ir1            ; return address
        bhid    return          ;   if so, return dividend
; 
; Normalize operands.  Use difference in exponents as shift count
; for divisor, and as repeat count for SUBC.
;
        float   r1, r1          ; normalize divisor
        pushf   r1              ; push as float 
        pop     rc              ; pop as int
        bzd     return          ; if (float)divisor was zero, return

        float   r0, r1          ; normalize dividend
        pushf   r1              ; push as float
        pop     r1              ; pop as int 

        lsh     -24, rc         ; get divisor exponent
        lsh     -24, r1         ; get dividend exponent
        subi    rc, r1          ; get difference in exponents
        lsh     r1, ar1         ; align divisor with dividend
; 
; Do COUNT+1 subtract & shifts.
;
        rpts    r1
                subc    ar1, r0
;
;  Remainder is in upper bits of R0
;
        addi    1, r1           ; shift count is -(r1+1)
        negi    r1, r1 
        lsh     r1, r0          ; shift right
;
;  Check sign and negate result if necessary.
;
return:
        bud     ir1             ; delayed branch to return
        negi    r0, r1          ; negate result
        cmpi    0, ar0          ; check sign
        ldin    r1, r0          ; if set, use negative result
        ; Return occurs here
;
; The following code handles cases of a full 32-bit dividend.  This occurs
; when R0 = abs(R0) = 080000000h.  Handle this by calling the unsigned mod
; function, then negating the result if necessary.
;
mod_32:
        push    ar0             ; remember sign
        call    umodqi3n        ; do divide

        brd     return          ; return
        pop     ar0             ; restore sign
        pop     ir1             ; return address
        nop
#endif

#ifdef L_unsfltconst
        .section .const
        .global ___unsfltconst
___unsfltconst:   .float 4294967296.0
#endif

#ifdef L_unsfltcompare
        .section .const
        .global ___unsfltcompare
___unsfltcompare: .float 2147483648.0
#endif

; Integer 32-bit signed multiplication
;
; The TMS320C3x MPYI instruction takes two 24-bit signed integers
; and produces a 48-bit signed result which is truncated to 32-bits.
;
; A 32-bit by 32-bit multiplication thus requires a number of steps.
;
; Consider the product of two 32-bit signed integers,
;
;       z = x * y
;
; where x = (b << 16) + a,  y = (d << 16) + c
;
; This can be expressed as
;
;       z = ((b << 16) + a) * ((d << 16) + c)
;
;          = ((b * d) << 32) + ((b * c + a * d) << 16) + a * c
;
; Let z = (f << 16) + e where f < (1 << 16).
;
; Since we are only interested in a 32-bit result, we can ignore the 
; (b * d) << 32 term, and thus
;
;       f = b * c + a * d,  e = a * c
;
; We can simplify things if we have some a priori knowledge of the
; operands, for example, if -32768 <= y <= 32767, then y = c and d = 0 and thus
;
;       f = b * c,  e = a * c
;
; ar2 multiplier, r2 multiplicand, r0 product
; clobbers r1, r2, r3
#ifdef L_mulsi3 
        .text
        .global ___mulqi3
___mulqi3:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldi     *-ar0(1), ar2
        ldi     *-ar0(2), r2
        .endif

        pop     ir1             ; return address
        ldi     ar2, r0         ;
        and     0ffffh, r0      ; a
        lsh     -16, ar2        ; b
        ldi     r2, r3          ; 
        and     0ffffh, r3      ; c
        mpyi    r3, ar2         ; c * b         
        lsh     -16, r2         ; d
        mpyi    r0, r2          ; a * d
        addi    ar2, r2         ; c * b + a * d
        bd      ir1             ; delayed branch to return
        lsh     16, r2          ; (c * b + a * d) << 16
        mpyi    r3, r0          ; a * c
        addi    r2, r0          ; a * c + (c * b + a * d) << 16
; branch occurs here

#endif  

;
; Integer 64 by 64 multiply
; long1 and long2 on stack
; result in r0,r1
;
#ifdef L_muldi3
        .text
        .global ___mulhi3
#ifdef _TMS320C4x
___mulhi3:
        pop     ar0
        ldi     sp,ar2
        ldi     *-ar2(1),r2
        ldi     *-ar2(3),r3
        mpyi3   r2,r3,r0
        mpyuhi3 r2,r3,r1
        mpyi    *-ar2(2),r2
        bd      ar0
        mpyi    *-ar2(0),r3
        addi    r2,r1
        addi    r3,r1
#else
___mulhi3:
        ldi     sp,ar2
        ldi     -16,rs
        ldi     *-ar2(2),ar0
        ldi     *-ar2(4),ar1
        ldi     ar0,r2
        and     0ffffh,r2
        ldi     ar1,r3
        and     0ffffh,r3
        lsh     rs,ar0
        lsh     rs,ar1

        mpyi    r2,r3,r0
        mpyi    ar0,ar1,r1
        mpyi    r2,ar1,rc
        lsh     rs,rc,re
        addi    re,r1
        lsh     16,rc
        addi    rc,r0
        addc    0,r1
        mpyi    r3,ar0,rc
        lsh     rs,rc,re
        addi    re,r1
        lsh     16,rc
        addi    rc,r0
        addc    0,r1

        ldi     *-ar2(1),ar0
        ldi     ar0,r2
        and     0ffffh,r2
        lsh     rs,ar0
        mpyi    r2,r3,rc
        addi    rc,r1
        mpyi    r2,ar1,rc
        mpyi    r3,ar0,re
        addi    re,rc
        lsh     16,rc
        addi    rc,r1

        ldi     *-ar2(2),ar0
        ldi     *-ar2(3),ar1
        ldi     ar0,r2
        and     0ffffh,r2
        ldi     ar1,r3
        and     0ffffh,r3
        lsh     rs,ar0
        lsh     rs,ar1
        mpyi    r2,r3,rc
        addi    rc,r1
        mpyi    r2,ar1,rc
        mpyi    r3,ar0,re
        pop     ar0
        bd      ar0
        addi    re,rc
        lsh     16,rc
        addi    rc,r1
#endif
#endif

;
; Integer 32 by 32 multiply highpart unsigned
; src1 in ar2
; src2 in r2
; result in r0
;
#ifdef L_umuldi3_high
        .text
        .global ___umulhi3_high
___umulhi3_high:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldi     *-ar0(1), ar2
        ldi     *-ar0(2), r2
        .endif

        ldi     -16,rs
        ldi     r2,r3
        and     0ffffh,r2
        ldi     ar2,ar1
        and     0ffffh,ar2
        lsh     rs,r3
        lsh     rs,ar1

        mpyi    ar2,r2,r1
        mpyi    ar1,r3,r0
        mpyi    ar2,r3,rc
        lsh     rs,rc,re
        addi    re,r0
        lsh     16,rc
        addi    rc,r1
        addc    0,r0
        mpyi    r2,ar1,rc
        lsh     rs,rc,re
        addi    re,r0
        pop     ar0
        bd      ar0
        lsh     16,rc
        addi    rc,r1
        addc    0,r0
#endif

;
; Integer 32 by 32 multiply highpart signed
; src1 in ar2
; src2 in r2
; result in r0
;
#ifdef L_smuldi3_high
        .text
        .global ___smulhi3_high
___smulhi3_high:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldi     *-ar0(1), ar2
        ldi     *-ar0(2), r2
        .endif

        ldi     -16,rs
        ldi     0,rc
        subi3   ar2,rc,r0
        ldi     r2,r3
        ldilt   r0,rc
        subi3   r2,rc,r0
        ldi     ar2,ar1
        tstb    ar1,ar1
        ldilt   r0,rc
        and     0ffffh,r2
        and     0ffffh,ar2
        lsh     rs,r3
        lsh     rs,ar1

        mpyi    ar2,r2,r1
        mpyi    ar1,r3,r0
        addi    rc,r0
        mpyi    ar2,r3,rc
        lsh     rs,rc,re
        addi    re,r0
        lsh     16,rc
        addi    rc,r1
        addc    0,r0
        mpyi    r2,ar1,rc
        lsh     rs,rc,re
        addi    re,r0
        pop     ar0
        bd      ar0
        lsh     16,rc
        addi    rc,r1
        addc    0,r0
#endif

;
; Integer 64 by 64 unsigned divide
; long1 and long2 on stack
; divide in r0,r1
; modulo in r2,r3
; routine takes a maximum of 64*8+23=535 cycles = 21.4 us @ 50Mhz
;
#ifdef L_udivdi3
        .text
        .global ___udivhi3
        .global ___udivide
        .global ___umodulo
        .ref udivqi3n
        .ref umodqi3n
___udivhi3:
        ldi     sp,ar2
        ldi     *-ar2(4),ar0
        ldi     *-ar2(3),ar1
        ldi     *-ar2(2),r0
        ldi     *-ar2(1),r1

___udivide:
        or      r1,ar1,r2
        bne     udiv0
        ldi     ar0,r2
        ldi     r0,ar2
        call    udivqi3n
        ldiu    0,r1
        rets

___umodulo:
        or      r1,ar1,r2
        bne     udiv0
        ldi     ar0,r2
        ldi     r0,ar2
        call    umodqi3n
        ldi     r0,r2
        ldiu    0,r3
        rets

udiv0:
        tstb    ar1,ar1
        bne     udiv1
        tstb    ar0,ar0
        bn      udiv1

        ldiu    63,rc
#ifdef _TMS320C4x
        rptbd   udivend0
        ldiu    0,r2
        addi    r0,r0
        rolc    r1
#else
        ldiu    0,r2
        addi    r0,r0
        rolc    r1
        rptb    udivend0
#endif

        rolc    r2
        subi3   ar0,r2,r3
        ldinc   r3,r2
        rolc    r0
udivend0:
        rolc    r1

        not     r0
        not     r1
        ldiu    0,r3
        rets
udiv1:
        push    r4
        push    r5
        ldiu    63,rc
        ldiu    0,r2
#ifdef _TMS320C4x
        rptbd   udivend1
        ldiu    0,r3
        addi    r0,r0
        rolc    r1
#else
        ldiu    0,r3
        addi    r0,r0
        rolc    r1
        rptb    udivend1
#endif

        rolc    r2
        rolc    r3
        subi3   ar0,r2,r4
        subb3   ar1,r3,r5
        ldinc   r4,r2
        ldinc   r5,r3
        rolc    r0
udivend1:
        rolc    r1

        not     r0
        not     r1
        pop     r5
        pop     r4
        rets
#endif

;
; Integer 64 by 64 unsigned modulo
; long1 and long2 on stack
; result in r0,r1
;
#ifdef L_umoddi3
        .text
        .global ___umodhi3
        .ref ___modulo
___umodhi3:
        ldi     sp,ar2
        ldi     *-ar2(4),ar0
        ldi     *-ar2(3),ar1
        ldi     *-ar2(2),r0
        ldi     *-ar2(1),r1
        call    ___umodulo
        pop     ar0
        bd      ar0
        ldi     r2,r0
        ldi     r3,r1
        nop
#endif

;
; Integer 64 by 64 signed divide
; long1 and long2 on stack
; result in r0,r1
;
#ifdef L_divdi3
        .text
        .global ___divhi3
        .ref ___udivide
___divhi3:
        ldi     0,ir0
        ldi     sp,ar2
        ldi     *-ar2(4),r0
        ldi     *-ar2(3),r1
        bge     div1
        not     ir0
        negi    r0
        negb    r1
div1:
        ldi     r0,ar0
        ldi     r1,ar1
        ldi     *-ar2(2),r0
        ldi     *-ar2(1),r1
        bge     div2
        not     ir0
        negi    r0
        negb    r1
div2:
        call    ___udivide
        tstb    ir0,ir0
        bge     div3
        negi    r0
        negb    r1
div3:   
        rets
#endif

;
; Integer 64 by 64 signed modulo
; long1 and long2 on stack
; result in r0,r1
;
#ifdef L_moddi3
        .text
        .global ___modhi3
        .ref ___umodulo
___modhi3:
        ldi     0,ir0
        ldi     sp,ar2
        ldi     *-ar2(4),r0
        ldi     *-ar2(3),r1
        bge     mod1
        not     ir0
        negi    r0
        negb    r1
mod1:
        ldi     r0,ar0
        ldi     r1,ar1
        ldi     *-ar2(2),r0
        ldi     *-ar2(1),r1
        bge     mod2
        not     ir0
        negi    r0
        negb    r1
mod2:
        call    ___umodulo
        ldi     r2,r0
        ldi     r3,r1
        tstb    ir0,ir0
        bge     mod3
        negi    r0
        negb    r1
mod3:   
        rets
#endif

;
; double to signed long long conversion
; input in r2
; result in r0,r1
;
#ifdef L_fix_truncsfdi2
        .text
        .global ___fix_truncqfhi2
        .ref ufix_truncqfhi2n
___fix_truncqfhi2:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldf     *-ar0(1), r2
        .endif

        cmpf    0.0,r2
        bge     ufix_truncqfhi2n
        negf    r2
        call    ufix_truncqfhi2n
        negi    r0
        negb    r1
        rets
#endif

;
; double to unsigned long long conversion
; input in r2
; result in r0,r1
;
#ifdef L_ufix_truncsfdi2
        .text
        .global ___ufix_truncqfhi2
        .global ufix_truncqfhi2n
___ufix_truncqfhi2:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldf     *-ar0(1), r2
        .endif

ufix_truncqfhi2n:
        cmpf    0.0,r2
        ble     ufix1
        pushf   r2
        pop     r3
        ash     -24,r3
        subi    31,r3
        cmpi    32,r3
        bgt     ufix1
        cmpi    -32,r3
        ble     ufix1
        ldi     1,r0
        ash     31,r0
        or3     r0,r2,r0
        ldi     r0,r1
        lsh3    r3,r0,r0
        subi    32,r3
        cmpi    -32,r3
        ldile   0,r1
        lsh3    r3,r1,r1
        rets
ufix1:
        ldi     0,r0
        ldi     0,r1
        rets
#endif

;
; signed long long to double conversion
; input on stack
; result in r0
;
#ifdef L_floatdisf2
        .text
        .global ___floathiqf2
        .ref ufloathiqf2n
___floathiqf2:
        ldi     sp,ar2
        ldi     *-ar2(2),r0
        ldi     *-ar2(1),r1
        bge     ufloathiqf2n
        negi    r0
        negb    r1
        call    ufloathiqf2n
        negf    r0
        rets
#endif

;
; unsigned long long to double conversion
; input on stack
; result in r0
;
#ifdef L_ufloatdisf2
        .text
        .global ___ufloathiqf2
        .global ufloathiqf2n
        .ref ___unsfltconst
___ufloathiqf2:
        ldi     sp,ar2
        ldi     *-ar2(2),r0
        ldi     *-ar2(1),r1
ufloathiqf2n:
        .if .BIGMODEL
#ifdef _TMS320C4x
        ldpk    @___unsfltconst
#else
        ldp     @___unsfltconst
#endif
        .endif
        ldf     @___unsfltconst,r2
        float   r0
        bge     uflt1
        addf    r2,r0
uflt1:
        float   r1
        bge     uflt2
        addf    r2,r1
uflt2:
#ifdef _TMS320C4x
        pop     r3
        bd      r3
        mpyf    r2,r1
        addf    r1,r0
        nop
#else
        ldf     r1,r3
        and     0ffh,r3
        norm    r3,r3
        mpyf    r2,r3
        pop     ar2
        bd      ar2
        addf    r3,r0
        mpyf    r2,r1
        addf    r1,r0
#endif
#endif

;
; long double to signed long long conversion
; input in r2
; result in r0,r1
;
#ifdef L_fix_truncdfdi2
        .text
        .global ___fix_trunchfhi2
        .ref ufix_trunchfhi2n
___fix_trunchfhi2:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldf     *-ar0(2), r2
        ldi     *-ar0(1), r2
        .endif

        cmpf    0.0,r2
        bge     ufix_trunchfhi2n
        negf    r2
        call    ufix_trunchfhi2n
        negi    r0
        negb    r1
        rets
#endif

;
; long double to unsigned long long conversion
; input in r2
; result in r0,r1
;
#ifdef L_ufix_truncdfdi2
        .text
        .global ___ufix_trunchfhi2
        .global ufix_trunchfhi2n
___ufix_trunchfhi2:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldf     *-ar0(2), r2
        ldi     *-ar0(1), r2
        .endif

ufix_trunchfhi2n:
        cmpf    0.0,r2
        ble     ufixh1
        pushf   r2
        pop     r3
        ash     -24,r3
        subi    31,r3
        cmpi    32,r3
        bgt     ufixh1
        cmpi    -32,r3
        ble     ufixh1
        ldi     1,r0
        ash     31,r0
        or3     r0,r2,r0
        ldi     r0,r1
        lsh3    r3,r0,r0
        subi    32,r3
        cmpi    -32,r3
        ldile   0,r1
        lsh3    r3,r1,r1
        rets
ufixh1:
        ldi     0,r0
        ldi     0,r1
        rets
#endif

;
; signed long long to long double conversion
; input on stack
; result in r0
;
#ifdef L_floatdidf2
        .text
        .global ___floathihf2
        .ref ufloathihf2n
___floathihf2:
        ldi     sp,ar2
        ldi     *-ar2(2),r0
        ldi     *-ar2(1),r1
        bge     ufloathihf2n
        negi    r0
        negb    r1
        call    ufloathihf2n
        negf    r0
        rets
#endif

;
; unsigned long long to double conversion
; input on stack
; result in r0
;
#ifdef L_ufloatdidf2
        .text
        .global ___ufloathihf2
        .global ufloathihf2n
        .ref ___unsfltconst
___ufloathihf2:
        ldi     sp,ar2
        ldi     *-ar2(2),r0
        ldi     *-ar2(1),r1
ufloathihf2n
        .if .BIGMODEL
#ifdef _TMS320C4x
        ldpk    @___unsfltconst
#else
        ldp     @___unsfltconst
#endif
        .endif
        ldf     @___unsfltconst,r2
        float   r0
        bge     uflth1
        addf    r2,r0
uflth1:
        float   r1
        bge     uflth2
        addf    r2,r1
uflth2:
#ifdef _TMS320C4x
        pop     r3
        bd      r3
        mpyf    r2,r1
        addf    r1,r0
        nop
#else
        ldf     r1,r3
        and     0ffh,r3
        norm    r3,r3
        mpyf    r2,r3
        pop     ar2
        bd      ar2
        addf    r3,r0
        mpyf    r2,r1
        addf    r1,r0
#endif
#endif

;
; calculate ffs
; input in ar2
; result in r0
;
#ifdef L_ffs
        .global ___ffs
        .ref ___unsfltconst
        .text
___ffs:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldi     *-ar0(1), ar2
        .endif

        negi    ar2,r0
        and     ar2,r0
        float   r0,r0
        ldfu    0.0,r1
        .if .BIGMODEL
#ifdef _TMS320C4x
        ldpk    @___unsfltconst
#else
        ldp     @___unsfltconst
#endif
        .endif
        ldflt   @___unsfltconst,r1
        addf    r1,r0
        pushf   r0
        pop     r0
        pop     ar0
        bd      ar0
        ash     -24,r0
        ldilt   -1,r0
        addi    1,r0
#endif

;
; calculate long double * long double
; input in r2, r3
; output in r0
;
#ifdef L_muldf3
        .global ___mulhf3
        .text
___mulhf3:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldf     *-ar0(2), r2
        ldi     *-ar0(1), r2
        ldf     *-ar0(4), r3
        ldi     *-ar0(3), r3
        .endif

        pop     ar2             ; return ad
        ldf     r2,r0           ; copy lsb0
        ldf     r3,r1           ; copy lsb1
        and     0ffh,r0         ; mask lsb0
        and     0ffh,r1         ; mask lsb1
        norm    r0,r0           ; correct lsb0
        norm    r1,r1           ; correct lsb1
        mpyf    r2,r1           ; arg0*lsb1
        mpyf    r3,r0           ; arg1*lsb0
        bd      ar2             ; return (delayed)
        addf    r0,r1           ; arg0*lsb1 + arg1*lsb0
        mpyf    r2,r3,r0        ; msb0*msb1
        addf    r1,r0           ; msb0*msb1 + arg0*lsb1 + arg1*lsb0
#endif

;
; calculate long double / long double
; r2 dividend, r3 divisor, r0 quotient
;
#ifdef L_divdf3
        .global ___divhf3
        .text
___divhf3:
        .if .REGPARM == 0
#ifdef _TMS320C4x
        lda     sp,ar0
#else
        ldiu    sp,ar0
#endif
        ldf     *-ar0(2), r2
        ldi     *-ar0(1), r2
        ldf     *-ar0(4), r3
        ldi     *-ar0(3), r3
        .endif

#ifdef _TMS320C4x
        pop     ar1
        rcpf    r3, r0
        mpyf3   r0, r3, r1
        subrf   2.0, r1         
        mpyf    r1, r0  
        mpyf3   r0, r3, r1
        bud     ar1
        subrf   2.0, r1 
        mpyf    r1, r0
        mpyf    r2, r0
#else
        pop     ar1
        pushf   r3
        pop     r0
        not     r0      
        push    r0
        popf    r0
        ldf     -1.0, r1
        xor     r1, r0

        mpyf3   r0, r3, r1      ; r1 = r[0] * v
        subrf   2.0, r1         ; r1 = 2.0 - r[0] * v
        mpyf    r1, r0          ; r0 = r[0] * (2.0 - r[0] * v) = r[1]
; End of 1st iteration

        mpyf3   r0, r3, r1      ; r1 = r[1] * v
        subrf   2.0, r1         ; r1 = 2.0 - r[1] * v
        mpyf    r1, r0          ; r0 = r[1] * (2.0 - r[1] * v) = r[2]
; End of 2nd iteration

        mpyf3   r0, r3, r1      ; r1 = r[2] * v
        subrf   2.0, r1         ; r1 = 2.0 - r[2] * v
        mpyf    r1, r0          ; r0 = r[2] * (2.0 - r[2] * v) = r[3]
; End of 3rd iteration

        or      080h, r0
        rnd     r0

;       mpyf3   r0, r3, r1      ; r1 = r[3] * v
        push    r4
        pushf   r4
        mpyf    r0, r3, r1

        ldf     r0, r4
        and     0ffh, r4
        norm    r4, r4
        mpyf    r3, r4
        addf    r4, r1

        ldf     r3, r4
        and     0ffh, r4
        norm    r4, r4
        mpyf    r0, r4
        addf    r4, r1
        
        subrf   2.0, r1         ; r1 = 2.0 - r[3] * v

        mpyf    r1, r0, r3      ; r3 = r[3] * (2.0 - r[3] * v) = r[5]

        ldf     r1, r4
        and     0ffh, r4
        norm    r4, r4
        mpyf    r0, r4
        addf    r4, r3

        ldf     r0, r4
        and     0ffh, r4
        norm    r4, r4
        mpyf    r1, r4
        addf    r4, r3

        mpyf    r2, r3, r0      ; Multiply by the dividend

        ldf     r2, r4
        and     0ffh, r4
        norm    r4, r4
        mpyf    r3, r4
        addf    r4, r0

        ldf     r3, r4
        and     0ffh, r4
        norm    r4, r4
        mpyf    r2, r4
        bd      ar1
        addf    r4, r0

        popf    r4
        pop     r4
#endif
#endif