619 lines · plain
1//===-- divsf3.S - single-precision floating point division ---------------===//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 division with the IEEE-75410// default rounding (to nearest, ties to even), in optimized AArch32 assembly11// 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(__divsf3)24 push {r4, lr}25 vmov r0, s026 vmov r1, s127 bl __aeabi_fdiv28 vmov s0, r029 pop {r4, pc}30#else31DEFINE_COMPILERRT_FUNCTION_ALIAS(__divsf3, __aeabi_fdiv)32#endif33 34DEFINE_COMPILERRT_FUNCTION(__aeabi_fdiv)35 // Extract the exponents of the inputs into r2 and r3, occupying bits 16-2336 // of each register so that there will be space lower down to store extra37 // data without exponent arithmetic carrying into it. In the process, check38 // both exponents for 00 or FF and branch out of line to handle all the39 // uncommon types of value (infinity, NaN, zero, denormals).40 //41 // Chaining conditional instructions like this means that the second42 // instruction (setting up r3) might not be executed at all, so fdiv_uncommon43 // will have to redo it just in case. That saves an instruction here,44 // executed for _all_ inputs, and moves it to the uncommon path run for only45 // some inputs.46 mov r12, #0xFF000047 ands r2, r12, r0, lsr #7 // r2 has exponent of numerator. (Is it 0?)48 andsne r3, r12, r1, lsr #7 // r3 has exponent of denominator. (Is it 0?)49 teqne r2, r12 // if neither was 0, is one FF?50 teqne r3, r12 // or the other?51 beq LOCAL_LABEL(uncommon) // branch out of line if any answer was yes52 53 // Calculate the output sign, which is always just the XOR of the input54 // signs. Store it in bit 8 of r2, below the numerator exponent.55 teq r0, r1 // is the output sign bit 1?56 orrmi r2, r2, #0x100 // if so, set bit 8 of r257 58 // Isolate the mantissas of both values, by setting bit 23 of each one and59 // clearing the 8 bits above that.60 //61 // In the process, swap the register allocations (which doesn't cost extra62 // instructions if we do it as part of this manipulation). We want the63 // numerator not to be in r0, because r0 is where we'll build up the quotient64 // while subtracting things from the numerator.65 orr r12, r0, #1 << 2366 orr r0, r1, #1 << 2367 bic r1, r12, #0xFF00000068 bic r0, r0, #0xFF00000069 70LOCAL_LABEL(div):71 // Start of the main division. We get here knowing that:72 //73 // r0 = mantissa of denominator, with the leading 1 at bit 2374 // r1 = mantissa of numerator, similarly75 // r2 = (exponent of numerator << 16) + (result sign << 8)76 // r3 = (exponent of denominator << 16)77 78 push {r14} // we'll need an extra register79 80 // Calculate the initial result exponent by just subtracting the two input81 // exponents. This doesn't affect the sign bit lower down in r2.82 sub r2, r2, r383 84 // That initial exponent might need to be adjusted by 1, depending on whether85 // dividing the mantissas gives a value >=1 or <1. We don't need to wait86 // until the division is finished to work that out: we can tell immediately87 // by just comparing the mantissas.88 //89 // The basic idea is to do the comparison in a way that sets the C flag if90 // numerator >= denominator. Then we recombine the sign and exponent by doing91 // "ADC r2, r2, r2, asr #16": the exponent in the top half of r2 is shifted92 // down to the low 8 bits, just below the sign bit, and using ADC rather than93 // ADD folds in the conditional increment from the mantissa comparison.94 //95 // If we're not incrementing the output exponent, we instead shift the96 // numerator mantissa left by 1, so that it _is_ greater than the denominator97 // mantissa. Otherwise we'd generate only a 22-bit quotient, instead of 23.98 //99 // The exponent also needs to be rebiased, so that dividing two numbers the100 // same gives an output exponent of 0x7F. If the two inputs have the same101 // exponent then we'll have computed an exponent of 0 via the SUB instruction102 // above; if the mantissas are the same as well then the ADC will increment103 // it; also, the leading bit of the quotient will increment the exponent104 // again when we recombine it with the output mantissa later. So we need to105 // add (0x7F - 2) to the mantissa now, to make an exponent of 0 from the SUB106 // come to 0x7F after both of those increments.107 //108 // Putting all of that together, what we _want_ to do is this:109 //110 // [#1] CMP r1, r0 // set C if num >= den111 // [#2] MOVLO r1, r1, lsl #1 // if num < den, shift num left112 // [#3] ADD r2, r2, #0x7D0000 // rebias exponent113 // [#4] ADC r2, r2, r2, asr #16 // combine sign + exp + adjustment114 //115 // However, we only do the first of those four instructions right here. The116 // other three are distributed through the code below, after unrelated load117 // or multiply instructions which will have a result delay slot on simple118 // CPUs. Each is labelled "exponent setup [#n]" in a comment.119 //120 // (Since instruction #4 depends on the flags set up by #2, we must avoid121 // clobbering the flags in _any_ of the instructions interleaved with this!)122 cmp r1, r0 // exponent setup [#1]123 124 // Start the mantissa division by making an approximation to the reciprocal125 // of the denominator. We first obtain an 8-bit approximation using a table126 // lookup indexed by the top 7 denominator bits (counting the leading 1, so127 // really there are only 6 bits in the table index).128 //129 // (r0 >> 17) is the table index, and its top bit is always set, so it ranges130 // from 64 to 127 inclusive. So we point the base register 64 bytes before131 // the actual table.132 adr r12, LOCAL_LABEL(tab) - 64133#if __thumb__134 // Thumb can't do this particular shift+add+load in one instruction - it only135 // supports left shifts of 0 to 3 bits, not right shifts of 17. So we must136 // calculate the load offset separately.137 add r14, r12, r0, lsr #17138 ldrb r14, [r14]139#else140 ldrb r14, [r12, r0, lsr #17]141#endif142 143 // Now do an iteration of Newton-Raphson to improve that 8-bit approximation144 // to have 15-16 accurate bits.145 //146 // Basics of Newton-Raphson for finding a reciprocal: if you want to find 1/d147 // and you have some approximation x, your next approximation is X = x(2-dx).148 // Looked at one way, this is the result of applying the N-R formula149 // X=x-f(x)/f'(x) to the function f(x) = 1/x - d. Another way to look at it150 // is to suppose that dx = 1 - e, for some e which is small (because dx is151 // already reasonably close to 1). Then you want to double the number of152 // correct bits in the next approximation, i.e. square the error. So you want153 // dX = 1-e^2 = (1-e)(1+e) = dx(2-dx). Cancelling d gives X = x(2-dx) again.154 //155 // In this situation, we're working in fixed-point integers rather than real156 // numbers, and all the scales are different:157 // * our input denominator d is in the range [2^23,2^24)158 // * our input approximation x is in the range [2^7,2^8)159 // * we want the output approximation to be in the range [2^15,2^16)160 // Those factors combine to mean that we want161 // x(2^32-dx) / 2^23162 // = (2^9 x) - (dx^2 / 2^23)163 //164 // But we also want to compute this using ordinary MUL, not a long multiply165 // instruction (those are slower). So we need to worry about the product166 // overflowing. dx fits in 32 bits, because it's the product of something167 // <2^24 with something <2^8; but we must shift it right before multiplying168 // by x again.169 170 mul r12, r0, r14 // r12 = dx171 movlo r1, r1, lsl #1 // exponent setup [#2] in the MUL delay slot172 mvn r12, r12, lsr #8 // r12 ~= -dx/2^8173 mul r3, r12, r14 // r3 ~= -dx^2/2^8174 mov r14, r14, lsl #9 // r14 = 2^9 x175 add r14, r14, r3, asr #15 // r14 ~= 2^9 x - dx^2 / 2^23176 177 // Now r14 is a 16-bit approximation to the reciprocal of the input mantissa,178 // scaled by 2^39 (so that the min mantissa 2^23 would have reciprocal 2^16179 // in principle, and the max mantissa 2^24-1 would have reciprocal just over180 // 2^15). The error is always negative (r14 is an underestimate of the true181 // value), and the maximum error is 6 and a bit ULP (that is, the true182 // reciprocal is strictly less than (r14+7)). Also, r14 is always strictly183 // less than 0x10000 (even in the case of the min mantissa, where the true184 // value would be _exactly_ 0x10000), which eliminates a case of integer185 // overflow.186 //187 // All of these properties of the reciprocal approximation are checked by188 // exhaustively iterating over all 2^23 possible input mantissas. (The nice189 // thing about doing this in single rather than double precision!)190 //191 // Now we extract most of the quotient by two steps of long division, using192 // the reciprocal estimate to identify a multiple of the denominator to193 // subtract from the numerator. To avoid integer overflow, the numerator194 // mantissa is shifted down 8 bits so that it's less than 0x10000. After we195 // calculate an approximate quotient, we shift the numerator left and196 // subtract that multiple of the denominator, moving the next portion of the197 // numerator into range for the next iteration.198 199 // First iteration of long division. We shift the numerator left 11 bits, and200 // since the quotient approximation is scaled by 2^31, we must shift that201 // right by 20 to make the right product to subtract from the numerator.202 mov r12, r1, lsr #8 // shift the numerator down203 mul r12, r14, r12 // make the quotient approximation204 mov r1, r1, lsl #11 // shift numerator left, ready for subtraction205 mov r3, r12, lsr #20 // make first 12-bit block of quotient bits206 mls r1, r0, r3, r1 // subtract that multiple of den from num207 208 add r2, r2, #0x7D0000 // exponent setup [#3] in the MLS delay slot209 210 // Second iteration of long division. Differences from the first step: this211 // time we shift the numerator 12 bits instead of 11, so that the total of212 // both steps is 23 bits, i.e. we've shifted up by exactly the full width of213 // the output mantissa. Also, the block of output quotient bits is left in a214 // different register: it was in r3 the first time, and this time it's in215 // r12, so that we still have both available at the end of the process.216 mov r12, r1, lsr #8 // shift the numerator down217 mul r12, r14, r12 // make the quotient approximation218 mov r1, r1, lsl #12 // shift numerator left, ready for subtraction219 mov r12, r12, lsr #19 // make second 11-bit block of quotient220 mls r1, r0, r12, r1 // subtract that multiple of den from num221 222 adc r2, r2, r2, asr #16 // exponent setup [#4] in the MLS delay slot223 224 // Now r1 contains the original numerator, shifted left 23, minus _some_225 // multiple of the original denominator (which is still in r0). The bounds on226 // the error in the above steps should make the error at most 1: that is, we227 // may have to subtract the denominator one more time to make r1 < r0, and228 // increment the quotient by one more.229 //230 // Our quotient is still in two pieces, computed separately in the above long231 // division steps. We fold the final increment into the same instruction that232 // recombines them, by doing the comparison in such a way that it sets the233 // carry flag if the increment is needed.234 235 cmp r1, r0 // Set carry flag if num >= den236 subhs r1, r1, r0 // If so, subtract den from num237 adc r3, r12, r3, lsl #12 // Recombine quotient halves, plus optional +1238 239 // We've finished with r14 as a temporary register, so we can unstack it now.240 pop {r14}241 242 // Now r3 contains the _rounded-down_ output quotient, and r1 contains the243 // remainder. That is, (denominator * r3 + r1) = (numerator << 23), and244 // 0 <= r1 < denominator.245 //246 // Next we must round to nearest, by checking if r1 is greater than half the247 // denominator. In division, it's not possible to hit an exact round-to-even248 // halfway case, so we don't need to spend any time checking for it.249 //250 // Proof of no round-to-even: define the 'width' of a dyadic rational to be251 // the distance between the lowest and highest 1 bits in its binary252 // representation, or equivalently, the index of its high bit if you scale it253 // by a power of 2 to make it an odd integer. E.g. any actual power of 2 has254 // width 0, and all of 0b11110, 0b1111, 0b11.11 and 0b0.01111 have width 3.255 // Then for any dyadic rationals a,b, width(ab) >= width(a)+width(b). Let w256 // be the maximum width that the input precision supports (so that for single257 // precision, w=23). Then if some division n/d were a round-to-even case, the258 // true quotient q=n/d would have width exactly w+1. But we have qd=n, so259 // width(n) >= width(q)+width(d) > w, which can't happen, because n is in the260 // input precision, hence had width <= w.)261 //262 // So we don't need to check for an exact _halfway_ case and clear the low263 // bit of the quotient after rounding up, as addition and multiplication both264 // need to do. But we do need to remember if the quotient itself was exact,265 // that is, if there was no remainder at all. That's needed in underflow266 // handling.267 268 // The rounding check wants to compare remainder with denominator/2. But of269 // course in integers it's easier to compare 2*remainder with denominator. So270 // we start by shifting the remainder left by 1, and in the process, set Z if271 // it's exactly 0 (i.e. the result needs no rounding at all).272 lsls r1, r1, #1273 // Now trial-subtract the denominator. We don't do this at all if the result274 // was exact. If we do do it, r1 goes negative precisely if we need to round275 // up, which sets the C flag. (The previous instruction will have left C276 // clear, since r1 had its top 8 bits all clear. So now C is set _only_ if277 // we're rounding up.)278 subsne r1, r1, r0279 // Recombine the quotient with the sign + exponent, and use the C flag from280 // the previous instruction to increment the quotient if we're rounding up.281 adc r0, r3, r2, lsl #23282 283 // If we haven't either overflowed or underflowed, we're done. We can284 // identify most of the safe cases by doing an unsigned comparison of the285 // initial output exponent (in the top half of r2) with 0xFC: if 0 <= r2 <286 // 0xFC0000 then we have neither underflow nor overflow.287 //288 // Rationale: the value in the top half of r2 had three chances to be289 // incremented before becoming the exponent field of the actual output float.290 // It was incremented if we found the numerator mantissa was >= the291 // denominator (producing the value in the _bottom_ half of r2, which we just292 // ADCed into the output). Then it gets unconditionally incremented again293 // when the ADC combines it with the leading mantissa bit. And finally,294 // round-up might increment it a third time. So 0xFC is the smallest value295 // that can possibly turn into the overflowed value 0xFF after all those296 // increments.297 //298 // On the underflow side, (top half of r2) = 0 corresponds to a value of 1 in299 // the final result's exponent field (and then rounding might increase it300 // further); if the exponent was less than that then r2 wraps round and looks301 // like a very large positive integer from the point of view of this unsigned302 // comparison.303 cmp r2, #0xFC0000304 bxlo lr305 306 // The same comparison will have set the N and V flags to reflect the result307 // of comparing r2 with 0xFC0000 as a _signed_ integer. That reliably308 // distinguishes potential underflow (r2 is negative) from potential overflow309 // (r2 is positive and at least 0xFC0000)310 bge LOCAL_LABEL(overflow)311 312 // Here we might or might not have underflow (but we know we don't have313 // overflow). To check more carefully, we look at the _bottom_ half of r2,314 // which contains the exponent after the first adjustment (for num >= denom),315 // That is, it's still off by 1 (compensating for the leading quotient bit),316 // and is also before rounding.317 //318 // We neglect the effect of rounding: division results that are tiny (less319 // than the smallest normalised number) before rounding, but then round up to320 // the smallest normal number, are an acceptable edge case to handle slowly.321 // We pass those to funder without worrying about them.322 //323 // So we want to check whether the bottom half of r2 was negative. It would324 // be nice to check bits 8-15 of it, but unfortunately, it's already been325 // combined with the sign (at bit 8), so those bits don't tell us anything326 // useful. Instead we look at the top 4 bits of the exponent field, i.e. the327 // 0xF0 bits. The largest _non_-overflowing exponent that might reach here is328 // less than 3, so it doesn't reach those bits; the smallest possible329 // underflow, obtained by dividing the smallest denormal by the largest330 // finite number, is -151 (before the leading bit increments it), which will331 // set the low 8 bits of r2 to 0x69. That is, the 0xF0 nibble of r2 will be332 // 0x60 or greater for a (pre-rounding) underflow, and zero for a333 // non-underflow.334 335 tst r2, #0xF0336 bxeq lr // no underflow after all; return337 338 // Rebias the exponent for funder, which also corrects the sign bit.339 add r0, r0, #192 << 23340 // Tell funder whether the true value is greater or less than the number in341 // r0. This is obtained from the sign of the remainder (still in r1), with342 // the only problem being that it's currently reversed. So negate r1 (leaving343 // 0 at 0 to indicate exactness).344 rsbs r1, r1, #0345 b SYMBOL_NAME(__compiler_rt_funder)346 347LOCAL_LABEL(overflow):348 // Here we might or might not have overflow (but we know we don't have349 // underflow). We must check whether we really have overflowed.350 //351 // For this it's easiest to check the exponent field in the actual output352 // value in r0, after _all_ the adjustments have been completed. The largest353 // overflowed exponent is 0x193, and the smallest exponent that can reach354 // this is 0xFD (we checked against 0xFC above, but then the leading quotient355 // bit incremented it). So it's enough to shift the output left by one356 // (moving the exponent field to the top), increment it once more (so that357 // the smallest overflowed exponent 0xFF wraps round to 0), and then compare358 // against 0xFE000000 as an unsigned integer.359 mov r12, r0, lsl #1360 add r12, r12, #1 << 24361 cmp r12, #0xFE << 24 // Check for exp = 253 or 254362 bxhs lr363 // We have actual overflow. Rebias r0 to bring the exponent back into range,364 // which ensures its sign is correct. Then make an infinity of that sign to365 // return.366 subs r0, r0, #0xC0 << 23367 movs r12, #0xFF // exponent of infinity368 orrs r12, r12, r0, lsr #23 // exponent and sign at bottom of r12369 movs r0, r12, lsl #23 // shift it up to the top of r0 to return370 bx lr371 372LOCAL_LABEL(uncommon):373 // We come here from the start of the function if either input is an uncommon374 // value: zero, denormal, infinity or NaN.375 //376 // We arrive here with r12 = 0xFF000000, and r2 containing the exponent of x377 // in bits 16..23. But r3 doesn't necessarily contain the exponent of y,378 // because the instruction that set it up was conditional. So first we379 // unconditionally repeat it.380 and r3, r12, r1, lsr #7381 382 // In all cases not involving a NaN as output, the sign of the output is made383 // in the same way as for finite numbers, as the XOR of the input signs. So384 // repeat the sign setup from the main branch.385 teq r0, r1 // is the output sign bit 1?386 orrmi r2, r2, #0x100 // if so, set bit 8 of r2387 388 // Detect infinities and NaNs, by checking if either of r2 or r3 is at least389 // 0xFF0000.390 cmp r2, #0xFF0000391 cmplo r3, #0xFF0000392 bhs LOCAL_LABEL(inf_NaN)393 394 // Now we know there are no infinities or NaNs, but there's at least one zero395 // or denormal.396 movs r12, r1, lsl #1 // is y zero?397 beq LOCAL_LABEL(divbyzero) // if so, go and handle division by zero398 movs r12, r0, lsl #1 // is x zero? (now we know that y is not)399 moveq r0, r2, lsl #23 // if so, 0/nonzero is just 0 (of right sign)400 bxeq lr401 402 // Now we've eliminated zeroes as well, leaving only denormals: either x or403 // y, or both, is a denormal. Call fnorm2 to convert both into a normalised404 // mantissa and a (potentially small) exponent.405 and r12, r2, #0x100 // save the result sign from r2406 lsr r2, #16 // shift extracted exponents down to bit 0407 lsr r3, #16 // where fnorm2 will expect them408 push {r0, r1, r2, r3, r12, lr}409 mov r0, sp // tell fnorm2 where to find its data410 bl SYMBOL_NAME(__compiler_rt_fnorm2)411 pop {r0, r1, r2, r3, r12, lr}412 lsl r3, #16 // shift exponents back up to bit 16413 orr r2, r12, r2, lsl #16 // and put the result sign back in r2414 415 // Now rejoin the main code path, having finished the setup it will expect:416 // swap x and y, and shift the fractions back down to the low 24 bits.417 mov r12, r0, lsr #8418 mov r0, r1, lsr #8419 mov r1, r12420 b LOCAL_LABEL(div)421 422LOCAL_LABEL(inf_NaN):423 // We come here if at least one input is a NaN or infinity. If either or both424 // inputs are NaN then we hand off to fnan2 to propagate a NaN from the425 // input.426 mov r12, #0xFF000000427 cmp r12, r0, lsl #1 // if (r0 << 1) > 0xFF000000, r0 is a NaN428 blo SYMBOL_NAME(__compiler_rt_fnan2)429 cmp r12, r1, lsl #1430 blo SYMBOL_NAME(__compiler_rt_fnan2)431 432 // No NaNs, so we have three options: inf/inf = NaN, inf/finite = inf, and433 // finite/inf = 0.434 435 // If both operands are infinity, we return a NaN. Since we know at436 // least _one_ is infinity, we can test this by checking if they're437 // equal apart from the sign bits.438 eor r3, r0, r1439 lsls r3, #1 // were all bits of XOR zero other than top?440 beq LOCAL_LABEL(invalid) // if so, both operands are infinity441 442 // See if x is infinite443 cmp r12, r0, lsl #1 // (r0 << 1) == 0xFF000000?444 beq LOCAL_LABEL(infret) // if so, infinity/finite = infinity445 446 // y is infinite and x is not, so we return a zero of the447 // combined sign.448 eor r0, r0, r1 // calculate the right sign449 and r0, r0, #0x80000000 // throw away everything else450 bx lr451 452LOCAL_LABEL(divbyzero):453 // Here, we know y is zero. But we don't know if x is zero or nonzero. So we454 // might be calculating 0/0 (invalid operation, generating a NaN), or455 // nonzero/0 (the IEEE "division by zero" exception, generating infinity).456 movs r12, r0, lsl #1 // is x zero too?457 beq LOCAL_LABEL(invalid) // if so, go and return a NaN458 459LOCAL_LABEL(infret):460 // Here, we're either dividing infinity by a finite number, or dividing a461 // nonzero number by 0. (Or both, if we're dividing infinity by 0.) In all462 // these cases we return infinity with the sign from r2.463 //464 // If we were implementing IEEE exceptions, we'd have to separate these465 // cases: infinity / finite is not an _exception_, it just returns infinity,466 // whereas (finite and nonzero) / 0 is a division-by-zero exception. But here467 // we're not implementing exceptions, so we can treat all three cases the468 // same.469 //470 // r2 contains the output sign in bit 8, which is a convenient place to find471 // it when making an infinity, because we can fill in the 8 exponent bits472 // below that and then shift it left.473 orr r2, r2, #0xff // sign + maximum exponent474 lsl r0, r2, #23 // shift up to the top475 bx lr476 477LOCAL_LABEL(invalid):478 // Return the default NaN, from an invalid operation (either dividing479 // infinity by infinity, or 0 by 0).480 ldr r0, =0x7FC00000481 bx lr482 483// Finally, the lookup table for the initial reciprocal approximation.484//485// The table index is made from the top 7 bits of the denominator mantissa. But486// the topmost bit is always 1, so only the other 6 bits vary. So it only has487// 64 entries, not 128.488//489// Each table entry is a single byte, with its top bit set. So the table490// entries correspond to the reciprocal of a 7-bit mantissa prefix scaled up by491// 2^14, or the reciprocal of a whole 24-bit mantissa scaled up by 2^31.492//493// Each of these 64 entries corresponds to a large interval of possible494// mantissas. For example, if the top 7 bits are 1000001 then the overall495// mantissa could be anything from 0x820000 to 0x83FFFF. And because the output496// of this table provides more bits than the input, there are several choices497// of 8-bit reciprocal approximation for a number in that interval. The498// reciprocal of 0x820000 starts with 0xFC plus a fraction, and the reciprocal499// of 0x83FFFF starts with 0xF9 minus a fraction, so there are four reasonable500// choices for that table entry: F9, FA, FB or FC. Which do we pick?501//502// The table below is generated by choosing whichever value minimises the503// maximum possible error _after_ the approximation is improved by the504// Newton-Raphson step. In the example above, we end up with FA.505//506// The Python code below will regenerate the table, complete with the per-entry507// comments.508 509/*510 511for prefix in range(64, 128):512 best = None513 514 # Max and min 23-bit mantissas with this 7-bit prefix515 mmin, mmax = prefix * 2**17, (prefix + 1) * 2**17 - 1516 517 # Max and min table entry corresponding to the reciprocal of something in518 # that range of mantissas: round up the reciprocal of mmax, and round down519 # the reciprocal of mmin. Also clamp to the range [0x80,0xff], because520 # 0x100 can't be used as a table entry due to not fitting in a byte, even521 # though it's the exact reciprocal of the overall-smallest mantissa522 # 0x800000.523 gmin = max(128, (2**31 + mmin - 1) // mmax)524 gmax = min(255, 2**31 // mmin)525 526 # For each of those table entries, compute the result of starting from that527 # value and doing a Newton-Raphson iteration, with the mantissa at each end528 # of the mantissa interval. One of these will be the worst possible error.529 # Choose the table entry whose worst error is as small as possible.530 #531 # (To find the extreme values of a more general function on an interval,532 # you must consider its values not only at the interval endpoints but also533 # any turning points within the interval. Here, the function has only one534 # turning point, and by construction it takes value 0 there, so we needn't535 # worry.)536 g = max(537 range(gmin, gmax + 1),538 key=lambda g: min(539 (g * (2**32 - d * g) / 2**23 - 2**39 / d) for d in [mmin, mmax]540 ),541 )542 543 print(f" .byte 0x{g:02x} // input [0x{mmin:06x},0x{mmax:06x}]"544 f", candidate outputs [0x{gmin:02x},0x{gmax:02x}]"545 )546 547*/548 549 .p2align 2 // make sure we start on a 4-byte boundary, even in Thumb550LOCAL_LABEL(tab):551 .byte 0xfe // input [0x800000,0x81ffff], candidate outputs [0xfd,0xff]552 .byte 0xfa // input [0x820000,0x83ffff], candidate outputs [0xf9,0xfc]553 .byte 0xf6 // input [0x840000,0x85ffff], candidate outputs [0xf5,0xf8]554 .byte 0xf3 // input [0x860000,0x87ffff], candidate outputs [0xf1,0xf4]555 .byte 0xef // input [0x880000,0x89ffff], candidate outputs [0xee,0xf0]556 .byte 0xec // input [0x8a0000,0x8bffff], candidate outputs [0xeb,0xed]557 .byte 0xe8 // input [0x8c0000,0x8dffff], candidate outputs [0xe7,0xea]558 .byte 0xe5 // input [0x8e0000,0x8fffff], candidate outputs [0xe4,0xe6]559 .byte 0xe2 // input [0x900000,0x91ffff], candidate outputs [0xe1,0xe3]560 .byte 0xdf // input [0x920000,0x93ffff], candidate outputs [0xde,0xe0]561 .byte 0xdc // input [0x940000,0x95ffff], candidate outputs [0xdb,0xdd]562 .byte 0xd9 // input [0x960000,0x97ffff], candidate outputs [0xd8,0xda]563 .byte 0xd6 // input [0x980000,0x99ffff], candidate outputs [0xd5,0xd7]564 .byte 0xd3 // input [0x9a0000,0x9bffff], candidate outputs [0xd3,0xd4]565 .byte 0xd1 // input [0x9c0000,0x9dffff], candidate outputs [0xd0,0xd2]566 .byte 0xce // input [0x9e0000,0x9fffff], candidate outputs [0xcd,0xcf]567 .byte 0xcc // input [0xa00000,0xa1ffff], candidate outputs [0xcb,0xcc]568 .byte 0xc9 // input [0xa20000,0xa3ffff], candidate outputs [0xc8,0xca]569 .byte 0xc7 // input [0xa40000,0xa5ffff], candidate outputs [0xc6,0xc7]570 .byte 0xc4 // input [0xa60000,0xa7ffff], candidate outputs [0xc4,0xc5]571 .byte 0xc2 // input [0xa80000,0xa9ffff], candidate outputs [0xc1,0xc3]572 .byte 0xc0 // input [0xaa0000,0xabffff], candidate outputs [0xbf,0xc0]573 .byte 0xbd // input [0xac0000,0xadffff], candidate outputs [0xbd,0xbe]574 .byte 0xbb // input [0xae0000,0xafffff], candidate outputs [0xbb,0xbc]575 .byte 0xb9 // input [0xb00000,0xb1ffff], candidate outputs [0xb9,0xba]576 .byte 0xb7 // input [0xb20000,0xb3ffff], candidate outputs [0xb7,0xb8]577 .byte 0xb5 // input [0xb40000,0xb5ffff], candidate outputs [0xb5,0xb6]578 .byte 0xb3 // input [0xb60000,0xb7ffff], candidate outputs [0xb3,0xb4]579 .byte 0xb1 // input [0xb80000,0xb9ffff], candidate outputs [0xb1,0xb2]580 .byte 0xaf // input [0xba0000,0xbbffff], candidate outputs [0xaf,0xb0]581 .byte 0xad // input [0xbc0000,0xbdffff], candidate outputs [0xad,0xae]582 .byte 0xac // input [0xbe0000,0xbfffff], candidate outputs [0xab,0xac]583 .byte 0xaa // input [0xc00000,0xc1ffff], candidate outputs [0xa9,0xaa]584 .byte 0xa8 // input [0xc20000,0xc3ffff], candidate outputs [0xa8,0xa8]585 .byte 0xa6 // input [0xc40000,0xc5ffff], candidate outputs [0xa6,0xa7]586 .byte 0xa5 // input [0xc60000,0xc7ffff], candidate outputs [0xa4,0xa5]587 .byte 0xa3 // input [0xc80000,0xc9ffff], candidate outputs [0xa3,0xa3]588 .byte 0xa1 // input [0xca0000,0xcbffff], candidate outputs [0xa1,0xa2]589 .byte 0xa0 // input [0xcc0000,0xcdffff], candidate outputs [0xa0,0xa0]590 .byte 0x9e // input [0xce0000,0xcfffff], candidate outputs [0x9e,0x9f]591 .byte 0x9d // input [0xd00000,0xd1ffff], candidate outputs [0x9d,0x9d]592 .byte 0x9b // input [0xd20000,0xd3ffff], candidate outputs [0x9b,0x9c]593 .byte 0x9a // input [0xd40000,0xd5ffff], candidate outputs [0x9a,0x9a]594 .byte 0x98 // input [0xd60000,0xd7ffff], candidate outputs [0x98,0x99]595 .byte 0x97 // input [0xd80000,0xd9ffff], candidate outputs [0x97,0x97]596 .byte 0x96 // input [0xda0000,0xdbffff], candidate outputs [0x95,0x96]597 .byte 0x94 // input [0xdc0000,0xddffff], candidate outputs [0x94,0x94]598 .byte 0x93 // input [0xde0000,0xdfffff], candidate outputs [0x93,0x93]599 .byte 0x92 // input [0xe00000,0xe1ffff], candidate outputs [0x91,0x92]600 .byte 0x90 // input [0xe20000,0xe3ffff], candidate outputs [0x90,0x90]601 .byte 0x8f // input [0xe40000,0xe5ffff], candidate outputs [0x8f,0x8f]602 .byte 0x8e // input [0xe60000,0xe7ffff], candidate outputs [0x8e,0x8e]603 .byte 0x8d // input [0xe80000,0xe9ffff], candidate outputs [0x8d,0x8d]604 .byte 0x8b // input [0xea0000,0xebffff], candidate outputs [0x8b,0x8c]605 .byte 0x8a // input [0xec0000,0xedffff], candidate outputs [0x8a,0x8a]606 .byte 0x89 // input [0xee0000,0xefffff], candidate outputs [0x89,0x89]607 .byte 0x88 // input [0xf00000,0xf1ffff], candidate outputs [0x88,0x88]608 .byte 0x87 // input [0xf20000,0xf3ffff], candidate outputs [0x87,0x87]609 .byte 0x86 // input [0xf40000,0xf5ffff], candidate outputs [0x86,0x86]610 .byte 0x85 // input [0xf60000,0xf7ffff], candidate outputs [0x85,0x85]611 .byte 0x84 // input [0xf80000,0xf9ffff], candidate outputs [0x84,0x84]612 .byte 0x83 // input [0xfa0000,0xfbffff], candidate outputs [0x83,0x83]613 .byte 0x82 // input [0xfc0000,0xfdffff], candidate outputs [0x82,0x82]614 .byte 0x81 // input [0xfe0000,0xffffff], candidate outputs [0x80,0x81]615 616END_COMPILERRT_FUNCTION(__aeabi_fdiv)617 618NO_EXEC_STACK_DIRECTIVE619