414 lines · plain
1//===----------------------Hexagon builtin routine ------------------------===//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// Double Precision Multiply10#define A r1:011#define AH r112#define AL r013#define B r3:214#define BH r315#define BL r216 17#define BTMP r5:418#define BTMPH r519#define BTMPL r420 21#define PP_ODD r7:622#define PP_ODD_H r723#define PP_ODD_L r624 25#define ONE r9:826#define S_ONE r827#define S_ZERO r928 29#define PP_HH r11:1030#define PP_HH_H r1131#define PP_HH_L r1032 33#define ATMP r13:1234#define ATMPH r1335#define ATMPL r1236 37#define PP_LL r15:1438#define PP_LL_H r1539#define PP_LL_L r1440 41#define TMP r2842 43#define MANTBITS 5244#define HI_MANTBITS 2045#define EXPBITS 1146#define BIAS 102447#define MANTISSA_TO_INT_BIAS 5248 49// Some constant to adjust normalization amount in error code50// Amount to right shift the partial product to get to a denorm51#define FUDGE 552 53#define Q6_ALIAS(TAG) .global __qdsp_##TAG ; .set __qdsp_##TAG, __hexagon_##TAG54#define FAST_ALIAS(TAG) .global __hexagon_fast_##TAG ; .set __hexagon_fast_##TAG, __hexagon_##TAG55#define FAST2_ALIAS(TAG) .global __hexagon_fast2_##TAG ; .set __hexagon_fast2_##TAG, __hexagon_##TAG56#define END(TAG) .size TAG,.-TAG57 58#define SR_ROUND_OFF 2259 .text60 .global __hexagon_muldf361 .type __hexagon_muldf3,@function62 Q6_ALIAS(muldf3)63 FAST_ALIAS(muldf3)64 FAST2_ALIAS(muldf3)65 .p2align 566__hexagon_muldf3:67 {68 p0 = dfclass(A,#2)69 p0 = dfclass(B,#2)70 ATMP = combine(##0x40000000,#0)71 }72 {73 ATMP = insert(A,#MANTBITS,#EXPBITS-1)74 BTMP = asl(B,#EXPBITS-1)75 TMP = #-BIAS76 ONE = #177 }78 {79 PP_ODD = mpyu(BTMPL,ATMPH)80 BTMP = insert(ONE,#2,#62)81 }82 // since we know that the MSB of the H registers is zero, we should never carry83 // H <= 2^31-1. L <= 2^32-1. Therefore, HL <= 2^63-2^32-2^31+184 // Adding 2 HLs, we get 2^64-3*2^32+2 maximum.85 // Therefore, we can add 3 2^32-1 values safely without carry. We only need one.86 {87 PP_LL = mpyu(ATMPL,BTMPL)88 PP_ODD += mpyu(ATMPL,BTMPH)89 }90 {91 PP_ODD += lsr(PP_LL,#32)92 PP_HH = mpyu(ATMPH,BTMPH)93 BTMP = combine(##BIAS+BIAS-4,#0)94 }95 {96 PP_HH += lsr(PP_ODD,#32)97 if (!p0) jump .Lmul_abnormal98 p1 = cmp.eq(PP_LL_L,#0) // 64 lsb's 0?99 p1 = cmp.eq(PP_ODD_L,#0) // 64 lsb's 0?100 }101 102 // PP_HH can have a maximum of 0x3FFF_FFFF_FFFF_FFFF or thereabouts103 // PP_HH can have a minimum of 0x1000_0000_0000_0000 or so104 105#undef PP_ODD106#undef PP_ODD_H107#undef PP_ODD_L108#define EXP10 r7:6109#define EXP1 r7110#define EXP0 r6111 {112 if (!p1) PP_HH_L = or(PP_HH_L,S_ONE)113 EXP0 = extractu(AH,#EXPBITS,#HI_MANTBITS)114 EXP1 = extractu(BH,#EXPBITS,#HI_MANTBITS)115 }116 {117 PP_LL = neg(PP_HH)118 EXP0 += add(TMP,EXP1)119 TMP = xor(AH,BH)120 }121 {122 if (!p2.new) PP_HH = PP_LL123 p2 = cmp.gt(TMP,#-1)124 p0 = !cmp.gt(EXP0,BTMPH)125 p0 = cmp.gt(EXP0,BTMPL)126 if (!p0.new) jump:nt .Lmul_ovf_unf127 }128 {129 A = convert_d2df(PP_HH)130 EXP0 = add(EXP0,#-BIAS-58)131 }132 {133 AH += asl(EXP0,#HI_MANTBITS)134 jumpr r31135 }136 137 .falign138.Lpossible_unf:139 // We end up with a positive exponent140 // But we may have rounded up to an exponent of 1.141 // If the exponent is 1, if we rounded up to it142 // we need to also raise underflow143 // Fortunately, this is pretty easy to detect, we must have +/- 0x0010_0000_0000_0000144 // And the PP should also have more than one bit set145 //146 // Note: ATMP should have abs(PP_HH)147 // Note: BTMPL should have 0x7FEFFFFF148 {149 p0 = cmp.eq(AL,#0)150 p0 = bitsclr(AH,BTMPL)151 if (!p0.new) jumpr:t r31152 BTMPH = #0x7fff153 }154 {155 p0 = bitsset(ATMPH,BTMPH)156 BTMPL = USR157 BTMPH = #0x030158 }159 {160 if (p0) BTMPL = or(BTMPL,BTMPH)161 }162 {163 USR = BTMPL164 }165 {166 p0 = dfcmp.eq(A,A)167 jumpr r31168 }169 .falign170.Lmul_ovf_unf:171 {172 A = convert_d2df(PP_HH)173 ATMP = abs(PP_HH) // take absolute value174 EXP1 = add(EXP0,#-BIAS-58)175 }176 {177 AH += asl(EXP1,#HI_MANTBITS)178 EXP1 = extractu(AH,#EXPBITS,#HI_MANTBITS)179 BTMPL = ##0x7FEFFFFF180 }181 {182 EXP1 += add(EXP0,##-BIAS-58)183 //BTMPH = add(clb(ATMP),#-2)184 BTMPH = #0185 }186 {187 p0 = cmp.gt(EXP1,##BIAS+BIAS-2) // overflow188 if (p0.new) jump:nt .Lmul_ovf189 }190 {191 p0 = cmp.gt(EXP1,#0)192 if (p0.new) jump:nt .Lpossible_unf193 BTMPH = sub(EXP0,BTMPH)194 TMP = #63 // max amount to shift195 }196 // Underflow197 //198 // PP_HH has the partial product with sticky LSB.199 // PP_HH can have a maximum of 0x3FFF_FFFF_FFFF_FFFF or thereabouts200 // PP_HH can have a minimum of 0x1000_0000_0000_0000 or so201 // The exponent of PP_HH is in EXP1, which is non-positive (0 or negative)202 // That's the exponent that happens after the normalization203 //204 // EXP0 has the exponent that, when added to the normalized value, is out of range.205 //206 // Strategy:207 //208 // * Shift down bits, with sticky bit, such that the bits are aligned according209 // to the LZ count and appropriate exponent, but not all the way to mantissa210 // field, keep around the last few bits.211 // * Put a 1 near the MSB212 // * Check the LSBs for inexact; if inexact also set underflow213 // * Convert [u]d2df -- will correctly round according to rounding mode214 // * Replace exponent field with zero215 216 {217 BTMPL = #0 // offset for extract218 BTMPH = sub(#FUDGE,BTMPH) // amount to right shift219 }220 {221 p3 = cmp.gt(PP_HH_H,#-1) // is it positive?222 BTMPH = min(BTMPH,TMP) // Don't shift more than 63223 PP_HH = ATMP224 }225 {226 TMP = USR227 PP_LL = extractu(PP_HH,BTMP)228 }229 {230 PP_HH = asr(PP_HH,BTMPH)231 BTMPL = #0x0030 // underflow flag232 AH = insert(S_ZERO,#EXPBITS,#HI_MANTBITS)233 }234 {235 p0 = cmp.gtu(ONE,PP_LL) // Did we extract all zeros?236 if (!p0.new) PP_HH_L = or(PP_HH_L,S_ONE) // add sticky bit237 PP_HH_H = setbit(PP_HH_H,#HI_MANTBITS+3) // Add back in a bit so we can use convert instruction238 }239 {240 PP_LL = neg(PP_HH)241 p1 = bitsclr(PP_HH_L,#0x7) // Are the LSB's clear?242 if (!p1.new) TMP = or(BTMPL,TMP) // If not, Inexact+Underflow243 }244 {245 if (!p3) PP_HH = PP_LL246 USR = TMP247 }248 {249 A = convert_d2df(PP_HH) // Do rounding250 p0 = dfcmp.eq(A,A) // realize exception251 }252 {253 AH = insert(S_ZERO,#EXPBITS-1,#HI_MANTBITS+1) // Insert correct exponent254 jumpr r31255 }256 .falign257.Lmul_ovf:258 // We get either max finite value or infinity. Either way, overflow+inexact259 {260 TMP = USR261 ATMP = combine(##0x7fefffff,#-1) // positive max finite262 A = PP_HH263 }264 {265 PP_LL_L = extractu(TMP,#2,#SR_ROUND_OFF) // rounding bits266 TMP = or(TMP,#0x28) // inexact + overflow267 BTMP = combine(##0x7ff00000,#0) // positive infinity268 }269 {270 USR = TMP271 PP_LL_L ^= lsr(AH,#31) // Does sign match rounding?272 TMP = PP_LL_L // unmodified rounding mode273 }274 {275 p0 = !cmp.eq(TMP,#1) // If not round-to-zero and276 p0 = !cmp.eq(PP_LL_L,#2) // Not rounding the other way,277 if (p0.new) ATMP = BTMP // we should get infinity278 p0 = dfcmp.eq(A,A) // Realize FP exception if enabled279 }280 {281 A = insert(ATMP,#63,#0) // insert inf/maxfinite, leave sign282 jumpr r31283 }284 285.Lmul_abnormal:286 {287 ATMP = extractu(A,#63,#0) // strip off sign288 BTMP = extractu(B,#63,#0) // strip off sign289 }290 {291 p3 = cmp.gtu(ATMP,BTMP)292 if (!p3.new) A = B // sort values293 if (!p3.new) B = A // sort values294 }295 {296 // Any NaN --> NaN, possibly raise invalid if sNaN297 p0 = dfclass(A,#0x0f) // A not NaN?298 if (!p0.new) jump:nt .Linvalid_nan299 if (!p3) ATMP = BTMP300 if (!p3) BTMP = ATMP301 }302 {303 // Infinity * nonzero number is infinity304 p1 = dfclass(A,#0x08) // A is infinity305 p1 = dfclass(B,#0x0e) // B is nonzero306 }307 {308 // Infinity * zero --> NaN, raise invalid309 // Other zeros return zero310 p0 = dfclass(A,#0x08) // A is infinity311 p0 = dfclass(B,#0x01) // B is zero312 }313 {314 if (p1) jump .Ltrue_inf315 p2 = dfclass(B,#0x01)316 }317 {318 if (p0) jump .Linvalid_zeroinf319 if (p2) jump .Ltrue_zero // so return zero320 TMP = ##0x7c000000321 }322 // We are left with a normal or subnormal times a subnormal. A > B323 // If A and B are both very small (exp(a) < BIAS-MANTBITS),324 // we go to a single sticky bit, which we can round easily.325 // If A and B might multiply to something bigger, decrease A exponent and increase326 // B exponent and try again327 {328 p0 = bitsclr(AH,TMP)329 if (p0.new) jump:nt .Lmul_tiny330 }331 {332 TMP = cl0(BTMP)333 }334 {335 TMP = add(TMP,#-EXPBITS)336 }337 {338 BTMP = asl(BTMP,TMP)339 }340 {341 B = insert(BTMP,#63,#0)342 AH -= asl(TMP,#HI_MANTBITS)343 }344 jump __hexagon_muldf3345.Lmul_tiny:346 {347 TMP = USR348 A = xor(A,B) // get sign bit349 }350 {351 TMP = or(TMP,#0x30) // Inexact + Underflow352 A = insert(ONE,#63,#0) // put in rounded up value353 BTMPH = extractu(TMP,#2,#SR_ROUND_OFF) // get rounding mode354 }355 {356 USR = TMP357 p0 = cmp.gt(BTMPH,#1) // Round towards pos/neg inf?358 if (!p0.new) AL = #0 // If not, zero359 BTMPH ^= lsr(AH,#31) // rounding my way --> set LSB360 }361 {362 p0 = cmp.eq(BTMPH,#3) // if rounding towards right inf363 if (!p0.new) AL = #0 // don't go to zero364 jumpr r31365 }366.Linvalid_zeroinf:367 {368 TMP = USR369 }370 {371 A = #-1372 TMP = or(TMP,#2)373 }374 {375 USR = TMP376 }377 {378 p0 = dfcmp.uo(A,A) // force exception if enabled379 jumpr r31380 }381.Linvalid_nan:382 {383 p0 = dfclass(B,#0x0f) // if B is not NaN384 TMP = convert_df2sf(A) // will generate invalid if sNaN385 if (p0.new) B = A // make it whatever A is386 }387 {388 BL = convert_df2sf(B) // will generate invalid if sNaN389 A = #-1390 jumpr r31391 }392 .falign393.Ltrue_zero:394 {395 A = B396 B = A397 }398.Ltrue_inf:399 {400 BH = extract(BH,#1,#31)401 }402 {403 AH ^= asl(BH,#31)404 jumpr r31405 }406END(__hexagon_muldf3)407 408#undef ATMP409#undef ATMPL410#undef ATMPH411#undef BTMP412#undef BTMPL413#undef BTMPH414