252 lines · plain
1//===-- mulsf3.S - single-precision floating point multiplication ---------===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8//9// This file implements single-precision soft-float multiplication with the10// IEEE-754 default rounding (to nearest, ties to even), in optimized Thumb111// assembly language.12//13//===----------------------------------------------------------------------===//14 15#include "../../assembly.h"16 17 .syntax unified18 .text19 .thumb20 .p2align 221 22DEFINE_AEABI_FUNCTION_ALIAS(__aeabi_fmul, __mulsf3)23 24DEFINE_COMPILERRT_THUMB_FUNCTION(__mulsf3)25 push {r4,r5,r6,lr}26 27 // Get exponents of the inputs, and check for uncommon values. In the process28 // of this we also compute the sign, because it's marginally quicker that29 // way.30 lsls r2, r0, #131 adcs r4, r4, r4 // set r4[0] to sign bit of x32 lsls r3, r1, #133 adcs r4, r4, r3 // set r4[0] to the output sign34 lsrs r2, r2, #2435 beq LOCAL_LABEL(zerodenorm0) // still do the next LSRS36 lsrs r3, r3, #2437 beq LOCAL_LABEL(zerodenorm)38 cmp r2, #25539 beq LOCAL_LABEL(naninf)40 cmp r3, #25541 beq LOCAL_LABEL(naninf)42 // Compute the output exponent. We'll be generating our product _without_ the43 // leading bit, so we subtract 0x7f rather than 0x80.44 adds r2, r2, r345 subs r2, r2, #0x7f46 // Blank off everything above the mantissas.47 lsls r0, r0, #948 lsls r1, r1, #949LOCAL_LABEL(normalised): // we may come back here from zerodenorm50 lsrs r0, r0, #951 lsrs r1, r1, #952 // Multiply. r0 and r1 are the mantissas of the inputs but without their53 // leading bits, so the product we want in principle is P=(r0+2^23)(r1+2^23).54 // P is at most (2^24-1)^2 < 2^48, so it fits in a word and a half.55 //56 // The technique below will actually compute P - 2^46, by not adding on the57 // term where the two 2^23 are multiplied. The 48-bit result will be58 // delivered in two output registers, one containing its bottom 32 bits and59 // the other containing the top 32, so they overlap in the middle 16 bits.60 // This is done using only two multiply instructions and some bookkeeping.61 //62 // In the comments I'll write X and Y for the original input mantissas (again63 // without their leading bits). I'll also decompose them as X = xh + xl and64 // Y = yh + yl, where xl and yl are in the range 0..2^8-1 and xh,yh are65 // multiples of 2^8.66 adds r5, r0, r167 lsls r5, r5, #7 // r5 = (X+Y) << 768 movs r6, r069 muls r6, r1, r6 // r6 is congruent mod 2^32 to X*Y70 lsrs r0, r0, #871 lsrs r1, r1, #872 muls r0, r1, r073 lsls r1, r0, #16 // r1 is congruent mod 2^32 to xh*yh74 subs r3, r6, r1 // now r3 is congruent mod 2^32 to75 // (X*Y) - (xh*yh) = xh*yl + xl*yh + xl*yl76 // and hence, since that is at most 0xfeff0001,77 // is _exactly_ equal to that78 adds r0, r0, r5 // r0 is now (xh*yh + (X+Y)<<23) >> 1679 lsrs r1, r3, #16 // r1 is the top 16 bits of r3, i.e.80 // (xh*yl + xl*yh + xl*yl) >> 1681 adds r3, r0, r1 // now r3 equals82 // (xh*yh + xh*yl + xl*yh + xl*yl + (X+Y)<<23) >> 1683 // i.e. (X*Y + (X+Y)<<23) >> 16,84 // i.e. (the right answer) >> 16.85 // Meanwhile, r6 is exactly the bottom 32 bits of the86 // right answer.87 // Renormalise if necessary.88 lsrs r1, r3, #3089 beq LOCAL_LABEL(norenorm)90 // Here we have to do something fiddly. Renormalisation would be a trivial91 // job if we had the leading mantissa bit - just note that it's one bit92 // position above where it should be, and shift right by one. But without93 // that bit, we currently have (2x - 2^30), and we want (x - 2^30); just94 // shifting right would of course give us (x - 2^29), so we must subtract an95 // extra 2^29 to fix this up.96 lsrs r3, r3, #197 movs r1, #198 lsls r1, r1, #2999 subs r3, r3, r1100 adds r2, r2, #1101LOCAL_LABEL(norenorm):102 // Round and shift down to the right bit position.103 lsrs r0, r3, #7 // round bit goes into the carry flag104 bcc LOCAL_LABEL(rounded)105 adds r0, r0, #1106 // In the round-up branch, we must also check if we have to round to even, by107 // testing all the bits below the round bit. We will normally not expect to,108 // so we do RTE by branching out of line and back again to avoid spending a109 // branch in the common case.110 lsls r5, r3, #32-7+1 // check the bits shifted out of r3 above111 bne LOCAL_LABEL(rounded) // if any is nonzero, we're not rounding to even112 lsls r5, r6, #15 // check the bottom 17 bits of the low-order 32113 // (enough to overlap r3 even if we renormalised)114 beq LOCAL_LABEL(rte) // if any is nonzero, fall through, else RTE115LOCAL_LABEL(rounded):116 // Put on the sign and exponent, check for underflow and overflow, and117 // return.118 //119 // Underflow occurs iff r2 (the output exponent) <= 0. Overflow occurs if120 // it's >= 0xFF. (Also if it's 0xFE and we rounded up to overflow, but since121 // this code doesn't report exceptions, we can ignore this case because it'll122 // happen to return the right answer regardless). So we handle most of this123 // via an unsigned comparison against 0xFF, which leaves the one case of a124 // zero exponent that we have to filter separately by testing the Z flag125 // after we shift the exponent back up into place.126 cmp r2, #0xFF // check for most over/underflows127 bhs LOCAL_LABEL(outflow) // ... and branch out of line for them128 lsls r5, r2, #23 // shift the exponent into its output location129 beq LOCAL_LABEL(outflow) // ... and branch again if it was 0130 lsls r4, r4, #31 // shift the output sign into place131 orrs r0, r0, r4 // and OR it in to the output132 adds r0, r0, r5 // OR in the mantissa133 pop {r4,r5,r6,pc} // and return134 135LOCAL_LABEL(rte):136 // Out-of-line handler for the round-to-even case. Clear the low mantissa bit137 // and go back to the post-rounding code.138 movs r5, #1139 bics r0, r0, r5140 b LOCAL_LABEL(rounded)141 142LOCAL_LABEL(outflow):143 cmp r2, #0144 bgt LOCAL_LABEL(overflow)145 // To handle underflow, we construct an intermediate value in the IEEE 754146 // style (using our existing full-length mantissa, and bias the exponent by147 // +0xC0), and indicate whether that intermediate was rounded up, down or not148 // at all. Then call the helper function funder, which will denormalise and149 // re-round correctly.150 lsls r1, r0, #7 // shift up the post-rounding mantissa151 subs r1, r3, r1 // and subtract it from the pre-rounding version152 lsls r6, r6, #15153 cmp r6, #1 // if the rest of the low bits are nonzero154 adcs r1, r1, r1 // then set an extra bit at the bottom155 156 lsls r4, r4, #31157 orrs r0, r0, r4 // put on the sign158 adds r2, r2, #192 // bias the exponent159 lsls r3, r2, #23160 adds r0, r0, r3 // put on the biased exponent161 162 bl SYMBOL_NAME(__compiler_rt_funder)163 pop {r4,r5,r6,pc}164 165LOCAL_LABEL(overflow):166 // Handle overflow by returning an infinity of the correct sign.167 lsls r4, r4, #8 // move the sign up to bit 8168 movs r0, #0xff169 orrs r0, r0, r4 // fill in an exponent just below it170 lsls r0, r0, #23 // and shift those 9 bits up to the top of the word171 pop {r4,r5,r6,pc}172 173 // We come here if there's at least one zero or denormal. On the fast path174 // above, it was convenient to check these before checking NaNs and175 // infinities, but NaNs take precedence, so now we're off the fast path, we176 // must still check for those.177 //178 // At the main entry point 'zerodenorm' we want r2 and r3 to be the two input179 // exponents. So if we branched after shifting-and-checking r2, we come to180 // this earlier entry point 'zerodenorm0' so that we still shift r3.181LOCAL_LABEL(zerodenorm0):182 lsrs r3, r3, #24183LOCAL_LABEL(zerodenorm):184 cmp r2, #255185 beq LOCAL_LABEL(naninf)186 cmp r3, #255187 beq LOCAL_LABEL(naninf)188 // Now we know we have at least one zero or denormal, and no NaN or infinity.189 // Check if either input is actually zero. We've ruled out 0 * infinity by190 // this point, so any zero input means we return zero of the correct sign.191 lsls r6, r0, #1 // is one input zero?192 beq LOCAL_LABEL(zero) // yes, go and return zero193 lsls r6, r1, #1 // is the other one zero?194 bne LOCAL_LABEL(denorm) // if not, one must have been a denormal195LOCAL_LABEL(zero):196 lsls r0, r4, #31 // shift up the output sign to make the return value197 pop {r4,r5,r6,pc}198 199 // Handle denormals via the helper function fnorm2, which will break both200 // inputs up into mantissa and exponent, renormalising and generating a201 // negative exponent if necessary.202LOCAL_LABEL(denorm):203 push {r0,r1,r2,r3}204 mov r0, sp205 bl SYMBOL_NAME(__compiler_rt_fnorm2)206 pop {r0,r1,r2,r3}207 // Convert fnorm2's return values into the right form to rejoin the main208 // code path.209 lsls r0, r0, #1210 lsls r1, r1, #1211 adds r2, r2, r3212 subs r2, r2, #0x7f213 b LOCAL_LABEL(normalised)214 215 // We come here if at least one input is a NaN or infinity. There may still216 // be zeroes (or denormals, though they make no difference at this stage).217LOCAL_LABEL(naninf):218 movs r6, #0xff219 lsls r6, r6, #24220 lsls r5, r0, #1221 cmp r5, r6222 bhi LOCAL_LABEL(nan) // first operand is a NaN223 lsls r5, r1, #1224 cmp r5, r6225 bhi LOCAL_LABEL(nan) // second operand is a NaN226 227 // We know we have at least one infinity, and no NaNs. We might also have a228 // zero, in which case we return the default quiet NaN.229 lsls r6, r0, #1230 beq LOCAL_LABEL(infzero) // if r0 is a zero, r1 must be inf231 lsls r6, r1, #1232 beq LOCAL_LABEL(infzero) // if r1 is a zero, r0 must be inf233 // Otherwise we have infinity * infinity, or infinity * finite. Just return234 // an appropriately signed infinity.235 b LOCAL_LABEL(overflow) // reuse the code there236 237 // We come here if at least one input is a NaN. Hand off to fnan2, which238 // propagates an appropriate NaN to the output, dealing with the special239 // cases of signalling/quiet NaNs.240LOCAL_LABEL(nan):241 bl SYMBOL_NAME(__compiler_rt_fnan2)242 pop {r4,r5,r6,pc}243 244 // Return a quiet NaN as the result of infinity * zero.245LOCAL_LABEL(infzero):246 ldr r0, =0x7fc00000247 pop {r4,r5,r6,pc}248 249END_COMPILERRT_FUNCTION(__mulsf3)250 251NO_EXEC_STACK_DIRECTIVE252