//===----------------------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 Multiply
#define A r1:0
#define AH r1
#define AL r0
#define B r3:2
#define BH r3
#define BL r2
#define BTMP r5:4
#define BTMPH r5
#define BTMPL r4
#define PP_ODD r7:6
#define PP_ODD_H r7
#define PP_ODD_L r6
#define ONE r9:8
#define S_ONE r8
#define S_ZERO r9
#define PP_HH r11:10
#define PP_HH_H r11
#define PP_HH_L r10
#define ATMP r13:12
#define ATMPH r13
#define ATMPL r12
#define PP_LL r15:14
#define PP_LL_H r15
#define PP_LL_L r14
#define TMP r28
#define MANTBITS 52
#define HI_MANTBITS 20
#define EXPBITS 11
#define BIAS 1024
#define MANTISSA_TO_INT_BIAS 52
// Some constant to adjust normalization amount in error code
// Amount to right shift the partial product to get to a denorm
#define FUDGE 5
#define Q6_ALIAS(TAG) .global __qdsp_##TAG ; .set __qdsp_##TAG, __hexagon_##TAG
#define FAST_ALIAS(TAG) .global __hexagon_fast_##TAG ; .set __hexagon_fast_##TAG, __hexagon_##TAG
#define FAST2_ALIAS(TAG) .global __hexagon_fast2_##TAG ; .set __hexagon_fast2_##TAG, __hexagon_##TAG
#define END(TAG) .size TAG,.-TAG
#define SR_ROUND_OFF 22
.text
.global __hexagon_muldf3
.type __hexagon_muldf3,@function
Q6_ALIAS(muldf3)
FAST_ALIAS(muldf3)
FAST2_ALIAS(muldf3)
.p2align 5
__hexagon_muldf3:
{
p0 = dfclass(A,#2)
p0 = dfclass(B,#2)
ATMP = combine(##0x40000000,#0)
}
{
ATMP = insert(A,#MANTBITS,#EXPBITS-1)
BTMP = asl(B,#EXPBITS-1)
TMP = #-BIAS
ONE = #1
}
{
PP_ODD = mpyu(BTMPL,ATMPH)
BTMP = insert(ONE,#2,#62)
}
// since we know that the MSB of the H registers is zero, we should never carry
// H <= 2^31-1. L <= 2^32-1. Therefore, HL <= 2^63-2^32-2^31+1
// Adding 2 HLs, we get 2^64-3*2^32+2 maximum.
// Therefore, we can add 3 2^32-1 values safely without carry. We only need one.
{
PP_LL = mpyu(ATMPL,BTMPL)
PP_ODD += mpyu(ATMPL,BTMPH)
}
{
PP_ODD += lsr(PP_LL,#32)
PP_HH = mpyu(ATMPH,BTMPH)
BTMP = combine(##BIAS+BIAS-4,#0)
}
{
PP_HH += lsr(PP_ODD,#32)
if (!p0) jump .Lmul_abnormal
p1 = cmp.eq(PP_LL_L,#0) // 64 lsb's 0?
p1 = cmp.eq(PP_ODD_L,#0) // 64 lsb's 0?
}
// PP_HH can have a maximum of 0x3FFF_FFFF_FFFF_FFFF or thereabouts
// PP_HH can have a minimum of 0x1000_0000_0000_0000 or so
#undef PP_ODD
#undef PP_ODD_H
#undef PP_ODD_L
#define EXP10 r7:6
#define EXP1 r7
#define EXP0 r6
{
if (!p1) PP_HH_L = or(PP_HH_L,S_ONE)
EXP0 = extractu(AH,#EXPBITS,#HI_MANTBITS)
EXP1 = extractu(BH,#EXPBITS,#HI_MANTBITS)
}
{
PP_LL = neg(PP_HH)
EXP0 += add(TMP,EXP1)
TMP = xor(AH,BH)
}
{
if (!p2.new) PP_HH = PP_LL
p2 = cmp.gt(TMP,#-1)
p0 = !cmp.gt(EXP0,BTMPH)
p0 = cmp.gt(EXP0,BTMPL)
if (!p0.new) jump:nt .Lmul_ovf_unf
}
{
A = convert_d2df(PP_HH)
EXP0 = add(EXP0,#-BIAS-58)
}
{
AH += asl(EXP0,#HI_MANTBITS)
jumpr r31
}
.falign
.Lpossible_unf:
// We end up with a positive exponent
// But we may have rounded up to an exponent of 1.
// If the exponent is 1, if we rounded up to it
// we need to also raise underflow
// Fortunately, this is pretty easy to detect, we must have +/- 0x0010_0000_0000_0000
// And the PP should also have more than one bit set
//
// Note: ATMP should have abs(PP_HH)
// Note: BTMPL should have 0x7FEFFFFF
{
p0 = cmp.eq(AL,#0)
p0 = bitsclr(AH,BTMPL)
if (!p0.new) jumpr:t r31
BTMPH = #0x7fff
}
{
p0 = bitsset(ATMPH,BTMPH)
BTMPL = USR
BTMPH = #0x030
}
{
if (p0) BTMPL = or(BTMPL,BTMPH)
}
{
USR = BTMPL
}
{
p0 = dfcmp.eq(A,A)
jumpr r31
}
.falign
.Lmul_ovf_unf:
{
A = convert_d2df(PP_HH)
ATMP = abs(PP_HH) // take absolute value
EXP1 = add(EXP0,#-BIAS-58)
}
{
AH += asl(EXP1,#HI_MANTBITS)
EXP1 = extractu(AH,#EXPBITS,#HI_MANTBITS)
BTMPL = ##0x7FEFFFFF
}
{
EXP1 += add(EXP0,##-BIAS-58)
//BTMPH = add(clb(ATMP),#-2)
BTMPH = #0
}
{
p0 = cmp.gt(EXP1,##BIAS+BIAS-2) // overflow
if (p0.new) jump:nt .Lmul_ovf
}
{
p0 = cmp.gt(EXP1,#0)
if (p0.new) jump:nt .Lpossible_unf
BTMPH = sub(EXP0,BTMPH)
TMP = #63 // max amount to shift
}
// Underflow
//
// PP_HH has the partial product with sticky LSB.
// PP_HH can have a maximum of 0x3FFF_FFFF_FFFF_FFFF or thereabouts
// PP_HH can have a minimum of 0x1000_0000_0000_0000 or so
// The exponent of PP_HH is in EXP1, which is non-positive (0 or negative)
// That's the exponent that happens after the normalization
//
// EXP0 has the exponent that, when added to the normalized value, is out of range.
//
// Strategy:
//
// * Shift down bits, with sticky bit, such that the bits are aligned according
// to the LZ count and appropriate exponent, but not all the way to mantissa
// field, keep around the last few bits.
// * Put a 1 near the MSB
// * Check the LSBs for inexact; if inexact also set underflow
// * Convert [u]d2df -- will correctly round according to rounding mode
// * Replace exponent field with zero
{
BTMPL = #0 // offset for extract
BTMPH = sub(#FUDGE,BTMPH) // amount to right shift
}
{
p3 = cmp.gt(PP_HH_H,#-1) // is it positive?
BTMPH = min(BTMPH,TMP) // Don't shift more than 63
PP_HH = ATMP
}
{
TMP = USR
PP_LL = extractu(PP_HH,BTMP)
}
{
PP_HH = asr(PP_HH,BTMPH)
BTMPL = #0x0030 // underflow flag
AH = insert(S_ZERO,#EXPBITS,#HI_MANTBITS)
}
{
p0 = cmp.gtu(ONE,PP_LL) // Did we extract all zeros?
if (!p0.new) PP_HH_L = or(PP_HH_L,S_ONE) // add sticky bit
PP_HH_H = setbit(PP_HH_H,#HI_MANTBITS+3) // Add back in a bit so we can use convert instruction
}
{
PP_LL = neg(PP_HH)
p1 = bitsclr(PP_HH_L,#0x7) // Are the LSB's clear?
if (!p1.new) TMP = or(BTMPL,TMP) // If not, Inexact+Underflow
}
{
if (!p3) PP_HH = PP_LL
USR = TMP
}
{
A = convert_d2df(PP_HH) // Do rounding
p0 = dfcmp.eq(A,A) // realize exception
}
{
AH = insert(S_ZERO,#EXPBITS-1,#HI_MANTBITS+1) // Insert correct exponent
jumpr r31
}
.falign
.Lmul_ovf:
// We get either max finite value or infinity. Either way, overflow+inexact
{
TMP = USR
ATMP = combine(##0x7fefffff,#-1) // positive max finite
A = PP_HH
}
{
PP_LL_L = extractu(TMP,#2,#SR_ROUND_OFF) // rounding bits
TMP = or(TMP,#0x28) // inexact + overflow
BTMP = combine(##0x7ff00000,#0) // positive infinity
}
{
USR = TMP
PP_LL_L ^= lsr(AH,#31) // Does sign match rounding?
TMP = PP_LL_L // unmodified rounding mode
}
{
p0 = !cmp.eq(TMP,#1) // If not round-to-zero and
p0 = !cmp.eq(PP_LL_L,#2) // Not rounding the other way,
if (p0.new) ATMP = BTMP // we should get infinity
p0 = dfcmp.eq(A,A) // Realize FP exception if enabled
}
{
A = insert(ATMP,#63,#0) // insert inf/maxfinite, leave sign
jumpr r31
}
.Lmul_abnormal:
{
ATMP = extractu(A,#63,#0) // strip off sign
BTMP = extractu(B,#63,#0) // strip off sign
}
{
p3 = cmp.gtu(ATMP,BTMP)
if (!p3.new) A = B // sort values
if (!p3.new) B = A // sort values
}
{
// Any NaN --> NaN, possibly raise invalid if sNaN
p0 = dfclass(A,#0x0f) // A not NaN?
if (!p0.new) jump:nt .Linvalid_nan
if (!p3) ATMP = BTMP
if (!p3) BTMP = ATMP
}
{
// Infinity * nonzero number is infinity
p1 = dfclass(A,#0x08) // A is infinity
p1 = dfclass(B,#0x0e) // B is nonzero
}
{
// Infinity * zero --> NaN, raise invalid
// Other zeros return zero
p0 = dfclass(A,#0x08) // A is infinity
p0 = dfclass(B,#0x01) // B is zero
}
{
if (p1) jump .Ltrue_inf
p2 = dfclass(B,#0x01)
}
{
if (p0) jump .Linvalid_zeroinf
if (p2) jump .Ltrue_zero // so return zero
TMP = ##0x7c000000
}
// We are left with a normal or subnormal times a subnormal. A > B
// If A and B are both very small (exp(a) < BIAS-MANTBITS),
// we go to a single sticky bit, which we can round easily.
// If A and B might multiply to something bigger, decrease A exponent and increase
// B exponent and try again
{
p0 = bitsclr(AH,TMP)
if (p0.new) jump:nt .Lmul_tiny
}
{
TMP = cl0(BTMP)
}
{
TMP = add(TMP,#-EXPBITS)
}
{
BTMP = asl(BTMP,TMP)
}
{
B = insert(BTMP,#63,#0)
AH -= asl(TMP,#HI_MANTBITS)
}
jump __hexagon_muldf3
.Lmul_tiny:
{
TMP = USR
A = xor(A,B) // get sign bit
}
{
TMP = or(TMP,#0x30) // Inexact + Underflow
A = insert(ONE,#63,#0) // put in rounded up value
BTMPH = extractu(TMP,#2,#SR_ROUND_OFF) // get rounding mode
}
{
USR = TMP
p0 = cmp.gt(BTMPH,#1) // Round towards pos/neg inf?
if (!p0.new) AL = #0 // If not, zero
BTMPH ^= lsr(AH,#31) // rounding my way --> set LSB
}
{
p0 = cmp.eq(BTMPH,#3) // if rounding towards right inf
if (!p0.new) AL = #0 // don't go to zero
jumpr r31
}
.Linvalid_zeroinf:
{
TMP = USR
}
{
A = #-1
TMP = or(TMP,#2)
}
{
USR = TMP
}
{
p0 = dfcmp.uo(A,A) // force exception if enabled
jumpr r31
}
.Linvalid_nan:
{
p0 = dfclass(B,#0x0f) // if B is not NaN
TMP = convert_df2sf(A) // will generate invalid if sNaN
if (p0.new) B = A // make it whatever A is
}
{
BL = convert_df2sf(B) // will generate invalid if sNaN
A = #-1
jumpr r31
}
.falign
.Ltrue_zero:
{
A = B
B = A
}
.Ltrue_inf:
{
BH = extract(BH,#1,#31)
}
{
AH ^= asl(BH,#31)
jumpr r31
}
END(__hexagon_muldf3)
#undef ATMP
#undef ATMPL
#undef ATMPH
#undef BTMP
#undef BTMPL
#undef BTMPH