← Back to branch summary

~ubuntu-branches/ubuntu/trusty/mysql-5.6/trusty

« back to all changes in this revision

Viewing changes to strings/dtoa.c

  • Committer: Package Import Robot
  • Author(s): James Page
  • Date: 2014-02-12 11:54:27 UTC
  • Revision ID: package-import@ubuntu.com-20140212115427-oq6tfsqxl1wuwehi
Tags: upstream-5.6.15
Import upstream version 5.6.15

Show diffs side-by-side

added added

removed removed

Lines of Context:
strings/dtoa.c
 
1
/* Copyright (c) 2007, 2010, Oracle and/or its affiliates. All rights reserved.
 
2
 
 
3
   This library is free software; you can redistribute it and/or
 
4
   modify it under the terms of the GNU Library General Public
 
5
   License as published by the Free Software Foundation; version 2
 
6
   of the License.
 
7
 
 
8
   This program is distributed in the hope that it will be useful,
 
9
   but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
   GNU General Public License for more details.
 
12
 
 
13
   You should have received a copy of the GNU General Public License
 
14
   along with this program; if not, write to the Free Software
 
15
   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301  USA */
 
16
 
 
17
/****************************************************************
 
18
 
 
19
  This file incorporates work covered by the following copyright and
 
20
  permission notice:
 
21
 
 
22
  The author of this software is David M. Gay.
 
23
 
 
24
  Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
 
25
 
 
26
  Permission to use, copy, modify, and distribute this software for any
 
27
  purpose without fee is hereby granted, provided that this entire notice
 
28
  is included in all copies of any software which is or includes a copy
 
29
  or modification of this software and in all copies of the supporting
 
30
  documentation for such software.
 
31
 
 
32
  THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
 
33
  WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
 
34
  REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
 
35
  OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
 
36
 
 
37
 ***************************************************************/
 
38
 
 
39
#include <my_base.h> /* for EOVERFLOW on Windows */
 
40
#include <my_global.h>
 
41
#include <m_string.h>  /* for memcpy and NOT_FIXED_DEC */
 
42
 
 
43
/**
 
44
   Appears to suffice to not call malloc() in most cases.
 
45
   @todo
 
46
     see if it is possible to get rid of malloc().
 
47
     this constant is sufficient to avoid malloc() on all inputs I have tried.
 
48
*/
 
49
#define DTOA_BUFF_SIZE (460 * sizeof(void *))
 
50
 
 
51
/* Magic value returned by dtoa() to indicate overflow */
 
52
#define DTOA_OVERFLOW 9999
 
53
 
 
54
static double my_strtod_int(const char *, char **, int *, char *, size_t);
 
55
static char *dtoa(double, int, int, int *, int *, char **, char *, size_t);
 
56
static void dtoa_free(char *, char *, size_t);
 
57
 
 
58
/**
 
59
   @brief
 
60
   Converts a given floating point number to a zero-terminated string
 
61
   representation using the 'f' format.
 
62
 
 
63
   @details
 
64
   This function is a wrapper around dtoa() to do the same as
 
65
   sprintf(to, "%-.*f", precision, x), though the conversion is usually more
 
66
   precise. The only difference is in handling [-,+]infinity and nan values,
 
67
   in which case we print '0\0' to the output string and indicate an overflow.
 
68
 
 
69
   @param x           the input floating point number.
 
70
   @param precision   the number of digits after the decimal point.
 
71
                      All properties of sprintf() apply:
 
72
                      - if the number of significant digits after the decimal
 
73
                        point is less than precision, the resulting string is
 
74
                        right-padded with zeros
 
75
                      - if the precision is 0, no decimal point appears
 
76
                      - if a decimal point appears, at least one digit appears
 
77
                        before it
 
78
   @param to          pointer to the output buffer. The longest string which
 
79
                      my_fcvt() can return is FLOATING_POINT_BUFFER bytes
 
80
                      (including the terminating '\0').
 
81
   @param error       if not NULL, points to a location where the status of
 
82
                      conversion is stored upon return.
 
83
                      FALSE  successful conversion
 
84
                      TRUE   the input number is [-,+]infinity or nan.
 
85
                             The output string in this case is always '0'.
 
86
   @return            number of written characters (excluding terminating '\0')
 
87
*/
 
88
 
 
89
size_t my_fcvt(double x, int precision, char *to, my_bool *error)
 
90
{
 
91
  int decpt, sign, len, i;
 
92
  char *res, *src, *end, *dst= to;
 
93
  char buf[DTOA_BUFF_SIZE];
 
94
  DBUG_ASSERT(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
 
95
  
 
96
  res= dtoa(x, 5, precision, &decpt, &sign, &end, buf, sizeof(buf));
 
97
 
 
98
  if (decpt == DTOA_OVERFLOW)
 
99
  {
 
100
    dtoa_free(res, buf, sizeof(buf));
 
101
    *to++= '0';
 
102
    *to= '\0';
 
103
    if (error != NULL)
 
104
      *error= TRUE;
 
105
    return 1;
 
106
  }
 
107
 
 
108
  src= res;
 
109
  len= end - src;
 
110
 
 
111
  if (sign)
 
112
    *dst++= '-';
 
113
 
 
114
  if (decpt <= 0)
 
115
  {
 
116
    *dst++= '0';
 
117
    *dst++= '.';
 
118
    for (i= decpt; i < 0; i++)
 
119
      *dst++= '0';
 
120
  }
 
121
 
 
122
  for (i= 1; i <= len; i++)
 
123
  {
 
124
    *dst++= *src++;
 
125
    if (i == decpt && i < len)
 
126
      *dst++= '.';
 
127
  }
 
128
  while (i++ <= decpt)
 
129
    *dst++= '0';
 
130
 
 
131
  if (precision > 0)
 
132
  {
 
133
    if (len <= decpt)
 
134
      *dst++= '.';
 
135
    
 
136
    for (i= precision - MY_MAX(0, (len - decpt)); i > 0; i--)
 
137
      *dst++= '0';
 
138
  }
 
139
  
 
140
  *dst= '\0';
 
141
  if (error != NULL)
 
142
    *error= FALSE;
 
143
 
 
144
  dtoa_free(res, buf, sizeof(buf));
 
145
 
 
146
  return dst - to;
 
147
}
 
148
 
 
149
/**
 
150
   @brief
 
151
   Converts a given floating point number to a zero-terminated string
 
152
   representation with a given field width using the 'e' format
 
153
   (aka scientific notation) or the 'f' one.
 
154
 
 
155
   @details
 
156
   The format is chosen automatically to provide the most number of significant
 
157
   digits (and thus, precision) with a given field width. In many cases, the
 
158
   result is similar to that of sprintf(to, "%g", x) with a few notable
 
159
   differences:
 
160
   - the conversion is usually more precise than C library functions.
 
161
   - there is no 'precision' argument. instead, we specify the number of
 
162
     characters available for conversion (i.e. a field width).
 
163
   - the result never exceeds the specified field width. If the field is too
 
164
     short to contain even a rounded decimal representation, my_gcvt()
 
165
     indicates overflow and truncates the output string to the specified width.
 
166
   - float-type arguments are handled differently than double ones. For a
 
167
     float input number (i.e. when the 'type' argument is MY_GCVT_ARG_FLOAT)
 
168
     we deliberately limit the precision of conversion by FLT_DIG digits to
 
169
     avoid garbage past the significant digits.
 
170
   - unlike sprintf(), in cases where the 'e' format is preferred,  we don't
 
171
     zero-pad the exponent to save space for significant digits. The '+' sign
 
172
     for a positive exponent does not appear for the same reason.
 
173
 
 
174
   @param x           the input floating point number.
 
175
   @param type        is either MY_GCVT_ARG_FLOAT or MY_GCVT_ARG_DOUBLE.
 
176
                      Specifies the type of the input number (see notes above).
 
177
   @param width       field width in characters. The minimal field width to
 
178
                      hold any number representation (albeit rounded) is 7
 
179
                      characters ("-Ne-NNN").
 
180
   @param to          pointer to the output buffer. The result is always
 
181
                      zero-terminated, and the longest returned string is thus
 
182
                      'width + 1' bytes.
 
183
   @param error       if not NULL, points to a location where the status of
 
184
                      conversion is stored upon return.
 
185
                      FALSE  successful conversion
 
186
                      TRUE   the input number is [-,+]infinity or nan.
 
187
                             The output string in this case is always '0'.
 
188
   @return            number of written characters (excluding terminating '\0')
 
189
 
 
190
   @todo
 
191
   Check if it is possible and  makes sense to do our own rounding on top of
 
192
   dtoa() instead of calling dtoa() twice in (rare) cases when the resulting
 
193
   string representation does not fit in the specified field width and we want
 
194
   to re-round the input number with fewer significant digits. Examples:
 
195
 
 
196
     my_gcvt(-9e-3, ..., 4, ...);
 
197
     my_gcvt(-9e-3, ..., 2, ...);
 
198
     my_gcvt(1.87e-3, ..., 4, ...);
 
199
     my_gcvt(55, ..., 1, ...);
 
200
 
 
201
   We do our best to minimize such cases by:
 
202
   
 
203
   - passing to dtoa() the field width as the number of significant digits
 
204
   
 
205
   - removing the sign of the number early (and decreasing the width before
 
206
     passing it to dtoa())
 
207
   
 
208
   - choosing the proper format to preserve the most number of significant
 
209
     digits.
 
210
*/
 
211
 
 
212
size_t my_gcvt(double x, my_gcvt_arg_type type, int width, char *to,
 
213
               my_bool *error)
 
214
{
 
215
  int decpt, sign, len, exp_len;
 
216
  char *res, *src, *end, *dst= to, *dend= dst + width;
 
217
  char buf[DTOA_BUFF_SIZE];
 
218
  my_bool have_space, force_e_format;
 
219
  DBUG_ASSERT(width > 0 && to != NULL);
 
220
  
 
221
  /* We want to remove '-' from equations early */
 
222
  if (x < 0.)
 
223
    width--;
 
224
 
 
225
  res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? width : MY_MIN(width, FLT_DIG),
 
226
            &decpt, &sign, &end, buf, sizeof(buf));
 
227
  if (decpt == DTOA_OVERFLOW)
 
228
  {
 
229
    dtoa_free(res, buf, sizeof(buf));
 
230
    *to++= '0';
 
231
    *to= '\0';
 
232
    if (error != NULL)
 
233
      *error= TRUE;
 
234
    return 1;
 
235
  }
 
236
 
 
237
  if (error != NULL)
 
238
    *error= FALSE;
 
239
 
 
240
  src= res;
 
241
  len= end - res;
 
242
 
 
243
  /*
 
244
    Number of digits in the exponent from the 'e' conversion.
 
245
     The sign of the exponent is taken into account separetely, we don't need
 
246
     to count it here.
 
247
   */
 
248
  exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
 
249
  
 
250
  /*
 
251
     Do we have enough space for all digits in the 'f' format?
 
252
     Let 'len' be the number of significant digits returned by dtoa,
 
253
     and F be the length of the resulting decimal representation.
 
254
     Consider the following cases:
 
255
     1. decpt <= 0, i.e. we have "0.NNN" => F = len - decpt + 2
 
256
     2. 0 < decpt < len, i.e. we have "NNN.NNN" => F = len + 1
 
257
     3. len <= decpt, i.e. we have "NNN00" => F = decpt
 
258
  */
 
259
  have_space= (decpt <= 0 ? len - decpt + 2 :
 
260
               decpt > 0 && decpt < len ? len + 1 :
 
261
               decpt) <= width;
 
262
  /*
 
263
    The following is true when no significant digits can be placed with the
 
264
    specified field width using the 'f' format, and the 'e' format
 
265
    will not be truncated.
 
266
  */
 
267
  force_e_format= (decpt <= 0 && width <= 2 - decpt && width >= 3 + exp_len);
 
268
  /*
 
269
    Assume that we don't have enough space to place all significant digits in
 
270
    the 'f' format. We have to choose between the 'e' format and the 'f' one
 
271
    to keep as many significant digits as possible.
 
272
    Let E and F be the lengths of decimal representaion in the 'e' and 'f'
 
273
    formats, respectively. We want to use the 'f' format if, and only if F <= E.
 
274
    Consider the following cases:
 
275
    1. decpt <= 0.
 
276
       F = len - decpt + 2 (see above)
 
277
       E = len + (len > 1) + 1 + 1 (decpt <= -99) + (decpt <= -9) + 1
 
278
       ("N.NNe-MMM")
 
279
       (F <= E) <=> (len == 1 && decpt >= -1) || (len > 1 && decpt >= -2)
 
280
       We also need to ensure that if the 'f' format is chosen,
 
281
       the field width allows us to place at least one significant digit
 
282
       (i.e. width > 2 - decpt). If not, we prefer the 'e' format.
 
283
    2. 0 < decpt < len
 
284
       F = len + 1 (see above)
 
285
       E = len + 1 + 1 + ... ("N.NNeMMM")
 
286
       F is always less than E.
 
287
    3. len <= decpt <= width
 
288
       In this case we have enough space to represent the number in the 'f'
 
289
       format, so we prefer it with some exceptions.
 
290
    4. width < decpt
 
291
       The number cannot be represented in the 'f' format at all, always use
 
292
       the 'e' 'one.
 
293
  */
 
294
  if ((have_space ||
 
295
      /*
 
296
        Not enough space, let's see if the 'f' format provides the most number
 
297
        of significant digits.
 
298
      */
 
299
       ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
 
300
                                            (len > 1 || !force_e_format)))) &&
 
301
         !force_e_format)) &&
 
302
      
 
303
       /*
 
304
         Use the 'e' format in some cases even if we have enough space for the
 
305
         'f' one. See comment for MAX_DECPT_FOR_F_FORMAT.
 
306
       */
 
307
      (!have_space || (decpt >= -MAX_DECPT_FOR_F_FORMAT + 1 &&
 
308
                       (decpt <= MAX_DECPT_FOR_F_FORMAT || len > decpt))))
 
309
  {
 
310
    /* 'f' format */
 
311
    int i;
 
312
 
 
313
    width-= (decpt < len) + (decpt <= 0 ? 1 - decpt : 0);
 
314
 
 
315
    /* Do we have to truncate any digits? */
 
316
    if (width < len)
 
317
    {
 
318
      if (width < decpt)
 
319
      {
 
320
        if (error != NULL)
 
321
          *error= TRUE;
 
322
        width= decpt;
 
323
      }
 
324
      
 
325
      /*
 
326
        We want to truncate (len - width) least significant digits after the
 
327
        decimal point. For this we are calling dtoa with mode=5, passing the
 
328
        number of significant digits = (len-decpt) - (len-width) = width-decpt
 
329
      */
 
330
      dtoa_free(res, buf, sizeof(buf));
 
331
      res= dtoa(x, 5, width - decpt, &decpt, &sign, &end, buf, sizeof(buf));
 
332
      src= res;
 
333
      len= end - res;
 
334
    }
 
335
 
 
336
    if (len == 0)
 
337
    {
 
338
      /* Underflow. Just print '0' and exit */
 
339
      *dst++= '0';
 
340
      goto end;
 
341
    }
 
342
    
 
343
    /*
 
344
      At this point we are sure we have enough space to put all digits
 
345
      returned by dtoa
 
346
    */
 
347
    if (sign && dst < dend)
 
348
      *dst++= '-';
 
349
    if (decpt <= 0)
 
350
    {
 
351
      if (dst < dend)
 
352
        *dst++= '0';
 
353
      if (len > 0 && dst < dend)
 
354
        *dst++= '.';
 
355
      for (; decpt < 0 && dst < dend; decpt++)
 
356
        *dst++= '0';
 
357
    }
 
358
 
 
359
    for (i= 1; i <= len && dst < dend; i++)
 
360
    {
 
361
      *dst++= *src++;
 
362
      if (i == decpt && i < len && dst < dend)
 
363
        *dst++= '.';
 
364
    }
 
365
    while (i++ <= decpt && dst < dend)
 
366
      *dst++= '0';
 
367
  }
 
368
  else
 
369
  {
 
370
    /* 'e' format */
 
371
    int decpt_sign= 0;
 
372
 
 
373
    if (--decpt < 0)
 
374
    {
 
375
      decpt= -decpt;
 
376
      width--;
 
377
      decpt_sign= 1;
 
378
    }
 
379
    width-= 1 + exp_len; /* eNNN */
 
380
 
 
381
    if (len > 1)
 
382
      width--;
 
383
 
 
384
    if (width <= 0)
 
385
    {
 
386
      /* Overflow */
 
387
      if (error != NULL)
 
388
        *error= TRUE;
 
389
      width= 0;
 
390
    }
 
391
      
 
392
    /* Do we have to truncate any digits? */
 
393
    if (width < len)
 
394
    {
 
395
      /* Yes, re-convert with a smaller width */
 
396
      dtoa_free(res, buf, sizeof(buf));
 
397
      res= dtoa(x, 4, width, &decpt, &sign, &end, buf, sizeof(buf));
 
398
      src= res;
 
399
      len= end - res;
 
400
      if (--decpt < 0)
 
401
        decpt= -decpt;
 
402
    }
 
403
    /*
 
404
      At this point we are sure we have enough space to put all digits
 
405
      returned by dtoa
 
406
    */
 
407
    if (sign && dst < dend)
 
408
      *dst++= '-';
 
409
    if (dst < dend)
 
410
      *dst++= *src++;
 
411
    if (len > 1 && dst < dend)
 
412
    {
 
413
      *dst++= '.';
 
414
      while (src < end && dst < dend)
 
415
        *dst++= *src++;
 
416
    }
 
417
    if (dst < dend)
 
418
      *dst++= 'e';
 
419
    if (decpt_sign && dst < dend)
 
420
      *dst++= '-';
 
421
 
 
422
    if (decpt >= 100 && dst < dend)
 
423
    {
 
424
      *dst++= decpt / 100 + '0';
 
425
      decpt%= 100;
 
426
      if (dst < dend)
 
427
        *dst++= decpt / 10 + '0';
 
428
    }
 
429
    else if (decpt >= 10 && dst < dend)
 
430
      *dst++= decpt / 10 + '0';
 
431
    if (dst < dend)
 
432
      *dst++= decpt % 10 + '0';
 
433
 
 
434
  }
 
435
 
 
436
end:
 
437
  dtoa_free(res, buf, sizeof(buf));
 
438
  *dst= '\0';
 
439
 
 
440
  return dst - to;
 
441
}
 
442
 
 
443
/**
 
444
   @brief
 
445
   Converts string to double (string does not have to be zero-terminated)
 
446
 
 
447
   @details
 
448
   This is a wrapper around dtoa's version of strtod().
 
449
 
 
450
   @param str     input string
 
451
   @param end     address of a pointer to the first character after the input
 
452
                  string. Upon return the pointer is set to point to the first
 
453
                  rejected character.
 
454
   @param error   Upon return is set to EOVERFLOW in case of underflow or
 
455
                  overflow.
 
456
   
 
457
   @return        The resulting double value. In case of underflow, 0.0 is
 
458
                  returned. In case overflow, signed DBL_MAX is returned.
 
459
*/
 
460
 
 
461
double my_strtod(const char *str, char **end, int *error)
 
462
{
 
463
  char buf[DTOA_BUFF_SIZE];
 
464
  double res;
 
465
  DBUG_ASSERT(end != NULL && ((str != NULL && *end != NULL) ||
 
466
                              (str == NULL && *end == NULL)) &&
 
467
              error != NULL);
 
468
 
 
469
  res= my_strtod_int(str, end, error, buf, sizeof(buf));
 
470
  return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX);
 
471
}
 
472
 
 
473
 
 
474
double my_atof(const char *nptr)
 
475
{
 
476
  int error;
 
477
  const char *end= nptr+65535;                  /* Should be enough */
 
478
  return (my_strtod(nptr, (char**) &end, &error));
 
479
}
 
480
 
 
481
 
 
482
/****************************************************************
 
483
 *
 
484
 * The author of this software is David M. Gay.
 
485
 *
 
486
 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
 
487
 *
 
488
 * Permission to use, copy, modify, and distribute this software for any
 
489
 * purpose without fee is hereby granted, provided that this entire notice
 
490
 * is included in all copies of any software which is or includes a copy
 
491
 * or modification of this software and in all copies of the supporting
 
492
 * documentation for such software.
 
493
 *
 
494
 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
 
495
 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
 
496
 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
 
497
 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
 
498
 *
 
499
 ***************************************************************/
 
500
/* Please send bug reports to David M. Gay (dmg at acm dot org,
 
501
 * with " at " changed at "@" and " dot " changed to ".").      */
 
502
 
 
503
/*
 
504
  Original copy of the software is located at http://www.netlib.org/fp/dtoa.c
 
505
  It was adjusted to serve MySQL server needs:
 
506
  * strtod() was modified to not expect a zero-terminated string.
 
507
    It now honors 'se' (end of string) argument as the input parameter,
 
508
    not just as the output one.
 
509
  * in dtoa(), in case of overflow/underflow/NaN result string now contains "0";
 
510
    decpt is set to DTOA_OVERFLOW to indicate overflow.
 
511
  * support for VAX, IBM mainframe and 16-bit hardware removed
 
512
  * we always assume that 64-bit integer type is available
 
513
  * support for Kernigan-Ritchie style headers (pre-ANSI compilers)
 
514
    removed
 
515
  * all gcc warnings ironed out
 
516
  * we always assume multithreaded environment, so we had to change
 
517
    memory allocation procedures to use stack in most cases;
 
518
    malloc is used as the last resort.
 
519
  * pow5mult rewritten to use pre-calculated pow5 list instead of
 
520
    the one generated on the fly.
 
521
*/
 
522
 
 
523
 
 
524
/*
 
525
  On a machine with IEEE extended-precision registers, it is
 
526
  necessary to specify double-precision (53-bit) rounding precision
 
527
  before invoking strtod or dtoa.  If the machine uses (the equivalent
 
528
  of) Intel 80x87 arithmetic, the call
 
529
       _control87(PC_53, MCW_PC);
 
530
  does this with many compilers.  Whether this or another call is
 
531
  appropriate depends on the compiler; for this to work, it may be
 
532
  necessary to #include "float.h" or another system-dependent header
 
533
  file.
 
534
*/
 
535
 
 
536
/*
 
537
  #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
 
538
       and dtoa should round accordingly.
 
539
  #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
 
540
       and Honor_FLT_ROUNDS is not #defined.
 
541
 
 
542
  TODO: check if we can get rid of the above two
 
543
*/
 
544
 
 
545
typedef int32 Long;
 
546
typedef uint32 ULong;
 
547
typedef int64 LLong;
 
548
typedef uint64 ULLong;
 
549
 
 
550
typedef union { double d; ULong L[2]; } U;
 
551
 
 
552
#if defined(WORDS_BIGENDIAN) || (defined(__FLOAT_WORD_ORDER) &&        \
 
553
                                 (__FLOAT_WORD_ORDER == __BIG_ENDIAN))
 
554
#define word0(x) (x)->L[0]
 
555
#define word1(x) (x)->L[1]
 
556
#else
 
557
#define word0(x) (x)->L[1]
 
558
#define word1(x) (x)->L[0]
 
559
#endif
 
560
 
 
561
#define dval(x) (x)->d
 
562
 
 
563
/* #define P DBL_MANT_DIG */
 
564
/* Ten_pmax= floor(P*log(2)/log(5)) */
 
565
/* Bletch= (highest power of 2 < DBL_MAX_10_EXP) / 16 */
 
566
/* Quick_max= floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
 
567
/* Int_max= floor(P*log(FLT_RADIX)/log(10) - 1) */
 
568
 
 
569
#define Exp_shift  20
 
570
#define Exp_shift1 20
 
571
#define Exp_msk1    0x100000
 
572
#define Exp_mask  0x7ff00000
 
573
#define P 53
 
574
#define Bias 1023
 
575
#define Emin (-1022)
 
576
#define Exp_1  0x3ff00000
 
577
#define Exp_11 0x3ff00000
 
578
#define Ebits 11
 
579
#define Frac_mask  0xfffff
 
580
#define Frac_mask1 0xfffff
 
581
#define Ten_pmax 22
 
582
#define Bletch 0x10
 
583
#define Bndry_mask  0xfffff
 
584
#define Bndry_mask1 0xfffff
 
585
#define LSB 1
 
586
#define Sign_bit 0x80000000
 
587
#define Log2P 1
 
588
#define Tiny1 1
 
589
#define Quick_max 14
 
590
#define Int_max 14
 
591
 
 
592
#ifndef Flt_Rounds
 
593
#ifdef FLT_ROUNDS
 
594
#define Flt_Rounds FLT_ROUNDS
 
595
#else
 
596
#define Flt_Rounds 1
 
597
#endif
 
598
#endif /*Flt_Rounds*/
 
599
 
 
600
#ifdef Honor_FLT_ROUNDS
 
601
#define Rounding rounding
 
602
#undef Check_FLT_ROUNDS
 
603
#define Check_FLT_ROUNDS
 
604
#else
 
605
#define Rounding Flt_Rounds
 
606
#endif
 
607
 
 
608
#define rounded_product(a,b) a*= b
 
609
#define rounded_quotient(a,b) a/= b
 
610
 
 
611
#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
 
612
#define Big1 0xffffffff
 
613
#define FFFFFFFF 0xffffffffUL
 
614
 
 
615
/* This is tested to be enough for dtoa */
 
616
 
 
617
#define Kmax 15
 
618
 
 
619
#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign,   \
 
620
                          2*sizeof(int) + y->wds*sizeof(ULong))
 
621
 
 
622
/* Arbitrary-length integer */
 
623
 
 
624
typedef struct Bigint
 
625
{
 
626
  union {
 
627
    ULong *x;              /* points right after this Bigint object */
 
628
    struct Bigint *next;   /* to maintain free lists */
 
629
  } p;
 
630
  int k;                   /* 2^k = maxwds */
 
631
  int maxwds;              /* maximum length in 32-bit words */
 
632
  int sign;                /* not zero if number is negative */
 
633
  int wds;                 /* current length in 32-bit words */
 
634
} Bigint;
 
635
 
 
636
 
 
637
/* A simple stack-memory based allocator for Bigints */
 
638
 
 
639
typedef struct Stack_alloc
 
640
{
 
641
  char *begin;
 
642
  char *free;
 
643
  char *end;
 
644
  /*
 
645
    Having list of free blocks lets us reduce maximum required amount
 
646
    of memory from ~4000 bytes to < 1680 (tested on x86).
 
647
  */
 
648
  Bigint *freelist[Kmax+1];
 
649
} Stack_alloc;
 
650
 
 
651
 
 
652
/*
 
653
  Try to allocate object on stack, and resort to malloc if all
 
654
  stack memory is used. Ensure allocated objects to be aligned by the pointer
 
655
  size in order to not break the alignment rules when storing a pointer to a
 
656
  Bigint.
 
657
*/
 
658
 
 
659
static Bigint *Balloc(int k, Stack_alloc *alloc)
 
660
{
 
661
  Bigint *rv;
 
662
  DBUG_ASSERT(k <= Kmax);
 
663
  if (k <= Kmax &&  alloc->freelist[k])
 
664
  {
 
665
    rv= alloc->freelist[k];
 
666
    alloc->freelist[k]= rv->p.next;
 
667
  }
 
668
  else
 
669
  {
 
670
    int x, len;
 
671
 
 
672
    x= 1 << k;
 
673
    len= MY_ALIGN(sizeof(Bigint) + x * sizeof(ULong), SIZEOF_CHARP);
 
674
 
 
675
    if (alloc->free + len <= alloc->end)
 
676
    {
 
677
      rv= (Bigint*) alloc->free;
 
678
      alloc->free+= len;
 
679
    }
 
680
    else
 
681
      rv= (Bigint*) malloc(len);
 
682
 
 
683
    rv->k= k;
 
684
    rv->maxwds= x;
 
685
  }
 
686
  rv->sign= rv->wds= 0;
 
687
  rv->p.x= (ULong*) (rv + 1);
 
688
  return rv;
 
689
}
 
690
 
 
691
 
 
692
/*
 
693
  If object was allocated on stack, try putting it to the free
 
694
  list. Otherwise call free().
 
695
*/
 
696
 
 
697
static void Bfree(Bigint *v, Stack_alloc *alloc)
 
698
{
 
699
  char *gptr= (char*) v;                       /* generic pointer */
 
700
  if (gptr < alloc->begin || gptr >= alloc->end)
 
701
    free(gptr);
 
702
  else if (v->k <= Kmax)
 
703
  {
 
704
    /*
 
705
      Maintain free lists only for stack objects: this way we don't
 
706
      have to bother with freeing lists in the end of dtoa;
 
707
      heap should not be used normally anyway.
 
708
    */
 
709
    v->p.next= alloc->freelist[v->k];
 
710
    alloc->freelist[v->k]= v;
 
711
  }
 
712
}
 
713
 
 
714
 
 
715
/*
 
716
  This is to place return value of dtoa in: tries to use stack
 
717
  as well, but passes by free lists management and just aligns len by
 
718
  the pointer size in order to not break the alignment rules when storing a
 
719
  pointer to a Bigint.
 
720
*/
 
721
 
 
722
static char *dtoa_alloc(int i, Stack_alloc *alloc)
 
723
{
 
724
  char *rv;
 
725
  int aligned_size= MY_ALIGN(i, SIZEOF_CHARP);
 
726
  if (alloc->free + aligned_size <= alloc->end)
 
727
  {
 
728
    rv= alloc->free;
 
729
    alloc->free+= aligned_size;
 
730
  }
 
731
  else
 
732
    rv= malloc(i);
 
733
  return rv;
 
734
}
 
735
 
 
736
 
 
737
/*
 
738
  dtoa_free() must be used to free values s returned by dtoa()
 
739
  This is the counterpart of dtoa_alloc()
 
740
*/
 
741
 
 
742
static void dtoa_free(char *gptr, char *buf, size_t buf_size)
 
743
{
 
744
  if (gptr < buf || gptr >= buf + buf_size)
 
745
    free(gptr);
 
746
}
 
747
 
 
748
 
 
749
/* Bigint arithmetic functions */
 
750
 
 
751
/* Multiply by m and add a */
 
752
 
 
753
static Bigint *multadd(Bigint *b, int m, int a, Stack_alloc *alloc)
 
754
{
 
755
  int i, wds;
 
756
  ULong *x;
 
757
  ULLong carry, y;
 
758
  Bigint *b1;
 
759
 
 
760
  wds= b->wds;
 
761
  x= b->p.x;
 
762
  i= 0;
 
763
  carry= a;
 
764
  do
 
765
  {
 
766
    y= *x * (ULLong)m + carry;
 
767
    carry= y >> 32;
 
768
    *x++= (ULong)(y & FFFFFFFF);
 
769
  }
 
770
  while (++i < wds);
 
771
  if (carry)
 
772
  {
 
773
    if (wds >= b->maxwds)
 
774
    {
 
775
      b1= Balloc(b->k+1, alloc);
 
776
      Bcopy(b1, b);
 
777
      Bfree(b, alloc);
 
778
      b= b1;
 
779
    }
 
780
    b->p.x[wds++]= (ULong) carry;
 
781
    b->wds= wds;
 
782
  }
 
783
  return b;
 
784
}
 
785
 
 
786
/**
 
787
  Converts a string to Bigint.
 
788
  
 
789
  Now we have nd0 digits, starting at s, followed by a
 
790
  decimal point, followed by nd-nd0 digits.  
 
791
  Unless nd0 == nd, in which case we have a number of the form:
 
792
     ".xxxxxx"    or    "xxxxxx."
 
793
 
 
794
  @param s     Input string, already partially parsed by my_strtod_int().
 
795
  @param nd0   Number of digits before decimal point.
 
796
  @param nd    Total number of digits.
 
797
  @param y9    Pre-computed value of the first nine digits.
 
798
  @param alloc Stack allocator for Bigints.
 
799
 */
 
800
static Bigint *s2b(const char *s, int nd0, int nd, ULong y9, Stack_alloc *alloc)
 
801
{
 
802
  Bigint *b;
 
803
  int i, k;
 
804
  Long x, y;
 
805
 
 
806
  x= (nd + 8) / 9;
 
807
  for (k= 0, y= 1; x > y; y <<= 1, k++) ;
 
808
  b= Balloc(k, alloc);
 
809
  b->p.x[0]= y9;
 
810
  b->wds= 1;
 
811
  
 
812
  i= 9;
 
813
  if (9 < nd0)
 
814
  {
 
815
    s+= 9;
 
816
    do
 
817
      b= multadd(b, 10, *s++ - '0', alloc);
 
818
    while (++i < nd0);
 
819
    s++;                                        /* skip '.' */
 
820
  }
 
821
  else
 
822
    s+= 10;
 
823
  /* now do the fractional part */
 
824
  for(; i < nd; i++)
 
825
    b= multadd(b, 10, *s++ - '0', alloc);
 
826
  return b;
 
827
}
 
828
 
 
829
 
 
830
static int hi0bits(register ULong x)
 
831
{
 
832
  register int k= 0;
 
833
 
 
834
  if (!(x & 0xffff0000))
 
835
  {
 
836
    k= 16;
 
837
    x<<= 16;
 
838
  }
 
839
  if (!(x & 0xff000000))
 
840
  {
 
841
    k+= 8;
 
842
    x<<= 8;
 
843
  }
 
844
  if (!(x & 0xf0000000))
 
845
  {
 
846
    k+= 4;
 
847
    x<<= 4;
 
848
  }
 
849
  if (!(x & 0xc0000000))
 
850
  {
 
851
    k+= 2;
 
852
    x<<= 2;
 
853
  }
 
854
  if (!(x & 0x80000000))
 
855
  {
 
856
    k++;
 
857
    if (!(x & 0x40000000))
 
858
      return 32;
 
859
  }
 
860
  return k;
 
861
}
 
862
 
 
863
 
 
864
static int lo0bits(ULong *y)
 
865
{
 
866
  register int k;
 
867
  register ULong x= *y;
 
868
 
 
869
  if (x & 7)
 
870
  {
 
871
    if (x & 1)
 
872
      return 0;
 
873
    if (x & 2)
 
874
    {
 
875
      *y= x >> 1;
 
876
      return 1;
 
877
    }
 
878
    *y= x >> 2;
 
879
    return 2;
 
880
  }
 
881
  k= 0;
 
882
  if (!(x & 0xffff))
 
883
  {
 
884
    k= 16;
 
885
    x>>= 16;
 
886
  }
 
887
  if (!(x & 0xff))
 
888
  {
 
889
    k+= 8;
 
890
    x>>= 8;
 
891
  }
 
892
  if (!(x & 0xf))
 
893
  {
 
894
    k+= 4;
 
895
    x>>= 4;
 
896
  }
 
897
  if (!(x & 0x3))
 
898
  {
 
899
    k+= 2;
 
900
    x>>= 2;
 
901
  }
 
902
  if (!(x & 1))
 
903
  {
 
904
    k++;
 
905
    x>>= 1;
 
906
    if (!x)
 
907
      return 32;
 
908
  }
 
909
  *y= x;
 
910
  return k;
 
911
}
 
912
 
 
913
 
 
914
/* Convert integer to Bigint number */
 
915
 
 
916
static Bigint *i2b(int i, Stack_alloc *alloc)
 
917
{
 
918
  Bigint *b;
 
919
 
 
920
  b= Balloc(1, alloc);
 
921
  b->p.x[0]= i;
 
922
  b->wds= 1;
 
923
  return b;
 
924
}
 
925
 
 
926
 
 
927
/* Multiply two Bigint numbers */
 
928
 
 
929
static Bigint *mult(Bigint *a, Bigint *b, Stack_alloc *alloc)
 
930
{
 
931
  Bigint *c;
 
932
  int k, wa, wb, wc;
 
933
  ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
 
934
  ULong y;
 
935
  ULLong carry, z;
 
936
 
 
937
  if (a->wds < b->wds)
 
938
  {
 
939
    c= a;
 
940
    a= b;
 
941
    b= c;
 
942
  }
 
943
  k= a->k;
 
944
  wa= a->wds;
 
945
  wb= b->wds;
 
946
  wc= wa + wb;
 
947
  if (wc > a->maxwds)
 
948
    k++;
 
949
  c= Balloc(k, alloc);
 
950
  for (x= c->p.x, xa= x + wc; x < xa; x++)
 
951
    *x= 0;
 
952
  xa= a->p.x;
 
953
  xae= xa + wa;
 
954
  xb= b->p.x;
 
955
  xbe= xb + wb;
 
956
  xc0= c->p.x;
 
957
  for (; xb < xbe; xc0++)
 
958
  {
 
959
    if ((y= *xb++))
 
960
    {
 
961
      x= xa;
 
962
      xc= xc0;
 
963
      carry= 0;
 
964
      do
 
965
      {
 
966
        z= *x++ * (ULLong)y + *xc + carry;
 
967
        carry= z >> 32;
 
968
        *xc++= (ULong) (z & FFFFFFFF);
 
969
      }
 
970
      while (x < xae);
 
971
      *xc= (ULong) carry;
 
972
    }
 
973
  }
 
974
  for (xc0= c->p.x, xc= xc0 + wc; wc > 0 && !*--xc; --wc) ;
 
975
  c->wds= wc;
 
976
  return c;
 
977
}
 
978
 
 
979
 
 
980
/*
 
981
  Precalculated array of powers of 5: tested to be enough for
 
982
  vasting majority of dtoa_r cases.
 
983
*/
 
984
 
 
985
static ULong powers5[]=
 
986
{
 
987
  625UL,
 
988
 
 
989
  390625UL,
 
990
 
 
991
  2264035265UL, 35UL,
 
992
 
 
993
  2242703233UL, 762134875UL,  1262UL,
 
994
 
 
995
  3211403009UL, 1849224548UL, 3668416493UL, 3913284084UL, 1593091UL,
 
996
 
 
997
  781532673UL,  64985353UL,   253049085UL,  594863151UL,  3553621484UL,
 
998
  3288652808UL, 3167596762UL, 2788392729UL, 3911132675UL, 590UL,
 
999
 
 
1000
  2553183233UL, 3201533787UL, 3638140786UL, 303378311UL, 1809731782UL,
 
1001
  3477761648UL, 3583367183UL, 649228654UL, 2915460784UL, 487929380UL,
 
1002
  1011012442UL, 1677677582UL, 3428152256UL, 1710878487UL, 1438394610UL,
 
1003
  2161952759UL, 4100910556UL, 1608314830UL, 349175UL
 
1004
};
 
1005
 
 
1006
 
 
1007
static Bigint p5_a[]=
 
1008
{
 
1009
  /*  { x } - k - maxwds - sign - wds */
 
1010
  { { powers5 }, 1, 1, 0, 1 },
 
1011
  { { powers5 + 1 }, 1, 1, 0, 1 },
 
1012
  { { powers5 + 2 }, 1, 2, 0, 2 },
 
1013
  { { powers5 + 4 }, 2, 3, 0, 3 },
 
1014
  { { powers5 + 7 }, 3, 5, 0, 5 },
 
1015
  { { powers5 + 12 }, 4, 10, 0, 10 },
 
1016
  { { powers5 + 22 }, 5, 19, 0, 19 }
 
1017
};
 
1018
 
 
1019
#define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
 
1020
 
 
1021
static Bigint *pow5mult(Bigint *b, int k, Stack_alloc *alloc)
 
1022
{
 
1023
  Bigint *b1, *p5, *p51=NULL;
 
1024
  int i;
 
1025
  static int p05[3]= { 5, 25, 125 };
 
1026
  my_bool overflow= FALSE;
 
1027
 
 
1028
  if ((i= k & 3))
 
1029
    b= multadd(b, p05[i-1], 0, alloc);
 
1030
 
 
1031
  if (!(k>>= 2))
 
1032
    return b;
 
1033
  p5= p5_a;
 
1034
  for (;;)
 
1035
  {
 
1036
    if (k & 1)
 
1037
    {
 
1038
      b1= mult(b, p5, alloc);
 
1039
      Bfree(b, alloc);
 
1040
      b= b1;
 
1041
    }
 
1042
    if (!(k>>= 1))
 
1043
      break;
 
1044
    /* Calculate next power of 5 */
 
1045
    if (overflow)
 
1046
    {
 
1047
      p51= mult(p5, p5, alloc);
 
1048
      Bfree(p5, alloc);
 
1049
      p5= p51;
 
1050
    }
 
1051
    else if (p5 < p5_a + P5A_MAX)
 
1052
      ++p5;
 
1053
    else if (p5 == p5_a + P5A_MAX)
 
1054
    {
 
1055
      p5= mult(p5, p5, alloc);
 
1056
      overflow= TRUE;
 
1057
    }
 
1058
  }
 
1059
  if (p51)
 
1060
    Bfree(p51, alloc);
 
1061
  return b;
 
1062
}
 
1063
 
 
1064
 
 
1065
static Bigint *lshift(Bigint *b, int k, Stack_alloc *alloc)
 
1066
{
 
1067
  int i, k1, n, n1;
 
1068
  Bigint *b1;
 
1069
  ULong *x, *x1, *xe, z;
 
1070
 
 
1071
  n= k >> 5;
 
1072
  k1= b->k;
 
1073
  n1= n + b->wds + 1;
 
1074
  for (i= b->maxwds; n1 > i; i<<= 1)
 
1075
    k1++;
 
1076
  b1= Balloc(k1, alloc);
 
1077
  x1= b1->p.x;
 
1078
  for (i= 0; i < n; i++)
 
1079
    *x1++= 0;
 
1080
  x= b->p.x;
 
1081
  xe= x + b->wds;
 
1082
  if (k&= 0x1f)
 
1083
  {
 
1084
    k1= 32 - k;
 
1085
    z= 0;
 
1086
    do
 
1087
    {
 
1088
      *x1++= *x << k | z;
 
1089
      z= *x++ >> k1;
 
1090
    }
 
1091
    while (x < xe);
 
1092
    if ((*x1= z))
 
1093
      ++n1;
 
1094
  }
 
1095
  else
 
1096
    do
 
1097
      *x1++= *x++;
 
1098
    while (x < xe);
 
1099
  b1->wds= n1 - 1;
 
1100
  Bfree(b, alloc);
 
1101
  return b1;
 
1102
}
 
1103
 
 
1104
 
 
1105
static int cmp(Bigint *a, Bigint *b)
 
1106
{
 
1107
  ULong *xa, *xa0, *xb, *xb0;
 
1108
  int i, j;
 
1109
 
 
1110
  i= a->wds;
 
1111
  j= b->wds;
 
1112
  if (i-= j)
 
1113
    return i;
 
1114
  xa0= a->p.x;
 
1115
  xa= xa0 + j;
 
1116
  xb0= b->p.x;
 
1117
  xb= xb0 + j;
 
1118
  for (;;)
 
1119
  {
 
1120
    if (*--xa != *--xb)
 
1121
      return *xa < *xb ? -1 : 1;
 
1122
    if (xa <= xa0)
 
1123
      break;
 
1124
  }
 
1125
  return 0;
 
1126
}
 
1127
 
 
1128
 
 
1129
static Bigint *diff(Bigint *a, Bigint *b, Stack_alloc *alloc)
 
1130
{
 
1131
  Bigint *c;
 
1132
  int i, wa, wb;
 
1133
  ULong *xa, *xae, *xb, *xbe, *xc;
 
1134
  ULLong borrow, y;
 
1135
 
 
1136
  i= cmp(a,b);
 
1137
  if (!i)
 
1138
  {
 
1139
    c= Balloc(0, alloc);
 
1140
    c->wds= 1;
 
1141
    c->p.x[0]= 0;
 
1142
    return c;
 
1143
  }
 
1144
  if (i < 0)
 
1145
  {
 
1146
    c= a;
 
1147
    a= b;
 
1148
    b= c;
 
1149
    i= 1;
 
1150
  }
 
1151
  else
 
1152
    i= 0;
 
1153
  c= Balloc(a->k, alloc);
 
1154
  c->sign= i;
 
1155
  wa= a->wds;
 
1156
  xa= a->p.x;
 
1157
  xae= xa + wa;
 
1158
  wb= b->wds;
 
1159
  xb= b->p.x;
 
1160
  xbe= xb + wb;
 
1161
  xc= c->p.x;
 
1162
  borrow= 0;
 
1163
  do
 
1164
  {
 
1165
    y= (ULLong)*xa++ - *xb++ - borrow;
 
1166
    borrow= y >> 32 & (ULong)1;
 
1167
    *xc++= (ULong) (y & FFFFFFFF);
 
1168
  }
 
1169
  while (xb < xbe);
 
1170
  while (xa < xae)
 
1171
  {
 
1172
    y= *xa++ - borrow;
 
1173
    borrow= y >> 32 & (ULong)1;
 
1174
    *xc++= (ULong) (y & FFFFFFFF);
 
1175
  }
 
1176
  while (!*--xc)
 
1177
    wa--;
 
1178
  c->wds= wa;
 
1179
  return c;
 
1180
}
 
1181
 
 
1182
 
 
1183
static double ulp(U *x)
 
1184
{
 
1185
  register Long L;
 
1186
  U u;
 
1187
 
 
1188
  L= (word0(x) & Exp_mask) - (P - 1)*Exp_msk1;
 
1189
  word0(&u) = L;
 
1190
  word1(&u) = 0;
 
1191
  return dval(&u);
 
1192
}
 
1193
 
 
1194
 
 
1195
static double b2d(Bigint *a, int *e)
 
1196
{
 
1197
  ULong *xa, *xa0, w, y, z;
 
1198
  int k;
 
1199
  U d;
 
1200
#define d0 word0(&d)
 
1201
#define d1 word1(&d)
 
1202
 
 
1203
  xa0= a->p.x;
 
1204
  xa= xa0 + a->wds;
 
1205
  y= *--xa;
 
1206
  k= hi0bits(y);
 
1207
  *e= 32 - k;
 
1208
  if (k < Ebits)
 
1209
  {
 
1210
    d0= Exp_1 | y >> (Ebits - k);
 
1211
    w= xa > xa0 ? *--xa : 0;
 
1212
    d1= y << ((32-Ebits) + k) | w >> (Ebits - k);
 
1213
    goto ret_d;
 
1214
  }
 
1215
  z= xa > xa0 ? *--xa : 0;
 
1216
  if (k-= Ebits)
 
1217
  {
 
1218
    d0= Exp_1 | y << k | z >> (32 - k);
 
1219
    y= xa > xa0 ? *--xa : 0;
 
1220
    d1= z << k | y >> (32 - k);
 
1221
  }
 
1222
  else
 
1223
  {
 
1224
    d0= Exp_1 | y;
 
1225
    d1= z;
 
1226
  }
 
1227
 ret_d:
 
1228
#undef d0
 
1229
#undef d1
 
1230
  return dval(&d);
 
1231
}
 
1232
 
 
1233
 
 
1234
static Bigint *d2b(U *d, int *e, int *bits, Stack_alloc *alloc)
 
1235
{
 
1236
  Bigint *b;
 
1237
  int de, k;
 
1238
  ULong *x, y, z;
 
1239
  int i;
 
1240
#define d0 word0(d)
 
1241
#define d1 word1(d)
 
1242
 
 
1243
  b= Balloc(1, alloc);
 
1244
  x= b->p.x;
 
1245
 
 
1246
  z= d0 & Frac_mask;
 
1247
  d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
 
1248
  if ((de= (int)(d0 >> Exp_shift)))
 
1249
    z|= Exp_msk1;
 
1250
  if ((y= d1))
 
1251
  {
 
1252
    if ((k= lo0bits(&y)))
 
1253
    {
 
1254
      x[0]= y | z << (32 - k);
 
1255
      z>>= k;
 
1256
    }
 
1257
    else
 
1258
      x[0]= y;
 
1259
    i= b->wds= (x[1]= z) ? 2 : 1;
 
1260
  }
 
1261
  else
 
1262
  {
 
1263
    k= lo0bits(&z);
 
1264
    x[0]= z;
 
1265
    i= b->wds= 1;
 
1266
    k+= 32;
 
1267
  }
 
1268
  if (de)
 
1269
  {
 
1270
    *e= de - Bias - (P-1) + k;
 
1271
    *bits= P - k;
 
1272
  }
 
1273
  else
 
1274
  {
 
1275
    *e= de - Bias - (P-1) + 1 + k;
 
1276
    *bits= 32*i - hi0bits(x[i-1]);
 
1277
  }
 
1278
  return b;
 
1279
#undef d0
 
1280
#undef d1
 
1281
}
 
1282
 
 
1283
 
 
1284
static double ratio(Bigint *a, Bigint *b)
 
1285
{
 
1286
  U da, db;
 
1287
  int k, ka, kb;
 
1288
 
 
1289
  dval(&da)= b2d(a, &ka);
 
1290
  dval(&db)= b2d(b, &kb);
 
1291
  k= ka - kb + 32*(a->wds - b->wds);
 
1292
  if (k > 0)
 
1293
    word0(&da)+= k*Exp_msk1;
 
1294
  else
 
1295
  {
 
1296
    k= -k;
 
1297
    word0(&db)+= k*Exp_msk1;
 
1298
  }
 
1299
  return dval(&da) / dval(&db);
 
1300
}
 
1301
 
 
1302
static const double tens[] =
 
1303
{
 
1304
  1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
 
1305
  1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
 
1306
  1e20, 1e21, 1e22
 
1307
};
 
1308
 
 
1309
static const double bigtens[]= { 1e16, 1e32, 1e64, 1e128, 1e256 };
 
1310
static const double tinytens[]=
 
1311
{ 1e-16, 1e-32, 1e-64, 1e-128,
 
1312
  9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */
 
1313
};
 
1314
/*
 
1315
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow 
 
1316
  flag unnecessarily.  It leads to a song and dance at the end of strtod.
 
1317
*/
 
1318
#define Scale_Bit 0x10
 
1319
#define n_bigtens 5
 
1320
 
 
1321
/*
 
1322
  strtod for IEEE--arithmetic machines.
 
1323
 
 
1324
  This strtod returns a nearest machine number to the input decimal
 
1325
  string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even
 
1326
  rule.
 
1327
 
 
1328
  Inspired loosely by William D. Clinger's paper "How to Read Floating
 
1329
  Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
 
1330
 
 
1331
  Modifications:
 
1332
 
 
1333
   1. We only require IEEE (not IEEE double-extended).
 
1334
   2. We get by with floating-point arithmetic in a case that
 
1335
     Clinger missed -- when we're computing d * 10^n
 
1336
     for a small integer d and the integer n is not too
 
1337
     much larger than 22 (the maximum integer k for which
 
1338
     we can represent 10^k exactly), we may be able to
 
1339
     compute (d*10^k) * 10^(e-k) with just one roundoff.
 
1340
   3. Rather than a bit-at-a-time adjustment of the binary
 
1341
     result in the hard case, we use floating-point
 
1342
     arithmetic to determine the adjustment to within
 
1343
     one bit; only in really hard cases do we need to
 
1344
     compute a second residual.
 
1345
   4. Because of 3., we don't need a large table of powers of 10
 
1346
     for ten-to-e (just some small tables, e.g. of 10^k
 
1347
     for 0 <= k <= 22).
 
1348
*/
 
1349
 
 
1350
static double my_strtod_int(const char *s00, char **se, int *error, char *buf, size_t buf_size)
 
1351
{
 
1352
  int scale;
 
1353
  int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, UNINIT_VAR(c), dsign,
 
1354
     e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
 
1355
  const char *s, *s0, *s1, *end = *se;
 
1356
  double aadj, aadj1;
 
1357
  U aadj2, adj, rv, rv0;
 
1358
  Long L;
 
1359
  ULong y, z;
 
1360
  Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
 
1361
#ifdef SET_INEXACT
 
1362
  int inexact, oldinexact;
 
1363
#endif
 
1364
#ifdef Honor_FLT_ROUNDS
 
1365
  int rounding;
 
1366
#endif
 
1367
  Stack_alloc alloc;
 
1368
 
 
1369
  *error= 0;
 
1370
 
 
1371
  alloc.begin= alloc.free= buf;
 
1372
  alloc.end= buf + buf_size;
 
1373
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
 
1374
 
 
1375
  sign= nz0= nz= 0;
 
1376
  dval(&rv)= 0.;
 
1377
  for (s= s00; s < end; s++)
 
1378
    switch (*s) {
 
1379
    case '-':
 
1380
      sign= 1;
 
1381
      /* no break */
 
1382
    case '+':
 
1383
      s++;
 
1384
      goto break2;
 
1385
    case '\t':
 
1386
    case '\n':
 
1387
    case '\v':
 
1388
    case '\f':
 
1389
    case '\r':
 
1390
    case ' ':
 
1391
      continue;
 
1392
    default:
 
1393
      goto break2;
 
1394
    }
 
1395
 break2:
 
1396
  if (s >= end)
 
1397
    goto ret0;
 
1398
  
 
1399
  if (*s == '0')
 
1400
  {
 
1401
    nz0= 1;
 
1402
    while (++s < end && *s == '0') ;
 
1403
    if (s >= end)
 
1404
      goto ret;
 
1405
  }
 
1406
  s0= s;
 
1407
  y= z= 0;
 
1408
  for (nd= nf= 0; s < end && (c= *s) >= '0' && c <= '9'; nd++, s++)
 
1409
    if (nd < 9)
 
1410
      y= 10*y + c - '0';
 
1411
    else if (nd < 16)
 
1412
      z= 10*z + c - '0';
 
1413
  nd0= nd;
 
1414
  if (s < end && c == '.')
 
1415
  {
 
1416
    c= *++s;
 
1417
    if (!nd)
 
1418
    {
 
1419
      for (; s < end; ++s)
 
1420
      {
 
1421
        c= *s;
 
1422
        if (c != '0')
 
1423
          break;
 
1424
        nz++;
 
1425
      }
 
1426
      if (s < end && c > '0' && c <= '9')
 
1427
      {
 
1428
        s0= s;
 
1429
        nf+= nz;
 
1430
        nz= 0;
 
1431
      }
 
1432
      else
 
1433
        goto dig_done;
 
1434
    }
 
1435
    for (; s < end; ++s)
 
1436
    {
 
1437
      c= *s;
 
1438
      if (c < '0' || c > '9')
 
1439
        break;
 
1440
      /*
 
1441
        Here we are parsing the fractional part.
 
1442
        We can stop counting digits after a while: the extra digits
 
1443
        will not contribute to the actual result produced by s2b().
 
1444
        We have to continue scanning, in case there is an exponent part.
 
1445
       */
 
1446
      if (nd < 2 * DBL_DIG)
 
1447
      {
 
1448
        nz++;
 
1449
        if (c-= '0')
 
1450
        {
 
1451
          nf+= nz;
 
1452
          for (i= 1; i < nz; i++)
 
1453
            if (nd++ < 9)
 
1454
              y*= 10;
 
1455
            else if (nd <= DBL_DIG + 1)
 
1456
              z*= 10;
 
1457
          if (nd++ < 9)
 
1458
            y= 10*y + c;
 
1459
          else if (nd <= DBL_DIG + 1)
 
1460
            z= 10*z + c;
 
1461
          nz= 0;
 
1462
        }
 
1463
      }
 
1464
    }
 
1465
  }
 
1466
 dig_done:
 
1467
  e= 0;
 
1468
  if (s < end && (c == 'e' || c == 'E'))
 
1469
  {
 
1470
    if (!nd && !nz && !nz0)
 
1471
      goto ret0;
 
1472
    s00= s;
 
1473
    esign= 0;
 
1474
    if (++s < end)
 
1475
      switch (c= *s) {
 
1476
      case '-':
 
1477
        esign= 1;
 
1478
      case '+':
 
1479
        c= *++s;
 
1480
      }
 
1481
    if (s < end && c >= '0' && c <= '9')
 
1482
    {
 
1483
      while (s < end && c == '0')
 
1484
        c= *++s;
 
1485
      if (s < end && c > '0' && c <= '9') {
 
1486
        L= c - '0';
 
1487
        s1= s;
 
1488
        while (++s < end && (c= *s) >= '0' && c <= '9')
 
1489
          L= 10*L + c - '0';
 
1490
        if (s - s1 > 8 || L > 19999)
 
1491
          /* Avoid confusion from exponents
 
1492
           * so large that e might overflow.
 
1493
           */
 
1494
          e= 19999; /* safe for 16 bit ints */
 
1495
        else
 
1496
          e= (int)L;
 
1497
        if (esign)
 
1498
          e= -e;
 
1499
      }
 
1500
      else
 
1501
        e= 0;
 
1502
    }
 
1503
    else
 
1504
      s= s00;
 
1505
  }
 
1506
  if (!nd)
 
1507
  {
 
1508
    if (!nz && !nz0)
 
1509
    {
 
1510
 ret0:
 
1511
      s= s00;
 
1512
      sign= 0;
 
1513
    }
 
1514
    goto ret;
 
1515
  }
 
1516
  e1= e -= nf;
 
1517
 
 
1518
  /*
 
1519
    Now we have nd0 digits, starting at s0, followed by a
 
1520
    decimal point, followed by nd-nd0 digits.  The number we're
 
1521
    after is the integer represented by those digits times
 
1522
    10**e
 
1523
   */
 
1524
 
 
1525
  if (!nd0)
 
1526
    nd0= nd;
 
1527
  k= nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
 
1528
  dval(&rv)= y;
 
1529
  if (k > 9)
 
1530
  {
 
1531
#ifdef SET_INEXACT
 
1532
    if (k > DBL_DIG)
 
1533
      oldinexact = get_inexact();
 
1534
#endif
 
1535
    dval(&rv)= tens[k - 9] * dval(&rv) + z;
 
1536
  }
 
1537
  bd0= 0;
 
1538
  if (nd <= DBL_DIG
 
1539
#ifndef Honor_FLT_ROUNDS
 
1540
    && Flt_Rounds == 1
 
1541
#endif
 
1542
      )
 
1543
  {
 
1544
    if (!e)
 
1545
      goto ret;
 
1546
    if (e > 0)
 
1547
    {
 
1548
      if (e <= Ten_pmax)
 
1549
      {
 
1550
#ifdef Honor_FLT_ROUNDS
 
1551
        /* round correctly FLT_ROUNDS = 2 or 3 */
 
1552
        if (sign)
 
1553
        {
 
1554
          rv.d= -rv.d;
 
1555
          sign= 0;
 
1556
        }
 
1557
#endif
 
1558
        /* rv = */ rounded_product(dval(&rv), tens[e]);
 
1559
        goto ret;
 
1560
      }
 
1561
      i= DBL_DIG - nd;
 
1562
      if (e <= Ten_pmax + i)
 
1563
      {
 
1564
        /*
 
1565
          A fancier test would sometimes let us do
 
1566
          this for larger i values.
 
1567
        */
 
1568
#ifdef Honor_FLT_ROUNDS
 
1569
        /* round correctly FLT_ROUNDS = 2 or 3 */
 
1570
        if (sign)
 
1571
        {
 
1572
          rv.d= -rv.d;
 
1573
          sign= 0;
 
1574
        }
 
1575
#endif
 
1576
        e-= i;
 
1577
        dval(&rv)*= tens[i];
 
1578
        /* rv = */ rounded_product(dval(&rv), tens[e]);
 
1579
        goto ret;
 
1580
      }
 
1581
    }
 
1582
#ifndef Inaccurate_Divide
 
1583
    else if (e >= -Ten_pmax)
 
1584
    {
 
1585
#ifdef Honor_FLT_ROUNDS
 
1586
      /* round correctly FLT_ROUNDS = 2 or 3 */
 
1587
      if (sign)
 
1588
      {
 
1589
        rv.d= -rv.d;
 
1590
        sign= 0;
 
1591
      }
 
1592
#endif
 
1593
      /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
 
1594
      goto ret;
 
1595
    }
 
1596
#endif
 
1597
  }
 
1598
  e1+= nd - k;
 
1599
 
 
1600
#ifdef SET_INEXACT
 
1601
  inexact= 1;
 
1602
  if (k <= DBL_DIG)
 
1603
    oldinexact= get_inexact();
 
1604
#endif
 
1605
  scale= 0;
 
1606
#ifdef Honor_FLT_ROUNDS
 
1607
  if ((rounding= Flt_Rounds) >= 2)
 
1608
  {
 
1609
    if (sign)
 
1610
      rounding= rounding == 2 ? 0 : 2;
 
1611
    else
 
1612
      if (rounding != 2)
 
1613
        rounding= 0;
 
1614
  }
 
1615
#endif
 
1616
 
 
1617
  /* Get starting approximation = rv * 10**e1 */
 
1618
 
 
1619
  if (e1 > 0)
 
1620
  {
 
1621
    if ((i= e1 & 15))
 
1622
      dval(&rv)*= tens[i];
 
1623
    if (e1&= ~15)
 
1624
    {
 
1625
      if (e1 > DBL_MAX_10_EXP)
 
1626
      {
 
1627
 ovfl:
 
1628
        *error= EOVERFLOW;
 
1629
        /* Can't trust HUGE_VAL */
 
1630
#ifdef Honor_FLT_ROUNDS
 
1631
        switch (rounding)
 
1632
        {
 
1633
        case 0: /* toward 0 */
 
1634
        case 3: /* toward -infinity */
 
1635
          word0(&rv)= Big0;
 
1636
          word1(&rv)= Big1;
 
1637
          break;
 
1638
        default:
 
1639
          word0(&rv)= Exp_mask;
 
1640
          word1(&rv)= 0;
 
1641
        }
 
1642
#else /*Honor_FLT_ROUNDS*/
 
1643
        word0(&rv)= Exp_mask;
 
1644
        word1(&rv)= 0;
 
1645
#endif /*Honor_FLT_ROUNDS*/
 
1646
#ifdef SET_INEXACT
 
1647
        /* set overflow bit */
 
1648
        dval(&rv0)= 1e300;
 
1649
        dval(&rv0)*= dval(&rv0);
 
1650
#endif
 
1651
        if (bd0)
 
1652
          goto retfree;
 
1653
        goto ret;
 
1654
      }
 
1655
      e1>>= 4;
 
1656
      for(j= 0; e1 > 1; j++, e1>>= 1)
 
1657
        if (e1 & 1)
 
1658
          dval(&rv)*= bigtens[j];
 
1659
    /* The last multiplication could overflow. */
 
1660
      word0(&rv)-= P*Exp_msk1;
 
1661
      dval(&rv)*= bigtens[j];
 
1662
      if ((z= word0(&rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - P))
 
1663
        goto ovfl;
 
1664
      if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P))
 
1665
      {
 
1666
        /* set to largest number (Can't trust DBL_MAX) */
 
1667
        word0(&rv)= Big0;
 
1668
        word1(&rv)= Big1;
 
1669
      }
 
1670
      else
 
1671
        word0(&rv)+= P*Exp_msk1;
 
1672
    }
 
1673
  }
 
1674
  else if (e1 < 0)
 
1675
  {
 
1676
    e1= -e1;
 
1677
    if ((i= e1 & 15))
 
1678
      dval(&rv)/= tens[i];
 
1679
    if ((e1>>= 4))
 
1680
    {
 
1681
      if (e1 >= 1 << n_bigtens)
 
1682
        goto undfl;
 
1683
      if (e1 & Scale_Bit)
 
1684
        scale= 2 * P;
 
1685
      for(j= 0; e1 > 0; j++, e1>>= 1)
 
1686
        if (e1 & 1)
 
1687
          dval(&rv)*= tinytens[j];
 
1688
      if (scale && (j = 2 * P + 1 - ((word0(&rv) & Exp_mask) >> Exp_shift)) > 0)
 
1689
      {
 
1690
        /* scaled rv is denormal; zap j low bits */
 
1691
        if (j >= 32)
 
1692
        {
 
1693
          word1(&rv)= 0;
 
1694
          if (j >= 53)
 
1695
            word0(&rv)= (P + 2) * Exp_msk1;
 
1696
          else
 
1697
            word0(&rv)&= 0xffffffff << (j - 32);
 
1698
        }
 
1699
        else
 
1700
          word1(&rv)&= 0xffffffff << j;
 
1701
      }
 
1702
      if (!dval(&rv))
 
1703
      {
 
1704
 undfl:
 
1705
          dval(&rv)= 0.;
 
1706
          if (bd0)
 
1707
            goto retfree;
 
1708
          goto ret;
 
1709
      }
 
1710
    }
 
1711
  }
 
1712
 
 
1713
  /* Now the hard part -- adjusting rv to the correct value.*/
 
1714
 
 
1715
  /* Put digits into bd: true value = bd * 10^e */
 
1716
 
 
1717
  bd0= s2b(s0, nd0, nd, y, &alloc);
 
1718
 
 
1719
  for(;;)
 
1720
  {
 
1721
    bd= Balloc(bd0->k, &alloc);
 
1722
    Bcopy(bd, bd0);
 
1723
    bb= d2b(&rv, &bbe, &bbbits, &alloc);  /* rv = bb * 2^bbe */
 
1724
    bs= i2b(1, &alloc);
 
1725
 
 
1726
    if (e >= 0)
 
1727
    {
 
1728
      bb2= bb5= 0;
 
1729
      bd2= bd5= e;
 
1730
    }
 
1731
    else
 
1732
    {
 
1733
      bb2= bb5= -e;
 
1734
      bd2= bd5= 0;
 
1735
    }
 
1736
    if (bbe >= 0)
 
1737
      bb2+= bbe;
 
1738
    else
 
1739
      bd2-= bbe;
 
1740
    bs2= bb2;
 
1741
#ifdef Honor_FLT_ROUNDS
 
1742
    if (rounding != 1)
 
1743
      bs2++;
 
1744
#endif
 
1745
    j= bbe - scale;
 
1746
    i= j + bbbits - 1;  /* logb(rv) */
 
1747
    if (i < Emin)  /* denormal */
 
1748
      j+= P - Emin;
 
1749
    else
 
1750
      j= P + 1 - bbbits;
 
1751
    bb2+= j;
 
1752
    bd2+= j;
 
1753
    bd2+= scale;
 
1754
    i= bb2 < bd2 ? bb2 : bd2;
 
1755
    if (i > bs2)
 
1756
      i= bs2;
 
1757
    if (i > 0)
 
1758
    {
 
1759
      bb2-= i;
 
1760
      bd2-= i;
 
1761
      bs2-= i;
 
1762
    }
 
1763
    if (bb5 > 0)
 
1764
    {
 
1765
      bs= pow5mult(bs, bb5, &alloc);
 
1766
      bb1= mult(bs, bb, &alloc);
 
1767
      Bfree(bb, &alloc);
 
1768
      bb= bb1;
 
1769
    }
 
1770
    if (bb2 > 0)
 
1771
      bb= lshift(bb, bb2, &alloc);
 
1772
    if (bd5 > 0)
 
1773
      bd= pow5mult(bd, bd5, &alloc);
 
1774
    if (bd2 > 0)
 
1775
      bd= lshift(bd, bd2, &alloc);
 
1776
    if (bs2 > 0)
 
1777
      bs= lshift(bs, bs2, &alloc);
 
1778
    delta= diff(bb, bd, &alloc);
 
1779
    dsign= delta->sign;
 
1780
    delta->sign= 0;
 
1781
    i= cmp(delta, bs);
 
1782
#ifdef Honor_FLT_ROUNDS
 
1783
    if (rounding != 1)
 
1784
    {
 
1785
      if (i < 0)
 
1786
      {
 
1787
        /* Error is less than an ulp */
 
1788
        if (!delta->p.x[0] && delta->wds <= 1)
 
1789
        {
 
1790
          /* exact */
 
1791
#ifdef SET_INEXACT
 
1792
          inexact= 0;
 
1793
#endif
 
1794
          break;
 
1795
        }
 
1796
        if (rounding)
 
1797
        {
 
1798
          if (dsign)
 
1799
          {
 
1800
            adj.d= 1.;
 
1801
            goto apply_adj;
 
1802
          }
 
1803
        }
 
1804
        else if (!dsign)
 
1805
        {
 
1806
          adj.d= -1.;
 
1807
          if (!word1(&rv) && !(word0(&rv) & Frac_mask))
 
1808
          {
 
1809
            y= word0(&rv) & Exp_mask;
 
1810
            if (!scale || y > 2*P*Exp_msk1)
 
1811
            {
 
1812
              delta= lshift(delta, Log2P, &alloc);
 
1813
              if (cmp(delta, bs) <= 0)
 
1814
              adj.d= -0.5;
 
1815
            }
 
1816
          }
 
1817
 apply_adj:
 
1818
          if (scale && (y= word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1)
 
1819
            word0(&adj)+= (2 * P + 1) * Exp_msk1 - y;
 
1820
          dval(&rv)+= adj.d * ulp(&rv);
 
1821
        }
 
1822
        break;
 
1823
      }
 
1824
      adj.d= ratio(delta, bs);
 
1825
      if (adj.d < 1.)
 
1826
        adj.d= 1.;
 
1827
      if (adj.d <= 0x7ffffffe)
 
1828
      {
 
1829
        /* adj = rounding ? ceil(adj) : floor(adj); */
 
1830
        y= adj.d;
 
1831
        if (y != adj.d)
 
1832
        {
 
1833
          if (!((rounding >> 1) ^ dsign))
 
1834
            y++;
 
1835
          adj.d= y;
 
1836
        }
 
1837
      }
 
1838
      if (scale && (y= word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1)
 
1839
        word0(&adj)+= (2 * P + 1) * Exp_msk1 - y;
 
1840
      adj.d*= ulp(&rv);
 
1841
      if (dsign)
 
1842
        dval(&rv)+= adj.d;
 
1843
      else
 
1844
        dval(&rv)-= adj.d;
 
1845
      goto cont;
 
1846
    }
 
1847
#endif /*Honor_FLT_ROUNDS*/
 
1848
 
 
1849
    if (i < 0)
 
1850
    {
 
1851
      /*
 
1852
        Error is less than half an ulp -- check for special case of mantissa
 
1853
        a power of two.
 
1854
      */
 
1855
      if (dsign || word1(&rv) || word0(&rv) & Bndry_mask ||
 
1856
          (word0(&rv) & Exp_mask) <= (2 * P + 1) * Exp_msk1)
 
1857
      {
 
1858
#ifdef SET_INEXACT
 
1859
        if (!delta->x[0] && delta->wds <= 1)
 
1860
          inexact= 0;
 
1861
#endif
 
1862
        break;
 
1863
      }
 
1864
      if (!delta->p.x[0] && delta->wds <= 1)
 
1865
      {
 
1866
        /* exact result */
 
1867
#ifdef SET_INEXACT
 
1868
        inexact= 0;
 
1869
#endif
 
1870
        break;
 
1871
      }
 
1872
      delta= lshift(delta, Log2P, &alloc);
 
1873
      if (cmp(delta, bs) > 0)
 
1874
        goto drop_down;
 
1875
      break;
 
1876
    }
 
1877
    if (i == 0)
 
1878
    {
 
1879
      /* exactly half-way between */
 
1880
      if (dsign)
 
1881
      {
 
1882
        if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 &&
 
1883
            word1(&rv) ==
 
1884
            ((scale && (y = word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1) ?
 
1885
             (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
 
1886
             0xffffffff))
 
1887
        {
 
1888
          /*boundary case -- increment exponent*/
 
1889
          word0(&rv)= (word0(&rv) & Exp_mask) + Exp_msk1;
 
1890
          word1(&rv) = 0;
 
1891
          dsign = 0;
 
1892
          break;
 
1893
        }
 
1894
      }
 
1895
      else if (!(word0(&rv) & Bndry_mask) && !word1(&rv))
 
1896
      {
 
1897
 drop_down:
 
1898
        /* boundary case -- decrement exponent */
 
1899
        if (scale)
 
1900
        {
 
1901
          L= word0(&rv) & Exp_mask;
 
1902
          if (L <= (2 *P + 1) * Exp_msk1)
 
1903
          {
 
1904
            if (L > (P + 2) * Exp_msk1)
 
1905
              /* round even ==> accept rv */
 
1906
              break;
 
1907
            /* rv = smallest denormal */
 
1908
            goto undfl;
 
1909
          }
 
1910
        }
 
1911
        L= (word0(&rv) & Exp_mask) - Exp_msk1;
 
1912
        word0(&rv)= L | Bndry_mask1;
 
1913
        word1(&rv)= 0xffffffff;
 
1914
        break;
 
1915
      }
 
1916
      if (!(word1(&rv) & LSB))
 
1917
        break;
 
1918
      if (dsign)
 
1919
        dval(&rv)+= ulp(&rv);
 
1920
      else
 
1921
      {
 
1922
        dval(&rv)-= ulp(&rv);
 
1923
        if (!dval(&rv))
 
1924
          goto undfl;
 
1925
      }
 
1926
      dsign= 1 - dsign;
 
1927
      break;
 
1928
    }
 
1929
    if ((aadj= ratio(delta, bs)) <= 2.)
 
1930
    {
 
1931
      if (dsign)
 
1932
        aadj= aadj1= 1.;
 
1933
      else if (word1(&rv) || word0(&rv) & Bndry_mask)
 
1934
      {
 
1935
        if (word1(&rv) == Tiny1 && !word0(&rv))
 
1936
          goto undfl;
 
1937
        aadj= 1.;
 
1938
        aadj1= -1.;
 
1939
      }
 
1940
      else
 
1941
      {
 
1942
        /* special case -- power of FLT_RADIX to be rounded down... */
 
1943
        if (aadj < 2. / FLT_RADIX)
 
1944
          aadj= 1. / FLT_RADIX;
 
1945
        else
 
1946
          aadj*= 0.5;
 
1947
        aadj1= -aadj;
 
1948
      }
 
1949
    }
 
1950
    else
 
1951
    {
 
1952
      aadj*= 0.5;
 
1953
      aadj1= dsign ? aadj : -aadj;
 
1954
#ifdef Check_FLT_ROUNDS
 
1955
      switch (Rounding)
 
1956
      {
 
1957
      case 2: /* towards +infinity */
 
1958
        aadj1-= 0.5;
 
1959
        break;
 
1960
      case 0: /* towards 0 */
 
1961
      case 3: /* towards -infinity */
 
1962
        aadj1+= 0.5;
 
1963
      }
 
1964
#else
 
1965
      if (Flt_Rounds == 0)
 
1966
        aadj1+= 0.5;
 
1967
#endif /*Check_FLT_ROUNDS*/
 
1968
    }
 
1969
    y= word0(&rv) & Exp_mask;
 
1970
 
 
1971
    /* Check for overflow */
 
1972
 
 
1973
    if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
 
1974
    {
 
1975
      dval(&rv0)= dval(&rv);
 
1976
      word0(&rv)-= P * Exp_msk1;
 
1977
      adj.d= aadj1 * ulp(&rv);
 
1978
      dval(&rv)+= adj.d;
 
1979
      if ((word0(&rv) & Exp_mask) >= Exp_msk1 * (DBL_MAX_EXP + Bias - P))
 
1980
      {
 
1981
        if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
 
1982
          goto ovfl;
 
1983
        word0(&rv)= Big0;
 
1984
        word1(&rv)= Big1;
 
1985
        goto cont;
 
1986
      }
 
1987
      else
 
1988
        word0(&rv)+= P * Exp_msk1;
 
1989
    }
 
1990
    else
 
1991
    {
 
1992
      if (scale && y <= 2 * P * Exp_msk1)
 
1993
      {
 
1994
        if (aadj <= 0x7fffffff)
 
1995
        {
 
1996
          if ((z= (ULong) aadj) <= 0)
 
1997
            z= 1;
 
1998
          aadj= z;
 
1999
          aadj1= dsign ? aadj : -aadj;
 
2000
        }
 
2001
        dval(&aadj2) = aadj1;
 
2002
        word0(&aadj2)+= (2 * P + 1) * Exp_msk1 - y;
 
2003
        aadj1= dval(&aadj2);
 
2004
        adj.d= aadj1 * ulp(&rv);
 
2005
        dval(&rv)+= adj.d;
 
2006
        if (rv.d == 0.)
 
2007
          goto undfl;
 
2008
      }
 
2009
      else
 
2010
      {
 
2011
        adj.d= aadj1 * ulp(&rv);
 
2012
        dval(&rv)+= adj.d;
 
2013
      }
 
2014
    }
 
2015
    z= word0(&rv) & Exp_mask;
 
2016
#ifndef SET_INEXACT
 
2017
    if (!scale)
 
2018
      if (y == z)
 
2019
      {
 
2020
        /* Can we stop now? */
 
2021
        L= (Long)aadj;
 
2022
        aadj-= L;
 
2023
        /* The tolerances below are conservative. */
 
2024
        if (dsign || word1(&rv) || word0(&rv) & Bndry_mask)
 
2025
        {
 
2026
          if (aadj < .4999999 || aadj > .5000001)
 
2027
            break;
 
2028
        }
 
2029
        else if (aadj < .4999999 / FLT_RADIX)
 
2030
          break;
 
2031
      }
 
2032
#endif
 
2033
 cont:
 
2034
    Bfree(bb, &alloc);
 
2035
    Bfree(bd, &alloc);
 
2036
    Bfree(bs, &alloc);
 
2037
    Bfree(delta, &alloc);
 
2038
  }
 
2039
#ifdef SET_INEXACT
 
2040
  if (inexact)
 
2041
  {
 
2042
    if (!oldinexact)
 
2043
    {
 
2044
      word0(&rv0)= Exp_1 + (70 << Exp_shift);
 
2045
      word1(&rv0)= 0;
 
2046
      dval(&rv0)+= 1.;
 
2047
    }
 
2048
  }
 
2049
  else if (!oldinexact)
 
2050
    clear_inexact();
 
2051
#endif
 
2052
  if (scale)
 
2053
  {
 
2054
    word0(&rv0)= Exp_1 - 2 * P * Exp_msk1;
 
2055
    word1(&rv0)= 0;
 
2056
    dval(&rv)*= dval(&rv0);
 
2057
  }
 
2058
#ifdef SET_INEXACT
 
2059
  if (inexact && !(word0(&rv) & Exp_mask))
 
2060
  {
 
2061
    /* set underflow bit */
 
2062
    dval(&rv0)= 1e-300;
 
2063
    dval(&rv0)*= dval(&rv0);
 
2064
  }
 
2065
#endif
 
2066
 retfree:
 
2067
  Bfree(bb, &alloc);
 
2068
  Bfree(bd, &alloc);
 
2069
  Bfree(bs, &alloc);
 
2070
  Bfree(bd0, &alloc);
 
2071
  Bfree(delta, &alloc);
 
2072
 ret:
 
2073
  *se= (char *)s;
 
2074
  return sign ? -dval(&rv) : dval(&rv);
 
2075
}
 
2076
 
 
2077
 
 
2078
static int quorem(Bigint *b, Bigint *S)
 
2079
{
 
2080
  int n;
 
2081
  ULong *bx, *bxe, q, *sx, *sxe;
 
2082
  ULLong borrow, carry, y, ys;
 
2083
 
 
2084
  n= S->wds;
 
2085
  if (b->wds < n)
 
2086
    return 0;
 
2087
  sx= S->p.x;
 
2088
  sxe= sx + --n;
 
2089
  bx= b->p.x;
 
2090
  bxe= bx + n;
 
2091
  q= *bxe / (*sxe + 1);  /* ensure q <= true quotient */
 
2092
  if (q)
 
2093
  {
 
2094
    borrow= 0;
 
2095
    carry= 0;
 
2096
    do
 
2097
    {
 
2098
      ys= *sx++ * (ULLong)q + carry;
 
2099
      carry= ys >> 32;
 
2100
      y= *bx - (ys & FFFFFFFF) - borrow;
 
2101
      borrow= y >> 32 & (ULong)1;
 
2102
      *bx++= (ULong) (y & FFFFFFFF);
 
2103
    }
 
2104
    while (sx <= sxe);
 
2105
    if (!*bxe)
 
2106
    {
 
2107
      bx= b->p.x;
 
2108
      while (--bxe > bx && !*bxe)
 
2109
        --n;
 
2110
      b->wds= n;
 
2111
    }
 
2112
  }
 
2113
  if (cmp(b, S) >= 0)
 
2114
  {
 
2115
    q++;
 
2116
    borrow= 0;
 
2117
    carry= 0;
 
2118
    bx= b->p.x;
 
2119
    sx= S->p.x;
 
2120
    do
 
2121
    {
 
2122
      ys= *sx++ + carry;
 
2123
      carry= ys >> 32;
 
2124
      y= *bx - (ys & FFFFFFFF) - borrow;
 
2125
      borrow= y >> 32 & (ULong)1;
 
2126
      *bx++= (ULong) (y & FFFFFFFF);
 
2127
    }
 
2128
    while (sx <= sxe);
 
2129
    bx= b->p.x;
 
2130
    bxe= bx + n;
 
2131
    if (!*bxe)
 
2132
    {
 
2133
      while (--bxe > bx && !*bxe)
 
2134
        --n;
 
2135
      b->wds= n;
 
2136
    }
 
2137
  }
 
2138
  return q;
 
2139
}
 
2140
 
 
2141
 
 
2142
/*
 
2143
   dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
 
2144
 
 
2145
   Inspired by "How to Print Floating-Point Numbers Accurately" by
 
2146
   Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
 
2147
 
 
2148
   Modifications:
 
2149
        1. Rather than iterating, we use a simple numeric overestimate
 
2150
           to determine k= floor(log10(d)).  We scale relevant
 
2151
           quantities using O(log2(k)) rather than O(k) multiplications.
 
2152
        2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
 
2153
           try to generate digits strictly left to right.  Instead, we
 
2154
           compute with fewer bits and propagate the carry if necessary
 
2155
           when rounding the final digit up.  This is often faster.
 
2156
        3. Under the assumption that input will be rounded nearest,
 
2157
           mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
 
2158
           That is, we allow equality in stopping tests when the
 
2159
           round-nearest rule will give the same floating-point value
 
2160
           as would satisfaction of the stopping test with strict
 
2161
           inequality.
 
2162
        4. We remove common factors of powers of 2 from relevant
 
2163
           quantities.
 
2164
        5. When converting floating-point integers less than 1e16,
 
2165
           we use floating-point arithmetic rather than resorting
 
2166
           to multiple-precision integers.
 
2167
        6. When asked to produce fewer than 15 digits, we first try
 
2168
           to get by with floating-point arithmetic; we resort to
 
2169
           multiple-precision integer arithmetic only if we cannot
 
2170
           guarantee that the floating-point calculation has given
 
2171
           the correctly rounded result.  For k requested digits and
 
2172
           "uniformly" distributed input, the probability is
 
2173
           something like 10^(k-15) that we must resort to the Long
 
2174
           calculation.
 
2175
 */
 
2176
 
 
2177
static char *dtoa(double dd, int mode, int ndigits, int *decpt, int *sign,
 
2178
                  char **rve, char *buf, size_t buf_size)
 
2179
{
 
2180
  /*
 
2181
    Arguments ndigits, decpt, sign are similar to those
 
2182
    of ecvt and fcvt; trailing zeros are suppressed from
 
2183
    the returned string.  If not null, *rve is set to point
 
2184
    to the end of the return value.  If d is +-Infinity or NaN,
 
2185
    then *decpt is set to DTOA_OVERFLOW.
 
2186
 
 
2187
    mode:
 
2188
          0 ==> shortest string that yields d when read in
 
2189
                and rounded to nearest.
 
2190
          1 ==> like 0, but with Steele & White stopping rule;
 
2191
                e.g. with IEEE P754 arithmetic , mode 0 gives
 
2192
                1e23 whereas mode 1 gives 9.999999999999999e22.
 
2193
          2 ==> max(1,ndigits) significant digits.  This gives a
 
2194
                return value similar to that of ecvt, except
 
2195
                that trailing zeros are suppressed.
 
2196
          3 ==> through ndigits past the decimal point.  This
 
2197
                gives a return value similar to that from fcvt,
 
2198
                except that trailing zeros are suppressed, and
 
2199
                ndigits can be negative.
 
2200
          4,5 ==> similar to 2 and 3, respectively, but (in
 
2201
                round-nearest mode) with the tests of mode 0 to
 
2202
                possibly return a shorter string that rounds to d.
 
2203
                With IEEE arithmetic and compilation with
 
2204
                -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
 
2205
                as modes 2 and 3 when FLT_ROUNDS != 1.
 
2206
          6-9 ==> Debugging modes similar to mode - 4:  don't try
 
2207
                fast floating-point estimate (if applicable).
 
2208
 
 
2209
      Values of mode other than 0-9 are treated as mode 0.
 
2210
 
 
2211
    Sufficient space is allocated to the return value
 
2212
    to hold the suppressed trailing zeros.
 
2213
  */
 
2214
 
 
2215
  int bbits, b2, b5, be, dig, i, ieps, UNINIT_VAR(ilim), ilim0, 
 
2216
    UNINIT_VAR(ilim1), j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
 
2217
    spec_case, try_quick;
 
2218
  Long L;
 
2219
  int denorm;
 
2220
  ULong x;
 
2221
  Bigint *b, *b1, *delta, *mlo, *mhi, *S;
 
2222
  U d2, eps, u;
 
2223
  double ds;
 
2224
  char *s, *s0;
 
2225
#ifdef Honor_FLT_ROUNDS
 
2226
  int rounding;
 
2227
#endif
 
2228
  Stack_alloc alloc;
 
2229
  
 
2230
  alloc.begin= alloc.free= buf;
 
2231
  alloc.end= buf + buf_size;
 
2232
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
 
2233
 
 
2234
  u.d= dd;
 
2235
  if (word0(&u) & Sign_bit)
 
2236
  {
 
2237
    /* set sign for everything, including 0's and NaNs */
 
2238
    *sign= 1;
 
2239
    word0(&u) &= ~Sign_bit;  /* clear sign bit */
 
2240
  }
 
2241
  else
 
2242
    *sign= 0;
 
2243
 
 
2244
  /* If infinity, set decpt to DTOA_OVERFLOW, if 0 set it to 1 */
 
2245
  if (((word0(&u) & Exp_mask) == Exp_mask && (*decpt= DTOA_OVERFLOW)) ||
 
2246
      (!dval(&u) && (*decpt= 1)))
 
2247
  {
 
2248
    /* Infinity, NaN, 0 */
 
2249
    char *res= (char*) dtoa_alloc(2, &alloc);
 
2250
    res[0]= '0';
 
2251
    res[1]= '\0';
 
2252
    if (rve)
 
2253
      *rve= res + 1;
 
2254
    return res;
 
2255
  }
 
2256
  
 
2257
#ifdef Honor_FLT_ROUNDS
 
2258
  if ((rounding= Flt_Rounds) >= 2)
 
2259
  {
 
2260
    if (*sign)
 
2261
      rounding= rounding == 2 ? 0 : 2;
 
2262
    else
 
2263
      if (rounding != 2)
 
2264
        rounding= 0;
 
2265
  }
 
2266
#endif
 
2267
 
 
2268
  b= d2b(&u, &be, &bbits, &alloc);
 
2269
  if ((i= (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
 
2270
  {
 
2271
    dval(&d2)= dval(&u);
 
2272
    word0(&d2) &= Frac_mask1;
 
2273
    word0(&d2) |= Exp_11;
 
2274
 
 
2275
    /*
 
2276
      log(x)       ~=~ log(1.5) + (x-1.5)/1.5
 
2277
      log10(x)      =  log(x) / log(10)
 
2278
                   ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
 
2279
      log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
 
2280
     
 
2281
      This suggests computing an approximation k to log10(d) by
 
2282
     
 
2283
      k= (i - Bias)*0.301029995663981
 
2284
           + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
 
2285
     
 
2286
      We want k to be too large rather than too small.
 
2287
      The error in the first-order Taylor series approximation
 
2288
      is in our favor, so we just round up the constant enough
 
2289
      to compensate for any error in the multiplication of
 
2290
      (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
 
2291
      and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
 
2292
      adding 1e-13 to the constant term more than suffices.
 
2293
      Hence we adjust the constant term to 0.1760912590558.
 
2294
      (We could get a more accurate k by invoking log10,
 
2295
       but this is probably not worthwhile.)
 
2296
    */
 
2297
 
 
2298
    i-= Bias;
 
2299
    denorm= 0;
 
2300
  }
 
2301
  else
 
2302
  {
 
2303
    /* d is denormalized */
 
2304
 
 
2305
    i= bbits + be + (Bias + (P-1) - 1);
 
2306
    x= i > 32  ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
 
2307
      : word1(&u) << (32 - i);
 
2308
    dval(&d2)= x;
 
2309
    word0(&d2)-= 31*Exp_msk1; /* adjust exponent */
 
2310
    i-= (Bias + (P-1) - 1) + 1;
 
2311
    denorm= 1;
 
2312
  }
 
2313
  ds= (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
 
2314
  k= (int)ds;
 
2315
  if (ds < 0. && ds != k)
 
2316
    k--;    /* want k= floor(ds) */
 
2317
  k_check= 1;
 
2318
  if (k >= 0 && k <= Ten_pmax)
 
2319
  {
 
2320
    if (dval(&u) < tens[k])
 
2321
      k--;
 
2322
    k_check= 0;
 
2323
  }
 
2324
  j= bbits - i - 1;
 
2325
  if (j >= 0)
 
2326
  {
 
2327
    b2= 0;
 
2328
    s2= j;
 
2329
  }
 
2330
  else
 
2331
  {
 
2332
    b2= -j;
 
2333
    s2= 0;
 
2334
  }
 
2335
  if (k >= 0)
 
2336
  {
 
2337
    b5= 0;
 
2338
    s5= k;
 
2339
    s2+= k;
 
2340
  }
 
2341
  else
 
2342
  {
 
2343
    b2-= k;
 
2344
    b5= -k;
 
2345
    s5= 0;
 
2346
  }
 
2347
  if (mode < 0 || mode > 9)
 
2348
    mode= 0;
 
2349
 
 
2350
#ifdef Check_FLT_ROUNDS
 
2351
  try_quick= Rounding == 1;
 
2352
#else
 
2353
  try_quick= 1;
 
2354
#endif
 
2355
 
 
2356
  if (mode > 5)
 
2357
  {
 
2358
    mode-= 4;
 
2359
    try_quick= 0;
 
2360
  }
 
2361
  leftright= 1;
 
2362
  switch (mode) {
 
2363
  case 0:
 
2364
  case 1:
 
2365
    ilim= ilim1= -1;
 
2366
    i= 18;
 
2367
    ndigits= 0;
 
2368
    break;
 
2369
  case 2:
 
2370
    leftright= 0;
 
2371
    /* no break */
 
2372
  case 4:
 
2373
    if (ndigits <= 0)
 
2374
      ndigits= 1;
 
2375
    ilim= ilim1= i= ndigits;
 
2376
    break;
 
2377
  case 3:
 
2378
    leftright= 0;
 
2379
    /* no break */
 
2380
  case 5:
 
2381
    i= ndigits + k + 1;
 
2382
    ilim= i;
 
2383
    ilim1= i - 1;
 
2384
    if (i <= 0)
 
2385
      i= 1;
 
2386
  }
 
2387
  s= s0= dtoa_alloc(i, &alloc);
 
2388
 
 
2389
#ifdef Honor_FLT_ROUNDS
 
2390
  if (mode > 1 && rounding != 1)
 
2391
    leftright= 0;
 
2392
#endif
 
2393
 
 
2394
  if (ilim >= 0 && ilim <= Quick_max && try_quick)
 
2395
  {
 
2396
    /* Try to get by with floating-point arithmetic. */
 
2397
    i= 0;
 
2398
    dval(&d2)= dval(&u);
 
2399
    k0= k;
 
2400
    ilim0= ilim;
 
2401
    ieps= 2; /* conservative */
 
2402
    if (k > 0)
 
2403
    {
 
2404
      ds= tens[k&0xf];
 
2405
      j= k >> 4;
 
2406
      if (j & Bletch)
 
2407
      {
 
2408
        /* prevent overflows */
 
2409
        j&= Bletch - 1;
 
2410
        dval(&u)/= bigtens[n_bigtens-1];
 
2411
        ieps++;
 
2412
      }
 
2413
      for (; j; j>>= 1, i++)
 
2414
      {
 
2415
        if (j & 1)
 
2416
        {
 
2417
          ieps++;
 
2418
          ds*= bigtens[i];
 
2419
        }
 
2420
      }
 
2421
      dval(&u)/= ds;
 
2422
    }
 
2423
    else if ((j1= -k))
 
2424
    {
 
2425
      dval(&u)*= tens[j1 & 0xf];
 
2426
      for (j= j1 >> 4; j; j>>= 1, i++)
 
2427
      {
 
2428
        if (j & 1)
 
2429
        {
 
2430
          ieps++;
 
2431
          dval(&u)*= bigtens[i];
 
2432
        }
 
2433
      }
 
2434
    }
 
2435
    if (k_check && dval(&u) < 1. && ilim > 0)
 
2436
    {
 
2437
      if (ilim1 <= 0)
 
2438
        goto fast_failed;
 
2439
      ilim= ilim1;
 
2440
      k--;
 
2441
      dval(&u)*= 10.;
 
2442
      ieps++;
 
2443
    }
 
2444
    dval(&eps)= ieps*dval(&u) + 7.;
 
2445
    word0(&eps)-= (P-1)*Exp_msk1;
 
2446
    if (ilim == 0)
 
2447
    {
 
2448
      S= mhi= 0;
 
2449
      dval(&u)-= 5.;
 
2450
      if (dval(&u) > dval(&eps))
 
2451
        goto one_digit;
 
2452
      if (dval(&u) < -dval(&eps))
 
2453
        goto no_digits;
 
2454
      goto fast_failed;
 
2455
    }
 
2456
    if (leftright)
 
2457
    {
 
2458
      /* Use Steele & White method of only generating digits needed. */
 
2459
      dval(&eps)= 0.5/tens[ilim-1] - dval(&eps);
 
2460
      for (i= 0;;)
 
2461
      {
 
2462
        L= (Long) dval(&u);
 
2463
        dval(&u)-= L;
 
2464
        *s++= '0' + (int)L;
 
2465
        if (dval(&u) < dval(&eps))
 
2466
          goto ret1;
 
2467
        if (1. - dval(&u) < dval(&eps))
 
2468
          goto bump_up;
 
2469
        if (++i >= ilim)
 
2470
          break;
 
2471
        dval(&eps)*= 10.;
 
2472
        dval(&u)*= 10.;
 
2473
      }
 
2474
    }
 
2475
    else
 
2476
    {
 
2477
      /* Generate ilim digits, then fix them up. */
 
2478
      dval(&eps)*= tens[ilim-1];
 
2479
      for (i= 1;; i++, dval(&u)*= 10.)
 
2480
      {
 
2481
        L= (Long)(dval(&u));
 
2482
        if (!(dval(&u)-= L))
 
2483
          ilim= i;
 
2484
        *s++= '0' + (int)L;
 
2485
        if (i == ilim)
 
2486
        {
 
2487
          if (dval(&u) > 0.5 + dval(&eps))
 
2488
            goto bump_up;
 
2489
          else if (dval(&u) < 0.5 - dval(&eps))
 
2490
          {
 
2491
            while (*--s == '0');
 
2492
            s++;
 
2493
            goto ret1;
 
2494
          }
 
2495
          break;
 
2496
        }
 
2497
      }
 
2498
    }
 
2499
  fast_failed:
 
2500
    s= s0;
 
2501
    dval(&u)= dval(&d2);
 
2502
    k= k0;
 
2503
    ilim= ilim0;
 
2504
  }
 
2505
 
 
2506
  /* Do we have a "small" integer? */
 
2507
 
 
2508
  if (be >= 0 && k <= Int_max)
 
2509
  {
 
2510
    /* Yes. */
 
2511
    ds= tens[k];
 
2512
    if (ndigits < 0 && ilim <= 0)
 
2513
    {
 
2514
      S= mhi= 0;
 
2515
      if (ilim < 0 || dval(&u) <= 5*ds)
 
2516
        goto no_digits;
 
2517
      goto one_digit;
 
2518
    }
 
2519
    for (i= 1;; i++, dval(&u)*= 10.)
 
2520
    {
 
2521
      L= (Long)(dval(&u) / ds);
 
2522
      dval(&u)-= L*ds;
 
2523
#ifdef Check_FLT_ROUNDS
 
2524
      /* If FLT_ROUNDS == 2, L will usually be high by 1 */
 
2525
      if (dval(&u) < 0)
 
2526
      {
 
2527
        L--;
 
2528
        dval(&u)+= ds;
 
2529
      }
 
2530
#endif
 
2531
      *s++= '0' + (int)L;
 
2532
      if (!dval(&u))
 
2533
      {
 
2534
        break;
 
2535
      }
 
2536
      if (i == ilim)
 
2537
      {
 
2538
#ifdef Honor_FLT_ROUNDS
 
2539
        if (mode > 1)
 
2540
        {
 
2541
          switch (rounding) {
 
2542
          case 0: goto ret1;
 
2543
          case 2: goto bump_up;
 
2544
          }
 
2545
        }
 
2546
#endif
 
2547
        dval(&u)+= dval(&u);
 
2548
        if (dval(&u) > ds || (dval(&u) == ds && L & 1))
 
2549
        {
 
2550
bump_up:
 
2551
          while (*--s == '9')
 
2552
            if (s == s0)
 
2553
            {
 
2554
              k++;
 
2555
              *s= '0';
 
2556
              break;
 
2557
            }
 
2558
          ++*s++;
 
2559
        }
 
2560
        break;
 
2561
      }
 
2562
    }
 
2563
    goto ret1;
 
2564
  }
 
2565
 
 
2566
  m2= b2;
 
2567
  m5= b5;
 
2568
  mhi= mlo= 0;
 
2569
  if (leftright)
 
2570
  {
 
2571
    i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
 
2572
    b2+= i;
 
2573
    s2+= i;
 
2574
    mhi= i2b(1, &alloc);
 
2575
  }
 
2576
  if (m2 > 0 && s2 > 0)
 
2577
  {
 
2578
    i= m2 < s2 ? m2 : s2;
 
2579
    b2-= i;
 
2580
    m2-= i;
 
2581
    s2-= i;
 
2582
  }
 
2583
  if (b5 > 0)
 
2584
  {
 
2585
    if (leftright)
 
2586
    {
 
2587
      if (m5 > 0)
 
2588
      {
 
2589
        mhi= pow5mult(mhi, m5, &alloc);
 
2590
        b1= mult(mhi, b, &alloc);
 
2591
        Bfree(b, &alloc);
 
2592
        b= b1;
 
2593
      }
 
2594
      if ((j= b5 - m5))
 
2595
        b= pow5mult(b, j, &alloc);
 
2596
    }
 
2597
    else
 
2598
      b= pow5mult(b, b5, &alloc);
 
2599
  }
 
2600
  S= i2b(1, &alloc);
 
2601
  if (s5 > 0)
 
2602
    S= pow5mult(S, s5, &alloc);
 
2603
 
 
2604
  /* Check for special case that d is a normalized power of 2. */
 
2605
 
 
2606
  spec_case= 0;
 
2607
  if ((mode < 2 || leftright)
 
2608
#ifdef Honor_FLT_ROUNDS
 
2609
      && rounding == 1
 
2610
#endif
 
2611
     )
 
2612
  {
 
2613
    if (!word1(&u) && !(word0(&u) & Bndry_mask) &&
 
2614
        word0(&u) & (Exp_mask & ~Exp_msk1)
 
2615
       )
 
2616
    {
 
2617
      /* The special case */
 
2618
      b2+= Log2P;
 
2619
      s2+= Log2P;
 
2620
      spec_case= 1;
 
2621
    }
 
2622
  }
 
2623
 
 
2624
  /*
 
2625
    Arrange for convenient computation of quotients:
 
2626
    shift left if necessary so divisor has 4 leading 0 bits.
 
2627
    
 
2628
    Perhaps we should just compute leading 28 bits of S once
 
2629
    a nd for all and pass them and a shift to quorem, so it
 
2630
    can do shifts and ors to compute the numerator for q.
 
2631
  */
 
2632
  if ((i= ((s5 ? 32 - hi0bits(S->p.x[S->wds-1]) : 1) + s2) & 0x1f))
 
2633
    i= 32 - i;
 
2634
  if (i > 4)
 
2635
  {
 
2636
    i-= 4;
 
2637
    b2+= i;
 
2638
    m2+= i;
 
2639
    s2+= i;
 
2640
  }
 
2641
  else if (i < 4)
 
2642
  {
 
2643
    i+= 28;
 
2644
    b2+= i;
 
2645
    m2+= i;
 
2646
    s2+= i;
 
2647
  }
 
2648
  if (b2 > 0)
 
2649
    b= lshift(b, b2, &alloc);
 
2650
  if (s2 > 0)
 
2651
    S= lshift(S, s2, &alloc);
 
2652
  if (k_check)
 
2653
  {
 
2654
    if (cmp(b,S) < 0)
 
2655
    {
 
2656
      k--;
 
2657
      /* we botched the k estimate */
 
2658
      b= multadd(b, 10, 0, &alloc);
 
2659
      if (leftright)
 
2660
        mhi= multadd(mhi, 10, 0, &alloc);
 
2661
      ilim= ilim1;
 
2662
    }
 
2663
  }
 
2664
  if (ilim <= 0 && (mode == 3 || mode == 5))
 
2665
  {
 
2666
    if (ilim < 0 || cmp(b,S= multadd(S,5,0, &alloc)) <= 0)
 
2667
    {
 
2668
      /* no digits, fcvt style */
 
2669
no_digits:
 
2670
      k= -1 - ndigits;
 
2671
      goto ret;
 
2672
    }
 
2673
one_digit:
 
2674
    *s++= '1';
 
2675
    k++;
 
2676
    goto ret;
 
2677
  }
 
2678
  if (leftright)
 
2679
  {
 
2680
    if (m2 > 0)
 
2681
      mhi= lshift(mhi, m2, &alloc);
 
2682
 
 
2683
    /*
 
2684
      Compute mlo -- check for special case that d is a normalized power of 2.
 
2685
    */
 
2686
 
 
2687
    mlo= mhi;
 
2688
    if (spec_case)
 
2689
    {
 
2690
      mhi= Balloc(mhi->k, &alloc);
 
2691
      Bcopy(mhi, mlo);
 
2692
      mhi= lshift(mhi, Log2P, &alloc);
 
2693
    }
 
2694
 
 
2695
    for (i= 1;;i++)
 
2696
    {
 
2697
      dig= quorem(b,S) + '0';
 
2698
      /* Do we yet have the shortest decimal string that will round to d? */
 
2699
      j= cmp(b, mlo);
 
2700
      delta= diff(S, mhi, &alloc);
 
2701
      j1= delta->sign ? 1 : cmp(b, delta);
 
2702
      Bfree(delta, &alloc);
 
2703
      if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
 
2704
#ifdef Honor_FLT_ROUNDS
 
2705
          && rounding >= 1
 
2706
#endif
 
2707
         )
 
2708
      {
 
2709
        if (dig == '9')
 
2710
          goto round_9_up;
 
2711
        if (j > 0)
 
2712
          dig++;
 
2713
        *s++= dig;
 
2714
        goto ret;
 
2715
      }
 
2716
      if (j < 0 || (j == 0 && mode != 1 && !(word1(&u) & 1)))
 
2717
      {
 
2718
        if (!b->p.x[0] && b->wds <= 1)
 
2719
        {
 
2720
          goto accept_dig;
 
2721
        }
 
2722
#ifdef Honor_FLT_ROUNDS
 
2723
        if (mode > 1)
 
2724
          switch (rounding) {
 
2725
          case 0: goto accept_dig;
 
2726
          case 2: goto keep_dig;
 
2727
          }
 
2728
#endif /*Honor_FLT_ROUNDS*/
 
2729
        if (j1 > 0)
 
2730
        {
 
2731
          b= lshift(b, 1, &alloc);
 
2732
          j1= cmp(b, S);
 
2733
          if ((j1 > 0 || (j1 == 0 && dig & 1))
 
2734
              && dig++ == '9')
 
2735
            goto round_9_up;
 
2736
        }
 
2737
accept_dig:
 
2738
        *s++= dig;
 
2739
        goto ret;
 
2740
      }
 
2741
      if (j1 > 0)
 
2742
      {
 
2743
#ifdef Honor_FLT_ROUNDS
 
2744
        if (!rounding)
 
2745
          goto accept_dig;
 
2746
#endif
 
2747
        if (dig == '9')
 
2748
        { /* possible if i == 1 */
 
2749
round_9_up:
 
2750
          *s++= '9';
 
2751
          goto roundoff;
 
2752
        }
 
2753
        *s++= dig + 1;
 
2754
        goto ret;
 
2755
      }
 
2756
#ifdef Honor_FLT_ROUNDS
 
2757
keep_dig:
 
2758
#endif
 
2759
      *s++= dig;
 
2760
      if (i == ilim)
 
2761
        break;
 
2762
      b= multadd(b, 10, 0, &alloc);
 
2763
      if (mlo == mhi)
 
2764
        mlo= mhi= multadd(mhi, 10, 0, &alloc);
 
2765
      else
 
2766
      {
 
2767
        mlo= multadd(mlo, 10, 0, &alloc);
 
2768
        mhi= multadd(mhi, 10, 0, &alloc);
 
2769
      }
 
2770
    }
 
2771
  }
 
2772
  else
 
2773
    for (i= 1;; i++)
 
2774
    {
 
2775
      *s++= dig= quorem(b,S) + '0';
 
2776
      if (!b->p.x[0] && b->wds <= 1)
 
2777
      {
 
2778
        goto ret;
 
2779
      }
 
2780
      if (i >= ilim)
 
2781
        break;
 
2782
      b= multadd(b, 10, 0, &alloc);
 
2783
    }
 
2784
 
 
2785
  /* Round off last digit */
 
2786
 
 
2787
#ifdef Honor_FLT_ROUNDS
 
2788
  switch (rounding) {
 
2789
  case 0: goto trimzeros;
 
2790
  case 2: goto roundoff;
 
2791
  }
 
2792
#endif
 
2793
  b= lshift(b, 1, &alloc);
 
2794
  j= cmp(b, S);
 
2795
  if (j > 0 || (j == 0 && dig & 1))
 
2796
  {
 
2797
roundoff:
 
2798
    while (*--s == '9')
 
2799
      if (s == s0)
 
2800
      {
 
2801
        k++;
 
2802
        *s++= '1';
 
2803
        goto ret;
 
2804
      }
 
2805
    ++*s++;
 
2806
  }
 
2807
  else
 
2808
  {
 
2809
#ifdef Honor_FLT_ROUNDS
 
2810
trimzeros:
 
2811
#endif
 
2812
    while (*--s == '0');
 
2813
    s++;
 
2814
  }
 
2815
ret:
 
2816
  Bfree(S, &alloc);
 
2817
  if (mhi)
 
2818
  {
 
2819
    if (mlo && mlo != mhi)
 
2820
      Bfree(mlo, &alloc);
 
2821
    Bfree(mhi, &alloc);
 
2822
  }
 
2823
ret1:
 
2824
  Bfree(b, &alloc);
 
2825
  *s= 0;
 
2826
  *decpt= k + 1;
 
2827
  if (rve)
 
2828
    *rve= s;
 
2829
  return s0;
 
2830
}