1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054 |
- /**
- @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);
- }
|