1/* $NetBSD: g_ddfmt.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */ 2 3/**************************************************************** 4 5The author of this software is David M. Gay. 6 7Copyright (C) 1998 by Lucent Technologies 8All Rights Reserved 9 10Permission to use, copy, modify, and distribute this software and 11its documentation for any purpose and without fee is hereby 12granted, provided that the above copyright notice appear in all 13copies and that both that the copyright notice and this 14permission notice and warranty disclaimer appear in supporting 15documentation, and that the name of Lucent or any of its entities 16not be used in advertising or publicity pertaining to 17distribution of the software without specific, written prior 18permission. 19 20LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 21INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 22IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 23SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 24WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 25IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 26ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 27THIS 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 38g_ddfmt(buf, dd0, ndig, bufsize) char *buf; double *dd0; int ndig; size_t bufsize; 39#else 40g_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