00001
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 #define _GNU_SOURCE
00050 #include "compat/fpclassify.h"
00051
00052 #include <sys/cdefs.h>
00053 #define __FBSDID(x)
00054 __FBSDID("$FreeBSD: src/lib/libc/gdtoa/_hdtoa.c,v 1.3 2005/01/18 18:44:07 das Exp $");
00055
00056 #include <float.h>
00057 #include <limits.h>
00058 #include <math.h>
00059 #include "fpmath.h"
00060 #include "gdtoaimp.h"
00061
00062
00063 #define INFSTR "Infinity"
00064 #define NANSTR "NaN"
00065
00066 #define DBL_ADJ (DBL_MAX_EXP - 2 + ((DBL_MANT_DIG - 1) % 4))
00067 #define LDBL_ADJ (LDBL_MAX_EXP - 2 + ((LDBL_MANT_DIG - 1) % 4))
00068
00069
00070
00071
00072
00073
00074 static int
00075 roundup(char *s0, int ndigits)
00076 {
00077 char *s;
00078
00079 for (s = s0 + ndigits - 1; *s == 0xf; s--) {
00080 if (s == s0) {
00081 *s = 1;
00082 return (1);
00083 }
00084 ++*s;
00085 }
00086 ++*s;
00087 return (0);
00088 }
00089
00090
00091
00092
00093
00094
00095
00096 static void
00097 dorounding(char *s0, int ndigits, int sign, int *decpt)
00098 {
00099 int adjust = 0;
00100
00101 switch (FLT_ROUNDS) {
00102 case 0:
00103 default:
00104 break;
00105 case 1:
00106 if ((s0[ndigits] > 8) ||
00107 (s0[ndigits] == 8 && s0[ndigits - 1] & 1))
00108 adjust = roundup(s0, ndigits);
00109 break;
00110 case 2:
00111 if (sign == 0)
00112 adjust = roundup(s0, ndigits);
00113 break;
00114 case 3:
00115 if (sign != 0)
00116 adjust = roundup(s0, ndigits);
00117 break;
00118 }
00119
00120 if (adjust)
00121 *decpt += 4;
00122 }
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148 char *
00149 __hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
00150 char **rve)
00151 {
00152 static const int sigfigs = (DBL_MANT_DIG + 3) / 4;
00153 union IEEEd2bits u;
00154 char *s, *s0;
00155 int bufsize;
00156
00157 u.d = d;
00158 *sign = u.bits.sign;
00159
00160 switch (fpclassify(d)) {
00161 case FP_NORMAL:
00162 *decpt = u.bits.exp - DBL_ADJ;
00163 break;
00164 case FP_ZERO:
00165 *decpt = 1;
00166 return (nrv_alloc("0", rve, 1));
00167 case FP_SUBNORMAL:
00168 u.d *= 0x1p514;
00169 *decpt = u.bits.exp - (514 + DBL_ADJ);
00170 break;
00171 case FP_INFINITE:
00172 *decpt = INT_MAX;
00173 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
00174 case FP_NAN:
00175 *decpt = INT_MAX;
00176 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
00177 default:
00178 abort();
00179 }
00180
00181
00182
00183 if (ndigits == 0)
00184 ndigits = 1;
00185
00186
00187
00188
00189
00190 bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
00191 s0 = rv_alloc(bufsize);
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
00202 *s = 0;
00203 for (; s > s0 + sigfigs - (DBL_MANL_SIZE / 4) - 1 && s > s0; s--) {
00204 *s = u.bits.manl & 0xf;
00205 u.bits.manl >>= 4;
00206 }
00207 for (; s > s0; s--) {
00208 *s = u.bits.manh & 0xf;
00209 u.bits.manh >>= 4;
00210 }
00211
00212
00213
00214
00215
00216
00217
00218 *s = u.bits.manh | (1U << ((DBL_MANT_DIG - 1) % 4));
00219
00220
00221 if (ndigits < 0) {
00222 for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
00223 ;
00224 }
00225
00226 if (sigfigs > ndigits && s0[ndigits] != 0)
00227 dorounding(s0, ndigits, u.bits.sign, decpt);
00228
00229 s = s0 + ndigits;
00230 if (rve != NULL)
00231 *rve = s;
00232 *s-- = '\0';
00233 for (; s >= s0; s--)
00234 *s = xdigs[(unsigned int)*s];
00235
00236 return (s0);
00237 }
00238
00239 #if (LDBL_MANT_DIG > DBL_MANT_DIG)
00240
00241
00242
00243
00244 char *
00245 __hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
00246 char **rve)
00247 {
00248 static const int sigfigs = (LDBL_MANT_DIG + 3) / 4;
00249 union IEEEl2bits u;
00250 char *s, *s0;
00251 int bufsize;
00252
00253 u.e = e;
00254 *sign = u.bits.sign;
00255
00256 switch (fpclassify(e)) {
00257 case FP_NORMAL:
00258 *decpt = u.bits.exp - LDBL_ADJ;
00259 break;
00260 case FP_ZERO:
00261 *decpt = 1;
00262 return (nrv_alloc("0", rve, 1));
00263 case FP_SUBNORMAL:
00264 u.e *= 0x1p514L;
00265 *decpt = u.bits.exp - (514 + LDBL_ADJ);
00266 break;
00267 case FP_INFINITE:
00268 *decpt = INT_MAX;
00269 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
00270 case FP_NAN:
00271 *decpt = INT_MAX;
00272 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
00273 default:
00274 abort();
00275 }
00276
00277
00278
00279 if (ndigits == 0)
00280 ndigits = 1;
00281
00282
00283
00284
00285
00286 bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
00287 s0 = rv_alloc(bufsize);
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
00298 *s = 0;
00299 for (; s > s0 + sigfigs - (LDBL_MANL_SIZE / 4) - 1 && s > s0; s--) {
00300 *s = u.bits.manl & 0xf;
00301 u.bits.manl >>= 4;
00302 }
00303 for (; s > s0; s--) {
00304 *s = u.bits.manh & 0xf;
00305 u.bits.manh >>= 4;
00306 }
00307
00308
00309
00310
00311
00312
00313
00314 *s = u.bits.manh | (1U << ((LDBL_MANT_DIG - 1) % 4));
00315
00316
00317 if (ndigits < 0) {
00318 for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
00319 ;
00320 }
00321
00322 if (sigfigs > ndigits && s0[ndigits] != 0)
00323 dorounding(s0, ndigits, u.bits.sign, decpt);
00324
00325 s = s0 + ndigits;
00326 if (rve != NULL)
00327 *rve = s;
00328 *s-- = '\0';
00329 for (; s >= s0; s--)
00330 *s = xdigs[(unsigned int)*s];
00331
00332 return (s0);
00333 }
00334
00335 #else
00336
00337 char *
00338 __hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
00339 char **rve)
00340 {
00341
00342 return (__hdtoa((double)e, xdigs, ndigits, decpt, sign, rve));
00343 }
00344
00345 #endif