00001
00024 #include <stdio.h>
00025 #include <stdlib.h>
00026 #include <string.h>
00027 #include <math.h>
00028 #include "util_hack.h"
00029 #include "epd.h"
00030
00031
00043 EpDouble *
00044 EpdAlloc()
00045 {
00046 EpDouble *epd;
00047
00048 epd = ALLOC(EpDouble, 1);
00049 return(epd);
00050 }
00051
00052
00064 int
00065 EpdCmp(const char *key1, const char *key2)
00066 {
00067 EpDouble *epd1 = (EpDouble *) key1;
00068 EpDouble *epd2 = (EpDouble *) key2;
00069 if (epd1->type.value != epd2->type.value ||
00070 epd1->exponent != epd2->exponent) {
00071 return(1);
00072 }
00073 return(0);
00074 }
00075
00076
00088 void
00089 EpdFree(EpDouble *epd)
00090 {
00091 FREE(epd);
00092 }
00093
00094
00106 void
00107 EpdGetString(EpDouble *epd, char *str)
00108 {
00109 double value;
00110 int exponent;
00111 char *pos;
00112
00113 if (IsNanDouble(epd->type.value)) {
00114 sprintf(str, "NaN");
00115 return;
00116 } else if (IsInfDouble(epd->type.value)) {
00117 if (epd->type.bits.sign == 1)
00118 sprintf(str, "-Inf");
00119 else
00120 sprintf(str, "Inf");
00121 return;
00122 }
00123
00124 assert(epd->type.bits.exponent == EPD_MAX_BIN ||
00125 epd->type.bits.exponent == 0);
00126
00127 EpdGetValueAndDecimalExponent(epd, &value, &exponent);
00128 sprintf(str, "%e", value);
00129 pos = strstr(str, "e");
00130 if (exponent >= 0) {
00131 if (exponent < 10)
00132 sprintf(pos + 1, "+0%d", exponent);
00133 else
00134 sprintf(pos + 1, "+%d", exponent);
00135 } else {
00136 exponent *= -1;
00137 if (exponent < 10)
00138 sprintf(pos + 1, "-0%d", exponent);
00139 else
00140 sprintf(pos + 1, "-%d", exponent);
00141 }
00142 }
00143
00144
00156 void
00157 EpdConvert(double value, EpDouble *epd)
00158 {
00159 epd->type.value = value;
00160 epd->exponent = 0;
00161 EpdNormalize(epd);
00162 }
00163
00164
00176 void
00177 EpdMultiply(EpDouble *epd1, double value)
00178 {
00179 EpDouble epd2;
00180 double tmp;
00181 int exponent;
00182
00183 if (EpdIsNan(epd1) || IsNanDouble(value)) {
00184 EpdMakeNan(epd1);
00185 return;
00186 } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00187 int sign;
00188
00189 EpdConvert(value, &epd2);
00190 sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00191 EpdMakeInf(epd1, sign);
00192 return;
00193 }
00194
00195 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00196
00197 EpdConvert(value, &epd2);
00198 tmp = epd1->type.value * epd2.type.value;
00199 exponent = epd1->exponent + epd2.exponent;
00200 epd1->type.value = tmp;
00201 epd1->exponent = exponent;
00202 EpdNormalize(epd1);
00203 }
00204
00205
00217 void
00218 EpdMultiply2(EpDouble *epd1, EpDouble *epd2)
00219 {
00220 double value;
00221 int exponent;
00222
00223 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00224 EpdMakeNan(epd1);
00225 return;
00226 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00227 int sign;
00228
00229 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00230 EpdMakeInf(epd1, sign);
00231 return;
00232 }
00233
00234 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00235 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00236
00237 value = epd1->type.value * epd2->type.value;
00238 exponent = epd1->exponent + epd2->exponent;
00239 epd1->type.value = value;
00240 epd1->exponent = exponent;
00241 EpdNormalize(epd1);
00242 }
00243
00244
00256 void
00257 EpdMultiply2Decimal(EpDouble *epd1, EpDouble *epd2)
00258 {
00259 double value;
00260 int exponent;
00261
00262 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00263 EpdMakeNan(epd1);
00264 return;
00265 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00266 int sign;
00267
00268 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00269 EpdMakeInf(epd1, sign);
00270 return;
00271 }
00272
00273 value = epd1->type.value * epd2->type.value;
00274 exponent = epd1->exponent + epd2->exponent;
00275 epd1->type.value = value;
00276 epd1->exponent = exponent;
00277 EpdNormalizeDecimal(epd1);
00278 }
00279
00280
00292 void
00293 EpdMultiply3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
00294 {
00295 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00296 EpdMakeNan(epd1);
00297 return;
00298 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00299 int sign;
00300
00301 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00302 EpdMakeInf(epd3, sign);
00303 return;
00304 }
00305
00306 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00307 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00308
00309 epd3->type.value = epd1->type.value * epd2->type.value;
00310 epd3->exponent = epd1->exponent + epd2->exponent;
00311 EpdNormalize(epd3);
00312 }
00313
00314
00326 void
00327 EpdMultiply3Decimal(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
00328 {
00329 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00330 EpdMakeNan(epd1);
00331 return;
00332 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00333 int sign;
00334
00335 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00336 EpdMakeInf(epd3, sign);
00337 return;
00338 }
00339
00340 epd3->type.value = epd1->type.value * epd2->type.value;
00341 epd3->exponent = epd1->exponent + epd2->exponent;
00342 EpdNormalizeDecimal(epd3);
00343 }
00344
00345
00357 void
00358 EpdDivide(EpDouble *epd1, double value)
00359 {
00360 EpDouble epd2;
00361 double tmp;
00362 int exponent;
00363
00364 if (EpdIsNan(epd1) || IsNanDouble(value)) {
00365 EpdMakeNan(epd1);
00366 return;
00367 } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00368 int sign;
00369
00370 EpdConvert(value, &epd2);
00371 if (EpdIsInf(epd1) && IsInfDouble(value)) {
00372 EpdMakeNan(epd1);
00373 } else if (EpdIsInf(epd1)) {
00374 sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00375 EpdMakeInf(epd1, sign);
00376 } else {
00377 sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00378 EpdMakeZero(epd1, sign);
00379 }
00380 return;
00381 }
00382
00383 if (value == 0.0) {
00384 EpdMakeNan(epd1);
00385 return;
00386 }
00387
00388 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00389
00390 EpdConvert(value, &epd2);
00391 tmp = epd1->type.value / epd2.type.value;
00392 exponent = epd1->exponent - epd2.exponent;
00393 epd1->type.value = tmp;
00394 epd1->exponent = exponent;
00395 EpdNormalize(epd1);
00396 }
00397
00398
00410 void
00411 EpdDivide2(EpDouble *epd1, EpDouble *epd2)
00412 {
00413 double value;
00414 int exponent;
00415
00416 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00417 EpdMakeNan(epd1);
00418 return;
00419 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00420 int sign;
00421
00422 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00423 EpdMakeNan(epd1);
00424 } else if (EpdIsInf(epd1)) {
00425 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00426 EpdMakeInf(epd1, sign);
00427 } else {
00428 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00429 EpdMakeZero(epd1, sign);
00430 }
00431 return;
00432 }
00433
00434 if (epd2->type.value == 0.0) {
00435 EpdMakeNan(epd1);
00436 return;
00437 }
00438
00439 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00440 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00441
00442 value = epd1->type.value / epd2->type.value;
00443 exponent = epd1->exponent - epd2->exponent;
00444 epd1->type.value = value;
00445 epd1->exponent = exponent;
00446 EpdNormalize(epd1);
00447 }
00448
00449
00461 void
00462 EpdDivide3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
00463 {
00464 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00465 EpdMakeNan(epd3);
00466 return;
00467 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00468 int sign;
00469
00470 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00471 EpdMakeNan(epd3);
00472 } else if (EpdIsInf(epd1)) {
00473 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00474 EpdMakeInf(epd3, sign);
00475 } else {
00476 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00477 EpdMakeZero(epd3, sign);
00478 }
00479 return;
00480 }
00481
00482 if (epd2->type.value == 0.0) {
00483 EpdMakeNan(epd3);
00484 return;
00485 }
00486
00487 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00488 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00489
00490 epd3->type.value = epd1->type.value / epd2->type.value;
00491 epd3->exponent = epd1->exponent - epd2->exponent;
00492 EpdNormalize(epd3);
00493 }
00494
00495
00507 void
00508 EpdAdd(EpDouble *epd1, double value)
00509 {
00510 EpDouble epd2;
00511 double tmp;
00512 int exponent, diff;
00513
00514 if (EpdIsNan(epd1) || IsNanDouble(value)) {
00515 EpdMakeNan(epd1);
00516 return;
00517 } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00518 int sign;
00519
00520 EpdConvert(value, &epd2);
00521 if (EpdIsInf(epd1) && IsInfDouble(value)) {
00522 sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00523 if (sign == 1)
00524 EpdMakeNan(epd1);
00525 } else if (EpdIsInf(&epd2)) {
00526 EpdCopy(&epd2, epd1);
00527 }
00528 return;
00529 }
00530
00531 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00532
00533 EpdConvert(value, &epd2);
00534 if (epd1->exponent > epd2.exponent) {
00535 diff = epd1->exponent - epd2.exponent;
00536 if (diff <= EPD_MAX_BIN)
00537 tmp = epd1->type.value + epd2.type.value / pow((double)2.0, (double)diff);
00538 else
00539 tmp = epd1->type.value;
00540 exponent = epd1->exponent;
00541 } else if (epd1->exponent < epd2.exponent) {
00542 diff = epd2.exponent - epd1->exponent;
00543 if (diff <= EPD_MAX_BIN)
00544 tmp = epd1->type.value / pow((double)2.0, (double)diff) + epd2.type.value;
00545 else
00546 tmp = epd2.type.value;
00547 exponent = epd2.exponent;
00548 } else {
00549 tmp = epd1->type.value + epd2.type.value;
00550 exponent = epd1->exponent;
00551 }
00552 epd1->type.value = tmp;
00553 epd1->exponent = exponent;
00554 EpdNormalize(epd1);
00555 }
00556
00557
00569 void
00570 EpdAdd2(EpDouble *epd1, EpDouble *epd2)
00571 {
00572 double value;
00573 int exponent, diff;
00574
00575 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00576 EpdMakeNan(epd1);
00577 return;
00578 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00579 int sign;
00580
00581 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00582 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00583 if (sign == 1)
00584 EpdMakeNan(epd1);
00585 } else if (EpdIsInf(epd2)) {
00586 EpdCopy(epd2, epd1);
00587 }
00588 return;
00589 }
00590
00591 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00592 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00593
00594 if (epd1->exponent > epd2->exponent) {
00595 diff = epd1->exponent - epd2->exponent;
00596 if (diff <= EPD_MAX_BIN) {
00597 value = epd1->type.value +
00598 epd2->type.value / pow((double)2.0, (double)diff);
00599 } else
00600 value = epd1->type.value;
00601 exponent = epd1->exponent;
00602 } else if (epd1->exponent < epd2->exponent) {
00603 diff = epd2->exponent - epd1->exponent;
00604 if (diff <= EPD_MAX_BIN) {
00605 value = epd1->type.value / pow((double)2.0, (double)diff) +
00606 epd2->type.value;
00607 } else
00608 value = epd2->type.value;
00609 exponent = epd2->exponent;
00610 } else {
00611 value = epd1->type.value + epd2->type.value;
00612 exponent = epd1->exponent;
00613 }
00614 epd1->type.value = value;
00615 epd1->exponent = exponent;
00616 EpdNormalize(epd1);
00617 }
00618
00619
00631 void
00632 EpdAdd3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
00633 {
00634 double value;
00635 int exponent, diff;
00636
00637 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00638 EpdMakeNan(epd3);
00639 return;
00640 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00641 int sign;
00642
00643 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00644 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00645 if (sign == 1)
00646 EpdMakeNan(epd3);
00647 else
00648 EpdCopy(epd1, epd3);
00649 } else if (EpdIsInf(epd1)) {
00650 EpdCopy(epd1, epd3);
00651 } else {
00652 EpdCopy(epd2, epd3);
00653 }
00654 return;
00655 }
00656
00657 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00658 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00659
00660 if (epd1->exponent > epd2->exponent) {
00661 diff = epd1->exponent - epd2->exponent;
00662 if (diff <= EPD_MAX_BIN) {
00663 value = epd1->type.value +
00664 epd2->type.value / pow((double)2.0, (double)diff);
00665 } else
00666 value = epd1->type.value;
00667 exponent = epd1->exponent;
00668 } else if (epd1->exponent < epd2->exponent) {
00669 diff = epd2->exponent - epd1->exponent;
00670 if (diff <= EPD_MAX_BIN) {
00671 value = epd1->type.value / pow((double)2.0, (double)diff) +
00672 epd2->type.value;
00673 } else
00674 value = epd2->type.value;
00675 exponent = epd2->exponent;
00676 } else {
00677 value = epd1->type.value + epd2->type.value;
00678 exponent = epd1->exponent;
00679 }
00680 epd3->type.value = value;
00681 epd3->exponent = exponent;
00682 EpdNormalize(epd3);
00683 }
00684
00685
00697 void
00698 EpdSubtract(EpDouble *epd1, double value)
00699 {
00700 EpDouble epd2;
00701 double tmp;
00702 int exponent, diff;
00703
00704 if (EpdIsNan(epd1) || IsNanDouble(value)) {
00705 EpdMakeNan(epd1);
00706 return;
00707 } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00708 int sign;
00709
00710 EpdConvert(value, &epd2);
00711 if (EpdIsInf(epd1) && IsInfDouble(value)) {
00712 sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00713 if (sign == 0)
00714 EpdMakeNan(epd1);
00715 } else if (EpdIsInf(&epd2)) {
00716 EpdCopy(&epd2, epd1);
00717 }
00718 return;
00719 }
00720
00721 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00722
00723 EpdConvert(value, &epd2);
00724 if (epd1->exponent > epd2.exponent) {
00725 diff = epd1->exponent - epd2.exponent;
00726 if (diff <= EPD_MAX_BIN)
00727 tmp = epd1->type.value - epd2.type.value / pow((double)2.0, (double)diff);
00728 else
00729 tmp = epd1->type.value;
00730 exponent = epd1->exponent;
00731 } else if (epd1->exponent < epd2.exponent) {
00732 diff = epd2.exponent - epd1->exponent;
00733 if (diff <= EPD_MAX_BIN)
00734 tmp = epd1->type.value / pow((double)2.0, (double)diff) - epd2.type.value;
00735 else
00736 tmp = epd2.type.value * (double)(-1.0);
00737 exponent = epd2.exponent;
00738 } else {
00739 tmp = epd1->type.value - epd2.type.value;
00740 exponent = epd1->exponent;
00741 }
00742 epd1->type.value = tmp;
00743 epd1->exponent = exponent;
00744 EpdNormalize(epd1);
00745 }
00746
00747
00759 void
00760 EpdSubtract2(EpDouble *epd1, EpDouble *epd2)
00761 {
00762 double value;
00763 int exponent, diff;
00764
00765 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00766 EpdMakeNan(epd1);
00767 return;
00768 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00769 int sign;
00770
00771 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00772 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00773 if (sign == 0)
00774 EpdMakeNan(epd1);
00775 } else if (EpdIsInf(epd2)) {
00776 EpdCopy(epd2, epd1);
00777 }
00778 return;
00779 }
00780
00781 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00782 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00783
00784 if (epd1->exponent > epd2->exponent) {
00785 diff = epd1->exponent - epd2->exponent;
00786 if (diff <= EPD_MAX_BIN) {
00787 value = epd1->type.value -
00788 epd2->type.value / pow((double)2.0, (double)diff);
00789 } else
00790 value = epd1->type.value;
00791 exponent = epd1->exponent;
00792 } else if (epd1->exponent < epd2->exponent) {
00793 diff = epd2->exponent - epd1->exponent;
00794 if (diff <= EPD_MAX_BIN) {
00795 value = epd1->type.value / pow((double)2.0, (double)diff) -
00796 epd2->type.value;
00797 } else
00798 value = epd2->type.value * (double)(-1.0);
00799 exponent = epd2->exponent;
00800 } else {
00801 value = epd1->type.value - epd2->type.value;
00802 exponent = epd1->exponent;
00803 }
00804 epd1->type.value = value;
00805 epd1->exponent = exponent;
00806 EpdNormalize(epd1);
00807 }
00808
00809
00821 void
00822 EpdSubtract3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
00823 {
00824 double value;
00825 int exponent, diff;
00826
00827 if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00828 EpdMakeNan(epd3);
00829 return;
00830 } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00831 int sign;
00832
00833 if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00834 sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00835 if (sign == 0)
00836 EpdCopy(epd1, epd3);
00837 else
00838 EpdMakeNan(epd3);
00839 } else if (EpdIsInf(epd1)) {
00840 EpdCopy(epd1, epd1);
00841 } else {
00842 sign = epd2->type.bits.sign ^ 0x1;
00843 EpdMakeInf(epd3, sign);
00844 }
00845 return;
00846 }
00847
00848 assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00849 assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00850
00851 if (epd1->exponent > epd2->exponent) {
00852 diff = epd1->exponent - epd2->exponent;
00853 if (diff <= EPD_MAX_BIN) {
00854 value = epd1->type.value -
00855 epd2->type.value / pow((double)2.0, (double)diff);
00856 } else
00857 value = epd1->type.value;
00858 exponent = epd1->exponent;
00859 } else if (epd1->exponent < epd2->exponent) {
00860 diff = epd2->exponent - epd1->exponent;
00861 if (diff <= EPD_MAX_BIN) {
00862 value = epd1->type.value / pow((double)2.0, (double)diff) -
00863 epd2->type.value;
00864 } else
00865 value = epd2->type.value * (double)(-1.0);
00866 exponent = epd2->exponent;
00867 } else {
00868 value = epd1->type.value - epd2->type.value;
00869 exponent = epd1->exponent;
00870 }
00871 epd3->type.value = value;
00872 epd3->exponent = exponent;
00873 EpdNormalize(epd3);
00874 }
00875
00876
00888 void
00889 EpdPow2(int n, EpDouble *epd)
00890 {
00891 if (n <= EPD_MAX_BIN) {
00892 EpdConvert(pow((double)2.0, (double)n), epd);
00893 } else {
00894 EpDouble epd1, epd2;
00895 int n1, n2;
00896
00897 n1 = n / 2;
00898 n2 = n - n1;
00899 EpdPow2(n1, &epd1);
00900 EpdPow2(n2, &epd2);
00901 EpdMultiply3(&epd1, &epd2, epd);
00902 }
00903 }
00904
00905
00917 void
00918 EpdPow2Decimal(int n, EpDouble *epd)
00919 {
00920 if (n <= EPD_MAX_BIN) {
00921 epd->type.value = pow((double)2.0, (double)n);
00922 epd->exponent = 0;
00923 EpdNormalizeDecimal(epd);
00924 } else {
00925 EpDouble epd1, epd2;
00926 int n1, n2;
00927
00928 n1 = n / 2;
00929 n2 = n - n1;
00930 EpdPow2Decimal(n1, &epd1);
00931 EpdPow2Decimal(n2, &epd2);
00932 EpdMultiply3Decimal(&epd1, &epd2, epd);
00933 }
00934 }
00935
00936
00948 void
00949 EpdNormalize(EpDouble *epd)
00950 {
00951 int exponent;
00952
00953 if (IsNanOrInfDouble(epd->type.value)) {
00954 epd->exponent = 0;
00955 return;
00956 }
00957
00958 exponent = EpdGetExponent(epd->type.value);
00959 if (exponent == EPD_MAX_BIN)
00960 return;
00961 exponent -= EPD_MAX_BIN;
00962 epd->type.bits.exponent = EPD_MAX_BIN;
00963 epd->exponent += exponent;
00964 }
00965
00966
00978 void
00979 EpdNormalizeDecimal(EpDouble *epd)
00980 {
00981 int exponent;
00982
00983 if (IsNanOrInfDouble(epd->type.value)) {
00984 epd->exponent = 0;
00985 return;
00986 }
00987
00988 exponent = EpdGetExponentDecimal(epd->type.value);
00989 epd->type.value /= pow((double)10.0, (double)exponent);
00990 epd->exponent += exponent;
00991 }
00992
00993
01005 void
01006 EpdGetValueAndDecimalExponent(EpDouble *epd, double *value, int *exponent)
01007 {
01008 EpDouble epd1, epd2;
01009
01010 if (EpdIsNanOrInf(epd))
01011 return;
01012
01013 if (EpdIsZero(epd)) {
01014 *value = 0.0;
01015 *exponent = 0;
01016 return;
01017 }
01018
01019 epd1.type.value = epd->type.value;
01020 epd1.exponent = 0;
01021 EpdPow2Decimal(epd->exponent, &epd2);
01022 EpdMultiply2Decimal(&epd1, &epd2);
01023
01024 *value = epd1.type.value;
01025 *exponent = epd1.exponent;
01026 }
01027
01039 int
01040 EpdGetExponent(double value)
01041 {
01042 int exponent;
01043 EpDouble epd;
01044
01045 epd.type.value = value;
01046 exponent = epd.type.bits.exponent;
01047 return(exponent);
01048 }
01049
01050
01062 int
01063 EpdGetExponentDecimal(double value)
01064 {
01065 char *pos, str[24];
01066 int exponent;
01067
01068 sprintf(str, "%E", value);
01069 pos = strstr(str, "E");
01070 sscanf(pos, "E%d", &exponent);
01071 return(exponent);
01072 }
01073
01074
01086 void
01087 EpdMakeInf(EpDouble *epd, int sign)
01088 {
01089 epd->type.bits.mantissa1 = 0;
01090 epd->type.bits.mantissa0 = 0;
01091 epd->type.bits.exponent = EPD_EXP_INF;
01092 epd->type.bits.sign = sign;
01093 epd->exponent = 0;
01094 }
01095
01096
01108 void
01109 EpdMakeZero(EpDouble *epd, int sign)
01110 {
01111 epd->type.bits.mantissa1 = 0;
01112 epd->type.bits.mantissa0 = 0;
01113 epd->type.bits.exponent = 0;
01114 epd->type.bits.sign = sign;
01115 epd->exponent = 0;
01116 }
01117
01118
01130 void
01131 EpdMakeNan(EpDouble *epd)
01132 {
01133 epd->type.nan.mantissa1 = 0;
01134 epd->type.nan.mantissa0 = 0;
01135 epd->type.nan.quiet_bit = 1;
01136 epd->type.nan.exponent = EPD_EXP_INF;
01137 epd->type.nan.sign = 1;
01138 epd->exponent = 0;
01139 }
01140
01141
01153 void
01154 EpdCopy(EpDouble *from, EpDouble *to)
01155 {
01156 to->type.value = from->type.value;
01157 to->exponent = from->exponent;
01158 }
01159
01160
01172 int
01173 EpdIsInf(EpDouble *epd)
01174 {
01175 return(IsInfDouble(epd->type.value));
01176 }
01177
01178
01190 int
01191 EpdIsZero(EpDouble *epd)
01192 {
01193 if (epd->type.value == 0.0)
01194 return(1);
01195 else
01196 return(0);
01197 }
01198
01199
01211 int
01212 EpdIsNan(EpDouble *epd)
01213 {
01214 return(IsNanDouble(epd->type.value));
01215 }
01216
01217
01229 int
01230 EpdIsNanOrInf(EpDouble *epd)
01231 {
01232 return(IsNanOrInfDouble(epd->type.value));
01233 }
01234
01235
01247 int
01248 IsInfDouble(double value)
01249 {
01250 IeeeDouble *ptr = (IeeeDouble *)(&value);
01251
01252 if (ptr->exponent == EPD_EXP_INF &&
01253 ptr->mantissa0 == 0 &&
01254 ptr->mantissa1 == 0) {
01255 if (ptr->sign == 0)
01256 return(1);
01257 else
01258 return(-1);
01259 }
01260 return(0);
01261 }
01262
01263
01275 int
01276 IsNanDouble(double value)
01277 {
01278 IeeeNan *ptr = (IeeeNan *)(&value);
01279
01280 if (ptr->exponent == EPD_EXP_INF &&
01281 ptr->sign == 1 &&
01282 ptr->quiet_bit == 1 &&
01283 ptr->mantissa0 == 0 &&
01284 ptr->mantissa1 == 0) {
01285 return(1);
01286 }
01287 return(0);
01288 }
01289
01290
01302 int
01303 IsNanOrInfDouble(double value)
01304 {
01305 IeeeNan *ptr = (IeeeNan *)(&value);
01306
01307 if (ptr->exponent == EPD_EXP_INF &&
01308 ptr->mantissa0 == 0 &&
01309 ptr->mantissa1 == 0 &&
01310 (ptr->sign == 1 || ptr->quiet_bit == 0)) {
01311 return(1);
01312 }
01313 return(0);
01314 }
