src/epd/epd.c File Reference

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "util.h"
#include "epd.h"
Include dependency graph for epd.c:

Go to the source code of this file.

Functions

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 535 of file epd.c.

00536 {
00537   EpDouble      epd2;
00538   double        tmp;
00539   int           exponent, diff;
00540 
00541   if (EpdIsNan(epd1) || IsNanDouble(value)) {
00542     EpdMakeNan(epd1);
00543     return;
00544   } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00545     int sign;
00546 
00547     EpdConvert(value, &epd2);
00548     if (EpdIsInf(epd1) && IsInfDouble(value)) {
00549       sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00550       if (sign == 1)
00551         EpdMakeNan(epd1);
00552     } else if (EpdIsInf(&epd2)) {
00553       EpdCopy(&epd2, epd1);
00554     }
00555     return;
00556   }
00557 
00558   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00559 
00560   EpdConvert(value, &epd2);
00561   if (epd1->exponent > epd2.exponent) {
00562     diff = epd1->exponent - epd2.exponent;
00563     if (diff <= EPD_MAX_BIN)
00564       tmp = epd1->type.value + epd2.type.value / pow((double)2.0, (double)diff);
00565     else
00566       tmp = epd1->type.value;
00567     exponent = epd1->exponent;
00568   } else if (epd1->exponent < epd2.exponent) {
00569     diff = epd2.exponent - epd1->exponent;
00570     if (diff <= EPD_MAX_BIN)
00571       tmp = epd1->type.value / pow((double)2.0, (double)diff) + epd2.type.value;
00572     else
00573       tmp = epd2.type.value;
00574     exponent = epd2.exponent;
00575   } else {
00576     tmp = epd1->type.value + epd2.type.value;
00577     exponent = epd1->exponent;
00578   }
00579   epd1->type.value = tmp;
00580   epd1->exponent = exponent;
00581   EpdNormalize(epd1);
00582 }

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 597 of file epd.c.

00598 {
00599   double        value;
00600   int           exponent, diff;
00601 
00602   if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00603     EpdMakeNan(epd1);
00604     return;
00605   } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00606     int sign;
00607 
00608     if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00609       sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00610       if (sign == 1)
00611         EpdMakeNan(epd1);
00612     } else if (EpdIsInf(epd2)) {
00613       EpdCopy(epd2, epd1);
00614     }
00615     return;
00616   }
00617 
00618   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00619   assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00620 
00621   if (epd1->exponent > epd2->exponent) {
00622     diff = epd1->exponent - epd2->exponent;
00623     if (diff <= EPD_MAX_BIN) {
00624       value = epd1->type.value +
00625                 epd2->type.value / pow((double)2.0, (double)diff);
00626     } else
00627       value = epd1->type.value;
00628     exponent = epd1->exponent;
00629   } else if (epd1->exponent < epd2->exponent) {
00630     diff = epd2->exponent - epd1->exponent;
00631     if (diff <= EPD_MAX_BIN) {
00632       value = epd1->type.value / pow((double)2.0, (double)diff) +
00633                 epd2->type.value;
00634     } else
00635       value = epd2->type.value;
00636     exponent = epd2->exponent;
00637   } else {
00638     value = epd1->type.value + epd2->type.value;
00639     exponent = epd1->exponent;
00640   }
00641   epd1->type.value = value;
00642   epd1->exponent = exponent;
00643   EpdNormalize(epd1);
00644 }

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 659 of file epd.c.

00660 {
00661   double        value;
00662   int           exponent, diff;
00663 
00664   if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00665     EpdMakeNan(epd3);
00666     return;
00667   } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00668     int sign;
00669 
00670     if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00671       sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00672       if (sign == 1)
00673         EpdMakeNan(epd3);
00674       else
00675         EpdCopy(epd1, epd3);
00676     } else if (EpdIsInf(epd1)) {
00677       EpdCopy(epd1, epd3);
00678     } else {
00679       EpdCopy(epd2, epd3);
00680     }
00681     return;
00682   }
00683 
00684   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00685   assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00686 
00687   if (epd1->exponent > epd2->exponent) {
00688     diff = epd1->exponent - epd2->exponent;
00689     if (diff <= EPD_MAX_BIN) {
00690       value = epd1->type.value +
00691                 epd2->type.value / pow((double)2.0, (double)diff);
00692     } else
00693       value = epd1->type.value;
00694     exponent = epd1->exponent;
00695   } else if (epd1->exponent < epd2->exponent) {
00696     diff = epd2->exponent - epd1->exponent;
00697     if (diff <= EPD_MAX_BIN) {
00698       value = epd1->type.value / pow((double)2.0, (double)diff) +
00699                 epd2->type.value;
00700     } else
00701       value = epd2->type.value;
00702     exponent = epd2->exponent;
00703   } else {
00704     value = epd1->type.value + epd2->type.value;
00705     exponent = epd1->exponent;
00706   }
00707   epd3->type.value = value;
00708   epd3->exponent = exponent;
00709   EpdNormalize(epd3);
00710 }

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 71 of file epd.c.

00072 {
00073   EpDouble      *epd;
00074 
00075   epd = ALLOC(EpDouble, 1);
00076   return(epd);
00077 }

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

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

Synopsis [Compares two EpDouble struct.]

Description [Compares two EpDouble struct.]

SideEffects []

SeeAlso []

Definition at line 92 of file epd.c.

00093 {
00094   EpDouble *epd1 = (EpDouble *) key1;
00095   EpDouble *epd2 = (EpDouble *) key2;
00096   if (epd1->type.value != epd2->type.value ||
00097       epd1->exponent != epd2->exponent) {
00098     return(1);
00099   }
00100   return(0);
00101 }

void EpdConvert ( double  value,
EpDouble *  epd 
)

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

Synopsis [Converts double to EpDouble struct.]

Description [Converts double to EpDouble struct.]

SideEffects []

SeeAlso []

Definition at line 184 of file epd.c.

00185 {
00186   epd->type.value = value;
00187   epd->exponent = 0;
00188   EpdNormalize(epd);
00189 }

void EpdCopy ( EpDouble *  from,
EpDouble *  to 
)

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

Synopsis [Copies a EpDouble struct.]

Description [Copies a EpDouble struct.]

SideEffects []

SeeAlso []

Definition at line 1181 of file epd.c.

01182 {
01183   to->type.value = from->type.value;
01184   to->exponent = from->exponent;
01185 }

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 385 of file epd.c.

00386 {
00387   EpDouble      epd2;
00388   double        tmp;
00389   int           exponent;
00390 
00391   if (EpdIsNan(epd1) || IsNanDouble(value)) {
00392     EpdMakeNan(epd1);
00393     return;
00394   } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00395     int sign;
00396 
00397     EpdConvert(value, &epd2);
00398     if (EpdIsInf(epd1) && IsInfDouble(value)) {
00399       EpdMakeNan(epd1);
00400     } else if (EpdIsInf(epd1)) {
00401       sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00402       EpdMakeInf(epd1, sign);
00403     } else {
00404       sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00405       EpdMakeZero(epd1, sign);
00406     }
00407     return;
00408   }
00409 
00410   if (value == 0.0) {
00411     EpdMakeNan(epd1);
00412     return;
00413   }
00414 
00415   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00416 
00417   EpdConvert(value, &epd2);
00418   tmp = epd1->type.value / epd2.type.value;
00419   exponent = epd1->exponent - epd2.exponent;
00420   epd1->type.value = tmp;
00421   epd1->exponent = exponent;
00422   EpdNormalize(epd1);
00423 }

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 438 of file epd.c.

00439 {
00440   double        value;
00441   int           exponent;
00442 
00443   if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00444     EpdMakeNan(epd1);
00445     return;
00446   } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00447     int sign;
00448 
00449     if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00450       EpdMakeNan(epd1);
00451     } else if (EpdIsInf(epd1)) {
00452       sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00453       EpdMakeInf(epd1, sign);
00454     } else {
00455       sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00456       EpdMakeZero(epd1, sign);
00457     }
00458     return;
00459   }
00460 
00461   if (epd2->type.value == 0.0) {
00462     EpdMakeNan(epd1);
00463     return;
00464   }
00465 
00466   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00467   assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00468 
00469   value = epd1->type.value / epd2->type.value;
00470   exponent = epd1->exponent - epd2->exponent;
00471   epd1->type.value = value;
00472   epd1->exponent = exponent;
00473   EpdNormalize(epd1);
00474 }

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 489 of file epd.c.

00490 {
00491   if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00492     EpdMakeNan(epd3);
00493     return;
00494   } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00495     int sign;
00496 
00497     if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00498       EpdMakeNan(epd3);
00499     } else if (EpdIsInf(epd1)) {
00500       sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00501       EpdMakeInf(epd3, sign);
00502     } else {
00503       sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00504       EpdMakeZero(epd3, sign);
00505     }
00506     return;
00507   }
00508 
00509   if (epd2->type.value == 0.0) {
00510     EpdMakeNan(epd3);
00511     return;
00512   }
00513 
00514   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00515   assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00516 
00517   epd3->type.value = epd1->type.value / epd2->type.value;
00518   epd3->exponent = epd1->exponent - epd2->exponent;
00519   EpdNormalize(epd3);
00520 }

void EpdFree ( EpDouble *  epd  ) 

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

Synopsis [Frees an EpDouble struct.]

Description [Frees an EpDouble struct.]

SideEffects []

SeeAlso []

Definition at line 116 of file epd.c.

00117 {
00118   FREE(epd);
00119 }

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 1067 of file epd.c.

01068 {
01069   int           exponent;
01070   EpDouble      epd;
01071 
01072   epd.type.value = value;
01073   exponent = epd.type.bits.exponent;
01074   return(exponent);
01075 }

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 1090 of file epd.c.

01091 {
01092   char  *pos, str[24];
01093   int   exponent;
01094 
01095   sprintf(str, "%E", value);
01096   pos = strstr(str, "E");
01097   sscanf(pos, "E%d", &exponent);
01098   return(exponent);
01099 }

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 134 of file epd.c.

00135 {
00136   double        value;
00137   int           exponent;
00138   char          *pos;
00139 
00140   if (IsNanDouble(epd->type.value)) {
00141     sprintf(str, "NaN");
00142     return;
00143   } else if (IsInfDouble(epd->type.value)) {
00144     if (epd->type.bits.sign == 1)
00145       sprintf(str, "-Inf");
00146     else
00147       sprintf(str, "Inf");
00148     return;
00149   }
00150 
00151   assert(epd->type.bits.exponent == EPD_MAX_BIN ||
00152          epd->type.bits.exponent == 0);
00153 
00154   EpdGetValueAndDecimalExponent(epd, &value, &exponent);
00155   sprintf(str, "%e", value);
00156   pos = strstr(str, "e");
00157   if (exponent >= 0) {
00158     if (exponent < 10)
00159       sprintf(pos + 1, "+0%d", exponent);
00160     else
00161       sprintf(pos + 1, "+%d", exponent);
00162   } else {
00163     exponent *= -1;
00164     if (exponent < 10)
00165       sprintf(pos + 1, "-0%d", exponent);
00166     else
00167       sprintf(pos + 1, "-%d", exponent);
00168   }
00169 }

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 1033 of file epd.c.

01034 {
01035   EpDouble      epd1, epd2;
01036 
01037   if (EpdIsNanOrInf(epd))
01038     return;
01039 
01040   if (EpdIsZero(epd)) {
01041     *value = 0.0;
01042     *exponent = 0;
01043     return;
01044   }
01045 
01046   epd1.type.value = epd->type.value;
01047   epd1.exponent = 0;
01048   EpdPow2Decimal(epd->exponent, &epd2);
01049   EpdMultiply2Decimal(&epd1, &epd2);
01050 
01051   *value = epd1.type.value;
01052   *exponent = epd1.exponent;
01053 }

int EpdIsInf ( EpDouble *  epd  ) 

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

Synopsis [Checks whether the value is Inf.]

Description [Checks whether the value is Inf.]

SideEffects []

SeeAlso []

Definition at line 1200 of file epd.c.

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

int EpdIsNan ( EpDouble *  epd  ) 

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

Synopsis [Checks whether the value is NaN.]

Description [Checks whether the value is NaN.]

SideEffects []

SeeAlso []

Definition at line 1239 of file epd.c.

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

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 1257 of file epd.c.

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

int EpdIsZero ( EpDouble *  epd  ) 

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

Synopsis [Checks whether the value is Zero.]

Description [Checks whether the value is Zero.]

SideEffects []

SeeAlso []

Definition at line 1218 of file epd.c.

01219 {
01220   if (epd->type.value == 0.0)
01221     return(1);
01222   else
01223     return(0);
01224 }

void EpdMakeInf ( EpDouble *  epd,
int  sign 
)

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

Synopsis [Makes EpDouble Inf.]

Description [Makes EpDouble Inf.]

SideEffects []

SeeAlso []

Definition at line 1114 of file epd.c.

01115 {
01116   epd->type.bits.mantissa1 = 0;
01117   epd->type.bits.mantissa0 = 0;
01118   epd->type.bits.exponent = EPD_EXP_INF;
01119   epd->type.bits.sign = sign;
01120   epd->exponent = 0;
01121 }

void EpdMakeNan ( EpDouble *  epd  ) 

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

Synopsis [Makes EpDouble NaN.]

Description [Makes EpDouble NaN.]

SideEffects []

SeeAlso []

Definition at line 1158 of file epd.c.

01159 {
01160   epd->type.nan.mantissa1 = 0;
01161   epd->type.nan.mantissa0 = 0;
01162   epd->type.nan.quiet_bit = 1;
01163   epd->type.nan.exponent = EPD_EXP_INF;
01164   epd->type.nan.sign = 1;
01165   epd->exponent = 0;
01166 }

void EpdMakeZero ( EpDouble *  epd,
int  sign 
)

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

Synopsis [Makes EpDouble Zero.]

Description [Makes EpDouble Zero.]

SideEffects []

SeeAlso []

Definition at line 1136 of file epd.c.

01137 {
01138   epd->type.bits.mantissa1 = 0;
01139   epd->type.bits.mantissa0 = 0;
01140   epd->type.bits.exponent = 0;
01141   epd->type.bits.sign = sign;
01142   epd->exponent = 0;
01143 }

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 204 of file epd.c.

00205 {
00206   EpDouble      epd2;
00207   double        tmp;
00208   int           exponent;
00209 
00210   if (EpdIsNan(epd1) || IsNanDouble(value)) {
00211     EpdMakeNan(epd1);
00212     return;
00213   } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00214     int sign;
00215 
00216     EpdConvert(value, &epd2);
00217     sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00218     EpdMakeInf(epd1, sign);
00219     return;
00220   }
00221 
00222   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00223 
00224   EpdConvert(value, &epd2);
00225   tmp = epd1->type.value * epd2.type.value;
00226   exponent = epd1->exponent + epd2.exponent;
00227   epd1->type.value = tmp;
00228   epd1->exponent = exponent;
00229   EpdNormalize(epd1);
00230 }

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 245 of file epd.c.

00246 {
00247   double        value;
00248   int           exponent;
00249 
00250   if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00251     EpdMakeNan(epd1);
00252     return;
00253   } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00254     int sign;
00255 
00256     sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00257     EpdMakeInf(epd1, sign);
00258     return;
00259   }
00260 
00261   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00262   assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00263 
00264   value = epd1->type.value * epd2->type.value;
00265   exponent = epd1->exponent + epd2->exponent;
00266   epd1->type.value = value;
00267   epd1->exponent = exponent;
00268   EpdNormalize(epd1);
00269 }

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 284 of file epd.c.

00285 {
00286   double        value;
00287   int           exponent;
00288 
00289   if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00290     EpdMakeNan(epd1);
00291     return;
00292   } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00293     int sign;
00294 
00295     sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00296     EpdMakeInf(epd1, sign);
00297     return;
00298   }
00299 
00300   value = epd1->type.value * epd2->type.value;
00301   exponent = epd1->exponent + epd2->exponent;
00302   epd1->type.value = value;
00303   epd1->exponent = exponent;
00304   EpdNormalizeDecimal(epd1);
00305 }

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 320 of file epd.c.

00321 {
00322   if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00323     EpdMakeNan(epd1);
00324     return;
00325   } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00326     int sign;
00327 
00328     sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00329     EpdMakeInf(epd3, sign);
00330     return;
00331   }
00332 
00333   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00334   assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00335 
00336   epd3->type.value = epd1->type.value * epd2->type.value;
00337   epd3->exponent = epd1->exponent + epd2->exponent;
00338   EpdNormalize(epd3);
00339 }

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 354 of file epd.c.

00355 {
00356   if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00357     EpdMakeNan(epd1);
00358     return;
00359   } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00360     int sign;
00361 
00362     sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00363     EpdMakeInf(epd3, sign);
00364     return;
00365   }
00366 
00367   epd3->type.value = epd1->type.value * epd2->type.value;
00368   epd3->exponent = epd1->exponent + epd2->exponent;
00369   EpdNormalizeDecimal(epd3);
00370 }

void EpdNormalize ( EpDouble *  epd  ) 

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

Synopsis [Normalize an arbitrary precision double value.]

Description [Normalize an arbitrary precision double value.]

SideEffects []

SeeAlso []

Definition at line 976 of file epd.c.

00977 {
00978   int           exponent;
00979 
00980   if (IsNanOrInfDouble(epd->type.value)) {
00981     epd->exponent = 0;
00982     return;
00983   }
00984 
00985   exponent = EpdGetExponent(epd->type.value);
00986   if (exponent == EPD_MAX_BIN)
00987     return;
00988   exponent -= EPD_MAX_BIN;
00989   epd->type.bits.exponent = EPD_MAX_BIN;
00990   epd->exponent += exponent;
00991 }

void EpdNormalizeDecimal ( EpDouble *  epd  ) 

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

Synopsis [Normalize an arbitrary precision double value.]

Description [Normalize an arbitrary precision double value.]

SideEffects []

SeeAlso []

Definition at line 1006 of file epd.c.

01007 {
01008   int           exponent;
01009 
01010   if (IsNanOrInfDouble(epd->type.value)) {
01011     epd->exponent = 0;
01012     return;
01013   }
01014 
01015   exponent = EpdGetExponentDecimal(epd->type.value);
01016   epd->type.value /= pow((double)10.0, (double)exponent);
01017   epd->exponent += exponent;
01018 }

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 916 of file epd.c.

00917 {
00918   if (n <= EPD_MAX_BIN) {
00919     EpdConvert(pow((double)2.0, (double)n), epd);
00920   } else {
00921     EpDouble    epd1, epd2;
00922     int         n1, n2;
00923 
00924     n1 = n / 2;
00925     n2 = n - n1;
00926     EpdPow2(n1, &epd1);
00927     EpdPow2(n2, &epd2);
00928     EpdMultiply3(&epd1, &epd2, epd);
00929   }
00930 }

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 945 of file epd.c.

00946 {
00947   if (n <= EPD_MAX_BIN) {
00948     epd->type.value = pow((double)2.0, (double)n);
00949     epd->exponent = 0;
00950     EpdNormalizeDecimal(epd);
00951   } else {
00952     EpDouble    epd1, epd2;
00953     int         n1, n2;
00954 
00955     n1 = n / 2;
00956     n2 = n - n1;
00957     EpdPow2Decimal(n1, &epd1);
00958     EpdPow2Decimal(n2, &epd2);
00959     EpdMultiply3Decimal(&epd1, &epd2, epd);
00960   }
00961 }

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 725 of file epd.c.

00726 {
00727   EpDouble      epd2;
00728   double        tmp;
00729   int           exponent, diff;
00730 
00731   if (EpdIsNan(epd1) || IsNanDouble(value)) {
00732     EpdMakeNan(epd1);
00733     return;
00734   } else if (EpdIsInf(epd1) || IsInfDouble(value)) {
00735     int sign;
00736 
00737     EpdConvert(value, &epd2);
00738     if (EpdIsInf(epd1) && IsInfDouble(value)) {
00739       sign = epd1->type.bits.sign ^ epd2.type.bits.sign;
00740       if (sign == 0)
00741         EpdMakeNan(epd1);
00742     } else if (EpdIsInf(&epd2)) {
00743       EpdCopy(&epd2, epd1);
00744     }
00745     return;
00746   }
00747 
00748   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00749 
00750   EpdConvert(value, &epd2);
00751   if (epd1->exponent > epd2.exponent) {
00752     diff = epd1->exponent - epd2.exponent;
00753     if (diff <= EPD_MAX_BIN)
00754       tmp = epd1->type.value - epd2.type.value / pow((double)2.0, (double)diff);
00755     else
00756       tmp = epd1->type.value;
00757     exponent = epd1->exponent;
00758   } else if (epd1->exponent < epd2.exponent) {
00759     diff = epd2.exponent - epd1->exponent;
00760     if (diff <= EPD_MAX_BIN)
00761       tmp = epd1->type.value / pow((double)2.0, (double)diff) - epd2.type.value;
00762     else
00763       tmp = epd2.type.value * (double)(-1.0);
00764     exponent = epd2.exponent;
00765   } else {
00766     tmp = epd1->type.value - epd2.type.value;
00767     exponent = epd1->exponent;
00768   }
00769   epd1->type.value = tmp;
00770   epd1->exponent = exponent;
00771   EpdNormalize(epd1);
00772 }

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 787 of file epd.c.

00788 {
00789   double        value;
00790   int           exponent, diff;
00791 
00792   if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00793     EpdMakeNan(epd1);
00794     return;
00795   } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00796     int sign;
00797 
00798     if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00799       sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00800       if (sign == 0)
00801         EpdMakeNan(epd1);
00802     } else if (EpdIsInf(epd2)) {
00803       EpdCopy(epd2, epd1);
00804     }
00805     return;
00806   }
00807 
00808   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00809   assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00810 
00811   if (epd1->exponent > epd2->exponent) {
00812     diff = epd1->exponent - epd2->exponent;
00813     if (diff <= EPD_MAX_BIN) {
00814       value = epd1->type.value -
00815                 epd2->type.value / pow((double)2.0, (double)diff);
00816     } else
00817       value = epd1->type.value;
00818     exponent = epd1->exponent;
00819   } else if (epd1->exponent < epd2->exponent) {
00820     diff = epd2->exponent - epd1->exponent;
00821     if (diff <= EPD_MAX_BIN) {
00822       value = epd1->type.value / pow((double)2.0, (double)diff) -
00823                 epd2->type.value;
00824     } else
00825       value = epd2->type.value * (double)(-1.0);
00826     exponent = epd2->exponent;
00827   } else {
00828     value = epd1->type.value - epd2->type.value;
00829     exponent = epd1->exponent;
00830   }
00831   epd1->type.value = value;
00832   epd1->exponent = exponent;
00833   EpdNormalize(epd1);
00834 }

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 849 of file epd.c.

00850 {
00851   double        value;
00852   int           exponent, diff;
00853 
00854   if (EpdIsNan(epd1) || EpdIsNan(epd2)) {
00855     EpdMakeNan(epd3);
00856     return;
00857   } else if (EpdIsInf(epd1) || EpdIsInf(epd2)) {
00858     int sign;
00859 
00860     if (EpdIsInf(epd1) && EpdIsInf(epd2)) {
00861       sign = epd1->type.bits.sign ^ epd2->type.bits.sign;
00862       if (sign == 0)
00863         EpdCopy(epd1, epd3);
00864       else
00865         EpdMakeNan(epd3);
00866     } else if (EpdIsInf(epd1)) {
00867       EpdCopy(epd1, epd1);
00868     } else {
00869       sign = epd2->type.bits.sign ^ 0x1;
00870       EpdMakeInf(epd3, sign);
00871     }
00872     return;
00873   }
00874 
00875   assert(epd1->type.bits.exponent == EPD_MAX_BIN);
00876   assert(epd2->type.bits.exponent == EPD_MAX_BIN);
00877 
00878   if (epd1->exponent > epd2->exponent) {
00879     diff = epd1->exponent - epd2->exponent;
00880     if (diff <= EPD_MAX_BIN) {
00881       value = epd1->type.value -
00882                 epd2->type.value / pow((double)2.0, (double)diff);
00883     } else
00884       value = epd1->type.value;
00885     exponent = epd1->exponent;
00886   } else if (epd1->exponent < epd2->exponent) {
00887     diff = epd2->exponent - epd1->exponent;
00888     if (diff <= EPD_MAX_BIN) {
00889       value = epd1->type.value / pow((double)2.0, (double)diff) -
00890                 epd2->type.value;
00891     } else
00892       value = epd2->type.value * (double)(-1.0);
00893     exponent = epd2->exponent;
00894   } else {
00895     value = epd1->type.value - epd2->type.value;
00896     exponent = epd1->exponent;
00897   }
00898   epd3->type.value = value;
00899   epd3->exponent = exponent;
00900   EpdNormalize(epd3);
00901 }

int IsInfDouble ( double  value  ) 

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

Synopsis [Checks whether the value is Inf.]

Description [Checks whether the value is Inf.]

SideEffects []

SeeAlso []

Definition at line 1275 of file epd.c.

01276 {
01277   EpType val;
01278 
01279   val.value = value;
01280   if (val.bits.exponent == EPD_EXP_INF &&
01281       val.bits.mantissa0 == 0 &&
01282       val.bits.mantissa1 == 0) {
01283     if (val.bits.sign == 0)
01284       return(1);
01285     else
01286       return(-1);
01287   }
01288   return(0);
01289 }

int IsNanDouble ( double  value  ) 

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

Synopsis [Checks whether the value is NaN.]

Description [Checks whether the value is NaN.]

SideEffects []

SeeAlso []

Definition at line 1304 of file epd.c.

01305 {
01306   EpType        val;
01307   
01308   val.value = value;
01309   if (val.nan.exponent == EPD_EXP_INF &&
01310       val.nan.sign == 1 &&
01311       val.nan.quiet_bit == 1 &&
01312       val.nan.mantissa0 == 0 &&
01313       val.nan.mantissa1 == 0) {
01314     return(1);
01315   }
01316   return(0);
01317 }

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 1332 of file epd.c.

01333 {
01334   EpType        val;
01335 
01336   val.value = value;
01337   if (val.nan.exponent == EPD_EXP_INF &&
01338       val.nan.mantissa0 == 0 &&
01339       val.nan.mantissa1 == 0 &&
01340       (val.nan.sign == 1 || val.nan.quiet_bit == 0)) {
01341     return(1);
01342   }
01343   return(0);
01344 }

Generated on Tue Jan 12 13:57:23 2010 for glu-2.2 by  doxygen 1.6.1