1/*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22/*
23 * Copyright 2008 Sun Microsystems, Inc.  All rights reserved.
24 * Use is subject to license terms.
25 */
26
27#pragma ident	"%Z%%M%	%I%	%E% SMI"
28
29/*
30 * Short cut for conversion from double precision to decimal
31 * floating point
32 */
33
34#include "lint.h"
35#include <sys/types.h>
36#include <sys/isa_defs.h>
37#include "base_conversion.h"
38
39/*
40 * Powers of ten rounded up.  If i is the largest index such that
41 * tbl_decade[i] <= x, then:
42 *
43 *  if i == 0 then x < 10^-49
44 *  else if i == TBL_DECADE_MAX then x >= 10^67
45 *  else 10^(i-TBL_DECADE_OFFSET) <= x < 10^(i-TBL_DECADE_OFFSET+1)
46 */
47
48#define	TBL_DECADE_OFFSET	50
49#define	TBL_DECADE_MAX		117
50
51static const double tbl_decade[TBL_DECADE_MAX + 1] = {
52	0.0,
53	1.00000000000000012631e-49, 1.00000000000000012631e-48,
54	1.00000000000000009593e-47, 1.00000000000000002300e-46,
55	1.00000000000000013968e-45, 1.00000000000000007745e-44,
56	1.00000000000000007745e-43, 1.00000000000000003762e-42,
57	1.00000000000000000576e-41, 1.00000000000000013321e-40,
58	1.00000000000000009243e-39, 1.00000000000000009243e-38,
59	1.00000000000000006632e-37, 1.00000000000000010809e-36,
60	1.00000000000000000786e-35, 1.00000000000000014150e-34,
61	1.00000000000000005597e-33, 1.00000000000000005597e-32,
62	1.00000000000000008334e-31, 1.00000000000000008334e-30,
63	1.00000000000000008334e-29, 1.00000000000000008334e-28,
64	1.00000000000000003849e-27, 1.00000000000000003849e-26,
65	1.00000000000000003849e-25, 1.00000000000000010737e-24,
66	1.00000000000000010737e-23, 1.00000000000000004860e-22,
67	1.00000000000000009562e-21, 1.00000000000000009562e-20,
68	1.00000000000000009562e-19, 1.00000000000000007154e-18,
69	1.00000000000000007154e-17, 1.00000000000000010236e-16,
70	1.00000000000000007771e-15, 1.00000000000000015659e-14,
71	1.00000000000000003037e-13, 1.00000000000000018184e-12,
72	1.00000000000000010106e-11, 1.00000000000000003643e-10,
73	1.00000000000000006228e-09, 1.00000000000000002092e-08,
74	1.00000000000000008710e-07, 1.00000000000000016651e-06,
75	1.00000000000000008180e-05, 1.00000000000000004792e-04,
76	1.00000000000000002082e-03, 1.00000000000000002082e-02,
77	1.00000000000000005551e-01, 1.00000000000000000000e+00,
78	1.00000000000000000000e+01, 1.00000000000000000000e+02,
79	1.00000000000000000000e+03, 1.00000000000000000000e+04,
80	1.00000000000000000000e+05, 1.00000000000000000000e+06,
81	1.00000000000000000000e+07, 1.00000000000000000000e+08,
82	1.00000000000000000000e+09, 1.00000000000000000000e+10,
83	1.00000000000000000000e+11, 1.00000000000000000000e+12,
84	1.00000000000000000000e+13, 1.00000000000000000000e+14,
85	1.00000000000000000000e+15, 1.00000000000000000000e+16,
86	1.00000000000000000000e+17, 1.00000000000000000000e+18,
87	1.00000000000000000000e+19, 1.00000000000000000000e+20,
88	1.00000000000000000000e+21, 1.00000000000000000000e+22,
89	1.00000000000000008389e+23, 1.00000000000000011744e+24,
90	1.00000000000000009060e+25, 1.00000000000000004765e+26,
91	1.00000000000000001329e+27, 1.00000000000000017821e+28,
92	1.00000000000000009025e+29, 1.00000000000000001988e+30,
93	1.00000000000000007618e+31, 1.00000000000000005366e+32,
94	1.00000000000000008969e+33, 1.00000000000000006087e+34,
95	1.00000000000000015310e+35, 1.00000000000000004242e+36,
96	1.00000000000000007194e+37, 1.00000000000000016638e+38,
97	1.00000000000000009082e+39, 1.00000000000000003038e+40,
98	1.00000000000000000620e+41, 1.00000000000000004489e+42,
99	1.00000000000000001394e+43, 1.00000000000000008821e+44,
100	1.00000000000000008821e+45, 1.00000000000000011990e+46,
101	1.00000000000000004385e+47, 1.00000000000000004385e+48,
102	1.00000000000000007630e+49, 1.00000000000000007630e+50,
103	1.00000000000000015937e+51, 1.00000000000000012614e+52,
104	1.00000000000000020590e+53, 1.00000000000000007829e+54,
105	1.00000000000000001024e+55, 1.00000000000000009190e+56,
106	1.00000000000000004835e+57, 1.00000000000000008319e+58,
107	1.00000000000000008319e+59, 1.00000000000000012779e+60,
108	1.00000000000000009211e+61, 1.00000000000000003502e+62,
109	1.00000000000000005786e+63, 1.00000000000000002132e+64,
110	1.00000000000000010901e+65, 1.00000000000000013239e+66,
111	1.00000000000000013239e+67
112};
113
114/*
115 * Convert a positive double precision integer x <= 2147483647999999744
116 * (the largest double less than 2^31 * 10^9; this implementation works
117 * up to the largest double less than 2^25 * 10^12) to a string of ASCII
118 * decimal digits, adding leading zeroes so that the result has at least
119 * n digits.  The string is terminated by a null byte, and its length
120 * is returned.
121 *
122 * This routine assumes round-to-nearest mode is in effect and any
123 * exceptions raised will be ignored.
124 */
125
126#define	tenm4	tbl_decade[TBL_DECADE_OFFSET - 4]
127#define	ten4	tbl_decade[TBL_DECADE_OFFSET + 4]
128#define	tenm12	tbl_decade[TBL_DECADE_OFFSET - 12]
129#define	ten12	tbl_decade[TBL_DECADE_OFFSET + 12]
130#define	one	tbl_decade[TBL_DECADE_OFFSET]
131
132static int
133__double_to_digits(double x, char *s, int n)
134{
135	double		y;
136	int		d[5], i, j;
137	char		*ss, tmp[4];
138
139	/* decompose x into four-digit chunks */
140	y = (int)(x * tenm12);
141	x -= y * ten12;
142	if (x < 0.0) {
143		y -= one;
144		x += ten12;
145	}
146	d[0] = (int)(y * tenm4);
147	d[1] = (int)(y - d[0] * ten4);
148	y = (int)(x * tenm4);
149	d[4] = (int)(x - y * ten4);
150	d[2] = (int)(y * tenm4);
151	d[3] = (int)(y - d[2] * ten4);
152
153	/*
154	 * Find the first nonzero chunk or the point at which to start
155	 * converting so we have n digits, whichever comes first.
156	 */
157	ss = s;
158	if (n > 20) {
159		for (j = 0; j < n - 20; j++)
160			*ss++ = '0';
161		i = 0;
162	} else {
163		for (i = 0; d[i] == 0 && n <= ((4 - i) << 2); i++)
164			;
165		__four_digits_quick(d[i], tmp);
166		for (j = 0; tmp[j] == '0' && n <= ((4 - i) << 2) + 3 - j; j++)
167			;
168		while (j < 4)
169			*ss++ = tmp[j++];
170		i++;
171	}
172
173	/* continue converting four-digit chunks */
174	while (i < 5) {
175		__four_digits_quick(d[i], ss);
176		ss += 4;
177		i++;
178	}
179
180	*ss = '\0';
181	return (ss - s);
182}
183
184/*
185 * Round a positive double precision number *x to the nearest integer,
186 * returning the result and passing back an indication of accuracy in
187 * *pe.  On entry, nrx is the number of rounding errors already com-
188 * mitted in forming *x.  On exit, *pe is 0 if *x was already integral
189 * and exact, 1 if the result is the correctly rounded integer value
190 * but not exact, and 2 if error in *x precludes determining the cor-
191 * rectly rounded integer value (i.e., the error might be larger than
192 * 1/2 - |*x - rx|, where rx is the nearest integer to *x).
193 */
194
195static union {
196	unsigned int	i[2];
197	double		d;
198} C[] = {
199#ifdef _LITTLE_ENDIAN
200	{ 0x00000000, 0x43300000 },
201	{ 0x00000000, 0x3ca00000 },
202	{ 0x00000000, 0x3fe00000 },
203	{ 0xffffffff, 0x3fdfffff },
204#else
205	{ 0x43300000, 0x00000000 },
206	{ 0x3ca00000, 0x00000000 },
207	{ 0x3fe00000, 0x00000000 },
208	{ 0x3fdfffff, 0xffffffff },	/* nextafter(1/2, 0) */
209#endif
210};
211
212#define	two52	C[0].d
213#define	twom53	C[1].d
214#define	half	C[2].d
215#define	halfdec	C[3].d
216
217static double
218__arint_set_n(double *x, int nrx, int *pe)
219{
220	int	hx;
221	double	rx, rmx;
222
223#ifdef _LITTLE_ENDIAN
224	hx = *(1+(int *)x);
225#else
226	hx = *(int *)x;
227#endif
228	if (hx >= 0x43300000) {
229		/* x >= 2^52, so it's already integral */
230		if (nrx == 0)
231			*pe = 0;
232		else if (nrx == 1 && hx < 0x43400000)
233			*pe = 1;
234		else
235			*pe = 2;
236		return (*x);
237	} else if (hx < 0x3fe00000) {
238		/* x < 1/2 */
239		if (nrx > 1 && hx == 0x3fdfffff)
240			*pe = (*x == halfdec)? 2 : 1;
241		else
242			*pe = 1;
243		return (0.0);
244	}
245
246	rx = (*x + two52) - two52;
247	if (nrx == 0) {
248		*pe = (rx == *x)? 0 : 1;
249	} else {
250		rmx = rx - *x;
251		if (rmx < 0.0)
252			rmx = -rmx;
253		*pe = (nrx * twom53 * *x < half - rmx)? 1 : 2;
254	}
255	return (rx);
256}
257
258/*
259 * Attempt to convert dd to a decimal record *pd according to the
260 * modes in *pm using double precision floating point.  Return zero
261 * and sets *ps to reflect any exceptions incurred if successful.
262 * Return a nonzero value if unsuccessful.
263 */
264int
265__fast_double_to_decimal(double *dd, decimal_mode *pm, decimal_record *pd,
266    fp_exception_field_type *ps)
267{
268	int			i, is, esum, eround, hd;
269	double			dds;
270	__ieee_flags_type	fb;
271
272	if (pm->rd != fp_nearest)
273		return (1);
274
275	if (pm->df == fixed_form) {
276		/* F format */
277		if (pm->ndigits < 0 || pm->ndigits > __TBL_TENS_MAX)
278			return (1);
279		__get_ieee_flags(&fb);
280		dds = __dabs(dd);
281		esum = 0;
282		if (pm->ndigits) {
283			/* scale by a positive power of ten */
284			if (pm->ndigits > __TBL_TENS_EXACT) {
285				dds *= __tbl_tens[pm->ndigits];
286				esum = 2;
287			} else {
288				dds = __mul_set(dds, __tbl_tens[pm->ndigits],
289				    &eround);
290				esum = eround;
291			}
292		}
293		if (dds > 2147483647999999744.0) {
294			__set_ieee_flags(&fb);
295			return (1);
296		}
297		dds = __arint_set_n(&dds, esum, &eround);
298		if (eround == 2) {
299			/* error is too large to round reliably; punt */
300			__set_ieee_flags(&fb);
301			return (1);
302		}
303		if (dds == 0.0) {
304			is = (pm->ndigits > 0)? pm->ndigits : 1;
305			for (i = 0; i < is; i++)
306				pd->ds[i] = '0';
307			pd->ds[is] = '\0';
308			eround++;
309		} else {
310			is = __double_to_digits(dds, pd->ds, pm->ndigits);
311		}
312		pd->ndigits = is;
313		pd->exponent = -pm->ndigits;
314	} else {
315		/* E format */
316		if (pm->ndigits < 1 || pm->ndigits > 18)
317			return (1);
318		__get_ieee_flags(&fb);
319		dds = __dabs(dd);
320		/* find the decade containing dds */
321#ifdef _LITTLE_ENDIAN
322		hd = *(1+(int *)dd);
323#else
324		hd = *(int *)dd;
325#endif
326		hd = (hd >> 20) & 0x7ff;
327		if (hd >= 0x400) {
328			if (hd > 0x4e0)
329				i = TBL_DECADE_MAX;
330			else
331				i = TBL_DECADE_MAX - ((0x4e0 - hd) >> 2);
332		} else {
333			if (hd < 0x358)
334				i = 0;
335			else
336				i = TBL_DECADE_OFFSET - ((0x3ff - hd) >> 2);
337		}
338		while (dds < tbl_decade[i])
339			i--;
340		/* determine the power of ten by which to scale */
341		i = pm->ndigits - 1 - (i - TBL_DECADE_OFFSET);
342		esum = 0;
343		if (i > 0) {
344			/* scale by a positive power of ten */
345			if (i > __TBL_TENS_EXACT) {
346				if (i > __TBL_TENS_MAX) {
347					__set_ieee_flags(&fb);
348					return (1);
349				}
350				dds *= __tbl_tens[i];
351				esum = 2;
352			} else {
353				dds = __mul_set(dds, __tbl_tens[i], &eround);
354				esum = eround;
355			}
356		} else if (i < 0) {
357			/* scale by a negative power of ten */
358			if (-i > __TBL_TENS_EXACT) {
359				if (-i > __TBL_TENS_MAX) {
360					__set_ieee_flags(&fb);
361					return (1);
362				}
363				dds /= __tbl_tens[-i];
364				esum = 2;
365			} else {
366				dds = __div_set(dds, __tbl_tens[-i], &eround);
367				esum = eround;
368			}
369		}
370		dds = __arint_set_n(&dds, esum, &eround);
371		if (eround == 2) {
372			/* error is too large to round reliably; punt */
373			__set_ieee_flags(&fb);
374			return (1);
375		}
376		is = __double_to_digits(dds, pd->ds, 1);
377		if (is > pm->ndigits) {
378			/*
379			 * The result rounded up to the next larger power
380			 * of ten; just discard the last zero and adjust
381			 * the exponent.
382			 */
383			pd->ds[--is] = '\0';
384			i--;
385		}
386		pd->ndigits = is;
387		pd->exponent = -i;
388	}
389	*ps = (eround == 0)? 0 : (1 << fp_inexact);
390	__set_ieee_flags(&fb);
391	return (0);
392}
393