406 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 square root10 11#define EXP r2812 13#define A r1:014#define AH r115#define AL r016 17#define SFSH r3:218#define SF_S r319#define SF_H r220 21#define SFHALF_SONE r5:422#define S_ONE r423#define SFHALF r524#define SF_D r625#define SF_E r726#define RECIPEST r827#define SFRAD r928 29#define FRACRAD r11:1030#define FRACRADH r1131#define FRACRADL r1032 33#define ROOT r13:1234#define ROOTHI r1335#define ROOTLO r1236 37#define PROD r15:1438#define PRODHI r1539#define PRODLO r1440 41#define P_TMP p042#define P_EXP1 p143#define NORMAL p244 45#define SF_EXPBITS 846#define SF_MANTBITS 2347 48#define DF_EXPBITS 1149#define DF_MANTBITS 5250 51#define DF_BIAS 0x3ff52 53#define DFCLASS_ZERO 0x0154#define DFCLASS_NORMAL 0x0255#define DFCLASS_DENORMAL 0x0256#define DFCLASS_INFINITE 0x0857#define DFCLASS_NAN 0x1058 59#define Q6_ALIAS(TAG) .global __qdsp_##TAG ; .set __qdsp_##TAG, __hexagon_##TAG; .type __qdsp_##TAG,@function60#define FAST_ALIAS(TAG) .global __hexagon_fast_##TAG ; .set __hexagon_fast_##TAG, __hexagon_##TAG; .type __hexagon_fast_##TAG,@function61#define FAST2_ALIAS(TAG) .global __hexagon_fast2_##TAG ; .set __hexagon_fast2_##TAG, __hexagon_##TAG; .type __hexagon_fast2_##TAG,@function62#define END(TAG) .size TAG,.-TAG63 64 .text65 .global __hexagon_sqrtdf266 .type __hexagon_sqrtdf2,@function67 .global __hexagon_sqrt68 .type __hexagon_sqrt,@function69 Q6_ALIAS(sqrtdf2)70 Q6_ALIAS(sqrt)71 FAST_ALIAS(sqrtdf2)72 FAST_ALIAS(sqrt)73 FAST2_ALIAS(sqrtdf2)74 FAST2_ALIAS(sqrt)75 .type sqrt,@function76 .p2align 577__hexagon_sqrtdf2:78__hexagon_sqrt:79 {80 PROD = extractu(A,#SF_MANTBITS+1,#DF_MANTBITS-SF_MANTBITS)81 EXP = extractu(AH,#DF_EXPBITS,#DF_MANTBITS-32)82 SFHALF_SONE = combine(##0x3f000004,#1)83 }84 {85 NORMAL = dfclass(A,#DFCLASS_NORMAL) // Is it normal86 NORMAL = cmp.gt(AH,#-1) // and positive?87 if (!NORMAL.new) jump:nt .Lsqrt_abnormal88 SFRAD = or(SFHALF,PRODLO)89 }90#undef NORMAL91.Ldenormal_restart:92 {93 FRACRAD = A94 SF_E,P_TMP = sfinvsqrta(SFRAD)95 SFHALF = and(SFHALF,#-16)96 SFSH = #097 }98#undef A99#undef AH100#undef AL101#define ERROR r1:0102#define ERRORHI r1103#define ERRORLO r0104 // SF_E : reciprocal square root105 // SF_H : half rsqrt106 // sf_S : square root107 // SF_D : error term108 // SFHALF: 0.5109 {110 SF_S += sfmpy(SF_E,SFRAD):lib // s0: root111 SF_H += sfmpy(SF_E,SFHALF):lib // h0: 0.5*y0. Could also decrement exponent...112 SF_D = SFHALF113#undef SFRAD114#define SHIFTAMT r9115 SHIFTAMT = and(EXP,#1)116 }117 {118 SF_D -= sfmpy(SF_S,SF_H):lib // d0: 0.5-H*S = 0.5-0.5*~1119 FRACRADH = insert(S_ONE,#DF_EXPBITS+1,#DF_MANTBITS-32) // replace upper bits with hidden120 P_EXP1 = cmp.gtu(SHIFTAMT,#0)121 }122 {123 SF_S += sfmpy(SF_S,SF_D):lib // s1: refine sqrt124 SF_H += sfmpy(SF_H,SF_D):lib // h1: refine half-recip125 SF_D = SFHALF126 SHIFTAMT = mux(P_EXP1,#8,#9)127 }128 {129 SF_D -= sfmpy(SF_S,SF_H):lib // d1: error term130 FRACRAD = asl(FRACRAD,SHIFTAMT) // Move fracrad bits to right place131 SHIFTAMT = mux(P_EXP1,#3,#2)132 }133 {134 SF_H += sfmpy(SF_H,SF_D):lib // d2: rsqrt135 // cool trick: half of 1/sqrt(x) has same mantissa as 1/sqrt(x).136 PROD = asl(FRACRAD,SHIFTAMT) // fracrad<<(2+exp1)137 }138 {139 SF_H = and(SF_H,##0x007fffff)140 }141 {142 SF_H = add(SF_H,##0x00800000 - 3)143 SHIFTAMT = mux(P_EXP1,#7,#8)144 }145 {146 RECIPEST = asl(SF_H,SHIFTAMT)147 SHIFTAMT = mux(P_EXP1,#15-(1+1),#15-(1+0))148 }149 {150 ROOT = mpyu(RECIPEST,PRODHI) // root = mpyu_full(recipest,hi(fracrad<<(2+exp1)))151 }152 153#undef SFSH // r3:2154#undef SF_H // r2155#undef SF_S // r3156#undef S_ONE // r4157#undef SFHALF // r5158#undef SFHALF_SONE // r5:4159#undef SF_D // r6160#undef SF_E // r7161 162#define HL r3:2163#define LL r5:4164#define HH r7:6165 166#undef P_EXP1167#define P_CARRY0 p1168#define P_CARRY1 p2169#define P_CARRY2 p3170 171 // Iteration 0172 // Maybe we can save a cycle by starting with ERROR=asl(fracrad), then as we multiply173 // We can shift and subtract instead of shift and add?174 {175 ERROR = asl(FRACRAD,#15)176 PROD = mpyu(ROOTHI,ROOTHI)177 P_CARRY0 = cmp.eq(r0,r0)178 }179 {180 ERROR -= asl(PROD,#15)181 PROD = mpyu(ROOTHI,ROOTLO)182 P_CARRY1 = cmp.eq(r0,r0)183 }184 {185 ERROR -= lsr(PROD,#16)186 P_CARRY2 = cmp.eq(r0,r0)187 }188 {189 ERROR = mpyu(ERRORHI,RECIPEST)190 }191 {192 ROOT += lsr(ERROR,SHIFTAMT)193 SHIFTAMT = add(SHIFTAMT,#16)194 ERROR = asl(FRACRAD,#31) // for next iter195 }196 // Iteration 1197 {198 PROD = mpyu(ROOTHI,ROOTHI)199 ERROR -= mpyu(ROOTHI,ROOTLO) // amount is 31, no shift needed200 }201 {202 ERROR -= asl(PROD,#31)203 PROD = mpyu(ROOTLO,ROOTLO)204 }205 {206 ERROR -= lsr(PROD,#33)207 }208 {209 ERROR = mpyu(ERRORHI,RECIPEST)210 }211 {212 ROOT += lsr(ERROR,SHIFTAMT)213 SHIFTAMT = add(SHIFTAMT,#16)214 ERROR = asl(FRACRAD,#47) // for next iter215 }216 // Iteration 2217 {218 PROD = mpyu(ROOTHI,ROOTHI)219 }220 {221 ERROR -= asl(PROD,#47)222 PROD = mpyu(ROOTHI,ROOTLO)223 }224 {225 ERROR -= asl(PROD,#16) // bidir shr 31-47226 PROD = mpyu(ROOTLO,ROOTLO)227 }228 {229 ERROR -= lsr(PROD,#17) // 64-47230 }231 {232 ERROR = mpyu(ERRORHI,RECIPEST)233 }234 {235 ROOT += lsr(ERROR,SHIFTAMT)236 }237#undef ERROR238#undef PROD239#undef PRODHI240#undef PRODLO241#define REM_HI r15:14242#define REM_HI_HI r15243#define REM_LO r1:0244#undef RECIPEST245#undef SHIFTAMT246#define TWOROOT_LO r9:8247 // Adjust Root248 {249 HL = mpyu(ROOTHI,ROOTLO)250 LL = mpyu(ROOTLO,ROOTLO)251 REM_HI = #0252 REM_LO = #0253 }254 {255 HL += lsr(LL,#33)256 LL += asl(HL,#33)257 P_CARRY0 = cmp.eq(r0,r0)258 }259 {260 HH = mpyu(ROOTHI,ROOTHI)261 REM_LO = sub(REM_LO,LL,P_CARRY0):carry262 TWOROOT_LO = #1263 }264 {265 HH += lsr(HL,#31)266 TWOROOT_LO += asl(ROOT,#1)267 }268#undef HL269#undef LL270#define REM_HI_TMP r3:2271#define REM_HI_TMP_HI r3272#define REM_LO_TMP r5:4273 {274 REM_HI = sub(FRACRAD,HH,P_CARRY0):carry275 REM_LO_TMP = sub(REM_LO,TWOROOT_LO,P_CARRY1):carry276#undef FRACRAD277#undef HH278#define ZERO r11:10279#define ONE r7:6280 ONE = #1281 ZERO = #0282 }283 {284 REM_HI_TMP = sub(REM_HI,ZERO,P_CARRY1):carry285 ONE = add(ROOT,ONE)286 EXP = add(EXP,#-DF_BIAS) // subtract bias --> signed exp287 }288 {289 // If carry set, no borrow: result was still positive290 if (P_CARRY1) ROOT = ONE291 if (P_CARRY1) REM_LO = REM_LO_TMP292 if (P_CARRY1) REM_HI = REM_HI_TMP293 }294 {295 REM_LO_TMP = sub(REM_LO,TWOROOT_LO,P_CARRY2):carry296 ONE = #1297 EXP = asr(EXP,#1) // divide signed exp by 2298 }299 {300 REM_HI_TMP = sub(REM_HI,ZERO,P_CARRY2):carry301 ONE = add(ROOT,ONE)302 }303 {304 if (P_CARRY2) ROOT = ONE305 if (P_CARRY2) REM_LO = REM_LO_TMP306 // since tworoot <= 2^32, remhi must be zero307#undef REM_HI_TMP308#undef REM_HI_TMP_HI309#define S_ONE r2310#define ADJ r3311 S_ONE = #1312 }313 {314 P_TMP = cmp.eq(REM_LO,ZERO) // is the low part zero315 if (!P_TMP.new) ROOTLO = or(ROOTLO,S_ONE) // if so, it's exact... hopefully316 ADJ = cl0(ROOT)317 EXP = add(EXP,#-63)318 }319#undef REM_LO320#define RET r1:0321#define RETHI r1322 {323 RET = convert_ud2df(ROOT) // set up mantissa, maybe set inexact flag324 EXP = add(EXP,ADJ) // add back bias325 }326 {327 RETHI += asl(EXP,#DF_MANTBITS-32) // add exponent adjust328 jumpr r31329 }330#undef REM_LO_TMP331#undef REM_HI_TMP332#undef REM_HI_TMP_HI333#undef REM_LO334#undef REM_HI335#undef TWOROOT_LO336 337#undef RET338#define A r1:0339#define AH r1340#define AL r1341#undef S_ONE342#define TMP r3:2343#define TMPHI r3344#define TMPLO r2345#undef P_CARRY0346#define P_NEG p1347 348 349#define SFHALF r5350#define SFRAD r9351.Lsqrt_abnormal:352 {353 P_TMP = dfclass(A,#DFCLASS_ZERO) // zero?354 if (P_TMP.new) jumpr:t r31355 }356 {357 P_TMP = dfclass(A,#DFCLASS_NAN)358 if (P_TMP.new) jump:nt .Lsqrt_nan359 }360 {361 P_TMP = cmp.gt(AH,#-1)362 if (!P_TMP.new) jump:nt .Lsqrt_invalid_neg363 if (!P_TMP.new) EXP = ##0x7F800001 // sNaN364 }365 {366 P_TMP = dfclass(A,#DFCLASS_INFINITE)367 if (P_TMP.new) jumpr:nt r31368 }369 // If we got here, we're denormal370 // prepare to restart371 {372 A = extractu(A,#DF_MANTBITS,#0) // Extract mantissa373 }374 {375 EXP = add(clb(A),#-DF_EXPBITS) // how much to normalize?376 }377 {378 A = asl(A,EXP) // Shift mantissa379 EXP = sub(#1,EXP) // Form exponent380 }381 {382 AH = insert(EXP,#1,#DF_MANTBITS-32) // insert lsb of exponent383 }384 {385 TMP = extractu(A,#SF_MANTBITS+1,#DF_MANTBITS-SF_MANTBITS) // get sf value (mant+exp1)386 SFHALF = ##0x3f000004 // form half constant387 }388 {389 SFRAD = or(SFHALF,TMPLO) // form sf value390 SFHALF = and(SFHALF,#-16)391 jump .Ldenormal_restart // restart392 }393.Lsqrt_nan:394 {395 EXP = convert_df2sf(A) // if sNaN, get invalid396 A = #-1 // qNaN397 jumpr r31398 }399.Lsqrt_invalid_neg:400 {401 A = convert_sf2df(EXP) // Invalid,NaNval402 jumpr r31403 }404END(__hexagon_sqrt)405END(__hexagon_sqrtdf2)406