320 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 AArch3211// assembly language suitable to be built as either Arm or Thumb2.12//13//===----------------------------------------------------------------------===//14 15#include "../assembly.h"16 17 18 .syntax unified19 .text20 .p2align 221 22#if __ARM_PCS_VFP23DEFINE_COMPILERRT_FUNCTION(__mulsf3)24 push {r4, lr}25 vmov r0, s026 vmov r1, s127 bl __aeabi_fmul28 vmov s0, r029 pop {r4, pc}30#else31DEFINE_COMPILERRT_FUNCTION_ALIAS(__mulsf3, __aeabi_fmul)32#endif33 34DEFINE_COMPILERRT_FUNCTION(__aeabi_fmul)35 36 // Check if either input exponent is 00 or FF (i.e. not a normalized number),37 // and if so, branch out of line. If we don't branch out of line, then we've38 // also extracted the exponents of the input values r0/r1 into bits 16..23 of39 // r2/r3. But if we do, then that hasn't necessarily been done (because the40 // second AND might have been skipped).41 mov r12, #0xFF000042 ands r2, r12, r0, lsr #7 // sets Z if exponent of x is 043 andsne r3, r12, r1, lsr #7 // otherwise, sets Z if exponent of y is 044 teqne r2, r12 // otherwise, sets Z if exponent of x is FF45 teqne r3, r12 // otherwise, sets Z if exponent of y is FF46 beq LOCAL_LABEL(uncommon) // branch out of line to handle inf/NaN/0/denorm47 48 // Calculate the sign of the result, and put it in an unused bit of r2.49 teq r0, r1 // sets N to the XOR of x and y's sign bits50 orrmi r2, r2, #0x100 // if N set, set bit 8 of r251 52 // Move the input mantissas to the high end of r0/r1, each with its leading53 // bit set explicitly, so that they're in the right form to be multiplied.54 mov r12, #0x8000000055 orr r0, r12, r0, lsl #856 orr r1, r12, r1, lsl #857 58 // Now we're ready to multiply mantissas. This is also the place we'll come59 // back to after decoding denormal inputs. The denormal decoding will also60 // have to set up the same register contents:61 // - decoded fractions at the top of r0 and r162 // - exponents in r2 and r3, starting at bit 1663 // - output sign in r2 bit 864LOCAL_LABEL(mul):65 66 // Here we multiply the mantissas, and compute the output exponent by adding67 // the input exponents and rebiasing. These operations are interleaved to68 // use a delay slot.69 //70 // The exponent is rebiased by subtracting 0x80, rather than the 0x7F you'd71 // expect. That compensates for the leading bit of the mantissa overlapping72 // it, when we recombine the exponent and mantissa by addition.73 add r2, r2, r3 // r2 has sum of exponents, freeing up r374 umull r1, r3, r0, r1 // r3:r1 has the double-width product75 sub r2, r2, #(0x80 << 16) // rebias the summed exponent76 77 // Compress the double-word product into just the high-order word r3, by78 // setting its bit 0 if any bit of the low-order word is nonzero. This79 // changes the represented value, but not by nearly enough to affect80 // rounding, because rounding only depends on the bit below the last output81 // bit, and the general question of whether _any_ nonzero bit exists below82 // that.83 cmp r1, #0 // if low word of full product is nonzero84 orrne r3, r3, #1 // then set LSB of high word85 86 // The two inputs to UMULL had their high bits set, that is, were at least87 // 0x80000000. So the 64-bit product was at least 0x4000000000000000, i.e.88 // the high bit of the product could be at the top of the word or one bit89 // below. Check which, by experimentally shifting left, and then undoing it90 // via RRX if we turned out to have shifted off a 1 bit.91 lsls r3, r3, #1 // shift left, setting C to the bit shifted off92 rrxcs r3, r3 // if that bit was 1, put it back again93 94 // That ensured the leading 1 bit of the product is now the top of r3, but95 // also, set C if the leading 1 was _already_ in the top bit. So now we know96 // whether to increment the exponent. The following instruction does the97 // conditional increment (because it's ADC), but also, copies the exponent98 // field from bit 16 of r2 into bit 0, so as to place it just below the99 // output sign bit.100 //101 // So, if the number hasn't overflowed or underflowed, the low 9 bits of r2102 // are exactly what we need to combine with the rounded mantissa. But the103 // full output exponent (with extra bits) is still available in the high half104 // of r2, so that we can check _whether_ we overflowed or underflowed.105 adc r2, r2, r2, asr #16106 107 // Recombine the exponent and mantissa, doing most of the rounding as a side108 // effect: we shift the mantissa right so as to put the round bit into C, and109 // then we recombine with the exponent using ADC, to increment the mantissa110 // if C was set.111 movs r12, r3, lsr #8112 adc r0, r12, r2, lsl #23113 114 // To complete the rounding, we must check for the round-to-even tiebreaking115 // case, by checking if we're in the exact halfway case, which occurs if and116 // only if we _did_ round up (we can tell this because C is still set from117 // the MOVS), and also, no bit of r3 is set _below_ the round bit.118 //119 // We combine this with an overflow check, so that C ends up set if anything120 // weird happened, and clear if we're completely finished and can return.121 //122 // The best instruction sequence for this part varies between Arm and Thumb.123#if !__thumb__124 // Arm state: if C was set then we check the low bits of r3, so that Z ends125 // up set if we need to round to even.126 //127 // (We rely here on Z reliably being clear to begin with, because shifting128 // down the output mantissa definitely gave a nonzero output. Also, the TST129 // doesn't change C, so if Z does end up set, then C was also set.)130 //131 // Then, if we're not rounding to even, we do a CMP which sets C if there's132 // been an overflow or an underflow. An overflow could occur for an output133 // exponent as low as 0xFC, because we might increment the exponent by 1 when134 // renormalizing, by another when recombining with the mantissa, and by one135 // more if rounding up causes a carry off the top of the mantissa. An136 // underflow occurs only if the output exponent is negative (because it's137 // offset by 1, so an exponent of 0 will be incremented to 1), in which case138 // the top 8 bits of r2 will all be set. Therefore, an unsigned comparison to139 // see if r2 > 0xFC0000 will catch all overflow and underflow cases. It also140 // catches a few very large cases that _don't_ quite overflow (exponents of141 // 0xFC and above that don't get maximally unlucky); those will also be142 // handled by the slow path.143 tstcs r3, #0x7F144 cmpne r2, #0xFC0000145#else146 // In Thumb, switching between different conditions has a higher cost due to147 // the (implicit in this code) IT instructions, so we prefer a strategy that148 // uses CC and CS conditions throughout, at the cost of requiring some extra149 // cleanup instructions on the slow path.150 //151 // If C is set (and hence round-to-even is a possibility), the basic idea is152 // to shift the full result word (r3) left by 25, leaving only its bottom 7153 // bits, which are now the top 7 bits; then we want to set C iff these are 0.154 //155 // The "CMP x,y" instruction sets C if y > x (as unsigned integers). So this156 // could be done in one instruction if only we had a register to use as x,157 // which has 0 in the top 7 bits and at least one nonzero. Then we could158 // compare that against the shifted-up value of r3, setting C precisely if159 // the top 7 bits of y are greater than 0. And happily, we _do_ have such a160 // register! r12 contains the shifted-down mantissa, which is guaranteed to161 // have a 1 in bit 23, and 0 above that.162 //163 // The shift of r3 happens only in the second operand of the compare, so we164 // don't lose the original value of r3 in this process.165 //166 // The check for over/underflow is exactly as in the Arm branch above, except167 // based on a different condition.168 cmpcs r12, r3, lsl #25 // now C is set iff we're rounding to even169 cmpcc r2, #0xFC0000 // and now it's also set if we've over/underflowed170#endif171 172 // That's all the checks for difficult cases done. If C is clear, we can173 // return.174 bxcc lr175 176 // Now the slower path begins. We have to recover enough information to177 // handle all of round-to-even, overflow and underflow.178 //179 // Round to even is the most likely of these, so we detect it first and180 // handle it as fast as possible.181 182#if __thumb__183 // First, Thumb-specific compensation code. The Arm branch of the #if above184 // will have set Z=0 to indicate round to even, but the Thumb branch didn't185 // leave any unambiguous indicator of RTE, so we must retest by checking all186 // the bits shifted off the bottom of the mantissa to see if they're exactly187 // the half-way value.188 lsl r12, r3, #24 // r12 = round bit and everything below189 cmp r12, #0x80000000 // set Z if that is exactly 0x80000000190#endif191 192 // Now Z is clear iff we have already rounded up and now must replace that193 // with rounding to even, which is done by just clearing the low bit of the194 // mantissa.195 biceq r0, r0, #1196 197 // Redo the over/underflow check (the same way as in both branches above),198 // and if it doesn't report a danger, we can return the rounded-to-even199 // answer.200 cmp r2, #0xFC0000 // check for over/underflow201 bxcc lr // and return if none.202 203 // Now we only have overflow and underflow left to handle. First, find out204 // which we're looking at. This is easy by testing the top bit of r2, but205 // even easier by using the fact that the possible positive and negative206 // values of r2 are widely enough separated that the 0xFC0000 subtracted by207 // the CMP above won't have made any difference. So the N flag output from208 // that comparison _already_ tells us which condition we have: if N is set we209 // have underflow, and if N is clear, overflow.210 bpl LOCAL_LABEL(overflow)211 212 // Here we're handling underflow.213 214 // Add the IEEE 754:1985 exponent bias which funder will expect. This also215 // brings the exponent back into a range where it can't possibly have carried216 // into the sign bit, so the output sign will now be right.217 add r0, r0, #(0xC0 << 23)218 219 // Determine whether we rounded up, down or not at all.220 lsls r2, r3, #1 // input mantissa, without its leading 1221 subs r1, r2, r0, lsl #9 // subtract the output mantissa (likewise)222 223 // And let funder handle the rest.224 b SYMBOL_NAME(__compiler_rt_funder)225 226LOCAL_LABEL(overflow):227 // We come here to handle overflow, but it's not guaranteed that an overflow228 // has actually happened: our check on the fast path erred on the side of229 // caution, by catching any output exponent that _could_ cause an overflow.230 // So first check whether this really is an overflow, by extracting the231 // output exponent. Exponent 0xFF, or anything that wrapped round to having232 // the high bit clear, are overflows; 0xFE down to 0xFC are not overflows.233 //234 // The value in r0 is correct to return, if there's no overflow.235 add r12, r0, #(1 << 23) // add 1 to the exponent so 0xFF wraps to 0236 movs r12, r12, lsl #1 // test the top bit of the modified value237 bxmi lr // if top bit is still 1, not an overflow238 239 // This is an overflow, so we need to replace it with an appropriately signed240 // infinity. First we correct the sign by applying a downward bias to the241 // exponent (the one suggested in IEEE 754:1985, which was chosen to bring242 // all possible overflowed results back into range).243 subs r0, r0, #(0xC0 << 23)244 245 // Now the sign bit of r0 is correct. Replace everything else with the246 // encoding of an infinity.247 mov r1, #0xFF248 and r0, r0, #0x80000000249 orr r0, r0, r1, lsl #23250 bx lr251 252LOCAL_LABEL(uncommon):253 // Handle zeros, denorms, infinities and NaNs. We arrive here knowing that254 // we've at least done the first _two_ instructions from the entry point,255 // even if all the rest were skipped. So r2 contains the sign and exponent of256 // x in bits 16..23, and r12 = 0xFF << 16.257 //258 // So, first repeat some instructions from the prologue, which were either259 // conditionally skipped in the sequence leading to the branch, or skipped260 // because they happened after the branch.261 and r3, r12, r1, lsr #7 // get exponent of y in r3 bits 16..23262 teq r0, r1 // calculate the sign of the result263 orrmi r2, r2, #0x100 // and put it in bit 8 of r2 as before264 265 // Check for infinities and NaNs, by testing each of r2,r3 to see if it's at266 // least 0xFF0000 (hence the exponent field is equal to 0xFF).267 cmp r2, r12268 cmplo r3, r12269 bhs LOCAL_LABEL(inf_NaN)270 271 // If we didn't take that branch, then we have only finite numbers, but at272 // least one is denormal or zero. A zero makes the result easy (and also is a273 // more likely input than a denormal), so check those first, as fast as274 // possible.275 movs r12, r0, lsl #1 // Z set if x == 0276 movsne r12, r1, lsl #1 // now Z set if either input is 0277 moveq r0, r2, lsl #23 // in either case, make 0 of the output sign278 bxeq lr // and return it279 280 // Now we know we only have denormals to deal with. Call fnorm2 to sort281 // them out, and rejoin the main code path above.282 and r12, r2, #0x100 // save the result sign from r2283 lsr r2, #16 // shift extracted exponents down to bit 0284 lsr r3, #16 // where fnorm2 will expect them285 push {r0, r1, r2, r3, r12, lr}286 mov r0, sp // tell fnorm2 where to find its data287 bl SYMBOL_NAME(__compiler_rt_fnorm2)288 pop {r0, r1, r2, r3, r12, lr}289 lsl r3, #16 // shift exponents back up to bit 16290 orr r2, r12, r2, lsl #16 // and put the result sign back in r2291 b LOCAL_LABEL(mul)292 293LOCAL_LABEL(inf_NaN):294 // We come here if at least one input is a NaN or infinity. If either or both295 // inputs are NaN then we hand off to fnan2 which will propagate a NaN from296 // the input; otherwise any multiplication involving infinity returns297 // infinity, unless it's infinity * 0 which is an invalid operation and298 // returns NaN again.299 mov r12, #0xFF000000300 cmp r12, r0, lsl #1 // if (r0 << 1) > 0xFF000000, r0 is a NaN301 blo SYMBOL_NAME(__compiler_rt_fnan2)302 cmp r12, r1, lsl #1303 blo SYMBOL_NAME(__compiler_rt_fnan2)304 305 // NaNs are dealt with, so now we have at least one infinity. Check if the306 // other operand is 0. This is conveniently done by XORing the two: because307 // we know that the low 31 bits of one operand are exactly 0x7F800000, we can308 // test if the low 31 bits of the other one are all 0 by checking whether the309 // low 31 bits of (x XOR y) equal 0x7F800000.310 eor r3, r0, r1311 cmp r12, r3, lsl #1 // if inf * 0, this sets Z312 lsr r0, r12, #1 // set up return value of +infinity313 orrne r0, r0, r2, lsl #23 // if not inf * 0, put on the output sign314 orreq r0, r0, #0x400000 // otherwise, set the 'quiet NaN' bit315 bx lr // and return316 317END_COMPILERRT_FUNCTION(__aeabi_fmul)318 319NO_EXEC_STACK_DIRECTIVE320