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