This is an archive of the discontinued LLVM Phabricator instance.

Differential D2799

[COMPILER-RT] Implement __multf3
ClosedPublic

Authored by koviankevin on Feb 13 2014, 9:48 PM.

Details

Reviewers
joerg
Commits
rG61d5d38e807b: Provide mul for IEEE quad. From GuanHong Liu.
rCRT211313: Provide mul for IEEE quad. From GuanHong Liu.
rL211313: Provide mul for IEEE quad. From GuanHong Liu.
Summary

Implement __multf3, which performs long double multiplication

Diff Detail

Event Timeline

koviankevin updated this revision to Unknown Object (????).Feb 19 2014, 6:12 PM

add test

koviankevin updated this revision to Unknown Object (????).Feb 24 2014, 7:50 PM

Use CRT_HAS_128BIT and modify the codes according to gribozavr's comments

koviankevin updated this revision to Unknown Object (????).Feb 27 2014, 8:43 PM

introduce CRT_LDBL_128BIT

koviankevin updated this revision to Unknown Object (????).Mar 10 2014, 1:12 AM

modify by gribozavr's comment and use COMPILER_RT_ABI

joerg added a subscriber: joerg.May 15 2014, 7:58 AM

Similar to D2798, please refactor first as one commit. Please make sure the changed "underflow to zero" case is covered by a test case.

koviankevin updated this revision to Diff 9734.May 23 2014, 1:06 AM

use UINT64_C and drop le64
muldf3.c and mulsf3.c are refactored in D3886

joerg accepted this revision.May 23 2014, 4:55 AM
joerg added a reviewer: joerg.
This revision is now accepted and ready to land.May 23 2014, 4:55 AM
joerg closed this revision.Jun 19 2014, 1:42 PM
joerg updated this revision to Diff 10657.

Closed by commit rL211313 (authored by @joerg).

Revision Contents

PathSize
lib/
132 lines
20 lines
lib/
compiler-rt/
trunk/
lib/
106 lines
96 lines

Diff 7108

lib/fp_mul.h

  • This file was added.
//===---- lib/fp_mul.h - floating point multiplication ------------*- C -*-===//
//
// The LLVM Compiler Infrastructure
//
// This file is dual licensed under the MIT and the University of Illinois Open
// Source Licenses. See LICENSE.TXT for details.
//
//===----------------------------------------------------------------------===//
//
// This file implements soft-float multiplication with the IEEE-754 default
// rounding (to nearest, ties to even).
//
//===----------------------------------------------------------------------===//
#ifndef FP_MUL_HEADER
#define FP_MUL_HEADER
#include "fp_lib.h"
#define __MULFP(precision) \
static inline fp_t __mul##precision##f3__(fp_t a, fp_t b) { \
const unsigned int aExponent = toRep(a) >> significandBits & maxExponent; \
const unsigned int bExponent = toRep(b) >> significandBits & maxExponent; \
const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit; \
\
rep_t aSignificand = toRep(a) & significandMask; \
rep_t bSignificand = toRep(b) & significandMask; \
int scale = 0; \
\
/* Detect if a or b is zero, denormal, infinity, or NaN.*/ \
if (aExponent-1U >= maxExponent-1U || bExponent-1U >= maxExponent-1U) { \
\
const rep_t aAbs = toRep(a) & absMask; \
const rep_t bAbs = toRep(b) & absMask; \
\
/* NaN * anything = qNaN*/ \
if (aAbs > infRep) return fromRep(toRep(a) | quietBit); \
/* anything * NaN = qNaN*/ \
if (bAbs > infRep) return fromRep(toRep(b) | quietBit); \
\
if (aAbs == infRep) { \
/* infinity * non-zero = +/- infinity*/ \
if (bAbs) return fromRep(aAbs | productSign); \
/* infinity * zero = NaN */ \
else return fromRep(qnanRep); \
} \
\
if (bAbs == infRep) { \
/*? non-zero * infinity = +/- infinity*/ \
if (aAbs) return fromRep(bAbs | productSign); \
/* zero * infinity = NaN */ \
else return fromRep(qnanRep); \
} \
\
/* zero * anything = +/- zero*/ \
if (!aAbs) return fromRep(productSign); \
/* anything * zero = +/- zero*/ \
if (!bAbs) return fromRep(productSign); \
\
/* one or both of a or b is denormal, the other (if applicable) is a
normal number. Renormalize one or both of a and b, and set scale to
include the necessary exponent adjustment. */ \
if (aAbs < implicitBit) scale += normalize(&aSignificand); \
if (bAbs < implicitBit) scale += normalize(&bSignificand); \
} \
\
/* Or in the implicit significand bit. (If we fell through from the
denormal path it was already set by normalize( ), but setting it twice
won't hurt anything.) */ \
aSignificand |= implicitBit; \
bSignificand |= implicitBit; \
\
/* Get the significand of a*b. Before multiplying the significands, shift
one of them left to left-align it in the field. Thus, the product will
have (exponentBits + 2) integral digits, all but two of which must be
zero. Normalizing this result is just a conditional left-shift by one
and bumping the exponent accordingly.*/ \
rep_t productHi, productLo; \
wideMultiply(aSignificand, bSignificand << exponentBits, \
&productHi, &productLo); \
\
int productExponent = aExponent + bExponent - exponentBias + scale; \
\
/* Normalize the significand, adjust exponent if needed.*/ \
if (productHi & implicitBit) productExponent++; \
else wideLeftShift(&productHi, &productLo, 1); \
\
/* If we have overflowed the type, return +/- infinity.*/ \
if (productExponent >= maxExponent) return fromRep(infRep | productSign); \
\
if (productExponent <= 0) { \
/* Result is denormal before rounding*/ \
\
/* If the result is so small that it just underflows to zero, return
a zero of the appropriate sign. Mathematically there is no need to
handle this case separately, but we make it a special case to
simplify the shift logic. */ \
const unsigned int shift = REP_C(1) - (unsigned int)productExponent; \
if (shift >= typeWidth) return fromRep(productSign); \
\
/* Otherwise, shift the significand of the result so that the round
bit is the high bit of productLo.*/ \
wideRightShiftWithSticky(&productHi, &productLo, shift); \
} \
else { \
/* Result is normal before rounding; insert the exponent.*/ \
productHi &= significandMask; \
productHi |= (rep_t)productExponent << significandBits; \
} \
\
/* Insert the sign of the result:*/ \
productHi |= productSign; \
\
/* Final rounding. The final result may overflow to infinity, or underflow
to zero, but those are the correct results in those cases. We use the
default IEEE-754 round-to-nearest, ties-to-even rounding mode.*/ \
if (productLo > signBit) productHi++; \
if (productLo == signBit) productHi += productHi & 1; \
return fromRep(productHi); \
}
#if defined SINGLE_PRECISION
COMPILER_RT_ABI __MULFP(s);
#elif defined DOUBLE_PRECISION
COMPILER_RT_ABI __MULFP(d);
#elif defined QUAD_PRECISION
COMPILER_RT_ABI __MULFP(t);
#else
#error SINGLE_PRECISION, DOUBLE_PRECISION or QUAD_PRECISION must be defined.
#endif
#endif // FP_MUL_HEADER

lib/muldf3.c

Context not available.
//===----------------------------------------------------------------------===// //===----------------------------------------------------------------------===//
#define DOUBLE_PRECISION #define DOUBLE_PRECISION
#include "fp_lib.h" #include "fp_mul.h"
ARM_EABI_FNALIAS(dmul, muldf3) ARM_EABI_FNALIAS(dmul, muldf3)
COMPILER_RT_ABI fp_t fp_t __muldf3(fp_t a, fp_t b) {
__muldf3(fp_t a, fp_t b) { return __muldf3__(a, b);
const unsigned int aExponent = toRep(a) >> significandBits & maxExponent;
const unsigned int bExponent = toRep(b) >> significandBits & maxExponent;
const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit;
rep_t aSignificand = toRep(a) & significandMask;
rep_t bSignificand = toRep(b) & significandMask;
int scale = 0;
// Detect if a or b is zero, denormal, infinity, or NaN.
if (aExponent-1U >= maxExponent-1U || bExponent-1U >= maxExponent-1U) {
const rep_t aAbs = toRep(a) & absMask;
const rep_t bAbs = toRep(b) & absMask;
// NaN * anything = qNaN
if (aAbs > infRep) return fromRep(toRep(a) | quietBit);
// anything * NaN = qNaN
if (bAbs > infRep) return fromRep(toRep(b) | quietBit);
if (aAbs == infRep) {
// infinity * non-zero = +/- infinity
if (bAbs) return fromRep(aAbs | productSign);
// infinity * zero = NaN
else return fromRep(qnanRep);
}
if (bAbs == infRep) {
// non-zero * infinity = +/- infinity
if (aAbs) return fromRep(bAbs | productSign);
// zero * infinity = NaN
else return fromRep(qnanRep);
}
// zero * anything = +/- zero
if (!aAbs) return fromRep(productSign);
// anything * zero = +/- zero
if (!bAbs) return fromRep(productSign);
// one or both of a or b is denormal, the other (if applicable) is a
// normal number. Renormalize one or both of a and b, and set scale to
// include the necessary exponent adjustment.
if (aAbs < implicitBit) scale += normalize(&aSignificand);
if (bAbs < implicitBit) scale += normalize(&bSignificand);
}
// Or in the implicit significand bit. (If we fell through from the
// denormal path it was already set by normalize( ), but setting it twice
// won't hurt anything.)
aSignificand |= implicitBit;
bSignificand |= implicitBit;
// Get the significand of a*b. Before multiplying the significands, shift
// one of them left to left-align it in the field. Thus, the product will
// have (exponentBits + 2) integral digits, all but two of which must be
// zero. Normalizing this result is just a conditional left-shift by one
// and bumping the exponent accordingly.
rep_t productHi, productLo;
wideMultiply(aSignificand, bSignificand << exponentBits,
&productHi, &productLo);
int productExponent = aExponent + bExponent - exponentBias + scale;
// Normalize the significand, adjust exponent if needed.
if (productHi & implicitBit) productExponent++;
else wideLeftShift(&productHi, &productLo, 1);
// If we have overflowed the type, return +/- infinity.
if (productExponent >= maxExponent) return fromRep(infRep | productSign);
if (productExponent <= 0) {
// Result is denormal before rounding
//
// If the result is so small that it just underflows to zero, return
// a zero of the appropriate sign. Mathematically there is no need to
// handle this case separately, but we make it a special case to
// simplify the shift logic.
const unsigned int shift = 1U - (unsigned int)productExponent;
if (shift >= typeWidth) return fromRep(productSign);
// Otherwise, shift the significand of the result so that the round
// bit is the high bit of productLo.
wideRightShiftWithSticky(&productHi, &productLo, shift);
}
else {
// Result is normal before rounding; insert the exponent.
productHi &= significandMask;
productHi |= (rep_t)productExponent << significandBits;
}
// Insert the sign of the result:
productHi |= productSign;
// Final rounding. The final result may overflow to infinity, or underflow
// to zero, but those are the correct results in those cases. We use the
// default IEEE-754 round-to-nearest, ties-to-even rounding mode.
if (productLo > signBit) productHi++;
if (productLo == signBit) productHi += productHi & 1;
return fromRep(productHi);
} }
Context not available.

lib/mulsf3.c

Context not available.
//===----------------------------------------------------------------------===// //===----------------------------------------------------------------------===//
#define SINGLE_PRECISION #define SINGLE_PRECISION
#include "fp_lib.h" #include "fp_mul.h"
ARM_EABI_FNALIAS(fmul, mulsf3) ARM_EABI_FNALIAS(fmul, mulsf3)
COMPILER_RT_ABI fp_t fp_t __mulsf3(fp_t a, fp_t b) {
__mulsf3(fp_t a, fp_t b) { return __mulsf3__(a, b);
const unsigned int aExponent = toRep(a) >> significandBits & maxExponent;
const unsigned int bExponent = toRep(b) >> significandBits & maxExponent;
const rep_t productSign = (toRep(a) ^ toRep(b)) & signBit;
rep_t aSignificand = toRep(a) & significandMask;
rep_t bSignificand = toRep(b) & significandMask;
int scale = 0;
// Detect if a or b is zero, denormal, infinity, or NaN.
if (aExponent-1U >= maxExponent-1U || bExponent-1U >= maxExponent-1U) {
const rep_t aAbs = toRep(a) & absMask;
const rep_t bAbs = toRep(b) & absMask;
// NaN * anything = qNaN
if (aAbs > infRep) return fromRep(toRep(a) | quietBit);
// anything * NaN = qNaN
if (bAbs > infRep) return fromRep(toRep(b) | quietBit);
if (aAbs == infRep) {
// infinity * non-zero = +/- infinity
if (bAbs) return fromRep(aAbs | productSign);
// infinity * zero = NaN
else return fromRep(qnanRep);
}
if (bAbs == infRep) {
// non-zero * infinity = +/- infinity
if (aAbs) return fromRep(bAbs | productSign);
// zero * infinity = NaN
else return fromRep(qnanRep);
}
// zero * anything = +/- zero
if (!aAbs) return fromRep(productSign);
// anything * zero = +/- zero
if (!bAbs) return fromRep(productSign);
// one or both of a or b is denormal, the other (if applicable) is a
// normal number. Renormalize one or both of a and b, and set scale to
// include the necessary exponent adjustment.
if (aAbs < implicitBit) scale += normalize(&aSignificand);
if (bAbs < implicitBit) scale += normalize(&bSignificand);
}
// Or in the implicit significand bit. (If we fell through from the
// denormal path it was already set by normalize( ), but setting it twice
// won't hurt anything.)
aSignificand |= implicitBit;
bSignificand |= implicitBit;
// Get the significand of a*b. Before multiplying the significands, shift
// one of them left to left-align it in the field. Thus, the product will
// have (exponentBits + 2) integral digits, all but two of which must be
// zero. Normalizing this result is just a conditional left-shift by one
// and bumping the exponent accordingly.
rep_t productHi, productLo;
wideMultiply(aSignificand, bSignificand << exponentBits,
&productHi, &productLo);
int productExponent = aExponent + bExponent - exponentBias + scale;
// Normalize the significand, adjust exponent if needed.
if (productHi & implicitBit) productExponent++;
else wideLeftShift(&productHi, &productLo, 1);
// If we have overflowed the type, return +/- infinity.
if (productExponent >= maxExponent) return fromRep(infRep | productSign);
if (productExponent <= 0) {
// Result is denormal before rounding, the exponent is zero and we
// need to shift the significand.
wideRightShiftWithSticky(&productHi, &productLo, 1U - (unsigned)productExponent);
}
else {
// Result is normal before rounding; insert the exponent.
productHi &= significandMask;
productHi |= (rep_t)productExponent << significandBits;
}
// Insert the sign of the result:
productHi |= productSign;
// Final rounding. The final result may overflow to infinity, or underflow
// to zero, but those are the correct results in those cases.
if (productLo > signBit) productHi++;
if (productLo == signBit) productHi += productHi & 1;
return fromRep(productHi);
} }
Context not available.

lib/multf3.c

  • This file was added.
//===-- lib/multf3.c - Quad-precision multiplication --------------*- C -*-===//
//
// The LLVM Compiler Infrastructure
//
// This file is dual licensed under the MIT and the University of Illinois Open
// Source Licenses. See LICENSE.TXT for details.
//
//===----------------------------------------------------------------------===//
//
// This file implements quad-precision soft-float multiplication
// with the IEEE-754 default rounding (to nearest, ties to even).
//
//===----------------------------------------------------------------------===//
#define QUAD_PRECISION
#include "fp_mul.h"
fp_t __multf3(fp_t a, fp_t b) {
return __multf3__(a, b);
}