← Back to branch summary

~drizzle-trunk/drizzle/development

« back to all changes in this revision

Viewing changes to strings/dtoa.c

  • Committer: brian
  • Date: 2008-06-25 05:29:13 UTC
  • Revision ID: brian@localhost.localdomain-20080625052913-6upwo0jsrl4lnapl
clean slate

Show diffs side-by-side

added added

removed removed

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