Subversion Repositories openrisc

/* ieee754-df.S double-precision floating point support for ARM

   Copyright (C) 2003, 2004, 2005  Free Software Foundation, Inc.
   Contributed by Nicolas Pitre (nico@cam.org)

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

/*
 * Notes: 
 * 
 * The goal of this code is to be as fast as possible.  This is
 * not meant to be easy to understand for the casual reader.
 * For slightly simpler code please see the single precision version
 * of this file.
 * 
 * Only the default rounding mode is intended for best performances.
 * Exceptions aren't supported yet, but that can be added quite easily
 * if necessary without impacting performances.
 */


@ For FPA, float words are always big-endian.
@ For VFP, floats words follow the memory system mode.
#if defined(__VFP_FP__) && !defined(__ARMEB__)
#define xl r0
#define xh r1
#define yl r2
#define yh r3
#else
#define xh r0
#define xl r1
#define yh r2
#define yl r3
#endif


#ifdef L_negdf2

ARM_FUNC_START negdf2
ARM_FUNC_ALIAS aeabi_dneg negdf2

        @ flip sign bit
        eor     xh, xh, #0x80000000
        RET

        FUNC_END aeabi_dneg
        FUNC_END negdf2

#endif

#ifdef L_addsubdf3

ARM_FUNC_START aeabi_drsub

        eor     xh, xh, #0x80000000     @ flip sign bit of first arg
        b       1f      

ARM_FUNC_START subdf3
ARM_FUNC_ALIAS aeabi_dsub subdf3

        eor     yh, yh, #0x80000000     @ flip sign bit of second arg
#if defined(__INTERWORKING_STUBS__)
        b       1f                      @ Skip Thumb-code prologue
#endif

ARM_FUNC_START adddf3
ARM_FUNC_ALIAS aeabi_dadd adddf3

1:      stmfd   sp!, {r4, r5, lr}

        @ Look for zeroes, equal values, INF, or NAN.
        mov     r4, xh, lsl #1
        mov     r5, yh, lsl #1
        teq     r4, r5
        teqeq   xl, yl
        orrnes  ip, r4, xl
        orrnes  ip, r5, yl
        mvnnes  ip, r4, asr #21
        mvnnes  ip, r5, asr #21
        beq     LSYM(Lad_s)

        @ Compute exponent difference.  Make largest exponent in r4,
        @ corresponding arg in xh-xl, and positive exponent difference in r5.
        mov     r4, r4, lsr #21
        rsbs    r5, r4, r5, lsr #21
        rsblt   r5, r5, #0
        ble     1f
        add     r4, r4, r5
        eor     yl, xl, yl
        eor     yh, xh, yh
        eor     xl, yl, xl
        eor     xh, yh, xh
        eor     yl, xl, yl
        eor     yh, xh, yh
1:
        @ If exponent difference is too large, return largest argument
        @ already in xh-xl.  We need up to 54 bit to handle proper rounding
        @ of 0x1p54 - 1.1.
        cmp     r5, #54
        RETLDM  "r4, r5" hi

        @ Convert mantissa to signed integer.
        tst     xh, #0x80000000
        mov     xh, xh, lsl #12
        mov     ip, #0x00100000
        orr     xh, ip, xh, lsr #12
        beq     1f
        rsbs    xl, xl, #0
        rsc     xh, xh, #0
1:
        tst     yh, #0x80000000
        mov     yh, yh, lsl #12
        orr     yh, ip, yh, lsr #12
        beq     1f
        rsbs    yl, yl, #0
        rsc     yh, yh, #0
1:
        @ If exponent == difference, one or both args were denormalized.
        @ Since this is not common case, rescale them off line.
        teq     r4, r5
        beq     LSYM(Lad_d)
LSYM(Lad_x):

        @ Compensate for the exponent overlapping the mantissa MSB added later
        sub     r4, r4, #1

        @ Shift yh-yl right per r5, add to xh-xl, keep leftover bits into ip.
        rsbs    lr, r5, #32
        blt     1f
        mov     ip, yl, lsl lr
        adds    xl, xl, yl, lsr r5
        adc     xh, xh, #0
        adds    xl, xl, yh, lsl lr
        adcs    xh, xh, yh, asr r5
        b       2f
1:      sub     r5, r5, #32
        add     lr, lr, #32
        cmp     yl, #1
        mov     ip, yh, lsl lr
        orrcs   ip, ip, #2              @ 2 not 1, to allow lsr #1 later
        adds    xl, xl, yh, asr r5
        adcs    xh, xh, yh, asr #31
2:
        @ We now have a result in xh-xl-ip.
        @ Keep absolute value in xh-xl-ip, sign in r5 (the n bit was set above)
        and     r5, xh, #0x80000000
        bpl     LSYM(Lad_p)
        rsbs    ip, ip, #0
        rscs    xl, xl, #0
        rsc     xh, xh, #0

        @ Determine how to normalize the result.
LSYM(Lad_p):
        cmp     xh, #0x00100000
        bcc     LSYM(Lad_a)
        cmp     xh, #0x00200000
        bcc     LSYM(Lad_e)

        @ Result needs to be shifted right.
        movs    xh, xh, lsr #1
        movs    xl, xl, rrx
        mov     ip, ip, rrx
        add     r4, r4, #1

        @ Make sure we did not bust our exponent.
        mov     r2, r4, lsl #21
        cmn     r2, #(2 << 21)
        bcs     LSYM(Lad_o)

        @ Our result is now properly aligned into xh-xl, remaining bits in ip.
        @ Round with MSB of ip. If halfway between two numbers, round towards
        @ LSB of xl = 0.
        @ Pack final result together.
LSYM(Lad_e):
        cmp     ip, #0x80000000
        moveqs  ip, xl, lsr #1
        adcs    xl, xl, #0
        adc     xh, xh, r4, lsl #20
        orr     xh, xh, r5
        RETLDM  "r4, r5"

        @ Result must be shifted left and exponent adjusted.
LSYM(Lad_a):
        movs    ip, ip, lsl #1
        adcs    xl, xl, xl
        adc     xh, xh, xh
        tst     xh, #0x00100000
        sub     r4, r4, #1
        bne     LSYM(Lad_e)

        @ No rounding necessary since ip will always be 0 at this point.
LSYM(Lad_l):

#if __ARM_ARCH__ < 5

        teq     xh, #0
        movne   r3, #20
        moveq   r3, #52
        moveq   xh, xl
        moveq   xl, #0
        mov     r2, xh
        cmp     r2, #(1 << 16)
        movhs   r2, r2, lsr #16
        subhs   r3, r3, #16
        cmp     r2, #(1 << 8)
        movhs   r2, r2, lsr #8
        subhs   r3, r3, #8
        cmp     r2, #(1 << 4)
        movhs   r2, r2, lsr #4
        subhs   r3, r3, #4
        cmp     r2, #(1 << 2)
        subhs   r3, r3, #2
        sublo   r3, r3, r2, lsr #1
        sub     r3, r3, r2, lsr #3

#else

        teq     xh, #0
        moveq   xh, xl
        moveq   xl, #0
        clz     r3, xh
        addeq   r3, r3, #32
        sub     r3, r3, #11

#endif

        @ determine how to shift the value.
        subs    r2, r3, #32
        bge     2f
        adds    r2, r2, #12
        ble     1f

        @ shift value left 21 to 31 bits, or actually right 11 to 1 bits
        @ since a register switch happened above.
        add     ip, r2, #20
        rsb     r2, r2, #12
        mov     xl, xh, lsl ip
        mov     xh, xh, lsr r2
        b       3f

        @ actually shift value left 1 to 20 bits, which might also represent
        @ 32 to 52 bits if counting the register switch that happened earlier.
1:      add     r2, r2, #20
2:      rsble   ip, r2, #32
        mov     xh, xh, lsl r2
        orrle   xh, xh, xl, lsr ip
        movle   xl, xl, lsl r2

        @ adjust exponent accordingly.
3:      subs    r4, r4, r3
        addge   xh, xh, r4, lsl #20
        orrge   xh, xh, r5
        RETLDM  "r4, r5" ge

        @ Exponent too small, denormalize result.
        @ Find out proper shift value.
        mvn     r4, r4
        subs    r4, r4, #31
        bge     2f
        adds    r4, r4, #12
        bgt     1f

        @ shift result right of 1 to 20 bits, sign is in r5.
        add     r4, r4, #20
        rsb     r2, r4, #32
        mov     xl, xl, lsr r4
        orr     xl, xl, xh, lsl r2
        orr     xh, r5, xh, lsr r4
        RETLDM  "r4, r5"

        @ shift result right of 21 to 31 bits, or left 11 to 1 bits after
        @ a register switch from xh to xl.
1:      rsb     r4, r4, #12
        rsb     r2, r4, #32
        mov     xl, xl, lsr r2
        orr     xl, xl, xh, lsl r4
        mov     xh, r5
        RETLDM  "r4, r5"

        @ Shift value right of 32 to 64 bits, or 0 to 32 bits after a switch
        @ from xh to xl.
2:      mov     xl, xh, lsr r4
        mov     xh, r5
        RETLDM  "r4, r5"

        @ Adjust exponents for denormalized arguments.
        @ Note that r4 must not remain equal to 0.
LSYM(Lad_d):
        teq     r4, #0
        eor     yh, yh, #0x00100000
        eoreq   xh, xh, #0x00100000
        addeq   r4, r4, #1
        subne   r5, r5, #1
        b       LSYM(Lad_x)


LSYM(Lad_s):
        mvns    ip, r4, asr #21
        mvnnes  ip, r5, asr #21
        beq     LSYM(Lad_i)

        teq     r4, r5
        teqeq   xl, yl
        beq     1f

        @ Result is x + 0.0 = x or 0.0 + y = y.
        orrs    ip, r4, xl
        moveq   xh, yh
        moveq   xl, yl
        RETLDM  "r4, r5"

1:      teq     xh, yh

        @ Result is x - x = 0.
        movne   xh, #0
        movne   xl, #0
        RETLDM  "r4, r5" ne

        @ Result is x + x = 2x.
        movs    ip, r4, lsr #21
        bne     2f
        movs    xl, xl, lsl #1
        adcs    xh, xh, xh
        orrcs   xh, xh, #0x80000000
        RETLDM  "r4, r5"
2:      adds    r4, r4, #(2 << 21)
        addcc   xh, xh, #(1 << 20)
        RETLDM  "r4, r5" cc
        and     r5, xh, #0x80000000

        @ Overflow: return INF.
LSYM(Lad_o):
        orr     xh, r5, #0x7f000000
        orr     xh, xh, #0x00f00000
        mov     xl, #0
        RETLDM  "r4, r5"

        @ At least one of x or y is INF/NAN.
        @   if xh-xl != INF/NAN: return yh-yl (which is INF/NAN)
        @   if yh-yl != INF/NAN: return xh-xl (which is INF/NAN)
        @   if either is NAN: return NAN
        @   if opposite sign: return NAN
        @   otherwise return xh-xl (which is INF or -INF)
LSYM(Lad_i):
        mvns    ip, r4, asr #21
        movne   xh, yh
        movne   xl, yl
        mvneqs  ip, r5, asr #21
        movne   yh, xh
        movne   yl, xl
        orrs    r4, xl, xh, lsl #12
        orreqs  r5, yl, yh, lsl #12
        teqeq   xh, yh
        orrne   xh, xh, #0x00080000     @ quiet NAN
        RETLDM  "r4, r5"

        FUNC_END aeabi_dsub
        FUNC_END subdf3
        FUNC_END aeabi_dadd
        FUNC_END adddf3

ARM_FUNC_START floatunsidf
ARM_FUNC_ALIAS aeabi_ui2d floatunsidf

        teq     r0, #0
        moveq   r1, #0
        RETc(eq)
        stmfd   sp!, {r4, r5, lr}
        mov     r4, #0x400              @ initial exponent
        add     r4, r4, #(52-1 - 1)
        mov     r5, #0                  @ sign bit is 0
        .ifnc   xl, r0
        mov     xl, r0
        .endif
        mov     xh, #0
        b       LSYM(Lad_l)

        FUNC_END aeabi_ui2d
        FUNC_END floatunsidf

ARM_FUNC_START floatsidf
ARM_FUNC_ALIAS aeabi_i2d floatsidf

        teq     r0, #0
        moveq   r1, #0
        RETc(eq)
        stmfd   sp!, {r4, r5, lr}
        mov     r4, #0x400              @ initial exponent
        add     r4, r4, #(52-1 - 1)
        ands    r5, r0, #0x80000000     @ sign bit in r5
        rsbmi   r0, r0, #0              @ absolute value
        .ifnc   xl, r0
        mov     xl, r0
        .endif
        mov     xh, #0
        b       LSYM(Lad_l)

        FUNC_END aeabi_i2d
        FUNC_END floatsidf

ARM_FUNC_START extendsfdf2
ARM_FUNC_ALIAS aeabi_f2d extendsfdf2

        movs    r2, r0, lsl #1          @ toss sign bit
        mov     xh, r2, asr #3          @ stretch exponent
        mov     xh, xh, rrx             @ retrieve sign bit
        mov     xl, r2, lsl #28         @ retrieve remaining bits
        andnes  r3, r2, #0xff000000     @ isolate exponent
        teqne   r3, #0xff000000         @ if not 0, check if INF or NAN
        eorne   xh, xh, #0x38000000     @ fixup exponent otherwise.
        RETc(ne)                        @ and return it.

        teq     r2, #0                  @ if actually 0
        teqne   r3, #0xff000000         @ or INF or NAN
        RETc(eq)                        @ we are done already.

        @ value was denormalized.  We can normalize it now.
        stmfd   sp!, {r4, r5, lr}
        mov     r4, #0x380              @ setup corresponding exponent
        and     r5, xh, #0x80000000     @ move sign bit in r5
        bic     xh, xh, #0x80000000
        b       LSYM(Lad_l)

        FUNC_END aeabi_f2d
        FUNC_END extendsfdf2

ARM_FUNC_START floatundidf
ARM_FUNC_ALIAS aeabi_ul2d floatundidf

        orrs    r2, r0, r1
#if !defined (__VFP_FP__) && !defined(__SOFTFP__)
        mvfeqd  f0, #0.0
#endif
        RETc(eq)

#if !defined (__VFP_FP__) && !defined(__SOFTFP__)
        @ For hard FPA code we want to return via the tail below so that
        @ we can return the result in f0 as well as in r0/r1 for backwards
        @ compatibility.
        adr     ip, LSYM(f0_ret)
        stmfd   sp!, {r4, r5, ip, lr}
#else
        stmfd   sp!, {r4, r5, lr}
#endif

        mov     r5, #0
        b       2f

ARM_FUNC_START floatdidf
ARM_FUNC_ALIAS aeabi_l2d floatdidf

        orrs    r2, r0, r1
#if !defined (__VFP_FP__) && !defined(__SOFTFP__)
        mvfeqd  f0, #0.0
#endif
        RETc(eq)

#if !defined (__VFP_FP__) && !defined(__SOFTFP__)
        @ For hard FPA code we want to return via the tail below so that
        @ we can return the result in f0 as well as in r0/r1 for backwards
        @ compatibility.
        adr     ip, LSYM(f0_ret)
        stmfd   sp!, {r4, r5, ip, lr}
#else
        stmfd   sp!, {r4, r5, lr}
#endif

        ands    r5, ah, #0x80000000     @ sign bit in r5
        bpl     2f
        rsbs    al, al, #0
        rsc     ah, ah, #0
2:
        mov     r4, #0x400              @ initial exponent
        add     r4, r4, #(52-1 - 1)

        @ FPA little-endian: must swap the word order.
        .ifnc   xh, ah
        mov     ip, al
        mov     xh, ah
        mov     xl, ip
        .endif

        movs    ip, xh, lsr #22
        beq     LSYM(Lad_p)

        @ The value is too big.  Scale it down a bit...
        mov     r2, #3
        movs    ip, ip, lsr #3
        addne   r2, r2, #3
        movs    ip, ip, lsr #3
        addne   r2, r2, #3
        add     r2, r2, ip, lsr #3

        rsb     r3, r2, #32
        mov     ip, xl, lsl r3
        mov     xl, xl, lsr r2
        orr     xl, xl, xh, lsl r3
        mov     xh, xh, lsr r2
        add     r4, r4, r2
        b       LSYM(Lad_p)

#if !defined (__VFP_FP__) && !defined(__SOFTFP__)

        @ Legacy code expects the result to be returned in f0.  Copy it
        @ there as well.
LSYM(f0_ret):
        stmfd   sp!, {r0, r1}
        ldfd    f0, [sp], #8
        RETLDM

#endif

        FUNC_END floatdidf
        FUNC_END aeabi_l2d
        FUNC_END floatundidf
        FUNC_END aeabi_ul2d

#endif /* L_addsubdf3 */

#ifdef L_muldivdf3

ARM_FUNC_START muldf3
ARM_FUNC_ALIAS aeabi_dmul muldf3
        stmfd   sp!, {r4, r5, r6, lr}

        @ Mask out exponents, trap any zero/denormal/INF/NAN.
        mov     ip, #0xff
        orr     ip, ip, #0x700
        ands    r4, ip, xh, lsr #20
        andnes  r5, ip, yh, lsr #20
        teqne   r4, ip
        teqne   r5, ip
        bleq    LSYM(Lml_s)

        @ Add exponents together
        add     r4, r4, r5

        @ Determine final sign.
        eor     r6, xh, yh

        @ Convert mantissa to unsigned integer.
        @ If power of two, branch to a separate path.
        bic     xh, xh, ip, lsl #21
        bic     yh, yh, ip, lsl #21
        orrs    r5, xl, xh, lsl #12
        orrnes  r5, yl, yh, lsl #12
        orr     xh, xh, #0x00100000
        orr     yh, yh, #0x00100000
        beq     LSYM(Lml_1)

#if __ARM_ARCH__ < 4

        @ Put sign bit in r6, which will be restored in yl later.
        and   r6, r6, #0x80000000

        @ Well, no way to make it shorter without the umull instruction.
        stmfd   sp!, {r6, r7, r8, r9, sl, fp}
        mov     r7, xl, lsr #16
        mov     r8, yl, lsr #16
        mov     r9, xh, lsr #16
        mov     sl, yh, lsr #16
        bic     xl, xl, r7, lsl #16
        bic     yl, yl, r8, lsl #16
        bic     xh, xh, r9, lsl #16
        bic     yh, yh, sl, lsl #16
        mul     ip, xl, yl
        mul     fp, xl, r8
        mov     lr, #0
        adds    ip, ip, fp, lsl #16
        adc     lr, lr, fp, lsr #16
        mul     fp, r7, yl
        adds    ip, ip, fp, lsl #16
        adc     lr, lr, fp, lsr #16
        mul     fp, xl, sl
        mov     r5, #0
        adds    lr, lr, fp, lsl #16
        adc     r5, r5, fp, lsr #16
        mul     fp, r7, yh
        adds    lr, lr, fp, lsl #16
        adc     r5, r5, fp, lsr #16
        mul     fp, xh, r8
        adds    lr, lr, fp, lsl #16
        adc     r5, r5, fp, lsr #16
        mul     fp, r9, yl
        adds    lr, lr, fp, lsl #16
        adc     r5, r5, fp, lsr #16
        mul     fp, xh, sl
        mul     r6, r9, sl
        adds    r5, r5, fp, lsl #16
        adc     r6, r6, fp, lsr #16
        mul     fp, r9, yh
        adds    r5, r5, fp, lsl #16
        adc     r6, r6, fp, lsr #16
        mul     fp, xl, yh
        adds    lr, lr, fp
        mul     fp, r7, sl
        adcs    r5, r5, fp
        mul     fp, xh, yl
        adc     r6, r6, #0
        adds    lr, lr, fp
        mul     fp, r9, r8
        adcs    r5, r5, fp
        mul     fp, r7, r8
        adc     r6, r6, #0
        adds    lr, lr, fp
        mul     fp, xh, yh
        adcs    r5, r5, fp
        adc     r6, r6, #0
        ldmfd   sp!, {yl, r7, r8, r9, sl, fp}

#else

        @ Here is the actual multiplication.
        umull   ip, lr, xl, yl
        mov     r5, #0
        umlal   lr, r5, xh, yl
        and     yl, r6, #0x80000000
        umlal   lr, r5, xl, yh
        mov     r6, #0
        umlal   r5, r6, xh, yh

#endif

        @ The LSBs in ip are only significant for the final rounding.
        @ Fold them into lr.
        teq     ip, #0
        orrne   lr, lr, #1

        @ Adjust result upon the MSB position.
        sub     r4, r4, #0xff
        cmp     r6, #(1 << (20-11))
        sbc     r4, r4, #0x300
        bcs     1f
        movs    lr, lr, lsl #1
        adcs    r5, r5, r5
        adc     r6, r6, r6
1:
        @ Shift to final position, add sign to result.
        orr     xh, yl, r6, lsl #11
        orr     xh, xh, r5, lsr #21
        mov     xl, r5, lsl #11
        orr     xl, xl, lr, lsr #21
        mov     lr, lr, lsl #11

        @ Check exponent range for under/overflow.
        subs    ip, r4, #(254 - 1)
        cmphi   ip, #0x700
        bhi     LSYM(Lml_u)

        @ Round the result, merge final exponent.
        cmp     lr, #0x80000000
        moveqs  lr, xl, lsr #1
        adcs    xl, xl, #0
        adc     xh, xh, r4, lsl #20
        RETLDM  "r4, r5, r6"

        @ Multiplication by 0x1p*: let''s shortcut a lot of code.
LSYM(Lml_1):
        and     r6, r6, #0x80000000
        orr     xh, r6, xh
        orr     xl, xl, yl
        eor     xh, xh, yh
        subs    r4, r4, ip, lsr #1
        rsbgts  r5, r4, ip
        orrgt   xh, xh, r4, lsl #20
        RETLDM  "r4, r5, r6" gt

        @ Under/overflow: fix things up for the code below.
        orr     xh, xh, #0x00100000
        mov     lr, #0
        subs    r4, r4, #1

LSYM(Lml_u):
        @ Overflow?
        bgt     LSYM(Lml_o)

        @ Check if denormalized result is possible, otherwise return signed 0.
        cmn     r4, #(53 + 1)
        movle   xl, #0
        bicle   xh, xh, #0x7fffffff
        RETLDM  "r4, r5, r6" le

        @ Find out proper shift value.
        rsb     r4, r4, #0
        subs    r4, r4, #32
        bge     2f
        adds    r4, r4, #12
        bgt     1f

        @ shift result right of 1 to 20 bits, preserve sign bit, round, etc.
        add     r4, r4, #20
        rsb     r5, r4, #32
        mov     r3, xl, lsl r5
        mov     xl, xl, lsr r4
        orr     xl, xl, xh, lsl r5
        and     r2, xh, #0x80000000
        bic     xh, xh, #0x80000000
        adds    xl, xl, r3, lsr #31
        adc     xh, r2, xh, lsr r4
        orrs    lr, lr, r3, lsl #1
        biceq   xl, xl, r3, lsr #31
        RETLDM  "r4, r5, r6"

        @ shift result right of 21 to 31 bits, or left 11 to 1 bits after
        @ a register switch from xh to xl. Then round.
1:      rsb     r4, r4, #12
        rsb     r5, r4, #32
        mov     r3, xl, lsl r4
        mov     xl, xl, lsr r5
        orr     xl, xl, xh, lsl r4
        bic     xh, xh, #0x7fffffff
        adds    xl, xl, r3, lsr #31
        adc     xh, xh, #0
        orrs    lr, lr, r3, lsl #1
        biceq   xl, xl, r3, lsr #31
        RETLDM  "r4, r5, r6"

        @ Shift value right of 32 to 64 bits, or 0 to 32 bits after a switch
        @ from xh to xl.  Leftover bits are in r3-r6-lr for rounding.
2:      rsb     r5, r4, #32
        orr     lr, lr, xl, lsl r5
        mov     r3, xl, lsr r4
        orr     r3, r3, xh, lsl r5
        mov     xl, xh, lsr r4
        bic     xh, xh, #0x7fffffff
        bic     xl, xl, xh, lsr r4
        add     xl, xl, r3, lsr #31
        orrs    lr, lr, r3, lsl #1
        biceq   xl, xl, r3, lsr #31
        RETLDM  "r4, r5, r6"

        @ One or both arguments are denormalized.
        @ Scale them leftwards and preserve sign bit.
LSYM(Lml_d):
        teq     r4, #0
        bne     2f
        and     r6, xh, #0x80000000
1:      movs    xl, xl, lsl #1
        adc     xh, xh, xh
        tst     xh, #0x00100000
        subeq   r4, r4, #1
        beq     1b
        orr     xh, xh, r6
        teq     r5, #0
        movne   pc, lr
2:      and     r6, yh, #0x80000000
3:      movs    yl, yl, lsl #1
        adc     yh, yh, yh
        tst     yh, #0x00100000
        subeq   r5, r5, #1
        beq     3b
        orr     yh, yh, r6
        mov     pc, lr

LSYM(Lml_s):
        @ Isolate the INF and NAN cases away
        teq     r4, ip
        and     r5, ip, yh, lsr #20
        teqne   r5, ip
        beq     1f

        @ Here, one or more arguments are either denormalized or zero.
        orrs    r6, xl, xh, lsl #1
        orrnes  r6, yl, yh, lsl #1
        bne     LSYM(Lml_d)

        @ Result is 0, but determine sign anyway.
LSYM(Lml_z):
        eor     xh, xh, yh
        bic     xh, xh, #0x7fffffff
        mov     xl, #0
        RETLDM  "r4, r5, r6"

1:      @ One or both args are INF or NAN.
        orrs    r6, xl, xh, lsl #1
        moveq   xl, yl
        moveq   xh, yh
        orrnes  r6, yl, yh, lsl #1
        beq     LSYM(Lml_n)             @ 0 * INF or INF * 0 -> NAN
        teq     r4, ip
        bne     1f
        orrs    r6, xl, xh, lsl #12
        bne     LSYM(Lml_n)             @ NAN * <anything> -> NAN
1:      teq     r5, ip
        bne     LSYM(Lml_i)
        orrs    r6, yl, yh, lsl #12
        movne   xl, yl
        movne   xh, yh
        bne     LSYM(Lml_n)             @ <anything> * NAN -> NAN

        @ Result is INF, but we need to determine its sign.
LSYM(Lml_i):
        eor     xh, xh, yh

        @ Overflow: return INF (sign already in xh).
LSYM(Lml_o):
        and     xh, xh, #0x80000000
        orr     xh, xh, #0x7f000000
        orr     xh, xh, #0x00f00000
        mov     xl, #0
        RETLDM  "r4, r5, r6"

        @ Return a quiet NAN.
LSYM(Lml_n):
        orr     xh, xh, #0x7f000000
        orr     xh, xh, #0x00f80000
        RETLDM  "r4, r5, r6"

        FUNC_END aeabi_dmul
        FUNC_END muldf3

ARM_FUNC_START divdf3
ARM_FUNC_ALIAS aeabi_ddiv divdf3
        
        stmfd   sp!, {r4, r5, r6, lr}

        @ Mask out exponents, trap any zero/denormal/INF/NAN.
        mov     ip, #0xff
        orr     ip, ip, #0x700
        ands    r4, ip, xh, lsr #20
        andnes  r5, ip, yh, lsr #20
        teqne   r4, ip
        teqne   r5, ip
        bleq    LSYM(Ldv_s)

        @ Substract divisor exponent from dividend''s.
        sub     r4, r4, r5

        @ Preserve final sign into lr.
        eor     lr, xh, yh

        @ Convert mantissa to unsigned integer.
        @ Dividend -> r5-r6, divisor -> yh-yl.
        orrs    r5, yl, yh, lsl #12
        mov     xh, xh, lsl #12
        beq     LSYM(Ldv_1)
        mov     yh, yh, lsl #12
        mov     r5, #0x10000000
        orr     yh, r5, yh, lsr #4
        orr     yh, yh, yl, lsr #24
        mov     yl, yl, lsl #8
        orr     r5, r5, xh, lsr #4
        orr     r5, r5, xl, lsr #24
        mov     r6, xl, lsl #8

        @ Initialize xh with final sign bit.
        and     xh, lr, #0x80000000

        @ Ensure result will land to known bit position.
        @ Apply exponent bias accordingly.
        cmp     r5, yh
        cmpeq   r6, yl
        adc     r4, r4, #(255 - 2)
        add     r4, r4, #0x300
        bcs     1f
        movs    yh, yh, lsr #1
        mov     yl, yl, rrx
1:
        @ Perform first substraction to align result to a nibble.
        subs    r6, r6, yl
        sbc     r5, r5, yh
        movs    yh, yh, lsr #1
        mov     yl, yl, rrx
        mov     xl, #0x00100000
        mov     ip, #0x00080000

        @ The actual division loop.
1:      subs    lr, r6, yl
        sbcs    lr, r5, yh
        subcs   r6, r6, yl
        movcs   r5, lr
        orrcs   xl, xl, ip
        movs    yh, yh, lsr #1
        mov     yl, yl, rrx
        subs    lr, r6, yl
        sbcs    lr, r5, yh
        subcs   r6, r6, yl
        movcs   r5, lr
        orrcs   xl, xl, ip, lsr #1
        movs    yh, yh, lsr #1
        mov     yl, yl, rrx
        subs    lr, r6, yl
        sbcs    lr, r5, yh
        subcs   r6, r6, yl
        movcs   r5, lr
        orrcs   xl, xl, ip, lsr #2
        movs    yh, yh, lsr #1
        mov     yl, yl, rrx
        subs    lr, r6, yl
        sbcs    lr, r5, yh
        subcs   r6, r6, yl
        movcs   r5, lr
        orrcs   xl, xl, ip, lsr #3

        orrs    lr, r5, r6
        beq     2f
        mov     r5, r5, lsl #4
        orr     r5, r5, r6, lsr #28
        mov     r6, r6, lsl #4
        mov     yh, yh, lsl #3
        orr     yh, yh, yl, lsr #29
        mov     yl, yl, lsl #3
        movs    ip, ip, lsr #4
        bne     1b

        @ We are done with a word of the result.
        @ Loop again for the low word if this pass was for the high word.
        tst     xh, #0x00100000
        bne     3f
        orr     xh, xh, xl
        mov     xl, #0
        mov     ip, #0x80000000
        b       1b
2:
        @ Be sure result starts in the high word.
        tst     xh, #0x00100000
        orreq   xh, xh, xl
        moveq   xl, #0
3:
        @ Check exponent range for under/overflow.
        subs    ip, r4, #(254 - 1)
        cmphi   ip, #0x700
        bhi     LSYM(Lml_u)

        @ Round the result, merge final exponent.
        subs    ip, r5, yh
        subeqs  ip, r6, yl
        moveqs  ip, xl, lsr #1
        adcs    xl, xl, #0
        adc     xh, xh, r4, lsl #20
        RETLDM  "r4, r5, r6"

        @ Division by 0x1p*: shortcut a lot of code.
LSYM(Ldv_1):
        and     lr, lr, #0x80000000
        orr     xh, lr, xh, lsr #12
        adds    r4, r4, ip, lsr #1
        rsbgts  r5, r4, ip
        orrgt   xh, xh, r4, lsl #20
        RETLDM  "r4, r5, r6" gt

        orr     xh, xh, #0x00100000
        mov     lr, #0
        subs    r4, r4, #1
        b       LSYM(Lml_u)

        @ Result mightt need to be denormalized: put remainder bits
        @ in lr for rounding considerations.
LSYM(Ldv_u):
        orr     lr, r5, r6
        b       LSYM(Lml_u)

        @ One or both arguments is either INF, NAN or zero.
LSYM(Ldv_s):
        and     r5, ip, yh, lsr #20
        teq     r4, ip
        teqeq   r5, ip
        beq     LSYM(Lml_n)             @ INF/NAN / INF/NAN -> NAN
        teq     r4, ip
        bne     1f
        orrs    r4, xl, xh, lsl #12
        bne     LSYM(Lml_n)             @ NAN / <anything> -> NAN
        teq     r5, ip
        bne     LSYM(Lml_i)             @ INF / <anything> -> INF
        mov     xl, yl
        mov     xh, yh
        b       LSYM(Lml_n)             @ INF / (INF or NAN) -> NAN
1:      teq     r5, ip
        bne     2f
        orrs    r5, yl, yh, lsl #12
        beq     LSYM(Lml_z)             @ <anything> / INF -> 0
        mov     xl, yl
        mov     xh, yh
        b       LSYM(Lml_n)             @ <anything> / NAN -> NAN
2:      @ If both are nonzero, we need to normalize and resume above.
        orrs    r6, xl, xh, lsl #1
        orrnes  r6, yl, yh, lsl #1
        bne     LSYM(Lml_d)
        @ One or both arguments are 0.
        orrs    r4, xl, xh, lsl #1
        bne     LSYM(Lml_i)             @ <non_zero> / 0 -> INF
        orrs    r5, yl, yh, lsl #1
        bne     LSYM(Lml_z)             @ 0 / <non_zero> -> 0
        b       LSYM(Lml_n)             @ 0 / 0 -> NAN

        FUNC_END aeabi_ddiv
        FUNC_END divdf3

#endif /* L_muldivdf3 */

#ifdef L_cmpdf2

@ Note: only r0 (return value) and ip are clobbered here.

ARM_FUNC_START gtdf2
ARM_FUNC_ALIAS gedf2 gtdf2
        mov     ip, #-1
        b       1f

ARM_FUNC_START ltdf2
ARM_FUNC_ALIAS ledf2 ltdf2
        mov     ip, #1
        b       1f

ARM_FUNC_START cmpdf2
ARM_FUNC_ALIAS nedf2 cmpdf2
ARM_FUNC_ALIAS eqdf2 cmpdf2
        mov     ip, #1                  @ how should we specify unordered here?

1:      str     ip, [sp, #-4]

        @ Trap any INF/NAN first.
        mov     ip, xh, lsl #1
        mvns    ip, ip, asr #21
        mov     ip, yh, lsl #1
        mvnnes  ip, ip, asr #21
        beq     3f

        @ Test for equality.
        @ Note that 0.0 is equal to -0.0.
2:      orrs    ip, xl, xh, lsl #1      @ if x == 0.0 or -0.0
        orreqs  ip, yl, yh, lsl #1      @ and y == 0.0 or -0.0
        teqne   xh, yh                  @ or xh == yh
        teqeq   xl, yl                  @ and xl == yl
        moveq   r0, #0                  @ then equal.
        RETc(eq)

        @ Clear C flag
        cmn     r0, #0

        @ Compare sign, 
        teq     xh, yh

        @ Compare values if same sign
        cmppl   xh, yh
        cmpeq   xl, yl

        @ Result:
        movcs   r0, yh, asr #31
        mvncc   r0, yh, asr #31
        orr     r0, r0, #1
        RET

        @ Look for a NAN.
3:      mov     ip, xh, lsl #1
        mvns    ip, ip, asr #21
        bne     4f
        orrs    ip, xl, xh, lsl #12
        bne     5f                      @ x is NAN
4:      mov     ip, yh, lsl #1
        mvns    ip, ip, asr #21
        bne     2b
        orrs    ip, yl, yh, lsl #12
        beq     2b                      @ y is not NAN
5:      ldr     r0, [sp, #-4]           @ unordered return code
        RET

        FUNC_END gedf2
        FUNC_END gtdf2
        FUNC_END ledf2
        FUNC_END ltdf2
        FUNC_END nedf2
        FUNC_END eqdf2
        FUNC_END cmpdf2

ARM_FUNC_START aeabi_cdrcmple

        mov     ip, r0
        mov     r0, r2
        mov     r2, ip
        mov     ip, r1
        mov     r1, r3
        mov     r3, ip
        b       6f
        
ARM_FUNC_START aeabi_cdcmpeq
ARM_FUNC_ALIAS aeabi_cdcmple aeabi_cdcmpeq

        @ The status-returning routines are required to preserve all
        @ registers except ip, lr, and cpsr.
6:      stmfd   sp!, {r0, lr}
        ARM_CALL cmpdf2
        @ Set the Z flag correctly, and the C flag unconditionally.
        cmp      r0, #0
        @ Clear the C flag if the return value was -1, indicating
        @ that the first operand was smaller than the second.
        cmnmi    r0, #0
        RETLDM   "r0"

        FUNC_END aeabi_cdcmple
        FUNC_END aeabi_cdcmpeq
        FUNC_END aeabi_cdrcmple
        
ARM_FUNC_START  aeabi_dcmpeq

        str     lr, [sp, #-8]!
        ARM_CALL aeabi_cdcmple
        moveq   r0, #1  @ Equal to.
        movne   r0, #0  @ Less than, greater than, or unordered.
        RETLDM

        FUNC_END aeabi_dcmpeq

ARM_FUNC_START  aeabi_dcmplt

        str     lr, [sp, #-8]!
        ARM_CALL aeabi_cdcmple
        movcc   r0, #1  @ Less than.
        movcs   r0, #0  @ Equal to, greater than, or unordered.
        RETLDM

        FUNC_END aeabi_dcmplt

ARM_FUNC_START  aeabi_dcmple

        str     lr, [sp, #-8]!
        ARM_CALL aeabi_cdcmple
        movls   r0, #1  @ Less than or equal to.
        movhi   r0, #0  @ Greater than or unordered.
        RETLDM

        FUNC_END aeabi_dcmple

ARM_FUNC_START  aeabi_dcmpge

        str     lr, [sp, #-8]!
        ARM_CALL aeabi_cdrcmple
        movls   r0, #1  @ Operand 2 is less than or equal to operand 1.
        movhi   r0, #0  @ Operand 2 greater than operand 1, or unordered.
        RETLDM

        FUNC_END aeabi_dcmpge

ARM_FUNC_START  aeabi_dcmpgt

        str     lr, [sp, #-8]!
        ARM_CALL aeabi_cdrcmple
        movcc   r0, #1  @ Operand 2 is less than operand 1.
        movcs   r0, #0  @ Operand 2 is greater than or equal to operand 1,
                        @ or they are unordered.
        RETLDM

        FUNC_END aeabi_dcmpgt

#endif /* L_cmpdf2 */

#ifdef L_unorddf2

ARM_FUNC_START unorddf2
ARM_FUNC_ALIAS aeabi_dcmpun unorddf2

        mov     ip, xh, lsl #1
        mvns    ip, ip, asr #21
        bne     1f
        orrs    ip, xl, xh, lsl #12
        bne     3f                      @ x is NAN
1:      mov     ip, yh, lsl #1
        mvns    ip, ip, asr #21
        bne     2f
        orrs    ip, yl, yh, lsl #12
        bne     3f                      @ y is NAN
2:      mov     r0, #0                  @ arguments are ordered.
        RET

3:      mov     r0, #1                  @ arguments are unordered.
        RET

        FUNC_END aeabi_dcmpun
        FUNC_END unorddf2

#endif /* L_unorddf2 */

#ifdef L_fixdfsi

ARM_FUNC_START fixdfsi
ARM_FUNC_ALIAS aeabi_d2iz fixdfsi

        @ check exponent range.
        mov     r2, xh, lsl #1
        adds    r2, r2, #(1 << 21)
        bcs     2f                      @ value is INF or NAN
        bpl     1f                      @ value is too small
        mov     r3, #(0xfffffc00 + 31)
        subs    r2, r3, r2, asr #21
        bls     3f                      @ value is too large

        @ scale value
        mov     r3, xh, lsl #11
        orr     r3, r3, #0x80000000
        orr     r3, r3, xl, lsr #21
        tst     xh, #0x80000000         @ the sign bit
        mov     r0, r3, lsr r2
        rsbne   r0, r0, #0
        RET

1:      mov     r0, #0
        RET

2:      orrs    xl, xl, xh, lsl #12
        bne     4f                      @ x is NAN.
3:      ands    r0, xh, #0x80000000     @ the sign bit
        moveq   r0, #0x7fffffff         @ maximum signed positive si
        RET

4:      mov     r0, #0                  @ How should we convert NAN?
        RET

        FUNC_END aeabi_d2iz
        FUNC_END fixdfsi

#endif /* L_fixdfsi */

#ifdef L_fixunsdfsi

ARM_FUNC_START fixunsdfsi
ARM_FUNC_ALIAS aeabi_d2uiz fixunsdfsi

        @ check exponent range.
        movs    r2, xh, lsl #1
        bcs     1f                      @ value is negative
        adds    r2, r2, #(1 << 21)
        bcs     2f                      @ value is INF or NAN
        bpl     1f                      @ value is too small
        mov     r3, #(0xfffffc00 + 31)
        subs    r2, r3, r2, asr #21
        bmi     3f                      @ value is too large

        @ scale value
        mov     r3, xh, lsl #11
        orr     r3, r3, #0x80000000
        orr     r3, r3, xl, lsr #21
        mov     r0, r3, lsr r2
        RET

1:      mov     r0, #0
        RET

2:      orrs    xl, xl, xh, lsl #12
        bne     4f                      @ value is NAN.
3:      mov     r0, #0xffffffff         @ maximum unsigned si
        RET

4:      mov     r0, #0                  @ How should we convert NAN?
        RET

        FUNC_END aeabi_d2uiz
        FUNC_END fixunsdfsi

#endif /* L_fixunsdfsi */

#ifdef L_truncdfsf2

ARM_FUNC_START truncdfsf2
ARM_FUNC_ALIAS aeabi_d2f truncdfsf2

        @ check exponent range.
        mov     r2, xh, lsl #1
        subs    r3, r2, #((1023 - 127) << 21)
        subcss  ip, r3, #(1 << 21)
        rsbcss  ip, ip, #(254 << 21)
        bls     2f                      @ value is out of range

1:      @ shift and round mantissa
        and     ip, xh, #0x80000000
        mov     r2, xl, lsl #3
        orr     xl, ip, xl, lsr #29
        cmp     r2, #0x80000000
        adc     r0, xl, r3, lsl #2
        biceq   r0, r0, #1
        RET

2:      @ either overflow or underflow
        tst     xh, #0x40000000
        bne     3f                      @ overflow

        @ check if denormalized value is possible
        adds    r2, r3, #(23 << 21)
        andlt   r0, xh, #0x80000000     @ too small, return signed 0.
        RETc(lt)

        @ denormalize value so we can resume with the code above afterwards.
        orr     xh, xh, #0x00100000
        mov     r2, r2, lsr #21
        rsb     r2, r2, #24
        rsb     ip, r2, #32
        movs    r3, xl, lsl ip
        mov     xl, xl, lsr r2
        orrne   xl, xl, #1              @ fold r3 for rounding considerations. 
        mov     r3, xh, lsl #11
        mov     r3, r3, lsr #11
        orr     xl, xl, r3, lsl ip
        mov     r3, r3, lsr r2
        mov     r3, r3, lsl #1
        b       1b

3:      @ chech for NAN
        mvns    r3, r2, asr #21
        bne     5f                      @ simple overflow
        orrs    r3, xl, xh, lsl #12
        movne   r0, #0x7f000000
        orrne   r0, r0, #0x00c00000
        RETc(ne)                        @ return NAN

5:      @ return INF with sign
        and     r0, xh, #0x80000000
        orr     r0, r0, #0x7f000000
        orr     r0, r0, #0x00800000
        RET

        FUNC_END aeabi_d2f
        FUNC_END truncdfsf2

#endif /* L_truncdfsf2 */