1 /**CFile***********************************************************************
7 Synopsis [Arithmetic functions with extended double precision.]
15 Copyright [Copyright (c) 1995-2004, Regents of the University of Colorado
19 Redistribution and use in source and binary forms, with or without
20 modification, are permitted provided that the following conditions
23 Redistributions of source code must retain the above copyright
24 notice, this list of conditions and the following disclaimer.
26 Redistributions in binary form must reproduce the above copyright
27 notice, this list of conditions and the following disclaimer in the
28 documentation and/or other materials provided with the distribution.
30 Neither the name of the University of Colorado nor the names of its
31 contributors may be used to endorse or promote products derived from
32 this software without specific prior written permission.
34 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
35 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
36 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
37 FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
38 COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
39 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
40 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
41 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
42 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
43 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
44 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
45 POSSIBILITY OF SUCH DAMAGE.]
47 Revision [$Id: epd.c,v 1.10 2004/08/13 18:20:30 fabio Exp $]
49 ******************************************************************************/
59 /**Function********************************************************************
61 Synopsis [Allocates an EpDouble struct.]
63 Description [Allocates an EpDouble struct.]
69 ******************************************************************************/
75 epd
= ALLOC(EpDouble
, 1);
80 /**Function********************************************************************
82 Synopsis [Compares two EpDouble struct.]
84 Description [Compares two EpDouble struct.]
90 ******************************************************************************/
92 EpdCmp(const char *key1
, const char *key2
)
94 EpDouble
*epd1
= (EpDouble
*) key1
;
95 EpDouble
*epd2
= (EpDouble
*) key2
;
96 if (epd1
->type
.value
!= epd2
->type
.value
||
97 epd1
->exponent
!= epd2
->exponent
) {
104 /**Function********************************************************************
106 Synopsis [Frees an EpDouble struct.]
108 Description [Frees an EpDouble struct.]
114 ******************************************************************************/
116 EpdFree(EpDouble
*epd
)
122 /**Function********************************************************************
124 Synopsis [Converts an arbitrary precision double value to a string.]
126 Description [Converts an arbitrary precision double value to a string.]
132 ******************************************************************************/
134 EpdGetString(EpDouble
*epd
, char *str
)
140 if (IsNanDouble(epd
->type
.value
)) {
143 } else if (IsInfDouble(epd
->type
.value
)) {
144 if (epd
->type
.bits
.sign
== 1)
145 sprintf(str
, "-Inf");
151 assert(epd
->type
.bits
.exponent
== EPD_MAX_BIN
||
152 epd
->type
.bits
.exponent
== 0);
154 EpdGetValueAndDecimalExponent(epd
, &value
, &exponent
);
155 sprintf(str
, "%e", value
);
156 pos
= strstr(str
, "e");
159 sprintf(pos
+ 1, "+0%d", exponent
);
161 sprintf(pos
+ 1, "+%d", exponent
);
165 sprintf(pos
+ 1, "-0%d", exponent
);
167 sprintf(pos
+ 1, "-%d", exponent
);
172 /**Function********************************************************************
174 Synopsis [Converts double to EpDouble struct.]
176 Description [Converts double to EpDouble struct.]
182 ******************************************************************************/
184 EpdConvert(double value
, EpDouble
*epd
)
186 epd
->type
.value
= value
;
192 /**Function********************************************************************
194 Synopsis [Multiplies two arbitrary precision double values.]
196 Description [Multiplies two arbitrary precision double values.]
202 ******************************************************************************/
204 EpdMultiply(EpDouble
*epd1
, double value
)
210 if (EpdIsNan(epd1
) || IsNanDouble(value
)) {
213 } else if (EpdIsInf(epd1
) || IsInfDouble(value
)) {
216 EpdConvert(value
, &epd2
);
217 sign
= epd1
->type
.bits
.sign
^ epd2
.type
.bits
.sign
;
218 EpdMakeInf(epd1
, sign
);
222 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
224 EpdConvert(value
, &epd2
);
225 tmp
= epd1
->type
.value
* epd2
.type
.value
;
226 exponent
= epd1
->exponent
+ epd2
.exponent
;
227 epd1
->type
.value
= tmp
;
228 epd1
->exponent
= exponent
;
233 /**Function********************************************************************
235 Synopsis [Multiplies two arbitrary precision double values.]
237 Description [Multiplies two arbitrary precision double values.]
243 ******************************************************************************/
245 EpdMultiply2(EpDouble
*epd1
, EpDouble
*epd2
)
250 if (EpdIsNan(epd1
) || EpdIsNan(epd2
)) {
253 } else if (EpdIsInf(epd1
) || EpdIsInf(epd2
)) {
256 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
257 EpdMakeInf(epd1
, sign
);
261 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
262 assert(epd2
->type
.bits
.exponent
== EPD_MAX_BIN
);
264 value
= epd1
->type
.value
* epd2
->type
.value
;
265 exponent
= epd1
->exponent
+ epd2
->exponent
;
266 epd1
->type
.value
= value
;
267 epd1
->exponent
= exponent
;
272 /**Function********************************************************************
274 Synopsis [Multiplies two arbitrary precision double values.]
276 Description [Multiplies two arbitrary precision double values.]
282 ******************************************************************************/
284 EpdMultiply2Decimal(EpDouble
*epd1
, EpDouble
*epd2
)
289 if (EpdIsNan(epd1
) || EpdIsNan(epd2
)) {
292 } else if (EpdIsInf(epd1
) || EpdIsInf(epd2
)) {
295 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
296 EpdMakeInf(epd1
, sign
);
300 value
= epd1
->type
.value
* epd2
->type
.value
;
301 exponent
= epd1
->exponent
+ epd2
->exponent
;
302 epd1
->type
.value
= value
;
303 epd1
->exponent
= exponent
;
304 EpdNormalizeDecimal(epd1
);
308 /**Function********************************************************************
310 Synopsis [Multiplies two arbitrary precision double values.]
312 Description [Multiplies two arbitrary precision double values.]
318 ******************************************************************************/
320 EpdMultiply3(EpDouble
*epd1
, EpDouble
*epd2
, EpDouble
*epd3
)
322 if (EpdIsNan(epd1
) || EpdIsNan(epd2
)) {
325 } else if (EpdIsInf(epd1
) || EpdIsInf(epd2
)) {
328 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
329 EpdMakeInf(epd3
, sign
);
333 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
334 assert(epd2
->type
.bits
.exponent
== EPD_MAX_BIN
);
336 epd3
->type
.value
= epd1
->type
.value
* epd2
->type
.value
;
337 epd3
->exponent
= epd1
->exponent
+ epd2
->exponent
;
342 /**Function********************************************************************
344 Synopsis [Multiplies two arbitrary precision double values.]
346 Description [Multiplies two arbitrary precision double values.]
352 ******************************************************************************/
354 EpdMultiply3Decimal(EpDouble
*epd1
, EpDouble
*epd2
, EpDouble
*epd3
)
356 if (EpdIsNan(epd1
) || EpdIsNan(epd2
)) {
359 } else if (EpdIsInf(epd1
) || EpdIsInf(epd2
)) {
362 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
363 EpdMakeInf(epd3
, sign
);
367 epd3
->type
.value
= epd1
->type
.value
* epd2
->type
.value
;
368 epd3
->exponent
= epd1
->exponent
+ epd2
->exponent
;
369 EpdNormalizeDecimal(epd3
);
373 /**Function********************************************************************
375 Synopsis [Divides two arbitrary precision double values.]
377 Description [Divides two arbitrary precision double values.]
383 ******************************************************************************/
385 EpdDivide(EpDouble
*epd1
, double value
)
391 if (EpdIsNan(epd1
) || IsNanDouble(value
)) {
394 } else if (EpdIsInf(epd1
) || IsInfDouble(value
)) {
397 EpdConvert(value
, &epd2
);
398 if (EpdIsInf(epd1
) && IsInfDouble(value
)) {
400 } else if (EpdIsInf(epd1
)) {
401 sign
= epd1
->type
.bits
.sign
^ epd2
.type
.bits
.sign
;
402 EpdMakeInf(epd1
, sign
);
404 sign
= epd1
->type
.bits
.sign
^ epd2
.type
.bits
.sign
;
405 EpdMakeZero(epd1
, sign
);
415 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
417 EpdConvert(value
, &epd2
);
418 tmp
= epd1
->type
.value
/ epd2
.type
.value
;
419 exponent
= epd1
->exponent
- epd2
.exponent
;
420 epd1
->type
.value
= tmp
;
421 epd1
->exponent
= exponent
;
426 /**Function********************************************************************
428 Synopsis [Divides two arbitrary precision double values.]
430 Description [Divides two arbitrary precision double values.]
436 ******************************************************************************/
438 EpdDivide2(EpDouble
*epd1
, EpDouble
*epd2
)
443 if (EpdIsNan(epd1
) || EpdIsNan(epd2
)) {
446 } else if (EpdIsInf(epd1
) || EpdIsInf(epd2
)) {
449 if (EpdIsInf(epd1
) && EpdIsInf(epd2
)) {
451 } else if (EpdIsInf(epd1
)) {
452 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
453 EpdMakeInf(epd1
, sign
);
455 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
456 EpdMakeZero(epd1
, sign
);
461 if (epd2
->type
.value
== 0.0) {
466 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
467 assert(epd2
->type
.bits
.exponent
== EPD_MAX_BIN
);
469 value
= epd1
->type
.value
/ epd2
->type
.value
;
470 exponent
= epd1
->exponent
- epd2
->exponent
;
471 epd1
->type
.value
= value
;
472 epd1
->exponent
= exponent
;
477 /**Function********************************************************************
479 Synopsis [Divides two arbitrary precision double values.]
481 Description [Divides two arbitrary precision double values.]
487 ******************************************************************************/
489 EpdDivide3(EpDouble
*epd1
, EpDouble
*epd2
, EpDouble
*epd3
)
491 if (EpdIsNan(epd1
) || EpdIsNan(epd2
)) {
494 } else if (EpdIsInf(epd1
) || EpdIsInf(epd2
)) {
497 if (EpdIsInf(epd1
) && EpdIsInf(epd2
)) {
499 } else if (EpdIsInf(epd1
)) {
500 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
501 EpdMakeInf(epd3
, sign
);
503 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
504 EpdMakeZero(epd3
, sign
);
509 if (epd2
->type
.value
== 0.0) {
514 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
515 assert(epd2
->type
.bits
.exponent
== EPD_MAX_BIN
);
517 epd3
->type
.value
= epd1
->type
.value
/ epd2
->type
.value
;
518 epd3
->exponent
= epd1
->exponent
- epd2
->exponent
;
523 /**Function********************************************************************
525 Synopsis [Adds two arbitrary precision double values.]
527 Description [Adds two arbitrary precision double values.]
533 ******************************************************************************/
535 EpdAdd(EpDouble
*epd1
, double value
)
541 if (EpdIsNan(epd1
) || IsNanDouble(value
)) {
544 } else if (EpdIsInf(epd1
) || IsInfDouble(value
)) {
547 EpdConvert(value
, &epd2
);
548 if (EpdIsInf(epd1
) && IsInfDouble(value
)) {
549 sign
= epd1
->type
.bits
.sign
^ epd2
.type
.bits
.sign
;
552 } else if (EpdIsInf(&epd2
)) {
553 EpdCopy(&epd2
, epd1
);
558 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
560 EpdConvert(value
, &epd2
);
561 if (epd1
->exponent
> epd2
.exponent
) {
562 diff
= epd1
->exponent
- epd2
.exponent
;
563 if (diff
<= EPD_MAX_BIN
)
564 tmp
= epd1
->type
.value
+ epd2
.type
.value
/ pow((double)2.0, (double)diff
);
566 tmp
= epd1
->type
.value
;
567 exponent
= epd1
->exponent
;
568 } else if (epd1
->exponent
< epd2
.exponent
) {
569 diff
= epd2
.exponent
- epd1
->exponent
;
570 if (diff
<= EPD_MAX_BIN
)
571 tmp
= epd1
->type
.value
/ pow((double)2.0, (double)diff
) + epd2
.type
.value
;
573 tmp
= epd2
.type
.value
;
574 exponent
= epd2
.exponent
;
576 tmp
= epd1
->type
.value
+ epd2
.type
.value
;
577 exponent
= epd1
->exponent
;
579 epd1
->type
.value
= tmp
;
580 epd1
->exponent
= exponent
;
585 /**Function********************************************************************
587 Synopsis [Adds two arbitrary precision double values.]
589 Description [Adds two arbitrary precision double values.]
595 ******************************************************************************/
597 EpdAdd2(EpDouble
*epd1
, EpDouble
*epd2
)
602 if (EpdIsNan(epd1
) || EpdIsNan(epd2
)) {
605 } else if (EpdIsInf(epd1
) || EpdIsInf(epd2
)) {
608 if (EpdIsInf(epd1
) && EpdIsInf(epd2
)) {
609 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
612 } else if (EpdIsInf(epd2
)) {
618 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
619 assert(epd2
->type
.bits
.exponent
== EPD_MAX_BIN
);
621 if (epd1
->exponent
> epd2
->exponent
) {
622 diff
= epd1
->exponent
- epd2
->exponent
;
623 if (diff
<= EPD_MAX_BIN
) {
624 value
= epd1
->type
.value
+
625 epd2
->type
.value
/ pow((double)2.0, (double)diff
);
627 value
= epd1
->type
.value
;
628 exponent
= epd1
->exponent
;
629 } else if (epd1
->exponent
< epd2
->exponent
) {
630 diff
= epd2
->exponent
- epd1
->exponent
;
631 if (diff
<= EPD_MAX_BIN
) {
632 value
= epd1
->type
.value
/ pow((double)2.0, (double)diff
) +
635 value
= epd2
->type
.value
;
636 exponent
= epd2
->exponent
;
638 value
= epd1
->type
.value
+ epd2
->type
.value
;
639 exponent
= epd1
->exponent
;
641 epd1
->type
.value
= value
;
642 epd1
->exponent
= exponent
;
647 /**Function********************************************************************
649 Synopsis [Adds two arbitrary precision double values.]
651 Description [Adds two arbitrary precision double values.]
657 ******************************************************************************/
659 EpdAdd3(EpDouble
*epd1
, EpDouble
*epd2
, EpDouble
*epd3
)
664 if (EpdIsNan(epd1
) || EpdIsNan(epd2
)) {
667 } else if (EpdIsInf(epd1
) || EpdIsInf(epd2
)) {
670 if (EpdIsInf(epd1
) && EpdIsInf(epd2
)) {
671 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
676 } else if (EpdIsInf(epd1
)) {
684 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
685 assert(epd2
->type
.bits
.exponent
== EPD_MAX_BIN
);
687 if (epd1
->exponent
> epd2
->exponent
) {
688 diff
= epd1
->exponent
- epd2
->exponent
;
689 if (diff
<= EPD_MAX_BIN
) {
690 value
= epd1
->type
.value
+
691 epd2
->type
.value
/ pow((double)2.0, (double)diff
);
693 value
= epd1
->type
.value
;
694 exponent
= epd1
->exponent
;
695 } else if (epd1
->exponent
< epd2
->exponent
) {
696 diff
= epd2
->exponent
- epd1
->exponent
;
697 if (diff
<= EPD_MAX_BIN
) {
698 value
= epd1
->type
.value
/ pow((double)2.0, (double)diff
) +
701 value
= epd2
->type
.value
;
702 exponent
= epd2
->exponent
;
704 value
= epd1
->type
.value
+ epd2
->type
.value
;
705 exponent
= epd1
->exponent
;
707 epd3
->type
.value
= value
;
708 epd3
->exponent
= exponent
;
713 /**Function********************************************************************
715 Synopsis [Subtracts two arbitrary precision double values.]
717 Description [Subtracts two arbitrary precision double values.]
723 ******************************************************************************/
725 EpdSubtract(EpDouble
*epd1
, double value
)
731 if (EpdIsNan(epd1
) || IsNanDouble(value
)) {
734 } else if (EpdIsInf(epd1
) || IsInfDouble(value
)) {
737 EpdConvert(value
, &epd2
);
738 if (EpdIsInf(epd1
) && IsInfDouble(value
)) {
739 sign
= epd1
->type
.bits
.sign
^ epd2
.type
.bits
.sign
;
742 } else if (EpdIsInf(&epd2
)) {
743 EpdCopy(&epd2
, epd1
);
748 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
750 EpdConvert(value
, &epd2
);
751 if (epd1
->exponent
> epd2
.exponent
) {
752 diff
= epd1
->exponent
- epd2
.exponent
;
753 if (diff
<= EPD_MAX_BIN
)
754 tmp
= epd1
->type
.value
- epd2
.type
.value
/ pow((double)2.0, (double)diff
);
756 tmp
= epd1
->type
.value
;
757 exponent
= epd1
->exponent
;
758 } else if (epd1
->exponent
< epd2
.exponent
) {
759 diff
= epd2
.exponent
- epd1
->exponent
;
760 if (diff
<= EPD_MAX_BIN
)
761 tmp
= epd1
->type
.value
/ pow((double)2.0, (double)diff
) - epd2
.type
.value
;
763 tmp
= epd2
.type
.value
* (double)(-1.0);
764 exponent
= epd2
.exponent
;
766 tmp
= epd1
->type
.value
- epd2
.type
.value
;
767 exponent
= epd1
->exponent
;
769 epd1
->type
.value
= tmp
;
770 epd1
->exponent
= exponent
;
775 /**Function********************************************************************
777 Synopsis [Subtracts two arbitrary precision double values.]
779 Description [Subtracts two arbitrary precision double values.]
785 ******************************************************************************/
787 EpdSubtract2(EpDouble
*epd1
, EpDouble
*epd2
)
792 if (EpdIsNan(epd1
) || EpdIsNan(epd2
)) {
795 } else if (EpdIsInf(epd1
) || EpdIsInf(epd2
)) {
798 if (EpdIsInf(epd1
) && EpdIsInf(epd2
)) {
799 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
802 } else if (EpdIsInf(epd2
)) {
808 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
809 assert(epd2
->type
.bits
.exponent
== EPD_MAX_BIN
);
811 if (epd1
->exponent
> epd2
->exponent
) {
812 diff
= epd1
->exponent
- epd2
->exponent
;
813 if (diff
<= EPD_MAX_BIN
) {
814 value
= epd1
->type
.value
-
815 epd2
->type
.value
/ pow((double)2.0, (double)diff
);
817 value
= epd1
->type
.value
;
818 exponent
= epd1
->exponent
;
819 } else if (epd1
->exponent
< epd2
->exponent
) {
820 diff
= epd2
->exponent
- epd1
->exponent
;
821 if (diff
<= EPD_MAX_BIN
) {
822 value
= epd1
->type
.value
/ pow((double)2.0, (double)diff
) -
825 value
= epd2
->type
.value
* (double)(-1.0);
826 exponent
= epd2
->exponent
;
828 value
= epd1
->type
.value
- epd2
->type
.value
;
829 exponent
= epd1
->exponent
;
831 epd1
->type
.value
= value
;
832 epd1
->exponent
= exponent
;
837 /**Function********************************************************************
839 Synopsis [Subtracts two arbitrary precision double values.]
841 Description [Subtracts two arbitrary precision double values.]
847 ******************************************************************************/
849 EpdSubtract3(EpDouble
*epd1
, EpDouble
*epd2
, EpDouble
*epd3
)
854 if (EpdIsNan(epd1
) || EpdIsNan(epd2
)) {
857 } else if (EpdIsInf(epd1
) || EpdIsInf(epd2
)) {
860 if (EpdIsInf(epd1
) && EpdIsInf(epd2
)) {
861 sign
= epd1
->type
.bits
.sign
^ epd2
->type
.bits
.sign
;
866 } else if (EpdIsInf(epd1
)) {
869 sign
= epd2
->type
.bits
.sign
^ 0x1;
870 EpdMakeInf(epd3
, sign
);
875 assert(epd1
->type
.bits
.exponent
== EPD_MAX_BIN
);
876 assert(epd2
->type
.bits
.exponent
== EPD_MAX_BIN
);
878 if (epd1
->exponent
> epd2
->exponent
) {
879 diff
= epd1
->exponent
- epd2
->exponent
;
880 if (diff
<= EPD_MAX_BIN
) {
881 value
= epd1
->type
.value
-
882 epd2
->type
.value
/ pow((double)2.0, (double)diff
);
884 value
= epd1
->type
.value
;
885 exponent
= epd1
->exponent
;
886 } else if (epd1
->exponent
< epd2
->exponent
) {
887 diff
= epd2
->exponent
- epd1
->exponent
;
888 if (diff
<= EPD_MAX_BIN
) {
889 value
= epd1
->type
.value
/ pow((double)2.0, (double)diff
) -
892 value
= epd2
->type
.value
* (double)(-1.0);
893 exponent
= epd2
->exponent
;
895 value
= epd1
->type
.value
- epd2
->type
.value
;
896 exponent
= epd1
->exponent
;
898 epd3
->type
.value
= value
;
899 epd3
->exponent
= exponent
;
904 /**Function********************************************************************
906 Synopsis [Computes arbitrary precision pow of base 2.]
908 Description [Computes arbitrary precision pow of base 2.]
914 ******************************************************************************/
916 EpdPow2(int n
, EpDouble
*epd
)
918 if (n
<= EPD_MAX_BIN
) {
919 EpdConvert(pow((double)2.0, (double)n
), epd
);
928 EpdMultiply3(&epd1
, &epd2
, epd
);
933 /**Function********************************************************************
935 Synopsis [Computes arbitrary precision pow of base 2.]
937 Description [Computes arbitrary precision pow of base 2.]
943 ******************************************************************************/
945 EpdPow2Decimal(int n
, EpDouble
*epd
)
947 if (n
<= EPD_MAX_BIN
) {
948 epd
->type
.value
= pow((double)2.0, (double)n
);
950 EpdNormalizeDecimal(epd
);
957 EpdPow2Decimal(n1
, &epd1
);
958 EpdPow2Decimal(n2
, &epd2
);
959 EpdMultiply3Decimal(&epd1
, &epd2
, epd
);
964 /**Function********************************************************************
966 Synopsis [Normalize an arbitrary precision double value.]
968 Description [Normalize an arbitrary precision double value.]
974 ******************************************************************************/
976 EpdNormalize(EpDouble
*epd
)
980 if (IsNanOrInfDouble(epd
->type
.value
)) {
985 exponent
= EpdGetExponent(epd
->type
.value
);
986 if (exponent
== EPD_MAX_BIN
)
988 exponent
-= EPD_MAX_BIN
;
989 epd
->type
.bits
.exponent
= EPD_MAX_BIN
;
990 epd
->exponent
+= exponent
;
994 /**Function********************************************************************
996 Synopsis [Normalize an arbitrary precision double value.]
998 Description [Normalize an arbitrary precision double value.]
1004 ******************************************************************************/
1006 EpdNormalizeDecimal(EpDouble
*epd
)
1010 if (IsNanOrInfDouble(epd
->type
.value
)) {
1015 exponent
= EpdGetExponentDecimal(epd
->type
.value
);
1016 epd
->type
.value
/= pow((double)10.0, (double)exponent
);
1017 epd
->exponent
+= exponent
;
1021 /**Function********************************************************************
1023 Synopsis [Returns value and decimal exponent of EpDouble.]
1025 Description [Returns value and decimal exponent of EpDouble.]
1031 ******************************************************************************/
1033 EpdGetValueAndDecimalExponent(EpDouble
*epd
, double *value
, int *exponent
)
1035 EpDouble epd1
, epd2
;
1037 if (EpdIsNanOrInf(epd
))
1040 if (EpdIsZero(epd
)) {
1046 epd1
.type
.value
= epd
->type
.value
;
1048 EpdPow2Decimal(epd
->exponent
, &epd2
);
1049 EpdMultiply2Decimal(&epd1
, &epd2
);
1051 *value
= epd1
.type
.value
;
1052 *exponent
= epd1
.exponent
;
1055 /**Function********************************************************************
1057 Synopsis [Returns the exponent value of a double.]
1059 Description [Returns the exponent value of a double.]
1065 ******************************************************************************/
1067 EpdGetExponent(double value
)
1072 epd
.type
.value
= value
;
1073 exponent
= epd
.type
.bits
.exponent
;
1078 /**Function********************************************************************
1080 Synopsis [Returns the decimal exponent value of a double.]
1082 Description [Returns the decimal exponent value of a double.]
1088 ******************************************************************************/
1090 EpdGetExponentDecimal(double value
)
1095 sprintf(str
, "%E", value
);
1096 pos
= strstr(str
, "E");
1097 sscanf(pos
, "E%d", &exponent
);
1102 /**Function********************************************************************
1104 Synopsis [Makes EpDouble Inf.]
1106 Description [Makes EpDouble Inf.]
1112 ******************************************************************************/
1114 EpdMakeInf(EpDouble
*epd
, int sign
)
1116 epd
->type
.bits
.mantissa1
= 0;
1117 epd
->type
.bits
.mantissa0
= 0;
1118 epd
->type
.bits
.exponent
= EPD_EXP_INF
;
1119 epd
->type
.bits
.sign
= sign
;
1124 /**Function********************************************************************
1126 Synopsis [Makes EpDouble Zero.]
1128 Description [Makes EpDouble Zero.]
1134 ******************************************************************************/
1136 EpdMakeZero(EpDouble
*epd
, int sign
)
1138 epd
->type
.bits
.mantissa1
= 0;
1139 epd
->type
.bits
.mantissa0
= 0;
1140 epd
->type
.bits
.exponent
= 0;
1141 epd
->type
.bits
.sign
= sign
;
1146 /**Function********************************************************************
1148 Synopsis [Makes EpDouble NaN.]
1150 Description [Makes EpDouble NaN.]
1156 ******************************************************************************/
1158 EpdMakeNan(EpDouble
*epd
)
1160 epd
->type
.nan
.mantissa1
= 0;
1161 epd
->type
.nan
.mantissa0
= 0;
1162 epd
->type
.nan
.quiet_bit
= 1;
1163 epd
->type
.nan
.exponent
= EPD_EXP_INF
;
1164 epd
->type
.nan
.sign
= 1;
1169 /**Function********************************************************************
1171 Synopsis [Copies a EpDouble struct.]
1173 Description [Copies a EpDouble struct.]
1179 ******************************************************************************/
1181 EpdCopy(EpDouble
*from
, EpDouble
*to
)
1183 to
->type
.value
= from
->type
.value
;
1184 to
->exponent
= from
->exponent
;
1188 /**Function********************************************************************
1190 Synopsis [Checks whether the value is Inf.]
1192 Description [Checks whether the value is Inf.]
1198 ******************************************************************************/
1200 EpdIsInf(EpDouble
*epd
)
1202 return(IsInfDouble(epd
->type
.value
));
1206 /**Function********************************************************************
1208 Synopsis [Checks whether the value is Zero.]
1210 Description [Checks whether the value is Zero.]
1216 ******************************************************************************/
1218 EpdIsZero(EpDouble
*epd
)
1220 if (epd
->type
.value
== 0.0)
1227 /**Function********************************************************************
1229 Synopsis [Checks whether the value is NaN.]
1231 Description [Checks whether the value is NaN.]
1237 ******************************************************************************/
1239 EpdIsNan(EpDouble
*epd
)
1241 return(IsNanDouble(epd
->type
.value
));
1245 /**Function********************************************************************
1247 Synopsis [Checks whether the value is NaN or Inf.]
1249 Description [Checks whether the value is NaN or Inf.]
1255 ******************************************************************************/
1257 EpdIsNanOrInf(EpDouble
*epd
)
1259 return(IsNanOrInfDouble(epd
->type
.value
));
1263 /**Function********************************************************************
1265 Synopsis [Checks whether the value is Inf.]
1267 Description [Checks whether the value is Inf.]
1273 ******************************************************************************/
1275 IsInfDouble(double value
)
1280 if (val
.bits
.exponent
== EPD_EXP_INF
&&
1281 val
.bits
.mantissa0
== 0 &&
1282 val
.bits
.mantissa1
== 0) {
1283 if (val
.bits
.sign
== 0)
1292 /**Function********************************************************************
1294 Synopsis [Checks whether the value is NaN.]
1296 Description [Checks whether the value is NaN.]
1302 ******************************************************************************/
1304 IsNanDouble(double value
)
1309 if (val
.nan
.exponent
== EPD_EXP_INF
&&
1310 val
.nan
.sign
== 1 &&
1311 val
.nan
.quiet_bit
== 1 &&
1312 val
.nan
.mantissa0
== 0 &&
1313 val
.nan
.mantissa1
== 0) {
1320 /**Function********************************************************************
1322 Synopsis [Checks whether the value is NaN or Inf.]
1324 Description [Checks whether the value is NaN or Inf.]
1330 ******************************************************************************/
1332 IsNanOrInfDouble(double value
)
1337 if (val
.nan
.exponent
== EPD_EXP_INF
&&
1338 val
.nan
.mantissa0
== 0 &&
1339 val
.nan
.mantissa1
== 0 &&
1340 (val
.nan
.sign
== 1 || val
.nan
.quiet_bit
== 0)) {