Replace our strtod with the gdtoa package.
[dragonfly.git] / lib / libc / gdtoa / _hldtoa.c
1 /*-
2  * Copyright (c) 2004-2008 David Schultz <das@FreeBSD.ORG>
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions
7  * are met:
8  * 1. Redistributions of source code must retain the above copyright
9  *    notice, this list of conditions and the following disclaimer.
10  * 2. Redistributions in binary form must reproduce the above copyright
11  *    notice, this list of conditions and the following disclaimer in the
12  *    documentation and/or other materials provided with the distribution.
13  *
14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24  * SUCH DAMAGE.
25  *
26  * $FreeBSD: src/lib/libc/gdtoa/_hldtoa.c,v 1.2 2008/04/12 14:53:52 das Exp $
27  */
28
29 #include <float.h>
30 #include <limits.h>
31 #include <math.h>
32 #include <stdint.h>
33
34 #ifdef __i386__
35 #include <ieeefp.h>
36 #endif
37
38 #include "../stdio/floatio.h"
39 #include "fpmath.h"
40 #include "gdtoaimp.h"
41
42 #if (LDBL_MANT_DIG > DBL_MANT_DIG)
43
44 /* Strings values used by dtoa() */
45 #define INFSTR  "Infinity"
46 #define NANSTR  "NaN"
47
48 #ifdef LDBL_IMPLICIT_NBIT
49 #define MANH_SIZE       LDBL_MANH_SIZE
50 #else
51 #define MANH_SIZE       (LDBL_MANH_SIZE - 1)
52 #endif
53
54 #if MANH_SIZE > 32
55 typedef uint64_t manh_t;
56 #else
57 typedef uint32_t manh_t;
58 #endif
59
60 #if LDBL_MANL_SIZE > 32
61 typedef uint64_t manl_t;
62 #else
63 typedef uint32_t manl_t;
64 #endif
65
66 #define LDBL_ADJ        (LDBL_MAX_EXP - 2)
67 #define SIGFIGS         ((LDBL_MANT_DIG + 3) / 4 + 1)
68
69 static const float one[] = { 1.0f, -1.0f };
70
71 /*
72  * This is the long double version of __hdtoa().
73  */
74 char *
75 __hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
76     char **rve)
77 {
78         union IEEEl2bits u;
79         char *s, *s0;
80         manh_t manh;
81         manl_t manl;
82         int bufsize;
83 #ifdef __i386__
84         fp_prec_t oldprec;
85 #endif
86
87         u.e = e;
88         *sign = u.bits.sign;
89
90         switch (fpclassify(e)) {
91         case FP_NORMAL:
92                 *decpt = u.bits.exp - LDBL_ADJ;
93                 break;
94         case FP_ZERO:
95                 *decpt = 1;
96                 return (nrv_alloc("0", rve, 1));
97         case FP_SUBNORMAL:
98 #ifdef __i386__
99                 oldprec = fpsetprec(FP_PE);
100 #endif
101                 u.e *= 0x1p514L;
102                 *decpt = u.bits.exp - (514 + LDBL_ADJ);
103 #ifdef __i386__
104                 fpsetprec(oldprec);
105 #endif
106                 break;
107         case FP_INFINITE:
108                 *decpt = INT_MAX;
109                 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
110         default:        /* FP_NAN or unrecognized */
111                 *decpt = INT_MAX;
112                 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
113         }
114
115         /* FP_NORMAL or FP_SUBNORMAL */
116
117         if (ndigits == 0)               /* dtoa() compatibility */
118                 ndigits = 1;
119
120         /*
121          * If ndigits < 0, we are expected to auto-size, so we allocate
122          * enough space for all the digits.
123          */
124         bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
125         s0 = rv_alloc(bufsize);
126
127         /* Round to the desired number of digits. */
128         if (SIGFIGS > ndigits && ndigits > 0) {
129                 float redux = one[u.bits.sign];
130                 int offset = 4 * ndigits + LDBL_MAX_EXP - 4 - LDBL_MANT_DIG;
131 #ifdef __i386__
132                 oldprec = fpsetprec(FP_PE);
133 #endif
134                 u.bits.exp = offset;
135                 u.e += redux;
136                 u.e -= redux;
137                 *decpt += u.bits.exp - offset;
138 #ifdef __i386__
139                 fpsetprec(oldprec);
140 #endif
141         }
142
143         mask_nbit_l(u);
144         manh = u.bits.manh;
145         manl = u.bits.manl;
146         *s0 = '1';
147         for (s = s0 + 1; s < s0 + bufsize; s++) {
148                 *s = xdigs[(manh >> (MANH_SIZE - 4)) & 0xf];
149                 manh = (manh << 4) | (manl >> (LDBL_MANL_SIZE - 4));
150                 manl <<= 4;
151         }
152
153         /* If ndigits < 0, we are expected to auto-size the precision. */
154         if (ndigits < 0) {
155                 for (ndigits = SIGFIGS; s0[ndigits - 1] == '0'; ndigits--)
156                         ;
157         }
158
159         s = s0 + ndigits;
160         *s = '\0';
161         if (rve != NULL)
162                 *rve = s;
163         return (s0);
164 }
165
166 #else   /* (LDBL_MANT_DIG == DBL_MANT_DIG) */
167
168 char *
169 __hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
170     char **rve)
171 {
172
173         return (__hdtoa((double)e, xdigs, ndigits, decpt, sign, rve));
174 }
175
176 #endif  /* (LDBL_MANT_DIG == DBL_MANT_DIG) */