compiler-rt/lib/builtins/hexagon/dfsqrt.S - llvm-project - Git at Google

 //===----------------------Hexagon builtin routine ------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//

 // Double Precision square root

 #define EXP r28

 #define A r1:0
 #define AH r1
 #define AL r0

 #define SFSH r3:2
 #define SF_S r3
 #define SF_H r2

 #define SFHALF_SONE r5:4
 #define S_ONE r4
 #define SFHALF r5
 #define SF_D r6
 #define SF_E r7
 #define RECIPEST r8
 #define SFRAD r9

 #define FRACRAD r11:10
 #define FRACRADH r11
 #define FRACRADL r10

 #define ROOT r13:12
 #define ROOTHI r13
 #define ROOTLO r12

 #define PROD r15:14
 #define PRODHI r15
 #define PRODLO r14

 #define P_TMP p0
 #define P_EXP1 p1
 #define NORMAL p2

 #define SF_EXPBITS 8
 #define SF_MANTBITS 23

 #define DF_EXPBITS 11
 #define DF_MANTBITS 52

 #define DF_BIAS 0x3ff

 #define DFCLASS_ZERO     0x01
 #define DFCLASS_NORMAL   0x02
 #define DFCLASS_DENORMAL 0x02
 #define DFCLASS_INFINITE 0x08
 #define DFCLASS_NAN      0x10

 #define Q6_ALIAS(TAG) .global __qdsp_##TAG ; .set __qdsp_##TAG, __hexagon_##TAG; .type __qdsp_##TAG,@function
 #define FAST_ALIAS(TAG) .global __hexagon_fast_##TAG ; .set __hexagon_fast_##TAG, __hexagon_##TAG; .type __hexagon_fast_##TAG,@function
 #define FAST2_ALIAS(TAG) .global __hexagon_fast2_##TAG ; .set __hexagon_fast2_##TAG, __hexagon_##TAG; .type __hexagon_fast2_##TAG,@function
 #define END(TAG) .size TAG,.-TAG

 	.text
 	.global __hexagon_sqrtdf2
 	.type __hexagon_sqrtdf2,@function
 	.global __hexagon_sqrt
 	.type __hexagon_sqrt,@function
 	Q6_ALIAS(sqrtdf2)
 	Q6_ALIAS(sqrt)
 	FAST_ALIAS(sqrtdf2)
 	FAST_ALIAS(sqrt)
 	FAST2_ALIAS(sqrtdf2)
 	FAST2_ALIAS(sqrt)
 	.type sqrt,@function
 	.p2align 5
 __hexagon_sqrtdf2:
 __hexagon_sqrt:
 	{
 		PROD = extractu(A,#SF_MANTBITS+1,#DF_MANTBITS-SF_MANTBITS)
 		EXP = extractu(AH,#DF_EXPBITS,#DF_MANTBITS-32)
 		SFHALF_SONE = combine(##0x3f000004,#1)
 	}
 	{
 		NORMAL = dfclass(A,#DFCLASS_NORMAL)		// Is it normal
 		NORMAL = cmp.gt(AH,#-1)				// and positive?
 		if (!NORMAL.new) jump:nt .Lsqrt_abnormal
 		SFRAD = or(SFHALF,PRODLO)
 	}
 #undef NORMAL
 .Ldenormal_restart:
 	{
 		FRACRAD = A
 		SF_E,P_TMP = sfinvsqrta(SFRAD)
 		SFHALF = and(SFHALF,#-16)
 		SFSH = #0
 	}
 #undef A
 #undef AH
 #undef AL
 #define ERROR r1:0
 #define ERRORHI r1
 #define ERRORLO r0
 	// SF_E : reciprocal square root
 	// SF_H : half rsqrt
 	// sf_S : square root
 	// SF_D : error term
 	// SFHALF: 0.5
 	{
 		SF_S += sfmpy(SF_E,SFRAD):lib		// s0: root
 		SF_H += sfmpy(SF_E,SFHALF):lib		// h0: 0.5*y0. Could also decrement exponent...
 		SF_D = SFHALF
 #undef SFRAD
 #define SHIFTAMT r9
 		SHIFTAMT = and(EXP,#1)
 	}
 	{
 		SF_D -= sfmpy(SF_S,SF_H):lib		// d0: 0.5-H*S = 0.5-0.5*~1
 		FRACRADH = insert(S_ONE,#DF_EXPBITS+1,#DF_MANTBITS-32)	// replace upper bits with hidden
 		P_EXP1 = cmp.gtu(SHIFTAMT,#0)
 	}
 	{
 		SF_S += sfmpy(SF_S,SF_D):lib		// s1: refine sqrt
 		SF_H += sfmpy(SF_H,SF_D):lib		// h1: refine half-recip
 		SF_D = SFHALF
 		SHIFTAMT = mux(P_EXP1,#8,#9)
 	}
 	{
 		SF_D -= sfmpy(SF_S,SF_H):lib		// d1: error term
 		FRACRAD = asl(FRACRAD,SHIFTAMT)		// Move fracrad bits to right place
 		SHIFTAMT = mux(P_EXP1,#3,#2)
 	}
 	{
 		SF_H += sfmpy(SF_H,SF_D):lib		// d2: rsqrt
 		// cool trick: half of 1/sqrt(x) has same mantissa as 1/sqrt(x).
 		PROD = asl(FRACRAD,SHIFTAMT)		// fracrad<<(2+exp1)
 	}
 	{
 		SF_H = and(SF_H,##0x007fffff)
 	}
 	{
 		SF_H = add(SF_H,##0x00800000 - 3)
 		SHIFTAMT = mux(P_EXP1,#7,#8)
 	}
 	{
 		RECIPEST = asl(SF_H,SHIFTAMT)
 		SHIFTAMT = mux(P_EXP1,#15-(1+1),#15-(1+0))
 	}
 	{
 		ROOT = mpyu(RECIPEST,PRODHI)		// root = mpyu_full(recipest,hi(fracrad<<(2+exp1)))
 	}

 #undef SFSH	// r3:2
 #undef SF_H	// r2
 #undef SF_S	// r3
 #undef S_ONE	// r4
 #undef SFHALF	// r5
 #undef SFHALF_SONE	// r5:4
 #undef SF_D	// r6
 #undef SF_E	// r7

 #define HL r3:2
 #define LL r5:4
 #define HH r7:6

 #undef P_EXP1
 #define P_CARRY0 p1
 #define P_CARRY1 p2
 #define P_CARRY2 p3

 	// Iteration 0
 	// Maybe we can save a cycle by starting with ERROR=asl(fracrad), then as we multiply
 	// We can shift and subtract instead of shift and add?
 	{
 		ERROR = asl(FRACRAD,#15)
 		PROD = mpyu(ROOTHI,ROOTHI)
 		P_CARRY0 = cmp.eq(r0,r0)
 	}
 	{
 		ERROR -= asl(PROD,#15)
 		PROD = mpyu(ROOTHI,ROOTLO)
 		P_CARRY1 = cmp.eq(r0,r0)
 	}
 	{
 		ERROR -= lsr(PROD,#16)
 		P_CARRY2 = cmp.eq(r0,r0)
 	}
 	{
 		ERROR = mpyu(ERRORHI,RECIPEST)
 	}
 	{
 		ROOT += lsr(ERROR,SHIFTAMT)
 		SHIFTAMT = add(SHIFTAMT,#16)
 		ERROR = asl(FRACRAD,#31)		// for next iter
 	}
 	// Iteration 1
 	{
 		PROD = mpyu(ROOTHI,ROOTHI)
 		ERROR -= mpyu(ROOTHI,ROOTLO)	// amount is 31, no shift needed
 	}
 	{
 		ERROR -= asl(PROD,#31)
 		PROD = mpyu(ROOTLO,ROOTLO)
 	}
 	{
 		ERROR -= lsr(PROD,#33)
 	}
 	{
 		ERROR = mpyu(ERRORHI,RECIPEST)
 	}
 	{
 		ROOT += lsr(ERROR,SHIFTAMT)
 		SHIFTAMT = add(SHIFTAMT,#16)
 		ERROR = asl(FRACRAD,#47)	// for next iter
 	}
 	// Iteration 2
 	{
 		PROD = mpyu(ROOTHI,ROOTHI)
 	}
 	{
 		ERROR -= asl(PROD,#47)
 		PROD = mpyu(ROOTHI,ROOTLO)
 	}
 	{
 		ERROR -= asl(PROD,#16)		// bidir shr 31-47
 		PROD = mpyu(ROOTLO,ROOTLO)
 	}
 	{
 		ERROR -= lsr(PROD,#17)		// 64-47
 	}
 	{
 		ERROR = mpyu(ERRORHI,RECIPEST)
 	}
 	{
 		ROOT += lsr(ERROR,SHIFTAMT)
 	}
 #undef ERROR
 #undef PROD
 #undef PRODHI
 #undef PRODLO
 #define REM_HI r15:14
 #define REM_HI_HI r15
 #define REM_LO r1:0
 #undef RECIPEST
 #undef SHIFTAMT
 #define TWOROOT_LO r9:8
 	// Adjust Root
 	{
 		HL = mpyu(ROOTHI,ROOTLO)
 		LL = mpyu(ROOTLO,ROOTLO)
 		REM_HI = #0
 		REM_LO = #0
 	}
 	{
 		HL += lsr(LL,#33)
 		LL += asl(HL,#33)
 		P_CARRY0 = cmp.eq(r0,r0)
 	}
 	{
 		HH = mpyu(ROOTHI,ROOTHI)
 		REM_LO = sub(REM_LO,LL,P_CARRY0):carry
 		TWOROOT_LO = #1
 	}
 	{
 		HH += lsr(HL,#31)
 		TWOROOT_LO += asl(ROOT,#1)
 	}
 #undef HL
 #undef LL
 #define REM_HI_TMP r3:2
 #define REM_HI_TMP_HI r3
 #define REM_LO_TMP r5:4
 	{
 		REM_HI = sub(FRACRAD,HH,P_CARRY0):carry
 		REM_LO_TMP = sub(REM_LO,TWOROOT_LO,P_CARRY1):carry
 #undef FRACRAD
 #undef HH
 #define ZERO r11:10
 #define ONE r7:6
 		ONE = #1
 		ZERO = #0
 	}
 	{
 		REM_HI_TMP = sub(REM_HI,ZERO,P_CARRY1):carry
 		ONE = add(ROOT,ONE)
 		EXP = add(EXP,#-DF_BIAS)			// subtract bias --> signed exp
 	}
 	{
 				// If carry set, no borrow: result was still positive
 		if (P_CARRY1) ROOT = ONE
 		if (P_CARRY1) REM_LO = REM_LO_TMP
 		if (P_CARRY1) REM_HI = REM_HI_TMP
 	}
 	{
 		REM_LO_TMP = sub(REM_LO,TWOROOT_LO,P_CARRY2):carry
 		ONE = #1
 		EXP = asr(EXP,#1)				// divide signed exp by 2
 	}
 	{
 		REM_HI_TMP = sub(REM_HI,ZERO,P_CARRY2):carry
 		ONE = add(ROOT,ONE)
 	}
 	{
 		if (P_CARRY2) ROOT = ONE
 		if (P_CARRY2) REM_LO = REM_LO_TMP
 								// since tworoot <= 2^32, remhi must be zero
 #undef REM_HI_TMP
 #undef REM_HI_TMP_HI
 #define S_ONE r2
 #define ADJ r3
 		S_ONE = #1
 	}
 	{
 		P_TMP = cmp.eq(REM_LO,ZERO)			// is the low part zero
 		if (!P_TMP.new) ROOTLO = or(ROOTLO,S_ONE)	// if so, it's exact... hopefully
 		ADJ = cl0(ROOT)
 		EXP = add(EXP,#-63)
 	}
 #undef REM_LO
 #define RET r1:0
 #define RETHI r1
 	{
 		RET = convert_ud2df(ROOT)			// set up mantissa, maybe set inexact flag
 		EXP = add(EXP,ADJ)				// add back bias
 	}
 	{
 		RETHI += asl(EXP,#DF_MANTBITS-32)		// add exponent adjust
 		jumpr r31
 	}
 #undef REM_LO_TMP
 #undef REM_HI_TMP
 #undef REM_HI_TMP_HI
 #undef REM_LO
 #undef REM_HI
 #undef TWOROOT_LO

 #undef RET
 #define A r1:0
 #define AH r1
 #define AL r1
 #undef S_ONE
 #define TMP r3:2
 #define TMPHI r3
 #define TMPLO r2
 #undef P_CARRY0
 #define P_NEG p1


 #define SFHALF r5
 #define SFRAD r9
 .Lsqrt_abnormal:
 	{
 		P_TMP = dfclass(A,#DFCLASS_ZERO)			// zero?
 		if (P_TMP.new) jumpr:t r31
 	}
 	{
 		P_TMP = dfclass(A,#DFCLASS_NAN)
 		if (P_TMP.new) jump:nt .Lsqrt_nan
 	}
 	{
 		P_TMP = cmp.gt(AH,#-1)
 		if (!P_TMP.new) jump:nt .Lsqrt_invalid_neg
 		if (!P_TMP.new) EXP = ##0x7F800001			// sNaN
 	}
 	{
 		P_TMP = dfclass(A,#DFCLASS_INFINITE)
 		if (P_TMP.new) jumpr:nt r31
 	}
 	// If we got here, we're denormal
 	// prepare to restart
 	{
 		A = extractu(A,#DF_MANTBITS,#0)		// Extract mantissa
 	}
 	{
 		EXP = add(clb(A),#-DF_EXPBITS)		// how much to normalize?
 	}
 	{
 		A = asl(A,EXP)				// Shift mantissa
 		EXP = sub(#1,EXP)			// Form exponent
 	}
 	{
 		AH = insert(EXP,#1,#DF_MANTBITS-32)		// insert lsb of exponent
 	}
 	{
 		TMP = extractu(A,#SF_MANTBITS+1,#DF_MANTBITS-SF_MANTBITS)	// get sf value (mant+exp1)
 		SFHALF = ##0x3f000004						// form half constant
 	}
 	{
 		SFRAD = or(SFHALF,TMPLO)			// form sf value
 		SFHALF = and(SFHALF,#-16)
 		jump .Ldenormal_restart				// restart
 	}
 .Lsqrt_nan:
 	{
 		EXP = convert_df2sf(A)				// if sNaN, get invalid
 		A = #-1						// qNaN
 		jumpr r31
 	}
 .Lsqrt_invalid_neg:
 	{
 		A = convert_sf2df(EXP)				// Invalid,NaNval
 		jumpr r31
 	}
 END(__hexagon_sqrt)
 END(__hexagon_sqrtdf2)
	//===----------------------Hexagon builtin routine ------------------------===//
	//
	// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
	// See https://llvm.org/LICENSE.txt for license information.
	// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
	//
	//===----------------------------------------------------------------------===//

	// Double Precision square root

	#define EXP r28

	#define A r1:0
	#define AH r1
	#define AL r0

	#define SFSH r3:2
	#define SF_S r3
	#define SF_H r2

	#define SFHALF_SONE r5:4
	#define S_ONE r4
	#define SFHALF r5
	#define SF_D r6
	#define SF_E r7
	#define RECIPEST r8
	#define SFRAD r9

	#define FRACRAD r11:10
	#define FRACRADH r11
	#define FRACRADL r10

	#define ROOT r13:12
	#define ROOTHI r13
	#define ROOTLO r12

	#define PROD r15:14
	#define PRODHI r15
	#define PRODLO r14

	#define P_TMP p0
	#define P_EXP1 p1
	#define NORMAL p2

	#define SF_EXPBITS 8
	#define SF_MANTBITS 23

	#define DF_EXPBITS 11
	#define DF_MANTBITS 52

	#define DF_BIAS 0x3ff

	#define DFCLASS_ZERO 0x01
	#define DFCLASS_NORMAL 0x02
	#define DFCLASS_DENORMAL 0x02
	#define DFCLASS_INFINITE 0x08
	#define DFCLASS_NAN 0x10

	#define Q6_ALIAS(TAG) .global __qdsp_##TAG ; .set __qdsp_##TAG, __hexagon_##TAG; .type __qdsp_##TAG,@function
	#define FAST_ALIAS(TAG) .global __hexagon_fast_##TAG ; .set __hexagon_fast_##TAG, __hexagon_##TAG; .type __hexagon_fast_##TAG,@function
	#define FAST2_ALIAS(TAG) .global __hexagon_fast2_##TAG ; .set __hexagon_fast2_##TAG, __hexagon_##TAG; .type __hexagon_fast2_##TAG,@function
	#define END(TAG) .size TAG,.-TAG

	.text
	.global __hexagon_sqrtdf2
	.type __hexagon_sqrtdf2,@function
	.global __hexagon_sqrt
	.type __hexagon_sqrt,@function
	Q6_ALIAS(sqrtdf2)
	Q6_ALIAS(sqrt)
	FAST_ALIAS(sqrtdf2)
	FAST_ALIAS(sqrt)
	FAST2_ALIAS(sqrtdf2)
	FAST2_ALIAS(sqrt)
	.type sqrt,@function
	.p2align 5
	__hexagon_sqrtdf2:
	__hexagon_sqrt:
	{
	PROD = extractu(A,#SF_MANTBITS+1,#DF_MANTBITS-SF_MANTBITS)
	EXP = extractu(AH,#DF_EXPBITS,#DF_MANTBITS-32)
	SFHALF_SONE = combine(##0x3f000004,#1)
	}
	{
	NORMAL = dfclass(A,#DFCLASS_NORMAL) // Is it normal
	NORMAL = cmp.gt(AH,#-1) // and positive?
	if (!NORMAL.new) jump:nt .Lsqrt_abnormal
	SFRAD = or(SFHALF,PRODLO)
	}
	#undef NORMAL
	.Ldenormal_restart:
	{
	FRACRAD = A
	SF_E,P_TMP = sfinvsqrta(SFRAD)
	SFHALF = and(SFHALF,#-16)
	SFSH = #0
	}
	#undef A
	#undef AH
	#undef AL
	#define ERROR r1:0
	#define ERRORHI r1
	#define ERRORLO r0
	// SF_E : reciprocal square root
	// SF_H : half rsqrt
	// sf_S : square root
	// SF_D : error term
	// SFHALF: 0.5
	{
	SF_S += sfmpy(SF_E,SFRAD):lib // s0: root
	SF_H += sfmpy(SF_E,SFHALF):lib // h0: 0.5*y0. Could also decrement exponent...
	SF_D = SFHALF
	#undef SFRAD
	#define SHIFTAMT r9
	SHIFTAMT = and(EXP,#1)
	}
	{
	SF_D -= sfmpy(SF_S,SF_H):lib // d0: 0.5-HS = 0.5-0.5~1
	FRACRADH = insert(S_ONE,#DF_EXPBITS+1,#DF_MANTBITS-32) // replace upper bits with hidden
	P_EXP1 = cmp.gtu(SHIFTAMT,#0)
	}
	{
	SF_S += sfmpy(SF_S,SF_D):lib // s1: refine sqrt
	SF_H += sfmpy(SF_H,SF_D):lib // h1: refine half-recip
	SF_D = SFHALF
	SHIFTAMT = mux(P_EXP1,#8,#9)
	}
	{
	SF_D -= sfmpy(SF_S,SF_H):lib // d1: error term
	FRACRAD = asl(FRACRAD,SHIFTAMT) // Move fracrad bits to right place
	SHIFTAMT = mux(P_EXP1,#3,#2)
	}
	{
	SF_H += sfmpy(SF_H,SF_D):lib // d2: rsqrt
	// cool trick: half of 1/sqrt(x) has same mantissa as 1/sqrt(x).
	PROD = asl(FRACRAD,SHIFTAMT) // fracrad<<(2+exp1)
	}
	{
	SF_H = and(SF_H,##0x007fffff)
	}
	{
	SF_H = add(SF_H,##0x00800000 - 3)
	SHIFTAMT = mux(P_EXP1,#7,#8)
	}
	{
	RECIPEST = asl(SF_H,SHIFTAMT)
	SHIFTAMT = mux(P_EXP1,#15-(1+1),#15-(1+0))
	}
	{
	ROOT = mpyu(RECIPEST,PRODHI) // root = mpyu_full(recipest,hi(fracrad<<(2+exp1)))
	}

	#undef SFSH // r3:2
	#undef SF_H // r2
	#undef SF_S // r3
	#undef S_ONE // r4
	#undef SFHALF // r5
	#undef SFHALF_SONE // r5:4
	#undef SF_D // r6
	#undef SF_E // r7

	#define HL r3:2
	#define LL r5:4
	#define HH r7:6

	#undef P_EXP1
	#define P_CARRY0 p1
	#define P_CARRY1 p2
	#define P_CARRY2 p3

	// Iteration 0
	// Maybe we can save a cycle by starting with ERROR=asl(fracrad), then as we multiply
	// We can shift and subtract instead of shift and add?
	{
	ERROR = asl(FRACRAD,#15)
	PROD = mpyu(ROOTHI,ROOTHI)
	P_CARRY0 = cmp.eq(r0,r0)
	}
	{
	ERROR -= asl(PROD,#15)
	PROD = mpyu(ROOTHI,ROOTLO)
	P_CARRY1 = cmp.eq(r0,r0)
	}
	{
	ERROR -= lsr(PROD,#16)
	P_CARRY2 = cmp.eq(r0,r0)
	}
	{
	ERROR = mpyu(ERRORHI,RECIPEST)
	}
	{
	ROOT += lsr(ERROR,SHIFTAMT)
	SHIFTAMT = add(SHIFTAMT,#16)
	ERROR = asl(FRACRAD,#31) // for next iter
	}
	// Iteration 1
	{
	PROD = mpyu(ROOTHI,ROOTHI)
	ERROR -= mpyu(ROOTHI,ROOTLO) // amount is 31, no shift needed
	}
	{
	ERROR -= asl(PROD,#31)
	PROD = mpyu(ROOTLO,ROOTLO)
	}
	{
	ERROR -= lsr(PROD,#33)
	}
	{
	ERROR = mpyu(ERRORHI,RECIPEST)
	}
	{
	ROOT += lsr(ERROR,SHIFTAMT)
	SHIFTAMT = add(SHIFTAMT,#16)
	ERROR = asl(FRACRAD,#47) // for next iter
	}
	// Iteration 2
	{
	PROD = mpyu(ROOTHI,ROOTHI)
	}
	{
	ERROR -= asl(PROD,#47)
	PROD = mpyu(ROOTHI,ROOTLO)
	}
	{
	ERROR -= asl(PROD,#16) // bidir shr 31-47
	PROD = mpyu(ROOTLO,ROOTLO)
	}
	{
	ERROR -= lsr(PROD,#17) // 64-47
	}
	{
	ERROR = mpyu(ERRORHI,RECIPEST)
	}
	{
	ROOT += lsr(ERROR,SHIFTAMT)
	}
	#undef ERROR
	#undef PROD
	#undef PRODHI
	#undef PRODLO
	#define REM_HI r15:14
	#define REM_HI_HI r15
	#define REM_LO r1:0
	#undef RECIPEST
	#undef SHIFTAMT
	#define TWOROOT_LO r9:8
	// Adjust Root
	{
	HL = mpyu(ROOTHI,ROOTLO)
	LL = mpyu(ROOTLO,ROOTLO)
	REM_HI = #0
	REM_LO = #0
	}
	{
	HL += lsr(LL,#33)
	LL += asl(HL,#33)
	P_CARRY0 = cmp.eq(r0,r0)
	}
	{
	HH = mpyu(ROOTHI,ROOTHI)
	REM_LO = sub(REM_LO,LL,P_CARRY0):carry
	TWOROOT_LO = #1
	}
	{
	HH += lsr(HL,#31)
	TWOROOT_LO += asl(ROOT,#1)
	}
	#undef HL
	#undef LL
	#define REM_HI_TMP r3:2
	#define REM_HI_TMP_HI r3
	#define REM_LO_TMP r5:4
	{
	REM_HI = sub(FRACRAD,HH,P_CARRY0):carry
	REM_LO_TMP = sub(REM_LO,TWOROOT_LO,P_CARRY1):carry
	#undef FRACRAD
	#undef HH
	#define ZERO r11:10
	#define ONE r7:6
	ONE = #1
	ZERO = #0
	}
	{
	REM_HI_TMP = sub(REM_HI,ZERO,P_CARRY1):carry
	ONE = add(ROOT,ONE)
	EXP = add(EXP,#-DF_BIAS) // subtract bias --> signed exp
	}
	{
	// If carry set, no borrow: result was still positive
	if (P_CARRY1) ROOT = ONE
	if (P_CARRY1) REM_LO = REM_LO_TMP
	if (P_CARRY1) REM_HI = REM_HI_TMP
	}
	{
	REM_LO_TMP = sub(REM_LO,TWOROOT_LO,P_CARRY2):carry
	ONE = #1
	EXP = asr(EXP,#1) // divide signed exp by 2
	}
	{
	REM_HI_TMP = sub(REM_HI,ZERO,P_CARRY2):carry
	ONE = add(ROOT,ONE)
	}
	{
	if (P_CARRY2) ROOT = ONE
	if (P_CARRY2) REM_LO = REM_LO_TMP
	// since tworoot <= 2^32, remhi must be zero
	#undef REM_HI_TMP
	#undef REM_HI_TMP_HI
	#define S_ONE r2
	#define ADJ r3
	S_ONE = #1
	}
	{
	P_TMP = cmp.eq(REM_LO,ZERO) // is the low part zero
	if (!P_TMP.new) ROOTLO = or(ROOTLO,S_ONE) // if so, it's exact... hopefully
	ADJ = cl0(ROOT)
	EXP = add(EXP,#-63)
	}
	#undef REM_LO
	#define RET r1:0
	#define RETHI r1
	{
	RET = convert_ud2df(ROOT) // set up mantissa, maybe set inexact flag
	EXP = add(EXP,ADJ) // add back bias
	}
	{
	RETHI += asl(EXP,#DF_MANTBITS-32) // add exponent adjust
	jumpr r31
	}
	#undef REM_LO_TMP
	#undef REM_HI_TMP
	#undef REM_HI_TMP_HI
	#undef REM_LO
	#undef REM_HI
	#undef TWOROOT_LO

	#undef RET
	#define A r1:0
	#define AH r1
	#define AL r1
	#undef S_ONE
	#define TMP r3:2
	#define TMPHI r3
	#define TMPLO r2
	#undef P_CARRY0
	#define P_NEG p1


	#define SFHALF r5
	#define SFRAD r9
	.Lsqrt_abnormal:
	{
	P_TMP = dfclass(A,#DFCLASS_ZERO) // zero?
	if (P_TMP.new) jumpr:t r31
	}
	{
	P_TMP = dfclass(A,#DFCLASS_NAN)
	if (P_TMP.new) jump:nt .Lsqrt_nan
	}
	{
	P_TMP = cmp.gt(AH,#-1)
	if (!P_TMP.new) jump:nt .Lsqrt_invalid_neg
	if (!P_TMP.new) EXP = ##0x7F800001 // sNaN
	}
	{
	P_TMP = dfclass(A,#DFCLASS_INFINITE)
	if (P_TMP.new) jumpr:nt r31
	}
	// If we got here, we're denormal
	// prepare to restart
	{
	A = extractu(A,#DF_MANTBITS,#0) // Extract mantissa
	}
	{
	EXP = add(clb(A),#-DF_EXPBITS) // how much to normalize?
	}
	{
	A = asl(A,EXP) // Shift mantissa
	EXP = sub(#1,EXP) // Form exponent
	}
	{
	AH = insert(EXP,#1,#DF_MANTBITS-32) // insert lsb of exponent
	}
	{
	TMP = extractu(A,#SF_MANTBITS+1,#DF_MANTBITS-SF_MANTBITS) // get sf value (mant+exp1)
	SFHALF = ##0x3f000004 // form half constant
	}
	{
	SFRAD = or(SFHALF,TMPLO) // form sf value
	SFHALF = and(SFHALF,#-16)
	jump .Ldenormal_restart // restart
	}
	.Lsqrt_nan:
	{
	EXP = convert_df2sf(A) // if sNaN, get invalid
	A = #-1 // qNaN
	jumpr r31
	}
	.Lsqrt_invalid_neg:
	{
	A = convert_sf2df(EXP) // Invalid,NaNval
	jumpr r31
	}
	END(__hexagon_sqrt)
	END(__hexagon_sqrtdf2)