distr/epd/epd.c

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