epd.c File Reference

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "misc/util/util_hack.h"
#include "epd.h"

Go to the source code of this file.

Functions
ABC_NAMESPACE_IMPL_START EpDouble *	EpdAlloc (void)
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)

Function Documentation

void EpdAdd	(	EpDouble *	epd1,
		double	value
	)

Function********************************************************************

Synopsis [Adds two arbitrary precision double values.]

Description [Adds two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 536 of file epd.c.

{
  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);
}

void EpdAdd2	(	EpDouble *	epd1,
		EpDouble *	epd2
	)

Function********************************************************************

Synopsis [Adds two arbitrary precision double values.]

Description [Adds two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 598 of file epd.c.

{
  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);
}

void EpdAdd3	(	EpDouble *	epd1,
		EpDouble *	epd2,
		EpDouble *	epd3
	)

Function********************************************************************

Synopsis [Adds two arbitrary precision double values.]

Description [Adds two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 660 of file epd.c.

{
  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);
}

ABC_NAMESPACE_IMPL_START EpDouble* EpdAlloc

(

void

)

CFile***********************************************************************

FileName [epd.c]

PackageName [epd]

Synopsis [Arithmetic functions with extended double precision.]

Description []

SeeAlso []

Author [In-Ho Moon]

Copyright [Copyright (c) 1995-2004, Regents of the University of Colorado

All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

Neither the name of the University of Colorado nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.]

Revision [

Id:

epd.c,v 1.10 2004/08/13 18:20:30 fabio Exp

] Function********************************************************************

Synopsis [Allocates an EpDouble struct.]

Description [Allocates an EpDouble struct.]

SideEffects []

SeeAlso []

Definition at line 72 of file epd.c.

{
  EpDouble      *epd;

  epd = ABC_ALLOC(EpDouble, 1);
  return(epd);
}

int EpdCmp	(	const char *	key1,
		const char *	key2
	)

Function********************************************************************

Synopsis [Compares two EpDouble struct.]

Description [Compares two EpDouble struct.]

SideEffects []

SeeAlso []

Definition at line 93 of file epd.c.

{
  EpDouble *epd1 = (EpDouble *) key1;
  EpDouble *epd2 = (EpDouble *) key2;
  if (epd1->type.value != epd2->type.value ||
      epd1->exponent != epd2->exponent) {
    return(1);
  }
  return(0);
}

void EpdConvert	(	double	value,
		EpDouble *	epd
	)

Function********************************************************************

Synopsis [Converts double to EpDouble struct.]

Description [Converts double to EpDouble struct.]

SideEffects []

SeeAlso []

Definition at line 185 of file epd.c.

{
  epd->type.value = value;
  epd->exponent = 0;
  EpdNormalize(epd);
}

void EpdCopy	(	EpDouble *	from,
		EpDouble *	to
	)

Function********************************************************************

Synopsis [Copies a EpDouble struct.]

Description [Copies a EpDouble struct.]

SideEffects []

SeeAlso []

Definition at line 1182 of file epd.c.

{
  to->type.value = from->type.value;
  to->exponent = from->exponent;
}

void EpdDivide	(	EpDouble *	epd1,
		double	value
	)

Function********************************************************************

Synopsis [Divides two arbitrary precision double values.]

Description [Divides two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 386 of file epd.c.

{
  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);
}

void EpdDivide2	(	EpDouble *	epd1,
		EpDouble *	epd2
	)

Function********************************************************************

Synopsis [Divides two arbitrary precision double values.]

Description [Divides two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 439 of file epd.c.

{
  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);
}

void EpdDivide3	(	EpDouble *	epd1,
		EpDouble *	epd2,
		EpDouble *	epd3
	)

Function********************************************************************

Synopsis [Divides two arbitrary precision double values.]

Description [Divides two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 490 of file epd.c.

{
  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);
}

void EpdFree

(

EpDouble *

epd

)

Function********************************************************************

Synopsis [Frees an EpDouble struct.]

Description [Frees an EpDouble struct.]

SideEffects []

SeeAlso []

Definition at line 117 of file epd.c.

{
  ABC_FREE(epd);
}

int EpdGetExponent

(

double

value

)

Function********************************************************************

Synopsis [Returns the exponent value of a double.]

Description [Returns the exponent value of a double.]

SideEffects []

SeeAlso []

Definition at line 1068 of file epd.c.

{
  int           exponent;
  EpDouble      epd;

  epd.type.value = value;
  exponent = epd.type.bits.exponent;
  return(exponent);
}

int EpdGetExponentDecimal

(

double

value

)

Function********************************************************************

Synopsis [Returns the decimal exponent value of a double.]

Description [Returns the decimal exponent value of a double.]

SideEffects []

SeeAlso []

Definition at line 1091 of file epd.c.

{
  char  *pos, str[24];
  int   exponent;

  sprintf(str, "%E", value);
  pos = strstr(str, "E");
  sscanf(pos, "E%d", &exponent);
  return(exponent);
}

void EpdGetString	(	EpDouble *	epd,
		char *	str
	)

Function********************************************************************

Synopsis [Converts an arbitrary precision double value to a string.]

Description [Converts an arbitrary precision double value to a string.]

SideEffects []

SeeAlso []

Definition at line 135 of file epd.c.

{
  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);
  }
}

void EpdGetValueAndDecimalExponent	(	EpDouble *	epd,
		double *	value,
		int *	exponent
	)

Function********************************************************************

Synopsis [Returns value and decimal exponent of EpDouble.]

Description [Returns value and decimal exponent of EpDouble.]

SideEffects []

SeeAlso []

Definition at line 1034 of file epd.c.

{
  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;
}

int EpdIsInf

(

EpDouble *

epd

)

Function********************************************************************

Synopsis [Checks whether the value is Inf.]

Description [Checks whether the value is Inf.]

SideEffects []

SeeAlso []

Definition at line 1201 of file epd.c.

{
  return(IsInfDouble(epd->type.value));
}

int EpdIsNan

(

EpDouble *

epd

)

Function********************************************************************

Synopsis [Checks whether the value is NaN.]

Description [Checks whether the value is NaN.]

SideEffects []

SeeAlso []

Definition at line 1240 of file epd.c.

{
  return(IsNanDouble(epd->type.value));
}

int EpdIsNanOrInf

(

EpDouble *

epd

)

Function********************************************************************

Synopsis [Checks whether the value is NaN or Inf.]

Description [Checks whether the value is NaN or Inf.]

SideEffects []

SeeAlso []

Definition at line 1258 of file epd.c.

{
  return(IsNanOrInfDouble(epd->type.value));
}

int EpdIsZero

(

EpDouble *

epd

)

Function********************************************************************

Synopsis [Checks whether the value is Zero.]

Description [Checks whether the value is Zero.]

SideEffects []

SeeAlso []

Definition at line 1219 of file epd.c.

{
  if (epd->type.value == 0.0)
    return(1);
  else
    return(0);
}

void EpdMakeInf	(	EpDouble *	epd,
		int	sign
	)

Function********************************************************************

Synopsis [Makes EpDouble Inf.]

Description [Makes EpDouble Inf.]

SideEffects []

SeeAlso []

Definition at line 1115 of file epd.c.

{
  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;
}

void EpdMakeNan

(

EpDouble *

epd

)

Function********************************************************************

Synopsis [Makes EpDouble NaN.]

Description [Makes EpDouble NaN.]

SideEffects []

SeeAlso []

Definition at line 1159 of file epd.c.

{
  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;
}

void EpdMakeZero	(	EpDouble *	epd,
		int	sign
	)

Function********************************************************************

Synopsis [Makes EpDouble Zero.]

Description [Makes EpDouble Zero.]

SideEffects []

SeeAlso []

Definition at line 1137 of file epd.c.

{
  epd->type.bits.mantissa1 = 0;
  epd->type.bits.mantissa0 = 0;
  epd->type.bits.exponent = 0;
  epd->type.bits.sign = sign;
  epd->exponent = 0;
}

void EpdMultiply	(	EpDouble *	epd1,
		double	value
	)

Function********************************************************************

Synopsis [Multiplies two arbitrary precision double values.]

Description [Multiplies two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 205 of file epd.c.

{
  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);
}

void EpdMultiply2	(	EpDouble *	epd1,
		EpDouble *	epd2
	)

Function********************************************************************

Synopsis [Multiplies two arbitrary precision double values.]

Description [Multiplies two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 246 of file epd.c.

{
  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);
}

void EpdMultiply2Decimal	(	EpDouble *	epd1,
		EpDouble *	epd2
	)

Function********************************************************************

Synopsis [Multiplies two arbitrary precision double values.]

Description [Multiplies two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 285 of file epd.c.

{
  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);
}

void EpdMultiply3	(	EpDouble *	epd1,
		EpDouble *	epd2,
		EpDouble *	epd3
	)

Function********************************************************************

Synopsis [Multiplies two arbitrary precision double values.]

Description [Multiplies two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 321 of file epd.c.

{
  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);
}

void EpdMultiply3Decimal	(	EpDouble *	epd1,
		EpDouble *	epd2,
		EpDouble *	epd3
	)

Function********************************************************************

Synopsis [Multiplies two arbitrary precision double values.]

Description [Multiplies two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 355 of file epd.c.

{
  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);
}

void EpdNormalize

(

EpDouble *

epd

)

Function********************************************************************

Synopsis [Normalize an arbitrary precision double value.]

Description [Normalize an arbitrary precision double value.]

SideEffects []

SeeAlso []

Definition at line 977 of file epd.c.

{
  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;
}

void EpdNormalizeDecimal

(

EpDouble *

epd

)

Function********************************************************************

Synopsis [Normalize an arbitrary precision double value.]

Description [Normalize an arbitrary precision double value.]

SideEffects []

SeeAlso []

Definition at line 1007 of file epd.c.

{
  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;
}

void EpdPow2	(	int	n,
		EpDouble *	epd
	)

Function********************************************************************

Synopsis [Computes arbitrary precision pow of base 2.]

Description [Computes arbitrary precision pow of base 2.]

SideEffects []

SeeAlso []

Definition at line 917 of file epd.c.

{
  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);
  }
}

void EpdPow2Decimal	(	int	n,
		EpDouble *	epd
	)

Function********************************************************************

Synopsis [Computes arbitrary precision pow of base 2.]

Description [Computes arbitrary precision pow of base 2.]

SideEffects []

SeeAlso []

Definition at line 946 of file epd.c.

{
  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);
  }
}

void EpdSubtract	(	EpDouble *	epd1,
		double	value
	)

Function********************************************************************

Synopsis [Subtracts two arbitrary precision double values.]

Description [Subtracts two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 726 of file epd.c.

{
  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);
}

void EpdSubtract2	(	EpDouble *	epd1,
		EpDouble *	epd2
	)

Function********************************************************************

Synopsis [Subtracts two arbitrary precision double values.]

Description [Subtracts two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 788 of file epd.c.

{
  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);
}

void EpdSubtract3	(	EpDouble *	epd1,
		EpDouble *	epd2,
		EpDouble *	epd3
	)

Function********************************************************************

Synopsis [Subtracts two arbitrary precision double values.]

Description [Subtracts two arbitrary precision double values.]

SideEffects []

SeeAlso []

Definition at line 850 of file epd.c.

{
  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);
}

int IsInfDouble

(

double

value

)

Function********************************************************************

Synopsis [Checks whether the value is Inf.]

Description [Checks whether the value is Inf.]

SideEffects []

SeeAlso []

Definition at line 1276 of file epd.c.

{
  EpType val;

  val.value = value;
  if (val.bits.exponent == EPD_EXP_INF &&
      val.bits.mantissa0 == 0 &&
      val.bits.mantissa1 == 0) {
    if (val.bits.sign == 0)
      return(1);
    else
      return(-1);
  }
  return(0);
}

int IsNanDouble

(

double

value

)

Function********************************************************************

Synopsis [Checks whether the value is NaN.]

Description [Checks whether the value is NaN.]

SideEffects []

SeeAlso []

Definition at line 1305 of file epd.c.

{
  EpType        val;
  
  val.value = value;
  if (val.nan.exponent == EPD_EXP_INF &&
      val.nan.sign == 1 &&
      val.nan.quiet_bit == 1 &&
      val.nan.mantissa0 == 0 &&
      val.nan.mantissa1 == 0) {
    return(1);
  }
  return(0);
}

int IsNanOrInfDouble

(

double

value

)

Function********************************************************************

Synopsis [Checks whether the value is NaN or Inf.]

Description [Checks whether the value is NaN or Inf.]

SideEffects []

SeeAlso []

Definition at line 1333 of file epd.c.

{
  EpType        val;

  val.value = value;
  if (val.nan.exponent == EPD_EXP_INF &&
      val.nan.mantissa0 == 0 &&
      val.nan.mantissa1 == 0 &&
      (val.nan.sign == 1 || val.nan.quiet_bit == 0)) {
    return(1);
  }
  return(0);
}