source: vis_dev/glu-2.1/src/epd/epd.c @ 8

Last change on this file since 8 was 8, checked in by cecile, 13 years ago

src glu
File size: 31.6 KB
Line 
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******************************************************************************/
70EpDouble *
71EpdAlloc(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******************************************************************************/
91int
92EpdCmp(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******************************************************************************/
115void
116EpdFree(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******************************************************************************/
133void
134EpdGetString(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******************************************************************************/
183void
184EpdConvert(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******************************************************************************/
203void
204EpdMultiply(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******************************************************************************/
244void
245EpdMultiply2(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******************************************************************************/
283void
284EpdMultiply2Decimal(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******************************************************************************/
319void
320EpdMultiply3(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******************************************************************************/
353void
354EpdMultiply3Decimal(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******************************************************************************/
384void
385EpdDivide(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******************************************************************************/
437void
438EpdDivide2(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******************************************************************************/
488void
489EpdDivide3(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******************************************************************************/
534void
535EpdAdd(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******************************************************************************/
596void
597EpdAdd2(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******************************************************************************/
658void
659EpdAdd3(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******************************************************************************/
724void
725EpdSubtract(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******************************************************************************/
786void
787EpdSubtract2(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******************************************************************************/
848void
849EpdSubtract3(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******************************************************************************/
915void
916EpdPow2(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******************************************************************************/
944void
945EpdPow2Decimal(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******************************************************************************/
975void
976EpdNormalize(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******************************************************************************/
1005void
1006EpdNormalizeDecimal(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******************************************************************************/
1032void
1033EpdGetValueAndDecimalExponent(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******************************************************************************/
1066int
1067EpdGetExponent(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******************************************************************************/
1089int
1090EpdGetExponentDecimal(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******************************************************************************/
1113void
1114EpdMakeInf(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******************************************************************************/
1135void
1136EpdMakeZero(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******************************************************************************/
1157void
1158EpdMakeNan(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******************************************************************************/
1180void
1181EpdCopy(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******************************************************************************/
1199int
1200EpdIsInf(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******************************************************************************/
1217int
1218EpdIsZero(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******************************************************************************/
1238int
1239EpdIsNan(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******************************************************************************/
1256int
1257EpdIsNanOrInf(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******************************************************************************/
1274int
1275IsInfDouble(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******************************************************************************/
1303int
1304IsNanDouble(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******************************************************************************/
1331int
1332IsNanOrInfDouble(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}
Note: See TracBrowser for help on using the repository browser.