00001
00051 #include <stdio.h>
00052 #include <stdlib.h>
00053 #include <string.h>
00054 #include <math.h>
00055 #include "util.h"
00056 #include "epd.h"
00057
00058
00070 EpDouble *
00071 EpdAlloc(void)
00072 {
00073 EpDouble *epd;
00074
00075 epd = ALLOC(EpDouble, 1);
00076 return(epd);
00077 }
00078
00079
00091 int
00092 EpdCmp(const char *key1, const char *key2)
00093 {
00094 EpDouble *epd1 = (EpDouble *) key1;
00095 EpDouble *epd2 = (EpDouble *) key2;
00096 if (epd1->type.value != epd2->type.value ||
00097 epd1->exponent != epd2->exponent) {
00098 return(1);
00099 }
00100 return(0);
00101 }
00102
00103
00115 void
00116 EpdFree(EpDouble *epd)
00117 {
00118 FREE(epd);
00119 }
00120
00121
00133 void
00134 EpdGetString(EpDouble *epd, char *str)
00135 {
00136 double value;
00137 int exponent;
00138 char *pos;
00139
00140 if (IsNanDouble(epd->type.value)) {
00141 sprintf(str, "NaN");
00142 return;
00143 } else if (IsInfDouble(epd->type.value)) {
00144 if (epd->type.bits.sign == 1)
00145 sprintf(str, "-Inf");
00146 else
00147 sprintf(str, "Inf");
00148 return;
00149 }
00150
00151 assert(epd->type.bits.exponent == EPD_MAX_BIN ||
00152 epd->type.bits.exponent == 0);
00153
00154 EpdGetValueAndDecimalExponent(epd, &value, &exponent);
00155 sprintf(str, "%e", value);
00156 pos = strstr(str, "e");
00157 if (exponent >= 0) {
00158 if (exponent < 10)
00159 sprintf(pos + 1, "+0%d", exponent);
00160 else
00161 sprintf(pos + 1, "+%d", exponent);
00162 } else {
00163 exponent *= -1;
00164 if (exponent < 10)
00165 sprintf(pos + 1, "-0%d", exponent);
00166 else
00167 sprintf(pos + 1, "-%d", exponent);
00168 }
00169 }
00170
00171
00183 void
00184 EpdConvert(double value, EpDouble *epd)
00185 {
00186 epd->type.value = value;
00187 epd->exponent = 0;
00188 EpdNormalize(epd);
00189 }
00190
00191
00203 void
00204 EpdMultiply(EpDouble *epd1, double value)
00205 {
00206 EpDouble epd2;
00207 double tmp;
00208 int exponent;
00209
00210 if (EpdIsNan(epd1) || IsNanDouble(value)) {
00211 EpdMakeNan(epd1);
00212 return;
00213 } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00214 int sign;
00215
00216 EpdConvert(value, &epd2);
00217 sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00218 EpdMakeInf(epd1, sign);
00219 return;
00220 }
00221
00222 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00223
00224 EpdConvert(value, &epd2);
00225 tmp = epd1->type.value * epd2.type.value;
00226 exponent = epd1->exponent + epd2.exponent;
00227 epd1->type.value = tmp;
00228 epd1->exponent = exponent;
00229 EpdNormalize(epd1);
00230 }
00231
00232
00244 void
00245 EpdMultiply2(EpDouble *epd1, EpDouble *epd2)
00246 {
00247 double value;
00248 int exponent;
00249
00250 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00251 EpdMakeNan(epd1);
00252 return;
00253 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00254 int sign;
00255
00256 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00257 EpdMakeInf(epd1, sign);
00258 return;
00259 }
00260
00261 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00262 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00263
00264 value = epd1->type.value * epd2->type.value;
00265 exponent = epd1->exponent + epd2->exponent;
00266 epd1->type.value = value;
00267 epd1->exponent = exponent;
00268 EpdNormalize(epd1);
00269 }
00270
00271
00283 void
00284 EpdMultiply2Decimal(EpDouble *epd1, EpDouble *epd2)
00285 {
00286 double value;
00287 int exponent;
00288
00289 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00290 EpdMakeNan(epd1);
00291 return;
00292 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00293 int sign;
00294
00295 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00296 EpdMakeInf(epd1, sign);
00297 return;
00298 }
00299
00300 value = epd1->type.value * epd2->type.value;
00301 exponent = epd1->exponent + epd2->exponent;
00302 epd1->type.value = value;
00303 epd1->exponent = exponent;
00304 EpdNormalizeDecimal(epd1);
00305 }
00306
00307
00319 void
00320 EpdMultiply3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
00321 {
00322 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00323 EpdMakeNan(epd1);
00324 return;
00325 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00326 int sign;
00327
00328 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00329 EpdMakeInf(epd3, sign);
00330 return;
00331 }
00332
00333 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00334 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00335
00336 epd3->type.value = epd1->type.value * epd2->type.value;
00337 epd3->exponent = epd1->exponent + epd2->exponent;
00338 EpdNormalize(epd3);
00339 }
00340
00341
00353 void
00354 EpdMultiply3Decimal(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
00355 {
00356 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00357 EpdMakeNan(epd1);
00358 return;
00359 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00360 int sign;
00361
00362 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00363 EpdMakeInf(epd3, sign);
00364 return;
00365 }
00366
00367 epd3->type.value = epd1->type.value * epd2->type.value;
00368 epd3->exponent = epd1->exponent + epd2->exponent;
00369 EpdNormalizeDecimal(epd3);
00370 }
00371
00372
00384 void
00385 EpdDivide(EpDouble *epd1, double value)
00386 {
00387 EpDouble epd2;
00388 double tmp;
00389 int exponent;
00390
00391 if (EpdIsNan(epd1) || IsNanDouble(value)) {
00392 EpdMakeNan(epd1);
00393 return;
00394 } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00395 int sign;
00396
00397 EpdConvert(value, &epd2);
00398 if (EpdIsInf(epd1) && IsInfDouble(value)) {
00399 EpdMakeNan(epd1);
00400 } else if (EpdIsInf(epd1)) {
00401 sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00402 EpdMakeInf(epd1, sign);
00403 } else {
00404 sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00405 EpdMakeZero(epd1, sign);
00406 }
00407 return;
00408 }
00409
00410 if (value == 0.0) {
00411 EpdMakeNan(epd1);
00412 return;
00413 }
00414
00415 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00416
00417 EpdConvert(value, &epd2);
00418 tmp = epd1->type.value / epd2.type.value;
00419 exponent = epd1->exponent - epd2.exponent;
00420 epd1->type.value = tmp;
00421 epd1->exponent = exponent;
00422 EpdNormalize(epd1);
00423 }
00424
00425
00437 void
00438 EpdDivide2(EpDouble *epd1, EpDouble *epd2)
00439 {
00440 double value;
00441 int exponent;
00442
00443 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00444 EpdMakeNan(epd1);
00445 return;
00446 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00447 int sign;
00448
00449 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00450 EpdMakeNan(epd1);
00451 } else if (EpdIsInf(epd1)) {
00452 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00453 EpdMakeInf(epd1, sign);
00454 } else {
00455 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00456 EpdMakeZero(epd1, sign);
00457 }
00458 return;
00459 }
00460
00461 if (epd2->type.value == 0.0) {
00462 EpdMakeNan(epd1);
00463 return;
00464 }
00465
00466 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00467 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00468
00469 value = epd1->type.value / epd2->type.value;
00470 exponent = epd1->exponent - epd2->exponent;
00471 epd1->type.value = value;
00472 epd1->exponent = exponent;
00473 EpdNormalize(epd1);
00474 }
00475
00476
00488 void
00489 EpdDivide3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
00490 {
00491 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00492 EpdMakeNan(epd3);
00493 return;
00494 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00495 int sign;
00496
00497 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00498 EpdMakeNan(epd3);
00499 } else if (EpdIsInf(epd1)) {
00500 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00501 EpdMakeInf(epd3, sign);
00502 } else {
00503 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00504 EpdMakeZero(epd3, sign);
00505 }
00506 return;
00507 }
00508
00509 if (epd2->type.value == 0.0) {
00510 EpdMakeNan(epd3);
00511 return;
00512 }
00513
00514 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00515 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00516
00517 epd3->type.value = epd1->type.value / epd2->type.value;
00518 epd3->exponent = epd1->exponent - epd2->exponent;
00519 EpdNormalize(epd3);
00520 }
00521
00522
00534 void
00535 EpdAdd(EpDouble *epd1, double value)
00536 {
00537 EpDouble epd2;
00538 double tmp;
00539 int exponent, diff;
00540
00541 if (EpdIsNan(epd1) || IsNanDouble(value)) {
00542 EpdMakeNan(epd1);
00543 return;
00544 } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00545 int sign;
00546
00547 EpdConvert(value, &epd2);
00548 if (EpdIsInf(epd1) && IsInfDouble(value)) {
00549 sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00550 if (sign == 1)
00551 EpdMakeNan(epd1);
00552 } else if (EpdIsInf(&epd2)) {
00553 EpdCopy(&epd2, epd1);
00554 }
00555 return;
00556 }
00557
00558 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00559
00560 EpdConvert(value, &epd2);
00561 if (epd1->exponent > epd2.exponent) {
00562 diff = epd1->exponent - epd2.exponent;
00563 if (diff <= EPD_MAX_BIN)
00564 tmp = epd1->type.value + epd2.type.value / pow((double)2.0, (double)diff);
00565 else
00566 tmp = epd1->type.value;
00567 exponent = epd1->exponent;
00568 } else if (epd1->exponent < epd2.exponent) {
00569 diff = epd2.exponent - epd1->exponent;
00570 if (diff <= EPD_MAX_BIN)
00571 tmp = epd1->type.value / pow((double)2.0, (double)diff) + epd2.type.value;
00572 else
00573 tmp = epd2.type.value;
00574 exponent = epd2.exponent;
00575 } else {
00576 tmp = epd1->type.value + epd2.type.value;
00577 exponent = epd1->exponent;
00578 }
00579 epd1->type.value = tmp;
00580 epd1->exponent = exponent;
00581 EpdNormalize(epd1);
00582 }
00583
00584
00596 void
00597 EpdAdd2(EpDouble *epd1, EpDouble *epd2)
00598 {
00599 double value;
00600 int exponent, diff;
00601
00602 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00603 EpdMakeNan(epd1);
00604 return;
00605 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00606 int sign;
00607
00608 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00609 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00610 if (sign == 1)
00611 EpdMakeNan(epd1);
00612 } else if (EpdIsInf(epd2)) {
00613 EpdCopy(epd2, epd1);
00614 }
00615 return;
00616 }
00617
00618 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00619 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00620
00621 if (epd1->exponent > epd2->exponent) {
00622 diff = epd1->exponent - epd2->exponent;
00623 if (diff <= EPD_MAX_BIN) {
00624 value = epd1->type.value +
00625 epd2->type.value / pow((double)2.0, (double)diff);
00626 } else
00627 value = epd1->type.value;
00628 exponent = epd1->exponent;
00629 } else if (epd1->exponent < epd2->exponent) {
00630 diff = epd2->exponent - epd1->exponent;
00631 if (diff <= EPD_MAX_BIN) {
00632 value = epd1->type.value / pow((double)2.0, (double)diff) +
00633 epd2->type.value;
00634 } else
00635 value = epd2->type.value;
00636 exponent = epd2->exponent;
00637 } else {
00638 value = epd1->type.value + epd2->type.value;
00639 exponent = epd1->exponent;
00640 }
00641 epd1->type.value = value;
00642 epd1->exponent = exponent;
00643 EpdNormalize(epd1);
00644 }
00645
00646
00658 void
00659 EpdAdd3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
00660 {
00661 double value;
00662 int exponent, diff;
00663
00664 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00665 EpdMakeNan(epd3);
00666 return;
00667 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00668 int sign;
00669
00670 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00671 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00672 if (sign == 1)
00673 EpdMakeNan(epd3);
00674 else
00675 EpdCopy(epd1, epd3);
00676 } else if (EpdIsInf(epd1)) {
00677 EpdCopy(epd1, epd3);
00678 } else {
00679 EpdCopy(epd2, epd3);
00680 }
00681 return;
00682 }
00683
00684 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00685 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00686
00687 if (epd1->exponent > epd2->exponent) {
00688 diff = epd1->exponent - epd2->exponent;
00689 if (diff <= EPD_MAX_BIN) {
00690 value = epd1->type.value +
00691 epd2->type.value / pow((double)2.0, (double)diff);
00692 } else
00693 value = epd1->type.value;
00694 exponent = epd1->exponent;
00695 } else if (epd1->exponent < epd2->exponent) {
00696 diff = epd2->exponent - epd1->exponent;
00697 if (diff <= EPD_MAX_BIN) {
00698 value = epd1->type.value / pow((double)2.0, (double)diff) +
00699 epd2->type.value;
00700 } else
00701 value = epd2->type.value;
00702 exponent = epd2->exponent;
00703 } else {
00704 value = epd1->type.value + epd2->type.value;
00705 exponent = epd1->exponent;
00706 }
00707 epd3->type.value = value;
00708 epd3->exponent = exponent;
00709 EpdNormalize(epd3);
00710 }
00711
00712
00724 void
00725 EpdSubtract(EpDouble *epd1, double value)
00726 {
00727 EpDouble epd2;
00728 double tmp;
00729 int exponent, diff;
00730
00731 if (EpdIsNan(epd1) || IsNanDouble(value)) {
00732 EpdMakeNan(epd1);
00733 return;
00734 } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00735 int sign;
00736
00737 EpdConvert(value, &epd2);
00738 if (EpdIsInf(epd1) && IsInfDouble(value)) {
00739 sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00740 if (sign == 0)
00741 EpdMakeNan(epd1);
00742 } else if (EpdIsInf(&epd2)) {
00743 EpdCopy(&epd2, epd1);
00744 }
00745 return;
00746 }
00747
00748 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00749
00750 EpdConvert(value, &epd2);
00751 if (epd1->exponent > epd2.exponent) {
00752 diff = epd1->exponent - epd2.exponent;
00753 if (diff <= EPD_MAX_BIN)
00754 tmp = epd1->type.value - epd2.type.value / pow((double)2.0, (double)diff);
00755 else
00756 tmp = epd1->type.value;
00757 exponent = epd1->exponent;
00758 } else if (epd1->exponent < epd2.exponent) {
00759 diff = epd2.exponent - epd1->exponent;
00760 if (diff <= EPD_MAX_BIN)
00761 tmp = epd1->type.value / pow((double)2.0, (double)diff) - epd2.type.value;
00762 else
00763 tmp = epd2.type.value * (double)(-1.0);
00764 exponent = epd2.exponent;
00765 } else {
00766 tmp = epd1->type.value - epd2.type.value;
00767 exponent = epd1->exponent;
00768 }
00769 epd1->type.value = tmp;
00770 epd1->exponent = exponent;
00771 EpdNormalize(epd1);
00772 }
00773
00774
00786 void
00787 EpdSubtract2(EpDouble *epd1, EpDouble *epd2)
00788 {
00789 double value;
00790 int exponent, diff;
00791
00792 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00793 EpdMakeNan(epd1);
00794 return;
00795 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00796 int sign;
00797
00798 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00799 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00800 if (sign == 0)
00801 EpdMakeNan(epd1);
00802 } else if (EpdIsInf(epd2)) {
00803 EpdCopy(epd2, epd1);
00804 }
00805 return;
00806 }
00807
00808 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00809 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00810
00811 if (epd1->exponent > epd2->exponent) {
00812 diff = epd1->exponent - epd2->exponent;
00813 if (diff <= EPD_MAX_BIN) {
00814 value = epd1->type.value -
00815 epd2->type.value / pow((double)2.0, (double)diff);
00816 } else
00817 value = epd1->type.value;
00818 exponent = epd1->exponent;
00819 } else if (epd1->exponent < epd2->exponent) {
00820 diff = epd2->exponent - epd1->exponent;
00821 if (diff <= EPD_MAX_BIN) {
00822 value = epd1->type.value / pow((double)2.0, (double)diff) -
00823 epd2->type.value;
00824 } else
00825 value = epd2->type.value * (double)(-1.0);
00826 exponent = epd2->exponent;
00827 } else {
00828 value = epd1->type.value - epd2->type.value;
00829 exponent = epd1->exponent;
00830 }
00831 epd1->type.value = value;
00832 epd1->exponent = exponent;
00833 EpdNormalize(epd1);
00834 }
00835
00836
00848 void
00849 EpdSubtract3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
00850 {
00851 double value;
00852 int exponent, diff;
00853
00854 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00855 EpdMakeNan(epd3);
00856 return;
00857 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00858 int sign;
00859
00860 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00861 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00862 if (sign == 0)
00863 EpdCopy(epd1, epd3);
00864 else
00865 EpdMakeNan(epd3);
00866 } else if (EpdIsInf(epd1)) {
00867 EpdCopy(epd1, epd1);
00868 } else {
00869 sign = epd2->type.bits.sign ^ 0x1;
00870 EpdMakeInf(epd3, sign);
00871 }
00872 return;
00873 }
00874
00875 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00876 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00877
00878 if (epd1->exponent > epd2->exponent) {
00879 diff = epd1->exponent - epd2->exponent;
00880 if (diff <= EPD_MAX_BIN) {
00881 value = epd1->type.value -
00882 epd2->type.value / pow((double)2.0, (double)diff);
00883 } else
00884 value = epd1->type.value;
00885 exponent = epd1->exponent;
00886 } else if (epd1->exponent < epd2->exponent) {
00887 diff = epd2->exponent - epd1->exponent;
00888 if (diff <= EPD_MAX_BIN) {
00889 value = epd1->type.value / pow((double)2.0, (double)diff) -
00890 epd2->type.value;
00891 } else
00892 value = epd2->type.value * (double)(-1.0);
00893 exponent = epd2->exponent;
00894 } else {
00895 value = epd1->type.value - epd2->type.value;
00896 exponent = epd1->exponent;
00897 }
00898 epd3->type.value = value;
00899 epd3->exponent = exponent;
00900 EpdNormalize(epd3);
00901 }
00902
00903
00915 void
00916 EpdPow2(int n, EpDouble *epd)
00917 {
00918 if (n <= EPD_MAX_BIN) {
00919 EpdConvert(pow((double)2.0, (double)n), epd);
00920 } else {
00921 EpDouble epd1, epd2;
00922 int n1, n2;
00923
00924 n1 = n / 2;
00925 n2 = n - n1;
00926 EpdPow2(n1, &epd1);
00927 EpdPow2(n2, &epd2);
00928 EpdMultiply3(&epd1, &epd2, epd);
00929 }
00930 }
00931
00932
00944 void
00945 EpdPow2Decimal(int n, EpDouble *epd)
00946 {
00947 if (n <= EPD_MAX_BIN) {
00948 epd->type.value = pow((double)2.0, (double)n);
00949 epd->exponent = 0;
00950 EpdNormalizeDecimal(epd);
00951 } else {
00952 EpDouble epd1, epd2;
00953 int n1, n2;
00954
00955 n1 = n / 2;
00956 n2 = n - n1;
00957 EpdPow2Decimal(n1, &epd1);
00958 EpdPow2Decimal(n2, &epd2);
00959 EpdMultiply3Decimal(&epd1, &epd2, epd);
00960 }
00961 }
00962
00963
00975 void
00976 EpdNormalize(EpDouble *epd)
00977 {
00978 int exponent;
00979
00980 if (IsNanOrInfDouble(epd->type.value)) {
00981 epd->exponent = 0;
00982 return;
00983 }
00984
00985 exponent = EpdGetExponent(epd->type.value);
00986 if (exponent == EPD_MAX_BIN)
00987 return;
00988 exponent -= EPD_MAX_BIN;
00989 epd->type.bits.exponent = EPD_MAX_BIN;
00990 epd->exponent += exponent;
00991 }
00992
00993
01005 void
01006 EpdNormalizeDecimal(EpDouble *epd)
01007 {
01008 int exponent;
01009
01010 if (IsNanOrInfDouble(epd->type.value)) {
01011 epd->exponent = 0;
01012 return;
01013 }
01014
01015 exponent = EpdGetExponentDecimal(epd->type.value);
01016 epd->type.value /= pow((double)10.0, (double)exponent);
01017 epd->exponent += exponent;
01018 }
01019
01020
01032 void
01033 EpdGetValueAndDecimalExponent(EpDouble *epd, double *value, int *exponent)
01034 {
01035 EpDouble epd1, epd2;
01036
01037 if (EpdIsNanOrInf(epd))
01038 return;
01039
01040 if (EpdIsZero(epd)) {
01041 *value = 0.0;
01042 *exponent = 0;
01043 return;
01044 }
01045
01046 epd1.type.value = epd->type.value;
01047 epd1.exponent = 0;
01048 EpdPow2Decimal(epd->exponent, &epd2);
01049 EpdMultiply2Decimal(&epd1, &epd2);
01050
01051 *value = epd1.type.value;
01052 *exponent = epd1.exponent;
01053 }
01054
01066 int
01067 EpdGetExponent(double value)
01068 {
01069 int exponent;
01070 EpDouble epd;
01071
01072 epd.type.value = value;
01073 exponent = epd.type.bits.exponent;
01074 return(exponent);
01075 }
01076
01077
01089 int
01090 EpdGetExponentDecimal(double value)
01091 {
01092 char *pos, str[24];
01093 int exponent;
01094
01095 sprintf(str, "%E", value);
01096 pos = strstr(str, "E");
01097 sscanf(pos, "E%d", &exponent);
01098 return(exponent);
01099 }
01100
01101
01113 void
01114 EpdMakeInf(EpDouble *epd, int sign)
01115 {
01116 epd->type.bits.mantissa1 = 0;
01117 epd->type.bits.mantissa0 = 0;
01118 epd->type.bits.exponent = EPD_EXP_INF;
01119 epd->type.bits.sign = sign;
01120 epd->exponent = 0;
01121 }
01122
01123
01135 void
01136 EpdMakeZero(EpDouble *epd, int sign)
01137 {
01138 epd->type.bits.mantissa1 = 0;
01139 epd->type.bits.mantissa0 = 0;
01140 epd->type.bits.exponent = 0;
01141 epd->type.bits.sign = sign;
01142 epd->exponent = 0;
01143 }
01144
01145
01157 void
01158 EpdMakeNan(EpDouble *epd)
01159 {
01160 epd->type.nan.mantissa1 = 0;
01161 epd->type.nan.mantissa0 = 0;
01162 epd->type.nan.quiet_bit = 1;
01163 epd->type.nan.exponent = EPD_EXP_INF;
01164 epd->type.nan.sign = 1;
01165 epd->exponent = 0;
01166 }
01167
01168
01180 void
01181 EpdCopy(EpDouble *from, EpDouble *to)
01182 {
01183 to->type.value = from->type.value;
01184 to->exponent = from->exponent;
01185 }
01186
01187
01199 int
01200 EpdIsInf(EpDouble *epd)
01201 {
01202 return(IsInfDouble(epd->type.value));
01203 }
01204
01205
01217 int
01218 EpdIsZero(EpDouble *epd)
01219 {
01220 if (epd->type.value == 0.0)
01221 return(1);
01222 else
01223 return(0);
01224 }
01225
01226
01238 int
01239 EpdIsNan(EpDouble *epd)
01240 {
01241 return(IsNanDouble(epd->type.value));
01242 }
01243
01244
01256 int
01257 EpdIsNanOrInf(EpDouble *epd)
01258 {
01259 return(IsNanOrInfDouble(epd->type.value));
01260 }
01261
01262
01274 int
01275 IsInfDouble(double value)
01276 {
01277 EpType val;
01278
01279 val.value = value;
01280 if (val.bits.exponent == EPD_EXP_INF &&
01281 val.bits.mantissa0 == 0 &&
01282 val.bits.mantissa1 == 0) {
01283 if (val.bits.sign == 0)
01284 return(1);
01285 else
01286 return(-1);
01287 }
01288 return(0);
01289 }
01290
01291
01303 int
01304 IsNanDouble(double value)
01305 {
01306 EpType val;
01307
01308 val.value = value;
01309 if (val.nan.exponent == EPD_EXP_INF &&
01310 val.nan.sign == 1 &&
01311 val.nan.quiet_bit == 1 &&
01312 val.nan.mantissa0 == 0 &&
01313 val.nan.mantissa1 == 0) {
01314 return(1);
01315 }
01316 return(0);
01317 }
01318
01319
01331 int
01332 IsNanOrInfDouble(double value)
01333 {
01334 EpType val;
01335
01336 val.value = value;
01337 if (val.nan.exponent == EPD_EXP_INF &&
01338 val.nan.mantissa0 == 0 &&
01339 val.nan.mantissa1 == 0 &&
01340 (val.nan.sign == 1 || val.nan.quiet_bit == 0)) {
01341 return(1);
01342 }
01343 return(0);
01344 }
