You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1053 lines
22 KiB
1053 lines
22 KiB
/**
|
|
@file
|
|
|
|
@ingroup epd
|
|
|
|
@brief Arithmetic functions with extended double precision.
|
|
|
|
@author In-Ho Moon
|
|
|
|
@copyright@parblock
|
|
Copyright (c) 1995-2015, Regents of the University of Colorado
|
|
|
|
All rights reserved.
|
|
|
|
Redistribution and use in source and binary forms, with or without
|
|
modification, are permitted provided that the following conditions
|
|
are met:
|
|
|
|
Redistributions of source code must retain the above copyright
|
|
notice, this list of conditions and the following disclaimer.
|
|
|
|
Redistributions in binary form must reproduce the above copyright
|
|
notice, this list of conditions and the following disclaimer in the
|
|
documentation and/or other materials provided with the distribution.
|
|
|
|
Neither the name of the University of Colorado nor the names of its
|
|
contributors may be used to endorse or promote products derived from
|
|
this software without specific prior written permission.
|
|
|
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
|
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
|
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
|
|
FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
|
|
COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
|
|
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
|
|
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
|
|
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
|
|
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
|
|
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
|
|
ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
|
POSSIBILITY OF SUCH DAMAGE.
|
|
@endparblock
|
|
|
|
*/
|
|
|
|
#include <math.h>
|
|
#include "util.h"
|
|
#include "epdInt.h"
|
|
|
|
|
|
/**
|
|
@brief Allocates an EpDouble struct.
|
|
*/
|
|
EpDouble *
|
|
EpdAlloc(void)
|
|
{
|
|
EpDouble *epd;
|
|
|
|
epd = ALLOC(EpDouble, 1);
|
|
return(epd);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Compares two EpDouble struct.
|
|
|
|
@return 0 if the two structures hold the same value; 1 otherwise.
|
|
*/
|
|
int
|
|
EpdCmp(const void *key1, const void *key2)
|
|
{
|
|
EpDouble const *epd1 = (EpDouble const *) key1;
|
|
EpDouble const *epd2 = (EpDouble const *) key2;
|
|
if (epd1->type.value != epd2->type.value ||
|
|
epd1->exponent != epd2->exponent) {
|
|
return(1);
|
|
}
|
|
return(0);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Frees an EpDouble struct.
|
|
*/
|
|
void
|
|
EpdFree(EpDouble *epd)
|
|
{
|
|
FREE(epd);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Converts an extended precision double value to a string.
|
|
|
|
@sideeffect The string is written at the address passed in `str`.
|
|
*/
|
|
void
|
|
EpdGetString(EpDouble const *epd, char *str)
|
|
{
|
|
double value;
|
|
int exponent;
|
|
char *pos;
|
|
|
|
if (!str) return;
|
|
|
|
if (IsNanDouble(epd->type.value)) {
|
|
sprintf(str, "NaN");
|
|
return;
|
|
} else if (IsInfDouble(epd->type.value)) {
|
|
if (epd->type.bits.sign == 1)
|
|
sprintf(str, "-inf");
|
|
else
|
|
sprintf(str, "inf");
|
|
return;
|
|
}
|
|
|
|
assert(epd->type.bits.exponent == EPD_MAX_BIN ||
|
|
epd->type.bits.exponent == 0);
|
|
|
|
EpdGetValueAndDecimalExponent(epd, &value, &exponent);
|
|
sprintf(str, "%e", value);
|
|
pos = strstr(str, "e");
|
|
if (exponent >= 0) {
|
|
if (exponent < 10)
|
|
sprintf(pos + 1, "+0%d", exponent);
|
|
else
|
|
sprintf(pos + 1, "+%d", exponent);
|
|
} else {
|
|
exponent *= -1;
|
|
if (exponent < 10)
|
|
sprintf(pos + 1, "-0%d", exponent);
|
|
else
|
|
sprintf(pos + 1, "-%d", exponent);
|
|
}
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Converts double to EpDouble struct.
|
|
*/
|
|
void
|
|
EpdConvert(double value, EpDouble *epd)
|
|
{
|
|
epd->type.value = value;
|
|
epd->exponent = 0;
|
|
EpdNormalize(epd);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Multiplies an extended precision double by a double.
|
|
*/
|
|
void
|
|
EpdMultiply(EpDouble *epd1, double value)
|
|
{
|
|
EpDouble epd2;
|
|
double tmp;
|
|
int exponent;
|
|
|
|
if (EpdIsNan(epd1) || IsNanDouble(value)) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || IsInfDouble(value)) {
|
|
int sign;
|
|
|
|
EpdConvert(value, &epd2);
|
|
sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
|
|
EpdMakeInf(epd1, sign);
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
EpdConvert(value, &epd2);
|
|
tmp = epd1->type.value * epd2.type.value;
|
|
exponent = epd1->exponent + epd2.exponent;
|
|
epd1->type.value = tmp;
|
|
epd1->exponent = exponent;
|
|
EpdNormalize(epd1);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Multiplies an extended precision double by another.
|
|
*/
|
|
void
|
|
EpdMultiply2(EpDouble *epd1, EpDouble const *epd2)
|
|
{
|
|
double value;
|
|
int exponent;
|
|
|
|
if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
|
|
int sign;
|
|
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
EpdMakeInf(epd1, sign);
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
assert(epd2->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
value = epd1->type.value * epd2->type.value;
|
|
exponent = epd1->exponent + epd2->exponent;
|
|
epd1->type.value = value;
|
|
epd1->exponent = exponent;
|
|
EpdNormalize(epd1);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Multiplies two extended precision double values.
|
|
*/
|
|
void
|
|
EpdMultiply2Decimal(EpDouble *epd1, EpDouble const *epd2)
|
|
{
|
|
double value;
|
|
int exponent;
|
|
|
|
if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
|
|
int sign;
|
|
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
EpdMakeInf(epd1, sign);
|
|
return;
|
|
}
|
|
|
|
value = epd1->type.value * epd2->type.value;
|
|
exponent = epd1->exponent + epd2->exponent;
|
|
epd1->type.value = value;
|
|
epd1->exponent = exponent;
|
|
EpdNormalizeDecimal(epd1);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Multiplies two extended precision double values.
|
|
|
|
@details The result goes in the third operand.
|
|
*/
|
|
void
|
|
EpdMultiply3(EpDouble const *epd1, EpDouble const *epd2, EpDouble *epd3)
|
|
{
|
|
if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
|
|
EpdMakeNan(epd3);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
|
|
int sign;
|
|
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
EpdMakeInf(epd3, sign);
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
assert(epd2->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
epd3->type.value = epd1->type.value * epd2->type.value;
|
|
epd3->exponent = epd1->exponent + epd2->exponent;
|
|
EpdNormalize(epd3);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Multiplies two extended precision double values.
|
|
*/
|
|
void
|
|
EpdMultiply3Decimal(EpDouble const *epd1, EpDouble const *epd2, EpDouble *epd3)
|
|
{
|
|
if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
|
|
EpdMakeNan(epd3);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
|
|
int sign;
|
|
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
EpdMakeInf(epd3, sign);
|
|
return;
|
|
}
|
|
|
|
epd3->type.value = epd1->type.value * epd2->type.value;
|
|
epd3->exponent = epd1->exponent + epd2->exponent;
|
|
EpdNormalizeDecimal(epd3);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Divides an extended precision double by a double.
|
|
*/
|
|
void
|
|
EpdDivide(EpDouble *epd1, double value)
|
|
{
|
|
EpDouble epd2;
|
|
double tmp;
|
|
int exponent;
|
|
|
|
if (EpdIsNan(epd1) || IsNanDouble(value)) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || IsInfDouble(value)) {
|
|
int sign;
|
|
|
|
EpdConvert(value, &epd2);
|
|
if (EpdIsInf(epd1) && IsInfDouble(value)) {
|
|
EpdMakeNan(epd1);
|
|
} else if (EpdIsInf(epd1)) {
|
|
sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
|
|
EpdMakeInf(epd1, sign);
|
|
} else {
|
|
sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
|
|
EpdMakeZero(epd1, sign);
|
|
}
|
|
return;
|
|
}
|
|
|
|
if (value == 0.0) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
EpdConvert(value, &epd2);
|
|
tmp = epd1->type.value / epd2.type.value;
|
|
exponent = epd1->exponent - epd2.exponent;
|
|
epd1->type.value = tmp;
|
|
epd1->exponent = exponent;
|
|
EpdNormalize(epd1);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Divides an extended precision double by another.
|
|
*/
|
|
void
|
|
EpdDivide2(EpDouble *epd1, EpDouble const *epd2)
|
|
{
|
|
double value;
|
|
int exponent;
|
|
|
|
if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
|
|
int sign;
|
|
|
|
if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
|
|
EpdMakeNan(epd1);
|
|
} else if (EpdIsInf(epd1)) {
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
EpdMakeInf(epd1, sign);
|
|
} else {
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
EpdMakeZero(epd1, sign);
|
|
}
|
|
return;
|
|
}
|
|
|
|
if (epd2->type.value == 0.0) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
assert(epd2->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
value = epd1->type.value / epd2->type.value;
|
|
exponent = epd1->exponent - epd2->exponent;
|
|
epd1->type.value = value;
|
|
epd1->exponent = exponent;
|
|
EpdNormalize(epd1);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Divides two extended precision double values.
|
|
*/
|
|
void
|
|
EpdDivide3(EpDouble const *epd1, EpDouble const *epd2, EpDouble *epd3)
|
|
{
|
|
if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
|
|
EpdMakeNan(epd3);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
|
|
int sign;
|
|
|
|
if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
|
|
EpdMakeNan(epd3);
|
|
} else if (EpdIsInf(epd1)) {
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
EpdMakeInf(epd3, sign);
|
|
} else {
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
EpdMakeZero(epd3, sign);
|
|
}
|
|
return;
|
|
}
|
|
|
|
if (epd2->type.value == 0.0) {
|
|
EpdMakeNan(epd3);
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
assert(epd2->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
epd3->type.value = epd1->type.value / epd2->type.value;
|
|
epd3->exponent = epd1->exponent - epd2->exponent;
|
|
EpdNormalize(epd3);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Adds a double to an extended precision double.
|
|
*/
|
|
void
|
|
EpdAdd(EpDouble *epd1, double value)
|
|
{
|
|
EpDouble epd2;
|
|
double tmp;
|
|
int exponent, diff;
|
|
|
|
if (EpdIsNan(epd1) || IsNanDouble(value)) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || IsInfDouble(value)) {
|
|
int sign;
|
|
|
|
EpdConvert(value, &epd2);
|
|
if (EpdIsInf(epd1) && IsInfDouble(value)) {
|
|
sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
|
|
if (sign == 1)
|
|
EpdMakeNan(epd1);
|
|
} else if (EpdIsInf(&epd2)) {
|
|
EpdCopy(&epd2, epd1);
|
|
}
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
EpdConvert(value, &epd2);
|
|
if (epd1->exponent > epd2.exponent) {
|
|
diff = epd1->exponent - epd2.exponent;
|
|
if (diff <= EPD_MAX_BIN)
|
|
tmp = epd1->type.value + epd2.type.value / pow((double)2.0, (double)diff);
|
|
else
|
|
tmp = epd1->type.value;
|
|
exponent = epd1->exponent;
|
|
} else if (epd1->exponent < epd2.exponent) {
|
|
diff = epd2.exponent - epd1->exponent;
|
|
if (diff <= EPD_MAX_BIN)
|
|
tmp = epd1->type.value / pow((double)2.0, (double)diff) + epd2.type.value;
|
|
else
|
|
tmp = epd2.type.value;
|
|
exponent = epd2.exponent;
|
|
} else {
|
|
tmp = epd1->type.value + epd2.type.value;
|
|
exponent = epd1->exponent;
|
|
}
|
|
epd1->type.value = tmp;
|
|
epd1->exponent = exponent;
|
|
EpdNormalize(epd1);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Adds an extended precision double to another.
|
|
|
|
@details The sum goes in the first argument.
|
|
*/
|
|
void
|
|
EpdAdd2(EpDouble *epd1, EpDouble const *epd2)
|
|
{
|
|
double value;
|
|
int exponent, diff;
|
|
|
|
if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
|
|
int sign;
|
|
|
|
if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
if (sign == 1)
|
|
EpdMakeNan(epd1);
|
|
} else if (EpdIsInf(epd2)) {
|
|
EpdCopy(epd2, epd1);
|
|
}
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
assert(epd2->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
if (epd1->exponent > epd2->exponent) {
|
|
diff = epd1->exponent - epd2->exponent;
|
|
if (diff <= EPD_MAX_BIN) {
|
|
value = epd1->type.value +
|
|
epd2->type.value / pow((double)2.0, (double)diff);
|
|
} else
|
|
value = epd1->type.value;
|
|
exponent = epd1->exponent;
|
|
} else if (epd1->exponent < epd2->exponent) {
|
|
diff = epd2->exponent - epd1->exponent;
|
|
if (diff <= EPD_MAX_BIN) {
|
|
value = epd1->type.value / pow((double)2.0, (double)diff) +
|
|
epd2->type.value;
|
|
} else
|
|
value = epd2->type.value;
|
|
exponent = epd2->exponent;
|
|
} else {
|
|
value = epd1->type.value + epd2->type.value;
|
|
exponent = epd1->exponent;
|
|
}
|
|
epd1->type.value = value;
|
|
epd1->exponent = exponent;
|
|
EpdNormalize(epd1);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Adds two extended precision double values.
|
|
*/
|
|
void
|
|
EpdAdd3(EpDouble const *epd1, EpDouble const *epd2, EpDouble *epd3)
|
|
{
|
|
double value;
|
|
int exponent, diff;
|
|
|
|
if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
|
|
EpdMakeNan(epd3);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
|
|
int sign;
|
|
|
|
if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
if (sign == 1)
|
|
EpdMakeNan(epd3);
|
|
else
|
|
EpdCopy(epd1, epd3);
|
|
} else if (EpdIsInf(epd1)) {
|
|
EpdCopy(epd1, epd3);
|
|
} else {
|
|
EpdCopy(epd2, epd3);
|
|
}
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
assert(epd2->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
if (epd1->exponent > epd2->exponent) {
|
|
diff = epd1->exponent - epd2->exponent;
|
|
if (diff <= EPD_MAX_BIN) {
|
|
value = epd1->type.value +
|
|
epd2->type.value / pow((double)2.0, (double)diff);
|
|
} else
|
|
value = epd1->type.value;
|
|
exponent = epd1->exponent;
|
|
} else if (epd1->exponent < epd2->exponent) {
|
|
diff = epd2->exponent - epd1->exponent;
|
|
if (diff <= EPD_MAX_BIN) {
|
|
value = epd1->type.value / pow((double)2.0, (double)diff) +
|
|
epd2->type.value;
|
|
} else
|
|
value = epd2->type.value;
|
|
exponent = epd2->exponent;
|
|
} else {
|
|
value = epd1->type.value + epd2->type.value;
|
|
exponent = epd1->exponent;
|
|
}
|
|
epd3->type.value = value;
|
|
epd3->exponent = exponent;
|
|
EpdNormalize(epd3);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Subtracts a double from an extended precision double.
|
|
*/
|
|
void
|
|
EpdSubtract(EpDouble *epd1, double value)
|
|
{
|
|
EpDouble epd2;
|
|
double tmp;
|
|
int exponent, diff;
|
|
|
|
if (EpdIsNan(epd1) || IsNanDouble(value)) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || IsInfDouble(value)) {
|
|
int sign;
|
|
|
|
EpdConvert(value, &epd2);
|
|
if (EpdIsInf(epd1) && IsInfDouble(value)) {
|
|
sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
|
|
if (sign == 0)
|
|
EpdMakeNan(epd1);
|
|
} else if (EpdIsInf(&epd2)) {
|
|
EpdCopy(&epd2, epd1);
|
|
}
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
EpdConvert(value, &epd2);
|
|
if (epd1->exponent > epd2.exponent) {
|
|
diff = epd1->exponent - epd2.exponent;
|
|
if (diff <= EPD_MAX_BIN)
|
|
tmp = epd1->type.value - epd2.type.value / pow((double)2.0, (double)diff);
|
|
else
|
|
tmp = epd1->type.value;
|
|
exponent = epd1->exponent;
|
|
} else if (epd1->exponent < epd2.exponent) {
|
|
diff = epd2.exponent - epd1->exponent;
|
|
if (diff <= EPD_MAX_BIN)
|
|
tmp = epd1->type.value / pow((double)2.0, (double)diff) - epd2.type.value;
|
|
else
|
|
tmp = epd2.type.value * (double)(-1.0);
|
|
exponent = epd2.exponent;
|
|
} else {
|
|
tmp = epd1->type.value - epd2.type.value;
|
|
exponent = epd1->exponent;
|
|
}
|
|
epd1->type.value = tmp;
|
|
epd1->exponent = exponent;
|
|
EpdNormalize(epd1);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Subtracts an extended precision double from another.
|
|
*/
|
|
void
|
|
EpdSubtract2(EpDouble *epd1, EpDouble const *epd2)
|
|
{
|
|
double value;
|
|
int exponent, diff;
|
|
|
|
if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
|
|
EpdMakeNan(epd1);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
|
|
int sign;
|
|
|
|
if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
if (sign == 0)
|
|
EpdMakeNan(epd1);
|
|
} else if (EpdIsInf(epd2)) {
|
|
EpdCopy(epd2, epd1);
|
|
}
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
assert(epd2->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
if (epd1->exponent > epd2->exponent) {
|
|
diff = epd1->exponent - epd2->exponent;
|
|
if (diff <= EPD_MAX_BIN) {
|
|
value = epd1->type.value -
|
|
epd2->type.value / pow((double)2.0, (double)diff);
|
|
} else
|
|
value = epd1->type.value;
|
|
exponent = epd1->exponent;
|
|
} else if (epd1->exponent < epd2->exponent) {
|
|
diff = epd2->exponent - epd1->exponent;
|
|
if (diff <= EPD_MAX_BIN) {
|
|
value = epd1->type.value / pow((double)2.0, (double)diff) -
|
|
epd2->type.value;
|
|
} else
|
|
value = epd2->type.value * (double)(-1.0);
|
|
exponent = epd2->exponent;
|
|
} else {
|
|
value = epd1->type.value - epd2->type.value;
|
|
exponent = epd1->exponent;
|
|
}
|
|
epd1->type.value = value;
|
|
epd1->exponent = exponent;
|
|
EpdNormalize(epd1);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Subtracts two extended precision double values.
|
|
*/
|
|
void
|
|
EpdSubtract3(EpDouble const *epd1, EpDouble const *epd2, EpDouble *epd3)
|
|
{
|
|
double value;
|
|
int exponent, diff;
|
|
|
|
if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
|
|
EpdMakeNan(epd3);
|
|
return;
|
|
} else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
|
|
int sign;
|
|
|
|
if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
|
|
sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
|
|
if (sign == 0)
|
|
EpdCopy(epd1, epd3);
|
|
else
|
|
EpdMakeNan(epd3);
|
|
} else if (EpdIsInf(epd1)) {
|
|
EpdCopy(epd1, epd3);
|
|
} else {
|
|
sign = epd2->type.bits.sign ^ 0x1;
|
|
EpdMakeInf(epd3, sign);
|
|
}
|
|
return;
|
|
}
|
|
|
|
assert(epd1->type.bits.exponent == EPD_MAX_BIN);
|
|
assert(epd2->type.bits.exponent == EPD_MAX_BIN);
|
|
|
|
if (epd1->exponent > epd2->exponent) {
|
|
diff = epd1->exponent - epd2->exponent;
|
|
if (diff <= EPD_MAX_BIN) {
|
|
value = epd1->type.value -
|
|
epd2->type.value / pow((double)2.0, (double)diff);
|
|
} else
|
|
value = epd1->type.value;
|
|
exponent = epd1->exponent;
|
|
} else if (epd1->exponent < epd2->exponent) {
|
|
diff = epd2->exponent - epd1->exponent;
|
|
if (diff <= EPD_MAX_BIN) {
|
|
value = epd1->type.value / pow((double)2.0, (double)diff) -
|
|
epd2->type.value;
|
|
} else
|
|
value = epd2->type.value * (double)(-1.0);
|
|
exponent = epd2->exponent;
|
|
} else {
|
|
value = epd1->type.value - epd2->type.value;
|
|
exponent = epd1->exponent;
|
|
}
|
|
epd3->type.value = value;
|
|
epd3->exponent = exponent;
|
|
EpdNormalize(epd3);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Computes extended precision pow of base 2.
|
|
*/
|
|
void
|
|
EpdPow2(int n, EpDouble *epd)
|
|
{
|
|
if (n <= EPD_MAX_BIN) {
|
|
EpdConvert(pow((double)2.0, (double)n), epd);
|
|
} else {
|
|
EpDouble epd1, epd2;
|
|
int n1, n2;
|
|
|
|
n1 = n / 2;
|
|
n2 = n - n1;
|
|
EpdPow2(n1, &epd1);
|
|
EpdPow2(n2, &epd2);
|
|
EpdMultiply3(&epd1, &epd2, epd);
|
|
}
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Computes extended precision pow of base 2.
|
|
*/
|
|
void
|
|
EpdPow2Decimal(int n, EpDouble *epd)
|
|
{
|
|
if (n <= EPD_MAX_BIN) {
|
|
epd->type.value = pow((double)2.0, (double)n);
|
|
epd->exponent = 0;
|
|
EpdNormalizeDecimal(epd);
|
|
} else {
|
|
EpDouble epd1, epd2;
|
|
int n1, n2;
|
|
|
|
n1 = n / 2;
|
|
n2 = n - n1;
|
|
EpdPow2Decimal(n1, &epd1);
|
|
EpdPow2Decimal(n2, &epd2);
|
|
EpdMultiply3Decimal(&epd1, &epd2, epd);
|
|
}
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Normalize an extended precision double value.
|
|
*/
|
|
void
|
|
EpdNormalize(EpDouble *epd)
|
|
{
|
|
int exponent;
|
|
|
|
if (IsNanOrInfDouble(epd->type.value)) {
|
|
epd->exponent = 0;
|
|
return;
|
|
}
|
|
|
|
exponent = EpdGetExponent(epd->type.value);
|
|
if (exponent == EPD_MAX_BIN)
|
|
return;
|
|
exponent -= EPD_MAX_BIN;
|
|
epd->type.bits.exponent = EPD_MAX_BIN;
|
|
epd->exponent += exponent;
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Normalize an extended precision double value.
|
|
*/
|
|
void
|
|
EpdNormalizeDecimal(EpDouble *epd)
|
|
{
|
|
int exponent;
|
|
|
|
if (IsNanOrInfDouble(epd->type.value)) {
|
|
epd->exponent = 0;
|
|
return;
|
|
}
|
|
|
|
exponent = EpdGetExponentDecimal(epd->type.value);
|
|
epd->type.value /= pow((double)10.0, (double)exponent);
|
|
epd->exponent += exponent;
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Returns value and decimal exponent of EpDouble.
|
|
*/
|
|
void
|
|
EpdGetValueAndDecimalExponent(EpDouble const *epd, double *value, int *exponent)
|
|
{
|
|
EpDouble epd1, epd2;
|
|
|
|
if (EpdIsNanOrInf(epd)) {
|
|
*exponent = EPD_EXP_INF;
|
|
*value = 0.0;
|
|
return;
|
|
}
|
|
|
|
if (EpdIsZero(epd)) {
|
|
*value = 0.0;
|
|
*exponent = 0;
|
|
return;
|
|
}
|
|
|
|
epd1.type.value = epd->type.value;
|
|
epd1.exponent = 0;
|
|
EpdPow2Decimal(epd->exponent, &epd2);
|
|
EpdMultiply2Decimal(&epd1, &epd2);
|
|
|
|
*value = epd1.type.value;
|
|
*exponent = epd1.exponent;
|
|
}
|
|
|
|
/**
|
|
@brief Returns the exponent value of a double.
|
|
*/
|
|
int
|
|
EpdGetExponent(double value)
|
|
{
|
|
int exponent;
|
|
EpDouble epd;
|
|
|
|
epd.type.value = value;
|
|
exponent = epd.type.bits.exponent;
|
|
return(exponent);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Returns the decimal exponent value of a double.
|
|
*/
|
|
int
|
|
EpdGetExponentDecimal(double value)
|
|
{
|
|
char *pos, str[24];
|
|
int exponent;
|
|
|
|
sprintf(str, "%E", value);
|
|
pos = strstr(str, "E");
|
|
sscanf(pos, "E%d", &exponent);
|
|
return(exponent);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Makes EpDouble Inf.
|
|
*/
|
|
void
|
|
EpdMakeInf(EpDouble *epd, int sign)
|
|
{
|
|
epd->type.bits.mantissa1 = 0;
|
|
epd->type.bits.mantissa0 = 0;
|
|
epd->type.bits.exponent = EPD_EXP_INF;
|
|
epd->type.bits.sign = sign;
|
|
epd->exponent = 0;
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Makes EpDouble Zero.
|
|
*/
|
|
void
|
|
EpdMakeZero(EpDouble *epd, int sign)
|
|
{
|
|
epd->type.bits.mantissa1 = 0;
|
|
epd->type.bits.mantissa0 = 0;
|
|
epd->type.bits.exponent = 0;
|
|
epd->type.bits.sign = sign;
|
|
epd->exponent = 0;
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Makes EpDouble NaN.
|
|
*/
|
|
void
|
|
EpdMakeNan(EpDouble *epd)
|
|
{
|
|
epd->type.nan.mantissa1 = 0;
|
|
epd->type.nan.mantissa0 = 0;
|
|
epd->type.nan.quiet_bit = 1;
|
|
epd->type.nan.exponent = EPD_EXP_INF;
|
|
epd->type.nan.sign = 1;
|
|
epd->exponent = 0;
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Copies an EpDouble struct.
|
|
*/
|
|
void
|
|
EpdCopy(EpDouble const *from, EpDouble *to)
|
|
{
|
|
to->type.value = from->type.value;
|
|
to->exponent = from->exponent;
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Checks whether the value is Inf.
|
|
*/
|
|
int
|
|
EpdIsInf(EpDouble const *epd)
|
|
{
|
|
return(IsInfDouble(epd->type.value));
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Checks whether the value is Zero.
|
|
*/
|
|
int
|
|
EpdIsZero(EpDouble const *epd)
|
|
{
|
|
if (epd->type.value == 0.0)
|
|
return(1);
|
|
else
|
|
return(0);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Checks whether the value is NaN.
|
|
*/
|
|
int
|
|
EpdIsNan(EpDouble const *epd)
|
|
{
|
|
return(IsNanDouble(epd->type.value));
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Checks whether the value is NaN or Inf.
|
|
*/
|
|
int
|
|
EpdIsNanOrInf(EpDouble const *epd)
|
|
{
|
|
return(IsNanOrInfDouble(epd->type.value));
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Checks whether the value is Inf.
|
|
*/
|
|
int
|
|
IsInfDouble(double value)
|
|
{
|
|
EpType val;
|
|
|
|
val.value = value;
|
|
if (val.bits.exponent == EPD_EXP_INF &&
|
|
val.bits.mantissa0 == 0 &&
|
|
val.bits.mantissa1 == 0) {
|
|
if (val.bits.sign == 0)
|
|
return(1);
|
|
else
|
|
return(-1);
|
|
}
|
|
return(0);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Checks whether the value is NaN.
|
|
*/
|
|
int
|
|
IsNanDouble(double value)
|
|
{
|
|
EpType val;
|
|
|
|
val.value = value;
|
|
if (val.nan.exponent == EPD_EXP_INF &&
|
|
val.nan.sign == 1 &&
|
|
val.nan.quiet_bit == 1 &&
|
|
val.nan.mantissa0 == 0 &&
|
|
val.nan.mantissa1 == 0) {
|
|
return(1);
|
|
}
|
|
return(0);
|
|
}
|
|
|
|
|
|
/**
|
|
@brief Checks whether the value is NaN or Inf.
|
|
*/
|
|
int
|
|
IsNanOrInfDouble(double value)
|
|
{
|
|
EpType val;
|
|
|
|
val.value = value;
|
|
if (val.nan.exponent == EPD_EXP_INF &&
|
|
val.nan.mantissa0 == 0 &&
|
|
val.nan.mantissa1 == 0 &&
|
|
(val.nan.sign == 1 || val.nan.quiet_bit == 0)) {
|
|
return(1);
|
|
}
|
|
return(0);
|
|
}
|