contrib/gdtoa/g_ddfmt.c

   1 /****************************************************************
   2
   3 The author of this software is David M. Gay.
   4
   5 Copyright (C) 1998 by Lucent Technologies
   6 All Rights Reserved
   7
   8 Permission to use, copy, modify, and distribute this software and
   9 its documentation for any purpose and without fee is hereby
  10 granted, provided that the above copyright notice appear in all
  11 copies and that both that the copyright notice and this
  12 permission notice and warranty disclaimer appear in supporting
  13 documentation, and that the name of Lucent or any of its entities
  14 not be used in advertising or publicity pertaining to
  15 distribution of the software without specific, written prior
  16 permission.
  17
  18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  25 THIS SOFTWARE.
  26
  27 ****************************************************************/
  28
  29 /* Please send bug reports to David M. Gay (dmg@acm.org). */
  30
  31 #include "gdtoaimp.h"
  32 #include <string.h>
  33
  34  char *
  35 #ifdef KR_headers
  36 g_ddfmt(buf, dd0, ndig, bufsize) char *buf; double *dd0; int ndig; size_t bufsize;
  37 #else
  38 g_ddfmt(char *buf, double *dd0, int ndig, size_t bufsize)
  39 #endif
  40 {
  41         FPI fpi;
  42         char *b, *s, *se;
  43         ULong *L, bits0[4], *bits, *zx;
  44         int bx, by, decpt, ex, ey, i, j, mode;
  45         Bigint *x, *y, *z;
  46         U *dd, ddx[2];
  47 #ifdef Honor_FLT_ROUNDS /*{{*/
  48         int Rounding;
  49 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
  50         Rounding = Flt_Rounds;
  51 #else /*}{*/
  52         Rounding = 1;
  53         switch(fegetround()) {
  54           case FE_TOWARDZERO:   Rounding = 0; break;
  55           case FE_UPWARD:       Rounding = 2; break;
  56           case FE_DOWNWARD:     Rounding = 3;
  57           }
  58 #endif /*}}*/
  59 #else /*}{*/
  60 #define Rounding FPI_Round_near
  61 #endif /*}}*/
  62
  63         if (bufsize < 10 || bufsize < ndig + 8)
  64                 return 0;
  65
  66         dd = (U*)dd0;
  67         L = dd->L;
  68         if ((L[_0] & 0x7ff00000L) == 0x7ff00000L) {
  69                 /* Infinity or NaN */
  70                 if (L[_0] & 0xfffff || L[_1]) {
  71  nanret:
  72                         return strcp(buf, "NaN");
  73                         }
  74                 if ((L[2+_0] & 0x7ff00000) == 0x7ff00000) {
  75                         if (L[2+_0] & 0xfffff || L[2+_1])
  76                                 goto nanret;
  77                         if ((L[_0] ^ L[2+_0]) & 0x80000000L)
  78                                 goto nanret;    /* Infinity - Infinity */
  79                         }
  80  infret:
  81                 b = buf;
  82                 if (L[_0] & 0x80000000L)
  83                         *b++ = '-';
  84                 return strcp(b, "Infinity");
  85                 }
  86         if ((L[2+_0] & 0x7ff00000) == 0x7ff00000) {
  87                 L += 2;
  88                 if (L[_0] & 0xfffff || L[_1])
  89                         goto nanret;
  90                 goto infret;
  91                 }
  92         if (dval(&dd[0]) + dval(&dd[1]) == 0.) {
  93                 b = buf;
  94 #ifndef IGNORE_ZERO_SIGN
  95                 if (L[_0] & L[2+_0] & 0x80000000L)
  96                         *b++ = '-';
  97 #endif
  98                 *b++ = '0';
  99                 *b = 0;
 100                 return b;
 101                 }
 102         if ((L[_0] & 0x7ff00000L) < (L[2+_0] & 0x7ff00000L)) {
 103                 dval(&ddx[1]) = dval(&dd[0]);
 104                 dval(&ddx[0]) = dval(&dd[1]);
 105                 dd = ddx;
 106                 L = dd->L;
 107                 }
 108         z = d2b(dval(&dd[0]), &ex, &bx);
 109         if (dval(&dd[1]) == 0.)
 110                 goto no_y;
 111         x = z;
 112         y = d2b(dval(&dd[1]), &ey, &by);
 113         if ( (i = ex - ey) !=0) {
 114                 if (i > 0) {
 115                         x = lshift(x, i);
 116                         ex = ey;
 117                         }
 118                 else
 119                         y = lshift(y, -i);
 120                 }
 121         if ((L[_0] ^ L[2+_0]) & 0x80000000L) {
 122                 z = diff(x, y);
 123                 if (L[_0] & 0x80000000L)
 124                         z->sign = 1 - z->sign;
 125                 }
 126         else {
 127                 z = sum(x, y);
 128                 if (L[_0] & 0x80000000L)
 129                         z->sign = 1;
 130                 }
 131         Bfree(x);
 132         Bfree(y);
 133  no_y:
 134         bits = zx = z->x;
 135         for(i = 0; !*zx; zx++)
 136                 i += 32;
 137         i += lo0bits(zx);
 138         if (i) {
 139                 rshift(z, i);
 140                 ex += i;
 141                 }
 142         fpi.nbits = z->wds * 32 - hi0bits(z->x[j = z->wds-1]);
 143         if (fpi.nbits < 106) {
 144                 fpi.nbits = 106;
 145                 if (j < 3) {
 146                         for(i = 0; i <= j; i++)
 147                                 bits0[i] = bits[i];
 148                         while(i < 4)
 149                                 bits0[i++] = 0;
 150                         bits = bits0;
 151                         }
 152                 }
 153         mode = 2;
 154         if (ndig <= 0) {
 155                 if (bufsize < (int)(fpi.nbits * .301029995664) + 10) {
 156                         Bfree(z);
 157                         return 0;
 158                         }
 159                 mode = 0;
 160                 }
 161         fpi.emin = 1-1023-53+1;
 162         fpi.emax = 2046-1023-106+1;
 163         fpi.rounding = Rounding;
 164         fpi.sudden_underflow = 0;
 165         i = STRTOG_Normal;
 166         s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
 167         b = g__fmt(buf, s, se, decpt, z->sign, bufsize);
 168         Bfree(z);
 169         return b;
 170         }