xref: /illumos-gate/usr/src/lib/libc/port/fp/__flt_decim.c (revision 7257d1b4)
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 
51 static 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 
132 static 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 
195 static 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 
217 static 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  */
264 int
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