Subversion Repositories openrisc_2011-10-31

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

   Copyright (C) 2003, 2004, 2005, 2007, 2008, 2009  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 3, or (at your option) any

   later version.

   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.

   Under Section 7 of GPL version 3, you are granted additional

   permissions described in the GCC Runtime Library Exception, version

   3.1, as published by the Free Software Foundation.

   You should have received a copy of the GNU General Public License and

   a copy of the GCC Runtime Library Exception along with this program;

   see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see

   .  */

/*

 * 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_arm_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_arm_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:      do_push {r4, r5, lr}

        @ Look for zeroes, equal values, INF, or NAN.

        shift1  lsl, r4, xh, #1

        shift1  lsl, r5, yh, #1

        teq     r4, r5

        do_it   eq

        teqeq   xl, yl

        do_it   ne, ttt

        COND(orr,s,ne)  ip, r4, xl

        COND(orr,s,ne)  ip, r5, yl

        COND(mvn,s,ne)  ip, r4, asr #21

        COND(mvn,s,ne)  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.

        shift1  lsr, r4, r4, #21

        rsbs    r5, r4, r5, lsr #21

        do_it   lt

        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

        do_it   hi

        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

#if defined(__thumb2__)

        negs    xl, xl

        sbc     xh, xh, xh, lsl #1

#else

        rsbs    xl, xl, #0

        rsc     xh, xh, #0

#endif

1:

        tst     yh, #0x80000000

        mov     yh, yh, lsl #12

        orr     yh, ip, yh, lsr #12

        beq     1f

#if defined(__thumb2__)

        negs    yl, yl

        sbc     yh, yh, yh, lsl #1

#else

        rsbs    yl, yl, #0

        rsc     yh, yh, #0

#endif

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

        shift1  lsl, ip, yl, lr

        shiftop adds xl xl yl lsr r5 yl

        adc     xh, xh, #0

        shiftop adds xl xl yh lsl lr yl

        shiftop adcs xh xh yh asr r5 yh

        b       2f

1:      sub     r5, r5, #32

        add     lr, lr, #32

        cmp     yl, #1

        shift1  lsl,ip, yh, lr

        do_it   cs

        orrcs   ip, ip, #2              @ 2 not 1, to allow lsr #1 later

        shiftop adds xl xl yh asr r5 yh

        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)

#if defined(__thumb2__)

        mov     lr, #0

        negs    ip, ip

        sbcs    xl, lr, xl

        sbc     xh, lr, xh

#else

        rsbs    ip, ip, #0

        rscs    xl, xl, #0

        rsc     xh, xh, #0

#endif

        @ 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

        do_it   eq

        COND(mov,s,eq)  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

        do_it   eq, t

        moveq   xh, xl

        moveq   xl, #0

        clz     r3, xh

        do_it   eq

        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

        shift1  lsl, xl, xh, ip

        shift1  lsr, xh, xh, 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:      do_it   le

        rsble   ip, r2, #32

        shift1  lsl, xh, xh, r2

#if defined(__thumb2__)

        lsr     ip, xl, ip

        itt     le

        orrle   xh, xh, ip

        lslle   xl, xl, r2

#else

        orrle   xh, xh, xl, lsr ip

        movle   xl, xl, lsl r2

#endif

        @ adjust exponent accordingly.

3:      subs    r4, r4, r3

        do_it   ge, tt

        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

        shift1  lsr, xl, xl, r4

        shiftop orr xl xl xh lsl r2 yh

        shiftop orr xh r5 xh lsr r4 yh

        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

        shift1  lsr, xl, xl, r2

        shiftop orr xl xl xh lsl r4 yh

        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:      shift1  lsr, xl, xh, 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

        do_it   eq, te

        eoreq   xh, xh, #0x00100000

        addeq   r4, r4, #1

        subne   r5, r5, #1

        b       LSYM(Lad_x)

LSYM(Lad_s):

        mvns    ip, r4, asr #21

        do_it   ne

        COND(mvn,s,ne)  ip, r5, asr #21

        beq     LSYM(Lad_i)

        teq     r4, r5

        do_it   eq

        teqeq   xl, yl

        beq     1f

        @ Result is x + 0.0 = x or 0.0 + y = y.

        orrs    ip, r4, xl

        do_it   eq, t

        moveq   xh, yh

        moveq   xl, yl

        RETLDM  "r4, r5"

1:      teq     xh, yh

        @ Result is x - x = 0.

        do_it   ne, tt

        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

        do_it   cs

        orrcs   xh, xh, #0x80000000

        RETLDM  "r4, r5"

2:      adds    r4, r4, #(2 << 21)

        do_it   cc, t

        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

        do_it   ne, te

        movne   xh, yh

        movne   xl, yl

        COND(mvn,s,eq)  ip, r5, asr #21

        do_it   ne, t

        movne   yh, xh

        movne   yl, xl

        orrs    r4, xl, xh, lsl #12

        do_it   eq, te

        COND(orr,s,eq)  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

        do_it   eq, t

        moveq   r1, #0

        RETc(eq)

        do_push {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

        do_it   eq, t

        moveq   r1, #0

        RETc(eq)

        do_push {r4, r5, lr}

        mov     r4, #0x400              @ initial exponent

        add     r4, r4, #(52-1 - 1)

        ands    r5, r0, #0x80000000     @ sign bit in r5

        do_it   mi

        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

        do_it   ne, ttt

        COND(and,s,ne)  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

        do_it   ne, e

        teqne   r3, #0xff000000         @ or INF or NAN

        RETc(eq)                        @ we are done already.

        @ value was denormalized.  We can normalize it now.

        do_push {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__)

        do_it   eq, t

        mvfeqd  f0, #0.0

#else

        do_it   eq

#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)

        @ Push pc as well so that RETLDM works correctly.

        do_push {r4, r5, ip, lr, pc}

#else

        do_push {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__)

        do_it   eq, t

        mvfeqd  f0, #0.0

#else

        do_it   eq

#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)

        @ Push pc as well so that RETLDM works correctly.

        do_push {r4, r5, ip, lr, pc}

#else

        do_push {r4, r5, lr}

#endif

        ands    r5, ah, #0x80000000     @ sign bit in r5

        bpl     2f

#if defined(__thumb2__)

        negs    al, al

        sbc     ah, ah, ah, lsl #1

#else

        rsbs    al, al, #0

        rsc     ah, ah, #0

#endif

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

        do_it   ne

        addne   r2, r2, #3

        movs    ip, ip, lsr #3

        do_it   ne

        addne   r2, r2, #3

        add     r2, r2, ip, lsr #3

        rsb     r3, r2, #32

        shift1  lsl, ip, xl, r3

        shift1  lsr, xl, xl, r2

        shiftop orr xl xl xh lsl r3 lr

        shift1  lsr, xh, xh, 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):

        do_push {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_arm_muldivdf3

ARM_FUNC_START muldf3

ARM_FUNC_ALIAS aeabi_dmul muldf3

        do_push {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

        do_it   ne, tte

        COND(and,s,ne)  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

        do_it   ne

        COND(orr,s,ne)  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

        do_it   ne

        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)

        do_it   hi

        cmphi   ip, #0x700

        bhi     LSYM(Lml_u)

        @ Round the result, merge final exponent.

        cmp     lr, #0x80000000

        do_it   eq

        COND(mov,s,eq)  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

        do_it   gt, tt

        COND(rsb,s,gt)  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)

        do_it   le, tt

        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

        shift1  lsl, r3, xl, r5

        shift1  lsr, xl, xl, r4

        shiftop orr xl xl xh lsl r5 r2

        and     r2, xh, #0x80000000

        bic     xh, xh, #0x80000000

        adds    xl, xl, r3, lsr #31

        shiftop adc xh r2 xh lsr r4 r6

        orrs    lr, lr, r3, lsl #1

        do_it   eq

        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

        shift1  lsl, r3, xl, r4

        shift1  lsr, xl, xl, r5

        shiftop orr xl xl xh lsl r4 r2

        bic     xh, xh, #0x7fffffff

        adds    xl, xl, r3, lsr #31

        adc     xh, xh, #0

        orrs    lr, lr, r3, lsl #1

        do_it   eq

        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

        shiftop orr lr lr xl lsl r5 r2

        shift1  lsr, r3, xl, r4

        shiftop orr r3 r3 xh lsl r5 r2

        shift1  lsr, xl, xh, r4

        bic     xh, xh, #0x7fffffff

        shiftop bic xl xl xh lsr r4 r2

        add     xl, xl, r3, lsr #31

        orrs    lr, lr, r3, lsl #1

        do_it   eq

        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

        do_it   eq

        subeq   r4, r4, #1

        beq     1b

        orr     xh, xh, r6

        teq     r5, #0

        do_it   ne

        RETc(ne)

2:      and     r6, yh, #0x80000000

3:      movs    yl, yl, lsl #1

        adc     yh, yh, yh

        tst     yh, #0x00100000

        do_it   eq

        subeq   r5, r5, #1

        beq     3b

        orr     yh, yh, r6

        RET

LSYM(Lml_s):

        @ Isolate the INF and NAN cases away

        teq     r4, ip

        and     r5, ip, yh, lsr #20

        do_it   ne

        teqne   r5, ip

        beq     1f

        @ Here, one or more arguments are either denormalized or zero.

        orrs    r6, xl, xh, lsl #1

        do_it   ne

        COND(orr,s,ne)  r6, yl, yh, lsl #1

        bne     LSYM(Lml_d)

        @ Result is 0, but determine sign anyway.

LSYM(Lml_z):

        eor     xh, xh, yh

        and     xh, xh, #0x80000000

        mov     xl, #0

        RETLDM  "r4, r5, r6"

1:      @ One or both args are INF or NAN.

        orrs    r6, xl, xh, lsl #1

        do_it   eq, te

        moveq   xl, yl

        moveq   xh, yh

        COND(orr,s,ne)  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 *  -> NAN

1:      teq     r5, ip

        bne     LSYM(Lml_i)

        orrs    r6, yl, yh, lsl #12

        do_it   ne, t

        movne   xl, yl

        movne   xh, yh

        bne     LSYM(Lml_n)             @  * 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

        do_push {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

        do_it   ne, tte

        COND(and,s,ne)  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