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