Version abc50729

3 files changed, 1475 insertions, 0 deletions

diff --git a/src/bdd/epd/epd.c b/src/bdd/epd/epd.c
new file mode 100644
index 00000000..a843b986
--- /dev/null
+++ b/src/bdd/epd/epd.c
@@ -0,0 +1,1314 @@
+/**CFile***********************************************************************
+
+  FileName    [epd.c]
+
+  PackageName [epd]
+
+  Synopsis    [Arithmetic functions with extended double precision.]
+
+  Description []
+
+  SeeAlso     []
+
+  Author      [In-Ho Moon]
+
+  Copyright [ This file was created at the University of Colorado at
+  Boulder.  The University of Colorado at Boulder makes no warranty
+  about the suitability of this software for any purpose.  It is
+  presented on an AS IS basis.]
+
+  Revision    [$Id: epd.c,v 1.1.1.1 2003/02/24 22:23:57 wjiang Exp $]
+
+******************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "util.h"
+#include "epd.h"
+
+
+/**Function********************************************************************
+
+  Synopsis    [Allocates an EpDouble struct.]
+
+  Description [Allocates an EpDouble struct.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+EpDouble *
+EpdAlloc()
+{
+  EpDouble    *epd;
+
+  epd = ALLOC(EpDouble, 1);
+  return(epd);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Compares two EpDouble struct.]
+
+  Description [Compares two EpDouble struct.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+int
+EpdCmp(const char *key1, const char *key2)
+{
+  EpDouble *epd1 = (EpDouble *) key1;
+  EpDouble *epd2 = (EpDouble *) key2;
+  if (epd1->type.value != epd2->type.value ||
+      epd1->exponent != epd2->exponent) {
+    return(1);
+  }
+  return(0);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Frees an EpDouble struct.]
+
+  Description [Frees an EpDouble struct.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdFree(EpDouble *epd)
+{
+  FREE(epd);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Multiplies two arbitrary precision double values.]
+
+  Description [Multiplies two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdGetString(EpDouble *epd, char *str)
+{
+  double    value;
+  int        exponent;
+  char        *pos;
+
+  if (IsNanDouble(epd->type.value)) {
+    sprintf(str, "NaN");
+    return;
+  } else if (IsInfDouble(epd->type.value)) {
+    if (epd->type.bits.sign == 1)
+      sprintf(str, "-Inf");
+    else
+      sprintf(str, "Inf");
+    return;
+  }
+
+  assert(epd->type.bits.exponent == EPD_MAX_BIN ||
+     epd->type.bits.exponent == 0);
+
+  EpdGetValueAndDecimalExponent(epd, &value, &exponent);
+  sprintf(str, "%e", value);
+  pos = strstr(str, "e");
+  if (exponent >= 0) {
+    if (exponent < 10)
+      sprintf(pos + 1, "+0%d", exponent);
+    else
+      sprintf(pos + 1, "+%d", exponent);
+  } else {
+    exponent *= -1;
+    if (exponent < 10)
+      sprintf(pos + 1, "-0%d", exponent);
+    else
+      sprintf(pos + 1, "-%d", exponent);
+  }
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Converts double to EpDouble struct.]
+
+  Description [Converts double to EpDouble struct.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdConvert(double value, EpDouble *epd)
+{
+  epd->type.value = value;
+  epd->exponent = 0;
+  EpdNormalize(epd);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Multiplies two arbitrary precision double values.]
+
+  Description [Multiplies two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdMultiply(EpDouble *epd1, double value)
+{
+  EpDouble    epd2;
+  double    tmp;
+  int        exponent;
+
+  if (EpdIsNan(epd1) || IsNanDouble(value)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
+    int    sign;
+
+    EpdConvert(value, &epd2);
+    sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
+    EpdMakeInf(epd1, sign);
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+
+  EpdConvert(value, &epd2);
+  tmp = epd1->type.value * epd2.type.value;
+  exponent = epd1->exponent + epd2.exponent;
+  epd1->type.value = tmp;
+  epd1->exponent = exponent;
+  EpdNormalize(epd1);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Multiplies two arbitrary precision double values.]
+
+  Description [Multiplies two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdMultiply2(EpDouble *epd1, EpDouble *epd2)
+{
+  double    value;
+  int        exponent;
+
+  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
+    int    sign;
+
+    sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+    EpdMakeInf(epd1, sign);
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
+
+  value = epd1->type.value * epd2->type.value;
+  exponent = epd1->exponent + epd2->exponent;
+  epd1->type.value = value;
+  epd1->exponent = exponent;
+  EpdNormalize(epd1);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Multiplies two arbitrary precision double values.]
+
+  Description [Multiplies two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdMultiply2Decimal(EpDouble *epd1, EpDouble *epd2)
+{
+  double    value;
+  int        exponent;
+
+  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
+    int    sign;
+
+    sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+    EpdMakeInf(epd1, sign);
+    return;
+  }
+
+  value = epd1->type.value * epd2->type.value;
+  exponent = epd1->exponent + epd2->exponent;
+  epd1->type.value = value;
+  epd1->exponent = exponent;
+  EpdNormalizeDecimal(epd1);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Multiplies two arbitrary precision double values.]
+
+  Description [Multiplies two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdMultiply3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
+{
+  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
+    int    sign;
+
+    sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+    EpdMakeInf(epd3, sign);
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
+
+  epd3->type.value = epd1->type.value * epd2->type.value;
+  epd3->exponent = epd1->exponent + epd2->exponent;
+  EpdNormalize(epd3);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Multiplies two arbitrary precision double values.]
+
+  Description [Multiplies two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdMultiply3Decimal(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
+{
+  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
+    int    sign;
+
+    sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+    EpdMakeInf(epd3, sign);
+    return;
+  }
+
+  epd3->type.value = epd1->type.value * epd2->type.value;
+  epd3->exponent = epd1->exponent + epd2->exponent;
+  EpdNormalizeDecimal(epd3);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Divides two arbitrary precision double values.]
+
+  Description [Divides two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdDivide(EpDouble *epd1, double value)
+{
+  EpDouble    epd2;
+  double    tmp;
+  int        exponent;
+
+  if (EpdIsNan(epd1) || IsNanDouble(value)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
+    int    sign;
+
+    EpdConvert(value, &epd2);
+    if (EpdIsInf(epd1) && IsInfDouble(value)) {
+      EpdMakeNan(epd1);
+    } else if (EpdIsInf(epd1)) {
+      sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
+      EpdMakeInf(epd1, sign);
+    } else {
+      sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
+      EpdMakeZero(epd1, sign);
+    }
+    return;
+  }
+
+  if (value == 0.0) {
+    EpdMakeNan(epd1);
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+
+  EpdConvert(value, &epd2);
+  tmp = epd1->type.value / epd2.type.value;
+  exponent = epd1->exponent - epd2.exponent;
+  epd1->type.value = tmp;
+  epd1->exponent = exponent;
+  EpdNormalize(epd1);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Divides two arbitrary precision double values.]
+
+  Description [Divides two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdDivide2(EpDouble *epd1, EpDouble *epd2)
+{
+  double    value;
+  int        exponent;
+
+  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
+    int    sign;
+
+    if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
+      EpdMakeNan(epd1);
+    } else if (EpdIsInf(epd1)) {
+      sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+      EpdMakeInf(epd1, sign);
+    } else {
+      sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+      EpdMakeZero(epd1, sign);
+    }
+    return;
+  }
+
+  if (epd2->type.value == 0.0) {
+    EpdMakeNan(epd1);
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
+
+  value = epd1->type.value / epd2->type.value;
+  exponent = epd1->exponent - epd2->exponent;
+  epd1->type.value = value;
+  epd1->exponent = exponent;
+  EpdNormalize(epd1);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Divides two arbitrary precision double values.]
+
+  Description [Divides two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdDivide3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
+{
+  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
+    EpdMakeNan(epd3);
+    return;
+  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
+    int    sign;
+
+    if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
+      EpdMakeNan(epd3);
+    } else if (EpdIsInf(epd1)) {
+      sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+      EpdMakeInf(epd3, sign);
+    } else {
+      sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+      EpdMakeZero(epd3, sign);
+    }
+    return;
+  }
+
+  if (epd2->type.value == 0.0) {
+    EpdMakeNan(epd3);
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
+
+  epd3->type.value = epd1->type.value / epd2->type.value;
+  epd3->exponent = epd1->exponent - epd2->exponent;
+  EpdNormalize(epd3);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Adds two arbitrary precision double values.]
+
+  Description [Adds two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdAdd(EpDouble *epd1, double value)
+{
+  EpDouble    epd2;
+  double    tmp;
+  int        exponent, diff;
+
+  if (EpdIsNan(epd1) || IsNanDouble(value)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
+    int    sign;
+
+    EpdConvert(value, &epd2);
+    if (EpdIsInf(epd1) && IsInfDouble(value)) {
+      sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
+      if (sign == 1)
+    EpdMakeNan(epd1);
+    } else if (EpdIsInf(&epd2)) {
+      EpdCopy(&epd2, epd1);
+    }
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+
+  EpdConvert(value, &epd2);
+  if (epd1->exponent > epd2.exponent) {
+    diff = epd1->exponent - epd2.exponent;
+    if (diff <= EPD_MAX_BIN)
+      tmp = epd1->type.value + epd2.type.value / pow((double)2.0, (double)diff);
+    else
+      tmp = epd1->type.value;
+    exponent = epd1->exponent;
+  } else if (epd1->exponent < epd2.exponent) {
+    diff = epd2.exponent - epd1->exponent;
+    if (diff <= EPD_MAX_BIN)
+      tmp = epd1->type.value / pow((double)2.0, (double)diff) + epd2.type.value;
+    else
+      tmp = epd2.type.value;
+    exponent = epd2.exponent;
+  } else {
+    tmp = epd1->type.value + epd2.type.value;
+    exponent = epd1->exponent;
+  }
+  epd1->type.value = tmp;
+  epd1->exponent = exponent;
+  EpdNormalize(epd1);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Adds two arbitrary precision double values.]
+
+  Description [Adds two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdAdd2(EpDouble *epd1, EpDouble *epd2)
+{
+  double    value;
+  int        exponent, diff;
+
+  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
+    int    sign;
+
+    if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
+      sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+      if (sign == 1)
+    EpdMakeNan(epd1);
+    } else if (EpdIsInf(epd2)) {
+      EpdCopy(epd2, epd1);
+    }
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
+
+  if (epd1->exponent > epd2->exponent) {
+    diff = epd1->exponent - epd2->exponent;
+    if (diff <= EPD_MAX_BIN) {
+      value = epd1->type.value +
+        epd2->type.value / pow((double)2.0, (double)diff);
+    } else
+      value = epd1->type.value;
+    exponent = epd1->exponent;
+  } else if (epd1->exponent < epd2->exponent) {
+    diff = epd2->exponent - epd1->exponent;
+    if (diff <= EPD_MAX_BIN) {
+      value = epd1->type.value / pow((double)2.0, (double)diff) +
+        epd2->type.value;
+    } else
+      value = epd2->type.value;
+    exponent = epd2->exponent;
+  } else {
+    value = epd1->type.value + epd2->type.value;
+    exponent = epd1->exponent;
+  }
+  epd1->type.value = value;
+  epd1->exponent = exponent;
+  EpdNormalize(epd1);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Adds two arbitrary precision double values.]
+
+  Description [Adds two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdAdd3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
+{
+  double    value;
+  int        exponent, diff;
+
+  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
+    EpdMakeNan(epd3);
+    return;
+  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
+    int    sign;
+
+    if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
+      sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+      if (sign == 1)
+    EpdMakeNan(epd3);
+      else
+    EpdCopy(epd1, epd3);
+    } else if (EpdIsInf(epd1)) {
+      EpdCopy(epd1, epd3);
+    } else {
+      EpdCopy(epd2, epd3);
+    }
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
+
+  if (epd1->exponent > epd2->exponent) {
+    diff = epd1->exponent - epd2->exponent;
+    if (diff <= EPD_MAX_BIN) {
+      value = epd1->type.value +
+        epd2->type.value / pow((double)2.0, (double)diff);
+    } else
+      value = epd1->type.value;
+    exponent = epd1->exponent;
+  } else if (epd1->exponent < epd2->exponent) {
+    diff = epd2->exponent - epd1->exponent;
+    if (diff <= EPD_MAX_BIN) {
+      value = epd1->type.value / pow((double)2.0, (double)diff) +
+        epd2->type.value;
+    } else
+      value = epd2->type.value;
+    exponent = epd2->exponent;
+  } else {
+    value = epd1->type.value + epd2->type.value;
+    exponent = epd1->exponent;
+  }
+  epd3->type.value = value;
+  epd3->exponent = exponent;
+  EpdNormalize(epd3);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Subtracts two arbitrary precision double values.]
+
+  Description [Subtracts two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdSubtract(EpDouble *epd1, double value)
+{
+  EpDouble    epd2;
+  double    tmp;
+  int        exponent, diff;
+
+  if (EpdIsNan(epd1) || IsNanDouble(value)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
+    int    sign;
+
+    EpdConvert(value, &epd2);
+    if (EpdIsInf(epd1) && IsInfDouble(value)) {
+      sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
+      if (sign == 0)
+    EpdMakeNan(epd1);
+    } else if (EpdIsInf(&epd2)) {
+      EpdCopy(&epd2, epd1);
+    }
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+
+  EpdConvert(value, &epd2);
+  if (epd1->exponent > epd2.exponent) {
+    diff = epd1->exponent - epd2.exponent;
+    if (diff <= EPD_MAX_BIN)
+      tmp = epd1->type.value - epd2.type.value / pow((double)2.0, (double)diff);
+    else
+      tmp = epd1->type.value;
+    exponent = epd1->exponent;
+  } else if (epd1->exponent < epd2.exponent) {
+    diff = epd2.exponent - epd1->exponent;
+    if (diff <= EPD_MAX_BIN)
+      tmp = epd1->type.value / pow((double)2.0, (double)diff) - epd2.type.value;
+    else
+      tmp = epd2.type.value * (double)(-1.0);
+    exponent = epd2.exponent;
+  } else {
+    tmp = epd1->type.value - epd2.type.value;
+    exponent = epd1->exponent;
+  }
+  epd1->type.value = tmp;
+  epd1->exponent = exponent;
+  EpdNormalize(epd1);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Subtracts two arbitrary precision double values.]
+
+  Description [Subtracts two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdSubtract2(EpDouble *epd1, EpDouble *epd2)
+{
+  double    value;
+  int        exponent, diff;
+
+  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
+    EpdMakeNan(epd1);
+    return;
+  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
+    int    sign;
+
+    if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
+      sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+      if (sign == 0)
+    EpdMakeNan(epd1);
+    } else if (EpdIsInf(epd2)) {
+      EpdCopy(epd2, epd1);
+    }
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
+
+  if (epd1->exponent > epd2->exponent) {
+    diff = epd1->exponent - epd2->exponent;
+    if (diff <= EPD_MAX_BIN) {
+      value = epd1->type.value -
+        epd2->type.value / pow((double)2.0, (double)diff);
+    } else
+      value = epd1->type.value;
+    exponent = epd1->exponent;
+  } else if (epd1->exponent < epd2->exponent) {
+    diff = epd2->exponent - epd1->exponent;
+    if (diff <= EPD_MAX_BIN) {
+      value = epd1->type.value / pow((double)2.0, (double)diff) -
+        epd2->type.value;
+    } else
+      value = epd2->type.value * (double)(-1.0);
+    exponent = epd2->exponent;
+  } else {
+    value = epd1->type.value - epd2->type.value;
+    exponent = epd1->exponent;
+  }
+  epd1->type.value = value;
+  epd1->exponent = exponent;
+  EpdNormalize(epd1);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Subtracts two arbitrary precision double values.]
+
+  Description [Subtracts two arbitrary precision double values.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdSubtract3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3)
+{
+  double    value;
+  int        exponent, diff;
+
+  if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
+    EpdMakeNan(epd3);
+    return;
+  } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
+    int    sign;
+
+    if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
+      sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
+      if (sign == 0)
+    EpdCopy(epd1, epd3);
+      else
+    EpdMakeNan(epd3);
+    } else if (EpdIsInf(epd1)) {
+      EpdCopy(epd1, epd1);
+    } else {
+      sign = epd2->type.bits.sign ^ 0x1;
+      EpdMakeInf(epd3, sign);
+    }
+    return;
+  }
+
+  assert(epd1->type.bits.exponent == EPD_MAX_BIN);
+  assert(epd2->type.bits.exponent == EPD_MAX_BIN);
+
+  if (epd1->exponent > epd2->exponent) {
+    diff = epd1->exponent - epd2->exponent;
+    if (diff <= EPD_MAX_BIN) {
+      value = epd1->type.value -
+        epd2->type.value / pow((double)2.0, (double)diff);
+    } else
+      value = epd1->type.value;
+    exponent = epd1->exponent;
+  } else if (epd1->exponent < epd2->exponent) {
+    diff = epd2->exponent - epd1->exponent;
+    if (diff <= EPD_MAX_BIN) {
+      value = epd1->type.value / pow((double)2.0, (double)diff) -
+        epd2->type.value;
+    } else
+      value = epd2->type.value * (double)(-1.0);
+    exponent = epd2->exponent;
+  } else {
+    value = epd1->type.value - epd2->type.value;
+    exponent = epd1->exponent;
+  }
+  epd3->type.value = value;
+  epd3->exponent = exponent;
+  EpdNormalize(epd3);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Computes arbitrary precision pow of base 2.]
+
+  Description [Computes arbitrary precision pow of base 2.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdPow2(int n, EpDouble *epd)
+{
+  if (n <= EPD_MAX_BIN) {
+    EpdConvert(pow((double)2.0, (double)n), epd);
+  } else {
+    EpDouble    epd1, epd2;
+    int        n1, n2;
+
+    n1 = n / 2;
+    n2 = n - n1;
+    EpdPow2(n1, &epd1);
+    EpdPow2(n2, &epd2);
+    EpdMultiply3(&epd1, &epd2, epd);
+  }
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Computes arbitrary precision pow of base 2.]
+
+  Description [Computes arbitrary precision pow of base 2.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdPow2Decimal(int n, EpDouble *epd)
+{
+  if (n <= EPD_MAX_BIN) {
+    epd->type.value = pow((double)2.0, (double)n);
+    epd->exponent = 0;
+    EpdNormalizeDecimal(epd);
+  } else {
+    EpDouble    epd1, epd2;
+    int        n1, n2;
+
+    n1 = n / 2;
+    n2 = n - n1;
+    EpdPow2Decimal(n1, &epd1);
+    EpdPow2Decimal(n2, &epd2);
+    EpdMultiply3Decimal(&epd1, &epd2, epd);
+  }
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Normalize an arbitrary precision double value.]
+
+  Description [Normalize an arbitrary precision double value.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdNormalize(EpDouble *epd)
+{
+  int        exponent;
+
+  if (IsNanOrInfDouble(epd->type.value)) {
+    epd->exponent = 0;
+    return;
+  }
+
+  exponent = EpdGetExponent(epd->type.value);
+  if (exponent == EPD_MAX_BIN)
+    return;
+  exponent -= EPD_MAX_BIN;
+  epd->type.bits.exponent = EPD_MAX_BIN;
+  epd->exponent += exponent;
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Normalize an arbitrary precision double value.]
+
+  Description [Normalize an arbitrary precision double value.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdNormalizeDecimal(EpDouble *epd)
+{
+  int        exponent;
+
+  if (IsNanOrInfDouble(epd->type.value)) {
+    epd->exponent = 0;
+    return;
+  }
+
+  exponent = EpdGetExponentDecimal(epd->type.value);
+  epd->type.value /= pow((double)10.0, (double)exponent);
+  epd->exponent += exponent;
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Returns value and decimal exponent of EpDouble.]
+
+  Description [Returns value and decimal exponent of EpDouble.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdGetValueAndDecimalExponent(EpDouble *epd, double *value, int *exponent)
+{
+  EpDouble    epd1, epd2;
+
+  if (EpdIsNanOrInf(epd))
+    return;
+
+  if (EpdIsZero(epd)) {
+    *value = 0.0;
+    *exponent = 0;
+    return;
+  }
+
+  epd1.type.value = epd->type.value;
+  epd1.exponent = 0;
+  EpdPow2Decimal(epd->exponent, &epd2);
+  EpdMultiply2Decimal(&epd1, &epd2);
+
+  *value = epd1.type.value;
+  *exponent = epd1.exponent;
+}
+
+/**Function********************************************************************
+
+  Synopsis    [Returns the exponent value of a double.]
+
+  Description [Returns the exponent value of a double.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+int
+EpdGetExponent(double value)
+{
+  int        exponent;
+  EpDouble    epd;
+
+  epd.type.value = value;
+  exponent = epd.type.bits.exponent;
+  return(exponent);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Returns the decimal exponent value of a double.]
+
+  Description [Returns the decimal exponent value of a double.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+int
+EpdGetExponentDecimal(double value)
+{
+  char    *pos, str[24];
+  int    exponent;
+
+  sprintf(str, "%E", value);
+  pos = strstr(str, "E");
+  sscanf(pos, "E%d", &exponent);
+  return(exponent);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Makes EpDouble Inf.]
+
+  Description [Makes EpDouble Inf.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdMakeInf(EpDouble *epd, int sign)
+{
+  epd->type.bits.mantissa1 = 0;
+  epd->type.bits.mantissa0 = 0;
+  epd->type.bits.exponent = EPD_EXP_INF;
+  epd->type.bits.sign = sign;
+  epd->exponent = 0;
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Makes EpDouble Zero.]
+
+  Description [Makes EpDouble Zero.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdMakeZero(EpDouble *epd, int sign)
+{
+  epd->type.bits.mantissa1 = 0;
+  epd->type.bits.mantissa0 = 0;
+  epd->type.bits.exponent = 0;
+  epd->type.bits.sign = sign;
+  epd->exponent = 0;
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Makes EpDouble NaN.]
+
+  Description [Makes EpDouble NaN.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdMakeNan(EpDouble *epd)
+{
+  epd->type.nan.mantissa1 = 0;
+  epd->type.nan.mantissa0 = 0;
+  epd->type.nan.quiet_bit = 1;
+  epd->type.nan.exponent = EPD_EXP_INF;
+  epd->type.nan.sign = 1;
+  epd->exponent = 0;
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Copies a EpDouble struct.]
+
+  Description [Copies a EpDouble struct.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+void
+EpdCopy(EpDouble *from, EpDouble *to)
+{
+  to->type.value = from->type.value;
+  to->exponent = from->exponent;
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Checks whether the value is Inf.]
+
+  Description [Checks whether the value is Inf.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+int
+EpdIsInf(EpDouble *epd)
+{
+  return(IsInfDouble(epd->type.value));
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Checks whether the value is Zero.]
+
+  Description [Checks whether the value is Zero.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+int
+EpdIsZero(EpDouble *epd)
+{
+  if (epd->type.value == 0.0)
+    return(1);
+  else
+    return(0);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Checks whether the value is NaN.]
+
+  Description [Checks whether the value is NaN.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+int
+EpdIsNan(EpDouble *epd)
+{
+  return(IsNanDouble(epd->type.value));
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Checks whether the value is NaN or Inf.]
+
+  Description [Checks whether the value is NaN or Inf.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+int
+EpdIsNanOrInf(EpDouble *epd)
+{
+  return(IsNanOrInfDouble(epd->type.value));
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Checks whether the value is Inf.]
+
+  Description [Checks whether the value is Inf.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+int
+IsInfDouble(double value)
+{
+  IeeeDouble    *ptr = (IeeeDouble *)(&value);
+
+  if (ptr->exponent == EPD_EXP_INF &&
+      ptr->mantissa0 == 0 &&
+      ptr->mantissa1 == 0) {
+    if (ptr->sign == 0)
+      return(1);
+    else
+      return(-1);
+  }
+  return(0);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Checks whether the value is NaN.]
+
+  Description [Checks whether the value is NaN.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+int
+IsNanDouble(double value)
+{
+  IeeeNan    *ptr = (IeeeNan *)(&value);
+
+  if (ptr->exponent == EPD_EXP_INF &&
+      ptr->sign == 1 &&
+      ptr->quiet_bit == 1 &&
+      ptr->mantissa0 == 0 &&
+      ptr->mantissa1 == 0) {
+    return(1);
+  }
+  return(0);
+}
+
+
+/**Function********************************************************************
+
+  Synopsis    [Checks whether the value is NaN or Inf.]
+
+  Description [Checks whether the value is NaN or Inf.]
+
+  SideEffects []
+
+  SeeAlso     []
+
+******************************************************************************/
+int
+IsNanOrInfDouble(double value)
+{
+  IeeeNan    *ptr = (IeeeNan *)(&value);
+
+  if (ptr->exponent == EPD_EXP_INF &&
+      ptr->mantissa0 == 0 &&
+      ptr->mantissa1 == 0 &&
+      (ptr->sign == 1 || ptr->quiet_bit == 0)) {
+    return(1);
+  }
+  return(0);
+}
diff --git a/src/bdd/epd/epd.h b/src/bdd/epd/epd.h
new file mode 100644
index 00000000..66db80e3
--- /dev/null
+++ b/src/bdd/epd/epd.h
@@ -0,0 +1,160 @@
+/**CHeaderFile*****************************************************************
+
+  FileName    [epd.h]
+
+  PackageName [epd]
+
+  Synopsis    [The University of Colorado extended double precision package.]
+
+  Description [arithmetic functions with extended double precision.]
+
+  SeeAlso     []
+
+  Author      [In-Ho Moon]
+
+  Copyright [This file was created at the University of Colorado at
+  Boulder.  The University of Colorado at Boulder makes no warranty
+  about the suitability of this software for any purpose.  It is
+  presented on an AS IS basis.]
+
+  Revision    [$Id: epd.h,v 1.1.1.1 2003/02/24 22:23:57 wjiang Exp $]
+
+******************************************************************************/
+
+#ifndef _EPD
+#define _EPD
+
+
+/*---------------------------------------------------------------------------*/
+/* Constant declarations                                                     */
+/*---------------------------------------------------------------------------*/
+
+#define    EPD_MAX_BIN    1023
+#define    EPD_MAX_DEC    308
+#define    EPD_EXP_INF    0x7ff
+
+/*---------------------------------------------------------------------------*/
+/* Structure declarations                                                    */
+/*---------------------------------------------------------------------------*/
+
+/**Struct**********************************************************************
+
+  Synopsis    [IEEE double struct.]
+
+  Description [IEEE double struct.]
+
+  SeeAlso     []
+
+******************************************************************************/
+#ifdef    EPD_BIG_ENDIAN
+struct IeeeDoubleStruct {    /* BIG_ENDIAN */
+  unsigned int sign: 1;
+  unsigned int exponent: 11;
+  unsigned int mantissa0: 20;
+  unsigned int mantissa1: 32;
+};
+#else
+struct IeeeDoubleStruct {    /* LITTLE_ENDIAN */
+  unsigned int mantissa1: 32;
+  unsigned int mantissa0: 20;
+  unsigned int exponent: 11;
+  unsigned int sign: 1;
+};
+#endif
+
+/**Struct**********************************************************************
+
+  Synopsis    [IEEE double NaN struct.]
+
+  Description [IEEE double NaN struct.]
+
+  SeeAlso     []
+
+******************************************************************************/
+#ifdef    EPD_BIG_ENDIAN
+struct IeeeNanStruct {    /* BIG_ENDIAN */
+  unsigned int sign: 1;
+  unsigned int exponent: 11;
+  unsigned int quiet_bit: 1;
+  unsigned int mantissa0: 19;
+  unsigned int mantissa1: 32;
+};
+#else
+struct IeeeNanStruct {    /* LITTLE_ENDIAN */
+  unsigned int mantissa1: 32;
+  unsigned int mantissa0: 19;
+  unsigned int quiet_bit: 1;
+  unsigned int exponent: 11;
+  unsigned int sign: 1;
+};
+#endif
+
+/**Struct**********************************************************************
+
+  Synopsis    [Extended precision double to keep very large value.]
+
+  Description [Extended precision double to keep very large value.]
+
+  SeeAlso     []
+
+******************************************************************************/
+struct EpDoubleStruct {
+  union {
+    double            value;
+    struct IeeeDoubleStruct    bits;
+    struct IeeeNanStruct    nan;
+  } type;
+  int        exponent;
+};
+
+/*---------------------------------------------------------------------------*/
+/* Type declarations                                                         */
+/*---------------------------------------------------------------------------*/
+typedef struct EpDoubleStruct EpDouble;
+typedef struct IeeeDoubleStruct IeeeDouble;
+typedef struct IeeeNanStruct IeeeNan;
+
+
+/*---------------------------------------------------------------------------*/
+/* Function prototypes                                                       */
+/*---------------------------------------------------------------------------*/
+
+EpDouble *EpdAlloc();
+int EpdCmp(const char *key1, const char *key2);
+void EpdFree(EpDouble *epd);
+void EpdGetString(EpDouble *epd, char *str);
+void EpdConvert(double value, EpDouble *epd);
+void EpdMultiply(EpDouble *epd1, double value);
+void EpdMultiply2(EpDouble *epd1, EpDouble *epd2);
+void EpdMultiply2Decimal(EpDouble *epd1, EpDouble *epd2);
+void EpdMultiply3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3);
+void EpdMultiply3Decimal(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3);
+void EpdDivide(EpDouble *epd1, double value);
+void EpdDivide2(EpDouble *epd1, EpDouble *epd2);
+void EpdDivide3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3);
+void EpdAdd(EpDouble *epd1, double value);
+void EpdAdd2(EpDouble *epd1, EpDouble *epd2);
+void EpdAdd3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3);
+void EpdSubtract(EpDouble *epd1, double value);
+void EpdSubtract2(EpDouble *epd1, EpDouble *epd2);
+void EpdSubtract3(EpDouble *epd1, EpDouble *epd2, EpDouble *epd3);
+void EpdPow2(int n, EpDouble *epd);
+void EpdPow2Decimal(int n, EpDouble *epd);
+void EpdNormalize(EpDouble *epd);
+void EpdNormalizeDecimal(EpDouble *epd);
+void EpdGetValueAndDecimalExponent(EpDouble *epd, double *value, int *exponent);
+int EpdGetExponent(double value);
+int EpdGetExponentDecimal(double value);
+void EpdMakeInf(EpDouble *epd, int sign);
+void EpdMakeZero(EpDouble *epd, int sign);
+void EpdMakeNan(EpDouble *epd);
+void EpdCopy(EpDouble *from, EpDouble *to);
+int EpdIsInf(EpDouble *epd);
+int EpdIsZero(EpDouble *epd);
+int EpdIsNan(EpDouble *epd);
+int EpdIsNanOrInf(EpDouble *epd);
+int IsInfDouble(double value);
+int IsNanDouble(double value);
+int IsNanOrInfDouble(double value);
+
+#endif /* _EPD */
diff --git a/src/bdd/epd/module.make b/src/bdd/epd/module.make
new file mode 100644
index 00000000..a8084db1
--- /dev/null
+++ b/src/bdd/epd/module.make
@@ -0,0 +1 @@
+SRC +=  src/bdd/epd/epd.c