usr/src/lib/libc/port/fp/__flt_decim.c

   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 }