abc-master
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
epd.c
Go to the documentation of this file.
1 /**CFile***********************************************************************
2 
3  FileName [epd.c]
4 
5  PackageName [epd]
6 
7  Synopsis [Arithmetic functions with extended double precision.]
8 
9  Description []
10 
11  SeeAlso []
12 
13  Author [In-Ho Moon]
14 
15  Copyright [Copyright (c) 1995-2004, Regents of the University of Colorado
16 
17  All rights reserved.
18 
19  Redistribution and use in source and binary forms, with or without
20  modification, are permitted provided that the following conditions
21  are met:
22 
23  Redistributions of source code must retain the above copyright
24  notice, this list of conditions and the following disclaimer.
25 
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.
29 
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.
33 
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.]
46 
47  Revision [$Id: epd.c,v 1.10 2004/08/13 18:20:30 fabio Exp $]
48 
49 ******************************************************************************/
50 
51 #include <stdio.h>
52 #include <stdlib.h>
53 #include <string.h>
54 #include <math.h>
55 #include "misc/util/util_hack.h"
56 #include "epd.h"
57 
58 ABC_NAMESPACE_IMPL_START
59 
60 /**Function********************************************************************
61 
62  Synopsis [Allocates an EpDouble struct.]
63 
64  Description [Allocates an EpDouble struct.]
65 
66  SideEffects []
67 
68  SeeAlso []
69 
70 ******************************************************************************/
71 EpDouble *
72 EpdAlloc(void)
73 {
74  EpDouble *epd;
75 
76  epd = ABC_ALLOC(EpDouble, 1);
77  return(epd);
78 }
79 
80 
81 /**Function********************************************************************
82 
83  Synopsis [Compares two EpDouble struct.]
84 
85  Description [Compares two EpDouble struct.]
86 
87  SideEffects []
88 
89  SeeAlso []
90 
91 ******************************************************************************/
92 int
93 EpdCmp(const char *key1, const char *key2)
94 {
95  EpDouble *epd1 = (EpDouble *) key1;
96  EpDouble *epd2 = (EpDouble *) key2;
97  if (epd1->type.value != epd2->type.value ||
98  epd1->exponent != epd2->exponent) {
99  return(1);
100  }
101  return(0);
102 }
103 
104 
105 /**Function********************************************************************
106 
107  Synopsis [Frees an EpDouble struct.]
108 
109  Description [Frees an EpDouble struct.]
110 
111  SideEffects []
112 
113  SeeAlso []
114 
115 ******************************************************************************/
116 void
117 EpdFree(EpDouble *epd)
118 {
119  ABC_FREE(epd);
120 }
121 
122 
123 /**Function********************************************************************
124 
125  Synopsis [Converts an arbitrary precision double value to a string.]
126 
127  Description [Converts an arbitrary precision double value to a string.]
128 
129  SideEffects []
130 
131  SeeAlso []
132 
133 ******************************************************************************/
134 void
135 EpdGetString(EpDouble *epd, char *str)
136 {
137  double value;
138  int exponent;
139  char *pos;
140 
141  if (IsNanDouble(epd->type.value)) {
142  sprintf(str, "NaN");
143  return;
144  } else if (IsInfDouble(epd->type.value)) {
145  if (epd->type.bits.sign == 1)
146  sprintf(str, "-Inf");
147  else
148  sprintf(str, "Inf");
149  return;
150  }
151 
152  assert(epd->type.bits.exponent == EPD_MAX_BIN ||
153  epd->type.bits.exponent == 0);
154 
155  EpdGetValueAndDecimalExponent(epd, &value, &exponent);
156  sprintf(str, "%e", value);
157  pos = strstr(str, "e");
158  if (exponent >= 0) {
159  if (exponent < 10)
160  sprintf(pos + 1, "+0%d", exponent);
161  else
162  sprintf(pos + 1, "+%d", exponent);
163  } else {
164  exponent *= -1;
165  if (exponent < 10)
166  sprintf(pos + 1, "-0%d", exponent);
167  else
168  sprintf(pos + 1, "-%d", exponent);
169  }
170 }
171 
172 
173 /**Function********************************************************************
174 
175  Synopsis [Converts double to EpDouble struct.]
176 
177  Description [Converts double to EpDouble struct.]
178 
179  SideEffects []
180 
181  SeeAlso []
182 
183 ******************************************************************************/
184 void
185 EpdConvert(double value, EpDouble *epd)
186 {
187  epd->type.value = value;
188  epd->exponent = 0;
189  EpdNormalize(epd);
190 }
191 
192 
193 /**Function********************************************************************
194 
195  Synopsis [Multiplies two arbitrary precision double values.]
196 
197  Description [Multiplies two arbitrary precision double values.]
198 
199  SideEffects []
200 
201  SeeAlso []
202 
203 ******************************************************************************/
204 void
205 EpdMultiply(EpDouble *epd1, double value)
206 {
207  EpDouble epd2;
208  double tmp;
209  int exponent;
210 
211  if (EpdIsNan(epd1) || IsNanDouble(value)) {
212  EpdMakeNan(epd1);
213  return;
214  } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
215  int sign;
216 
217  EpdConvert(value, &epd2);
218  sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
219  EpdMakeInf(epd1, sign);
220  return;
221  }
222 
223  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
224 
225  EpdConvert(value, &epd2);
226  tmp = epd1->type.value * epd2.type.value;
227  exponent = epd1->exponent + epd2.exponent;
228  epd1->type.value = tmp;
229  epd1->exponent = exponent;
230  EpdNormalize(epd1);
231 }
232 
233 
234 /**Function********************************************************************
235 
236  Synopsis [Multiplies two arbitrary precision double values.]
237 
238  Description [Multiplies two arbitrary precision double values.]
239 
240  SideEffects []
241 
242  SeeAlso []
243 
244 ******************************************************************************/
245 void
246 EpdMultiply2(EpDouble *epd1, EpDouble *epd2)
247 {
248  double value;
249  int exponent;
250 
251  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
252  EpdMakeNan(epd1);
253  return;
254  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
255  int sign;
256 
257  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
258  EpdMakeInf(epd1, sign);
259  return;
260  }
261 
262  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
263  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
264 
265  value = epd1->type.value * epd2->type.value;
266  exponent = epd1->exponent + epd2->exponent;
267  epd1->type.value = value;
268  epd1->exponent = exponent;
269  EpdNormalize(epd1);
270 }
271 
272 
273 /**Function********************************************************************
274 
275  Synopsis [Multiplies two arbitrary precision double values.]
276 
277  Description [Multiplies two arbitrary precision double values.]
278 
279  SideEffects []
280 
281  SeeAlso []
282 
283 ******************************************************************************/
284 void
285 EpdMultiply2Decimal(EpDouble *epd1, EpDouble *epd2)
286 {
287  double value;
288  int exponent;
289 
290  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
291  EpdMakeNan(epd1);
292  return;
293  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
294  int sign;
295 
296  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
297  EpdMakeInf(epd1, sign);
298  return;
299  }
300 
301  value = epd1->type.value * epd2->type.value;
302  exponent = epd1->exponent + epd2->exponent;
303  epd1->type.value = value;
304  epd1->exponent = exponent;
305  EpdNormalizeDecimal(epd1);
306 }
307 
308 
309 /**Function********************************************************************
310 
311  Synopsis [Multiplies two arbitrary precision double values.]
312 
313  Description [Multiplies two arbitrary precision double values.]
314 
315  SideEffects []
316 
317  SeeAlso []
318 
319 ******************************************************************************/
320 void
321 EpdMultiply3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
322 {
323  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
324  EpdMakeNan(epd1);
325  return;
326  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
327  int sign;
328 
329  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
330  EpdMakeInf(epd3, sign);
331  return;
332  }
333 
334  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
335  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
336 
337  epd3->type.value = epd1->type.value * epd2->type.value;
338  epd3->exponent = epd1->exponent + epd2->exponent;
339  EpdNormalize(epd3);
340 }
341 
342 
343 /**Function********************************************************************
344 
345  Synopsis [Multiplies two arbitrary precision double values.]
346 
347  Description [Multiplies two arbitrary precision double values.]
348 
349  SideEffects []
350 
351  SeeAlso []
352 
353 ******************************************************************************/
354 void
355 EpdMultiply3Decimal(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
356 {
357  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
358  EpdMakeNan(epd1);
359  return;
360  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
361  int sign;
362 
363  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
364  EpdMakeInf(epd3, sign);
365  return;
366  }
367 
368  epd3->type.value = epd1->type.value * epd2->type.value;
369  epd3->exponent = epd1->exponent + epd2->exponent;
370  EpdNormalizeDecimal(epd3);
371 }
372 
373 
374 /**Function********************************************************************
375 
376  Synopsis [Divides two arbitrary precision double values.]
377 
378  Description [Divides two arbitrary precision double values.]
379 
380  SideEffects []
381 
382  SeeAlso []
383 
384 ******************************************************************************/
385 void
386 EpdDivide(EpDouble *epd1, double value)
387 {
388  EpDouble epd2;
389  double tmp;
390  int exponent;
391 
392  if (EpdIsNan(epd1) || IsNanDouble(value)) {
393  EpdMakeNan(epd1);
394  return;
395  } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
396  int sign;
397 
398  EpdConvert(value, &epd2);
399  if (EpdIsInf(epd1) && IsInfDouble(value)) {
400  EpdMakeNan(epd1);
401  } else if (EpdIsInf(epd1)) {
402  sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
403  EpdMakeInf(epd1, sign);
404  } else {
405  sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
406  EpdMakeZero(epd1, sign);
407  }
408  return;
409  }
410 
411  if (value == 0.0) {
412  EpdMakeNan(epd1);
413  return;
414  }
415 
416  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
417 
418  EpdConvert(value, &epd2);
419  tmp = epd1->type.value / epd2.type.value;
420  exponent = epd1->exponent - epd2.exponent;
421  epd1->type.value = tmp;
422  epd1->exponent = exponent;
423  EpdNormalize(epd1);
424 }
425 
426 
427 /**Function********************************************************************
428 
429  Synopsis [Divides two arbitrary precision double values.]
430 
431  Description [Divides two arbitrary precision double values.]
432 
433  SideEffects []
434 
435  SeeAlso []
436 
437 ******************************************************************************/
438 void
439 EpdDivide2(EpDouble *epd1, EpDouble *epd2)
440 {
441  double value;
442  int exponent;
443 
444  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
445  EpdMakeNan(epd1);
446  return;
447  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
448  int sign;
449 
450  if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
451  EpdMakeNan(epd1);
452  } else if (EpdIsInf(epd1)) {
453  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
454  EpdMakeInf(epd1, sign);
455  } else {
456  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
457  EpdMakeZero(epd1, sign);
458  }
459  return;
460  }
461 
462  if (epd2->type.value == 0.0) {
463  EpdMakeNan(epd1);
464  return;
465  }
466 
467  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
468  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
469 
470  value = epd1->type.value / epd2->type.value;
471  exponent = epd1->exponent - epd2->exponent;
472  epd1->type.value = value;
473  epd1->exponent = exponent;
474  EpdNormalize(epd1);
475 }
476 
477 
478 /**Function********************************************************************
479 
480  Synopsis [Divides two arbitrary precision double values.]
481 
482  Description [Divides two arbitrary precision double values.]
483 
484  SideEffects []
485 
486  SeeAlso []
487 
488 ******************************************************************************/
489 void
490 EpdDivide3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
491 {
492  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
493  EpdMakeNan(epd3);
494  return;
495  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
496  int sign;
497 
498  if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
499  EpdMakeNan(epd3);
500  } else if (EpdIsInf(epd1)) {
501  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
502  EpdMakeInf(epd3, sign);
503  } else {
504  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
505  EpdMakeZero(epd3, sign);
506  }
507  return;
508  }
509 
510  if (epd2->type.value == 0.0) {
511  EpdMakeNan(epd3);
512  return;
513  }
514 
515  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
516  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
517 
518  epd3->type.value = epd1->type.value / epd2->type.value;
519  epd3->exponent = epd1->exponent - epd2->exponent;
520  EpdNormalize(epd3);
521 }
522 
523 
524 /**Function********************************************************************
525 
526  Synopsis [Adds two arbitrary precision double values.]
527 
528  Description [Adds two arbitrary precision double values.]
529 
530  SideEffects []
531 
532  SeeAlso []
533 
534 ******************************************************************************/
535 void
536 EpdAdd(EpDouble *epd1, double value)
537 {
538  EpDouble epd2;
539  double tmp;
540  int exponent, diff;
541 
542  if (EpdIsNan(epd1) || IsNanDouble(value)) {
543  EpdMakeNan(epd1);
544  return;
545  } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
546  int sign;
547 
548  EpdConvert(value, &epd2);
549  if (EpdIsInf(epd1) && IsInfDouble(value)) {
550  sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
551  if (sign == 1)
552  EpdMakeNan(epd1);
553  } else if (EpdIsInf(&epd2)) {
554  EpdCopy(&epd2, epd1);
555  }
556  return;
557  }
558 
559  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
560 
561  EpdConvert(value, &epd2);
562  if (epd1->exponent > epd2.exponent) {
563  diff = epd1->exponent - epd2.exponent;
564  if (diff <= EPD_MAX_BIN)
565  tmp = epd1->type.value + epd2.type.value / pow((double)2.0, (double)diff);
566  else
567  tmp = epd1->type.value;
568  exponent = epd1->exponent;
569  } else if (epd1->exponent < epd2.exponent) {
570  diff = epd2.exponent - epd1->exponent;
571  if (diff <= EPD_MAX_BIN)
572  tmp = epd1->type.value / pow((double)2.0, (double)diff) + epd2.type.value;
573  else
574  tmp = epd2.type.value;
575  exponent = epd2.exponent;
576  } else {
577  tmp = epd1->type.value + epd2.type.value;
578  exponent = epd1->exponent;
579  }
580  epd1->type.value = tmp;
581  epd1->exponent = exponent;
582  EpdNormalize(epd1);
583 }
584 
585 
586 /**Function********************************************************************
587 
588  Synopsis [Adds two arbitrary precision double values.]
589 
590  Description [Adds two arbitrary precision double values.]
591 
592  SideEffects []
593 
594  SeeAlso []
595 
596 ******************************************************************************/
597 void
598 EpdAdd2(EpDouble *epd1, EpDouble *epd2)
599 {
600  double value;
601  int exponent, diff;
602 
603  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
604  EpdMakeNan(epd1);
605  return;
606  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
607  int sign;
608 
609  if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
610  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
611  if (sign == 1)
612  EpdMakeNan(epd1);
613  } else if (EpdIsInf(epd2)) {
614  EpdCopy(epd2, epd1);
615  }
616  return;
617  }
618 
619  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
620  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
621 
622  if (epd1->exponent > epd2->exponent) {
623  diff = epd1->exponent - epd2->exponent;
624  if (diff <= EPD_MAX_BIN) {
625  value = epd1->type.value +
626  epd2->type.value / pow((double)2.0, (double)diff);
627  } else
628  value = epd1->type.value;
629  exponent = epd1->exponent;
630  } else if (epd1->exponent < epd2->exponent) {
631  diff = epd2->exponent - epd1->exponent;
632  if (diff <= EPD_MAX_BIN) {
633  value = epd1->type.value / pow((double)2.0, (double)diff) +
634  epd2->type.value;
635  } else
636  value = epd2->type.value;
637  exponent = epd2->exponent;
638  } else {
639  value = epd1->type.value + epd2->type.value;
640  exponent = epd1->exponent;
641  }
642  epd1->type.value = value;
643  epd1->exponent = exponent;
644  EpdNormalize(epd1);
645 }
646 
647 
648 /**Function********************************************************************
649 
650  Synopsis [Adds two arbitrary precision double values.]
651 
652  Description [Adds two arbitrary precision double values.]
653 
654  SideEffects []
655 
656  SeeAlso []
657 
658 ******************************************************************************/
659 void
660 EpdAdd3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
661 {
662  double value;
663  int exponent, diff;
664 
665  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
666  EpdMakeNan(epd3);
667  return;
668  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
669  int sign;
670 
671  if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
672  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
673  if (sign == 1)
674  EpdMakeNan(epd3);
675  else
676  EpdCopy(epd1, epd3);
677  } else if (EpdIsInf(epd1)) {
678  EpdCopy(epd1, epd3);
679  } else {
680  EpdCopy(epd2, epd3);
681  }
682  return;
683  }
684 
685  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
686  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
687 
688  if (epd1->exponent > epd2->exponent) {
689  diff = epd1->exponent - epd2->exponent;
690  if (diff <= EPD_MAX_BIN) {
691  value = epd1->type.value +
692  epd2->type.value / pow((double)2.0, (double)diff);
693  } else
694  value = epd1->type.value;
695  exponent = epd1->exponent;
696  } else if (epd1->exponent < epd2->exponent) {
697  diff = epd2->exponent - epd1->exponent;
698  if (diff <= EPD_MAX_BIN) {
699  value = epd1->type.value / pow((double)2.0, (double)diff) +
700  epd2->type.value;
701  } else
702  value = epd2->type.value;
703  exponent = epd2->exponent;
704  } else {
705  value = epd1->type.value + epd2->type.value;
706  exponent = epd1->exponent;
707  }
708  epd3->type.value = value;
709  epd3->exponent = exponent;
710  EpdNormalize(epd3);
711 }
712 
713 
714 /**Function********************************************************************
715 
716  Synopsis [Subtracts two arbitrary precision double values.]
717 
718  Description [Subtracts two arbitrary precision double values.]
719 
720  SideEffects []
721 
722  SeeAlso []
723 
724 ******************************************************************************/
725 void
726 EpdSubtract(EpDouble *epd1, double value)
727 {
728  EpDouble epd2;
729  double tmp;
730  int exponent, diff;
731 
732  if (EpdIsNan(epd1) || IsNanDouble(value)) {
733  EpdMakeNan(epd1);
734  return;
735  } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
736  int sign;
737 
738  EpdConvert(value, &epd2);
739  if (EpdIsInf(epd1) && IsInfDouble(value)) {
740  sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
741  if (sign == 0)
742  EpdMakeNan(epd1);
743  } else if (EpdIsInf(&epd2)) {
744  EpdCopy(&epd2, epd1);
745  }
746  return;
747  }
748 
749  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
750 
751  EpdConvert(value, &epd2);
752  if (epd1->exponent > epd2.exponent) {
753  diff = epd1->exponent - epd2.exponent;
754  if (diff <= EPD_MAX_BIN)
755  tmp = epd1->type.value - epd2.type.value / pow((double)2.0, (double)diff);
756  else
757  tmp = epd1->type.value;
758  exponent = epd1->exponent;
759  } else if (epd1->exponent < epd2.exponent) {
760  diff = epd2.exponent - epd1->exponent;
761  if (diff <= EPD_MAX_BIN)
762  tmp = epd1->type.value / pow((double)2.0, (double)diff) - epd2.type.value;
763  else
764  tmp = epd2.type.value * (double)(-1.0);
765  exponent = epd2.exponent;
766  } else {
767  tmp = epd1->type.value - epd2.type.value;
768  exponent = epd1->exponent;
769  }
770  epd1->type.value = tmp;
771  epd1->exponent = exponent;
772  EpdNormalize(epd1);
773 }
774 
775 
776 /**Function********************************************************************
777 
778  Synopsis [Subtracts two arbitrary precision double values.]
779 
780  Description [Subtracts two arbitrary precision double values.]
781 
782  SideEffects []
783 
784  SeeAlso []
785 
786 ******************************************************************************/
787 void
788 EpdSubtract2(EpDouble *epd1, EpDouble *epd2)
789 {
790  double value;
791  int exponent, diff;
792 
793  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
794  EpdMakeNan(epd1);
795  return;
796  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
797  int sign;
798 
799  if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
800  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
801  if (sign == 0)
802  EpdMakeNan(epd1);
803  } else if (EpdIsInf(epd2)) {
804  EpdCopy(epd2, epd1);
805  }
806  return;
807  }
808 
809  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
810  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
811 
812  if (epd1->exponent > epd2->exponent) {
813  diff = epd1->exponent - epd2->exponent;
814  if (diff <= EPD_MAX_BIN) {
815  value = epd1->type.value -
816  epd2->type.value / pow((double)2.0, (double)diff);
817  } else
818  value = epd1->type.value;
819  exponent = epd1->exponent;
820  } else if (epd1->exponent < epd2->exponent) {
821  diff = epd2->exponent - epd1->exponent;
822  if (diff <= EPD_MAX_BIN) {
823  value = epd1->type.value / pow((double)2.0, (double)diff) -
824  epd2->type.value;
825  } else
826  value = epd2->type.value * (double)(-1.0);
827  exponent = epd2->exponent;
828  } else {
829  value = epd1->type.value - epd2->type.value;
830  exponent = epd1->exponent;
831  }
832  epd1->type.value = value;
833  epd1->exponent = exponent;
834  EpdNormalize(epd1);
835 }
836 
837 
838 /**Function********************************************************************
839 
840  Synopsis [Subtracts two arbitrary precision double values.]
841 
842  Description [Subtracts two arbitrary precision double values.]
843 
844  SideEffects []
845 
846  SeeAlso []
847 
848 ******************************************************************************/
849 void
850 EpdSubtract3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
851 {
852  double value;
853  int exponent, diff;
854 
855  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
856  EpdMakeNan(epd3);
857  return;
858  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
859  int sign;
860 
861  if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
862  sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
863  if (sign == 0)
864  EpdCopy(epd1, epd3);
865  else
866  EpdMakeNan(epd3);
867  } else if (EpdIsInf(epd1)) {
868  EpdCopy(epd1, epd1);
869  } else {
870  sign = epd2->type.bits.sign ^ 0x1;
871  EpdMakeInf(epd3, sign);
872  }
873  return;
874  }
875 
876  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
877  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
878 
879  if (epd1->exponent > epd2->exponent) {
880  diff = epd1->exponent - epd2->exponent;
881  if (diff <= EPD_MAX_BIN) {
882  value = epd1->type.value -
883  epd2->type.value / pow((double)2.0, (double)diff);
884  } else
885  value = epd1->type.value;
886  exponent = epd1->exponent;
887  } else if (epd1->exponent < epd2->exponent) {
888  diff = epd2->exponent - epd1->exponent;
889  if (diff <= EPD_MAX_BIN) {
890  value = epd1->type.value / pow((double)2.0, (double)diff) -
891  epd2->type.value;
892  } else
893  value = epd2->type.value * (double)(-1.0);
894  exponent = epd2->exponent;
895  } else {
896  value = epd1->type.value - epd2->type.value;
897  exponent = epd1->exponent;
898  }
899  epd3->type.value = value;
900  epd3->exponent = exponent;
901  EpdNormalize(epd3);
902 }
903 
904 
905 /**Function********************************************************************
906 
907  Synopsis [Computes arbitrary precision pow of base 2.]
908 
909  Description [Computes arbitrary precision pow of base 2.]
910 
911  SideEffects []
912 
913  SeeAlso []
914 
915 ******************************************************************************/
916 void
917 EpdPow2(int n, EpDouble *epd)
918 {
919  if (n <= EPD_MAX_BIN) {
920  EpdConvert(pow((double)2.0, (double)n), epd);
921  } else {
922  EpDouble epd1, epd2;
923  int n1, n2;
924 
925  n1 = n / 2;
926  n2 = n - n1;
927  EpdPow2(n1, &epd1);
928  EpdPow2(n2, &epd2);
929  EpdMultiply3(&epd1, &epd2, epd);
930  }
931 }
932 
933 
934 /**Function********************************************************************
935 
936  Synopsis [Computes arbitrary precision pow of base 2.]
937 
938  Description [Computes arbitrary precision pow of base 2.]
939 
940  SideEffects []
941 
942  SeeAlso []
943 
944 ******************************************************************************/
945 void
946 EpdPow2Decimal(int n, EpDouble *epd)
947 {
948  if (n <= EPD_MAX_BIN) {
949  epd->type.value = pow((double)2.0, (double)n);
950  epd->exponent = 0;
951  EpdNormalizeDecimal(epd);
952  } else {
953  EpDouble epd1, epd2;
954  int n1, n2;
955 
956  n1 = n / 2;
957  n2 = n - n1;
958  EpdPow2Decimal(n1, &epd1);
959  EpdPow2Decimal(n2, &epd2);
960  EpdMultiply3Decimal(&epd1, &epd2, epd);
961  }
962 }
963 
964 
965 /**Function********************************************************************
966 
967  Synopsis [Normalize an arbitrary precision double value.]
968 
969  Description [Normalize an arbitrary precision double value.]
970 
971  SideEffects []
972 
973  SeeAlso []
974 
975 ******************************************************************************/
976 void
977 EpdNormalize(EpDouble *epd)
978 {
979  int exponent;
980 
981  if (IsNanOrInfDouble(epd->type.value)) {
982  epd->exponent = 0;
983  return;
984  }
985 
986  exponent = EpdGetExponent(epd->type.value);
987  if (exponent == EPD_MAX_BIN)
988  return;
989  exponent -= EPD_MAX_BIN;
990  epd->type.bits.exponent = EPD_MAX_BIN;
991  epd->exponent += exponent;
992 }
993 
994 
995 /**Function********************************************************************
996 
997  Synopsis [Normalize an arbitrary precision double value.]
998 
999  Description [Normalize an arbitrary precision double value.]
1000 
1001  SideEffects []
1002 
1003  SeeAlso []
1004 
1005 ******************************************************************************/
1006 void
1007 EpdNormalizeDecimal(EpDouble *epd)
1008 {
1009  int exponent;
1010 
1011  if (IsNanOrInfDouble(epd->type.value)) {
1012  epd->exponent = 0;
1013  return;
1014  }
1015 
1016  exponent = EpdGetExponentDecimal(epd->type.value);
1017  epd->type.value /= pow((double)10.0, (double)exponent);
1018  epd->exponent += exponent;
1019 }
1020 
1021 
1022 /**Function********************************************************************
1023 
1024  Synopsis [Returns value and decimal exponent of EpDouble.]
1025 
1026  Description [Returns value and decimal exponent of EpDouble.]
1027 
1028  SideEffects []
1029 
1030  SeeAlso []
1031 
1032 ******************************************************************************/
1033 void
1034 EpdGetValueAndDecimalExponent(EpDouble *epd, double *value, int *exponent)
1035 {
1036  EpDouble epd1, epd2;
1037 
1038  if (EpdIsNanOrInf(epd))
1039  return;
1040 
1041  if (EpdIsZero(epd)) {
1042  *value = 0.0;
1043  *exponent = 0;
1044  return;
1045  }
1046 
1047  epd1.type.value = epd->type.value;
1048  epd1.exponent = 0;
1049  EpdPow2Decimal(epd->exponent, &epd2);
1050  EpdMultiply2Decimal(&epd1, &epd2);
1051 
1052  *value = epd1.type.value;
1053  *exponent = epd1.exponent;
1054 }
1055 
1056 /**Function********************************************************************
1057 
1058  Synopsis [Returns the exponent value of a double.]
1059 
1060  Description [Returns the exponent value of a double.]
1061 
1062  SideEffects []
1063 
1064  SeeAlso []
1065 
1066 ******************************************************************************/
1067 int
1068 EpdGetExponent(double value)
1069 {
1070  int exponent;
1071  EpDouble epd;
1072 
1073  epd.type.value = value;
1074  exponent = epd.type.bits.exponent;
1075  return(exponent);
1076 }
1077 
1078 
1079 /**Function********************************************************************
1080 
1081  Synopsis [Returns the decimal exponent value of a double.]
1082 
1083  Description [Returns the decimal exponent value of a double.]
1084 
1085  SideEffects []
1086 
1087  SeeAlso []
1088 
1089 ******************************************************************************/
1090 int
1091 EpdGetExponentDecimal(double value)
1092 {
1093  char *pos, str[24];
1094  int exponent;
1095 
1096  sprintf(str, "%E", value);
1097  pos = strstr(str, "E");
1098  sscanf(pos, "E%d", &exponent);
1099  return(exponent);
1100 }
1101 
1102 
1103 /**Function********************************************************************
1104 
1105  Synopsis [Makes EpDouble Inf.]
1106 
1107  Description [Makes EpDouble Inf.]
1108 
1109  SideEffects []
1110 
1111  SeeAlso []
1112 
1113 ******************************************************************************/
1114 void
1115 EpdMakeInf(EpDouble *epd, int sign)
1116 {
1117  epd->type.bits.mantissa1 = 0;
1118  epd->type.bits.mantissa0 = 0;
1119  epd->type.bits.exponent = EPD_EXP_INF;
1120  epd->type.bits.sign = sign;
1121  epd->exponent = 0;
1122 }
1123 
1124 
1125 /**Function********************************************************************
1126 
1127  Synopsis [Makes EpDouble Zero.]
1128 
1129  Description [Makes EpDouble Zero.]
1130 
1131  SideEffects []
1132 
1133  SeeAlso []
1134 
1135 ******************************************************************************/
1136 void
1137 EpdMakeZero(EpDouble *epd, int sign)
1138 {
1139  epd->type.bits.mantissa1 = 0;
1140  epd->type.bits.mantissa0 = 0;
1141  epd->type.bits.exponent = 0;
1142  epd->type.bits.sign = sign;
1143  epd->exponent = 0;
1144 }
1145 
1146 
1147 /**Function********************************************************************
1148 
1149  Synopsis [Makes EpDouble NaN.]
1150 
1151  Description [Makes EpDouble NaN.]
1152 
1153  SideEffects []
1154 
1155  SeeAlso []
1156 
1157 ******************************************************************************/
1158 void
1159 EpdMakeNan(EpDouble *epd)
1160 {
1161  epd->type.nan.mantissa1 = 0;
1162  epd->type.nan.mantissa0 = 0;
1163  epd->type.nan.quiet_bit = 1;
1164  epd->type.nan.exponent = EPD_EXP_INF;
1165  epd->type.nan.sign = 1;
1166  epd->exponent = 0;
1167 }
1168 
1169 
1170 /**Function********************************************************************
1171 
1172  Synopsis [Copies a EpDouble struct.]
1173 
1174  Description [Copies a EpDouble struct.]
1175 
1176  SideEffects []
1177 
1178  SeeAlso []
1179 
1180 ******************************************************************************/
1181 void
1182 EpdCopy(EpDouble *from, EpDouble *to)
1183 {
1184  to->type.value = from->type.value;
1185  to->exponent = from->exponent;
1186 }
1187 
1188 
1189 /**Function********************************************************************
1190 
1191  Synopsis [Checks whether the value is Inf.]
1192 
1193  Description [Checks whether the value is Inf.]
1194 
1195  SideEffects []
1196 
1197  SeeAlso []
1198 
1199 ******************************************************************************/
1200 int
1201 EpdIsInf(EpDouble *epd)
1202 {
1203  return(IsInfDouble(epd->type.value));
1204 }
1205 
1206 
1207 /**Function********************************************************************
1208 
1209  Synopsis [Checks whether the value is Zero.]
1210 
1211  Description [Checks whether the value is Zero.]
1212 
1213  SideEffects []
1214 
1215  SeeAlso []
1216 
1217 ******************************************************************************/
1218 int
1219 EpdIsZero(EpDouble *epd)
1220 {
1221  if (epd->type.value == 0.0)
1222  return(1);
1223  else
1224  return(0);
1225 }
1226 
1227 
1228 /**Function********************************************************************
1229 
1230  Synopsis [Checks whether the value is NaN.]
1231 
1232  Description [Checks whether the value is NaN.]
1233 
1234  SideEffects []
1235 
1236  SeeAlso []
1237 
1238 ******************************************************************************/
1239 int
1240 EpdIsNan(EpDouble *epd)
1241 {
1242  return(IsNanDouble(epd->type.value));
1243 }
1244 
1245 
1246 /**Function********************************************************************
1247 
1248  Synopsis [Checks whether the value is NaN or Inf.]
1249 
1250  Description [Checks whether the value is NaN or Inf.]
1251 
1252  SideEffects []
1253 
1254  SeeAlso []
1255 
1256 ******************************************************************************/
1257 int
1258 EpdIsNanOrInf(EpDouble *epd)
1259 {
1260  return(IsNanOrInfDouble(epd->type.value));
1261 }
1262 
1263 
1264 /**Function********************************************************************
1265 
1266  Synopsis [Checks whether the value is Inf.]
1267 
1268  Description [Checks whether the value is Inf.]
1269 
1270  SideEffects []
1271 
1272  SeeAlso []
1273 
1274 ******************************************************************************/
1275 int
1276 IsInfDouble(double value)
1277 {
1278  EpType val;
1279 
1280  val.value = value;
1281  if (val.bits.exponent == EPD_EXP_INF &&
1282  val.bits.mantissa0 == 0 &&
1283  val.bits.mantissa1 == 0) {
1284  if (val.bits.sign == 0)
1285  return(1);
1286  else
1287  return(-1);
1288  }
1289  return(0);
1290 }
1291 
1292 
1293 /**Function********************************************************************
1294 
1295  Synopsis [Checks whether the value is NaN.]
1296 
1297  Description [Checks whether the value is NaN.]
1298 
1299  SideEffects []
1300 
1301  SeeAlso []
1302 
1303 ******************************************************************************/
1304 int
1305 IsNanDouble(double value)
1306 {
1307  EpType val;
1308 
1309  val.value = value;
1310  if (val.nan.exponent == EPD_EXP_INF &&
1311  val.nan.sign == 1 &&
1312  val.nan.quiet_bit == 1 &&
1313  val.nan.mantissa0 == 0 &&
1314  val.nan.mantissa1 == 0) {
1315  return(1);
1316  }
1317  return(0);
1318 }
1319 
1320 
1321 /**Function********************************************************************
1322 
1323  Synopsis [Checks whether the value is NaN or Inf.]
1324 
1325  Description [Checks whether the value is NaN or Inf.]
1326 
1327  SideEffects []
1328 
1329  SeeAlso []
1330 
1331 ******************************************************************************/
1332 int
1333 IsNanOrInfDouble(double value)
1334 {
1335  EpType val;
1336 
1337  val.value = value;
1338  if (val.nan.exponent == EPD_EXP_INF &&
1339  val.nan.mantissa0 == 0 &&
1340  val.nan.mantissa1 == 0 &&
1341  (val.nan.sign == 1 || val.nan.quiet_bit == 0)) {
1342  return(1);
1343  }
1344  return(0);
1345 }
1346 
1347 ABC_NAMESPACE_IMPL_END
IsNanDouble
int IsNanDouble(double value)
Definition: epd.c:1305
EpTypeUnion::value
double value
Definition: epd.h:130
EpdSubtract
void EpdSubtract(EpDouble *epd1, double value)
Definition: epd.c:726
EpdGetExponentDecimal
int EpdGetExponentDecimal(double value)
Definition: epd.c:1091
EpdCmp
int EpdCmp(const char *key1, const char *key2)
Definition: epd.c:93
EpdConvert
void EpdConvert(double value, EpDouble *epd)
Definition: epd.c:185
IeeeDoubleStruct::sign
unsigned int sign
Definition: epd.h:89
EpdDivide3
void EpdDivide3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
Definition: epd.c:490
EpdGetValueAndDecimalExponent
void EpdGetValueAndDecimalExponent(EpDouble *epd, double *value, int *exponent)
Definition: epd.c:1034
EpdCopy
void EpdCopy(EpDouble *from, EpDouble *to)
Definition: epd.c:1182
IeeeDoubleStruct::mantissa1
unsigned int mantissa1
Definition: epd.h:86
pos
bool pos
Definition: globals.c:30
ABC_ALLOC
#define ABC_ALLOC(type, num)
Definition: abc_global.h:229
EpTypeUnion::nan
struct IeeeNanStruct nan
Definition: epd.h:132
Minisat::sign
bool sign(Lit p)
Definition: SolverTypes.h:61
EpdIsNan
int EpdIsNan(EpDouble *epd)
Definition: epd.c:1240
IeeeDoubleStruct::mantissa0
unsigned int mantissa0
Definition: epd.h:87
EpdPow2
void EpdPow2(int n, EpDouble *epd)
Definition: epd.c:917
EpdSubtract3
void EpdSubtract3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
Definition: epd.c:850
EpdNormalize
void EpdNormalize(EpDouble *epd)
Definition: epd.c:977
util_hack.h
IsInfDouble
int IsInfDouble(double value)
Definition: epd.c:1276
EpdNormalizeDecimal
void EpdNormalizeDecimal(EpDouble *epd)
Definition: epd.c:1007
EPD_EXP_INF
#define EPD_EXP_INF
Definition: epd.h:62
IeeeNanStruct::exponent
unsigned int exponent
Definition: epd.h:115
EpDoubleStruct::type
union EpTypeUnion type
Definition: epd.h:136
IeeeNanStruct::sign
unsigned int sign
Definition: epd.h:116
EpdAdd3
void EpdAdd3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
Definition: epd.c:660
EpdDivide
void EpdDivide(EpDouble *epd1, double value)
Definition: epd.c:386
IeeeNanStruct::mantissa0
unsigned int mantissa0
Definition: epd.h:113
strstr
char * strstr()
IeeeNanStruct::quiet_bit
unsigned int quiet_bit
Definition: epd.h:114
EpdPow2Decimal
void EpdPow2Decimal(int n, EpDouble *epd)
Definition: epd.c:946
ABC_NAMESPACE_IMPL_END
#define ABC_NAMESPACE_IMPL_END
Definition: abc_global.h:108
EpdDivide2
void EpdDivide2(EpDouble *epd1, EpDouble *epd2)
Definition: epd.c:439
EpdMakeZero
void EpdMakeZero(EpDouble *epd, int sign)
Definition: epd.c:1137
EpdAdd
void EpdAdd(EpDouble *epd1, double value)
Definition: epd.c:536
sprintf
char * sprintf()
EpdAlloc
ABC_NAMESPACE_IMPL_START EpDouble * EpdAlloc(void)
Definition: epd.c:72
EpTypeUnion::bits
struct IeeeDoubleStruct bits
Definition: epd.h:131
EpdMakeInf
void EpdMakeInf(EpDouble *epd, int sign)
Definition: epd.c:1115
ABC_NAMESPACE_IMPL_START
#define ABC_NAMESPACE_IMPL_START
Definition: abc_global.h:107
IeeeDoubleStruct::exponent
unsigned int exponent
Definition: epd.h:88
EpdIsZero
int EpdIsZero(EpDouble *epd)
Definition: epd.c:1219
EpdMakeNan
void EpdMakeNan(EpDouble *epd)
Definition: epd.c:1159
IeeeNanStruct::mantissa1
unsigned int mantissa1
Definition: epd.h:112
EpdMultiply2Decimal
void EpdMultiply2Decimal(EpDouble *epd1, EpDouble *epd2)
Definition: epd.c:285
ABC_FREE
#define ABC_FREE(obj)
Definition: abc_global.h:232
IsNanOrInfDouble
int IsNanOrInfDouble(double value)
Definition: epd.c:1333
EpdMultiply2
void EpdMultiply2(EpDouble *epd1, EpDouble *epd2)
Definition: epd.c:246
EpdMultiply3Decimal
void EpdMultiply3Decimal(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
Definition: epd.c:355
EpdMultiply
void EpdMultiply(EpDouble *epd1, double value)
Definition: epd.c:205
value
int value
Definition: misc/espresso/main.h:90
EpdFree
void EpdFree(EpDouble *epd)
Definition: epd.c:117
assert
#define assert(ex)
Definition: util_old.h:213
EpDoubleStruct::exponent
int exponent
Definition: epd.h:137
EpdIsInf
int EpdIsInf(EpDouble *epd)
Definition: epd.c:1201
EpdGetString
void EpdGetString(EpDouble *epd, char *str)
Definition: epd.c:135
EpdMultiply3
void EpdMultiply3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
Definition: epd.c:321
EpdAdd2
void EpdAdd2(EpDouble *epd1, EpDouble *epd2)
Definition: epd.c:598
EPD_MAX_BIN
#define EPD_MAX_BIN
Definition: epd.h:60
EpdGetExponent
int EpdGetExponent(double value)
Definition: epd.c:1068
EpdSubtract2
void EpdSubtract2(EpDouble *epd1, EpDouble *epd2)
Definition: epd.c:788
EpdIsNanOrInf
int EpdIsNanOrInf(EpDouble *epd)
Definition: epd.c:1258