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