1/****************************************************************************
2**
3** Copyright (C) 2014 Digia Plc and/or its subsidiary(-ies).
4** Contact: http://www.qt-project.org/legal
5**
6** This file is part of the QtCore module of the Qt Toolkit.
7**
8** $QT_BEGIN_LICENSE:LGPL$
9** Commercial License Usage
10** Licensees holding valid commercial Qt licenses may use this file in
11** accordance with the commercial license agreement provided with the
12** Software or, alternatively, in accordance with the terms contained in
13** a written agreement between you and Digia. For licensing terms and
14** conditions see http://qt.digia.com/licensing. For further information
15** use the contact form at http://qt.digia.com/contact-us.
16**
17** GNU Lesser General Public License Usage
18** Alternatively, this file may be used under the terms of the GNU Lesser
19** General Public License version 2.1 as published by the Free Software
20** Foundation and appearing in the file LICENSE.LGPL included in the
21** packaging of this file. Please review the following information to
22** ensure the GNU Lesser General Public License version 2.1 requirements
23** will be met: http://www.gnu.org/licenses/old-licenses/lgpl-2.1.html.
24**
25** In addition, as a special exception, Digia gives you certain additional
26** rights. These rights are described in the Digia Qt LGPL Exception
27** version 1.1, included in the file LGPL_EXCEPTION.txt in this package.
28**
29** GNU General Public License Usage
30** Alternatively, this file may be used under the terms of the GNU
31** General Public License version 3.0 as published by the Free Software
32** Foundation and appearing in the file LICENSE.GPL included in the
33** packaging of this file. Please review the following information to
34** ensure the GNU General Public License version 3.0 requirements will be
35** met: http://www.gnu.org/copyleft/gpl.html.
36**
37**
38** $QT_END_LICENSE$
39**
40****************************************************************************/
41
42#include "qlocale_tools_p.h"
43#include "qlocale_p.h"
44#include "qstring.h"
45
46#include <ctype.h>
47#include <float.h>
48#include <limits.h>
49#include <math.h>
50#include <stdlib.h>
51#include <time.h>
52
53#ifdef Q_OS_WINCE
54# include "qfunctions_wince.h" // for _control87
55#endif
56
57#if defined(Q_OS_LINUX) && !defined(__UCLIBC__)
58# include <fenv.h>
59#endif
60
61// Sizes as defined by the ISO C99 standard - fallback
62#ifndef LLONG_MAX
63# define LLONG_MAX Q_INT64_C(0x7fffffffffffffff)
64#endif
65#ifndef LLONG_MIN
66# define LLONG_MIN (-LLONG_MAX - Q_INT64_C(1))
67#endif
68#ifndef ULLONG_MAX
69# define ULLONG_MAX Q_UINT64_C(0xffffffffffffffff)
70#endif
71
72QT_BEGIN_NAMESPACE
73
74#ifndef QT_QLOCALE_USES_FCVT
75static char *_qdtoa( NEEDS_VOLATILE double d, int mode, int ndigits, int *decpt,
76 int *sign, char **rve, char **digits_str);
77#endif
78
79QString qulltoa(qulonglong l, int base, const QChar _zero)
80{
81 ushort buff[65]; // length of MAX_ULLONG in base 2
82 ushort *p = buff + 65;
83
84 if (base != 10 || _zero.unicode() == '0') {
85 while (l != 0) {
86 int c = l % base;
87
88 --p;
89
90 if (c < 10)
91 *p = '0' + c;
92 else
93 *p = c - 10 + 'a';
94
95 l /= base;
96 }
97 }
98 else {
99 while (l != 0) {
100 int c = l % base;
101
102 *(--p) = _zero.unicode() + c;
103
104 l /= base;
105 }
106 }
107
108 return QString(reinterpret_cast<QChar *>(p), 65 - (p - buff));
109}
110
111QString qlltoa(qlonglong l, int base, const QChar zero)
112{
113 return qulltoa(l < 0 ? -l : l, base, zero);
114}
115
116QString &decimalForm(QChar zero, QChar decimal, QChar group,
117 QString &digits, int decpt, uint precision,
118 PrecisionMode pm,
119 bool always_show_decpt,
120 bool thousands_group)
121{
122 if (decpt < 0) {
123 for (int i = 0; i < -decpt; ++i)
124 digits.prepend(zero);
125 decpt = 0;
126 }
127 else if (decpt > digits.length()) {
128 for (int i = digits.length(); i < decpt; ++i)
129 digits.append(zero);
130 }
131
132 if (pm == PMDecimalDigits) {
133 uint decimal_digits = digits.length() - decpt;
134 for (uint i = decimal_digits; i < precision; ++i)
135 digits.append(zero);
136 }
137 else if (pm == PMSignificantDigits) {
138 for (uint i = digits.length(); i < precision; ++i)
139 digits.append(zero);
140 }
141 else { // pm == PMChopTrailingZeros
142 }
143
144 if (always_show_decpt || decpt < digits.length())
145 digits.insert(decpt, decimal);
146
147 if (thousands_group) {
148 for (int i = decpt - 3; i > 0; i -= 3)
149 digits.insert(i, group);
150 }
151
152 if (decpt == 0)
153 digits.prepend(zero);
154
155 return digits;
156}
157
158QString &exponentForm(QChar zero, QChar decimal, QChar exponential,
159 QChar group, QChar plus, QChar minus,
160 QString &digits, int decpt, uint precision,
161 PrecisionMode pm,
162 bool always_show_decpt)
163{
164 int exp = decpt - 1;
165
166 if (pm == PMDecimalDigits) {
167 for (uint i = digits.length(); i < precision + 1; ++i)
168 digits.append(zero);
169 }
170 else if (pm == PMSignificantDigits) {
171 for (uint i = digits.length(); i < precision; ++i)
172 digits.append(zero);
173 }
174 else { // pm == PMChopTrailingZeros
175 }
176
177 if (always_show_decpt || digits.length() > 1)
178 digits.insert(1, decimal);
179
180 digits.append(exponential);
181 digits.append(QLocalePrivate::longLongToString(zero, group, plus, minus,
182 exp, 2, 10, -1, QLocalePrivate::AlwaysShowSign));
183
184 return digits;
185}
186
187// Removes thousand-group separators in "C" locale.
188bool removeGroupSeparators(QLocalePrivate::CharBuff *num)
189{
190 int group_cnt = 0; // counts number of group chars
191 int decpt_idx = -1;
192
193 char *data = num->data();
194 int l = qstrlen(data);
195
196 // Find the decimal point and check if there are any group chars
197 int i = 0;
198 for (; i < l; ++i) {
199 char c = data[i];
200
201 if (c == ',') {
202 if (i == 0 || data[i - 1] < '0' || data[i - 1] > '9')
203 return false;
204 if (i == l - 1 || data[i + 1] < '0' || data[i + 1] > '9')
205 return false;
206 ++group_cnt;
207 }
208 else if (c == '.') {
209 // Fail if more than one decimal points
210 if (decpt_idx != -1)
211 return false;
212 decpt_idx = i;
213 } else if (c == 'e' || c == 'E') {
214 // an 'e' or 'E' - if we have not encountered a decimal
215 // point, this is where it "is".
216 if (decpt_idx == -1)
217 decpt_idx = i;
218 }
219 }
220
221 // If no group chars, we're done
222 if (group_cnt == 0)
223 return true;
224
225 // No decimal point means that it "is" at the end of the string
226 if (decpt_idx == -1)
227 decpt_idx = l;
228
229 i = 0;
230 while (i < l && group_cnt > 0) {
231 char c = data[i];
232
233 if (c == ',') {
234 // Don't allow group chars after the decimal point
235 if (i > decpt_idx)
236 return false;
237
238 // Check that it is placed correctly relative to the decpt
239 if ((decpt_idx - i) % 4 != 0)
240 return false;
241
242 // Remove it
243 memmove(data + i, data + i + 1, l - i - 1);
244 data[--l] = '\0';
245
246 --group_cnt;
247 --decpt_idx;
248 } else {
249 // Check that we are not missing a separator
250 if (i < decpt_idx
251 && (decpt_idx - i) % 4 == 0
252 && !(i == 0 && c == '-')) // check for negative sign at start of string
253 return false;
254 ++i;
255 }
256 }
257
258 return true;
259}
260
261#if defined(Q_CC_MWERKS) && defined(Q_OS_WIN32)
262inline bool isascii(int c)
263{
264 return (c >= 0 && c <=127);
265}
266#endif
267
268/*-
269 * Copyright (c) 1992, 1993
270 * The Regents of the University of California. All rights reserved.
271 *
272 * Redistribution and use in source and binary forms, with or without
273 * modification, are permitted provided that the following conditions
274 * are met:
275 * 1. Redistributions of source code must retain the above copyright
276 * notice, this list of conditions and the following disclaimer.
277 * 2. Redistributions in binary form must reproduce the above copyright
278 * notice, this list of conditions and the following disclaimer in the
279 * documentation and/or other materials provided with the distribution.
280 * 3. All advertising materials mentioning features or use of this software
281 * must display the following acknowledgment:
282 * This product includes software developed by the University of
283 * California, Berkeley and its contributors.
284 * 4. Neither the name of the University nor the names of its contributors
285 * may be used to endorse or promote products derived from this software
286 * without specific prior written permission.
287 *
288 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
289 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
290 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
291 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
292 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
293 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
294 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
295 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
296 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
297 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
298 * SUCH DAMAGE.
299 */
300
301// static char sccsid[] = "@(#)strtouq.c 8.1 (Berkeley) 6/4/93";
302// "$FreeBSD: src/lib/libc/stdlib/strtoull.c,v 1.5.2.1 2001/03/02 09:45:20 obrien Exp $";
303
304/*
305 * Convert a string to an unsigned long long integer.
306 *
307 * Ignores `locale' stuff. Assumes that the upper and lower case
308 * alphabets and digits are each contiguous.
309 */
310qulonglong qstrtoull(const char *nptr, const char **endptr, register int base, bool *ok)
311{
312 register const char *s = nptr;
313 register qulonglong acc;
314 register unsigned char c;
315 register qulonglong qbase, cutoff;
316 register int any, cutlim;
317
318 if (ok != 0)
319 *ok = true;
320
321 /*
322 * See strtoq for comments as to the logic used.
323 */
324 s = nptr;
325 do {
326 c = *s++;
327 } while (isspace(c));
328 if (c == '-') {
329 if (ok != 0)
330 *ok = false;
331 if (endptr != 0)
332 *endptr = s - 1;
333 return 0;
334 } else {
335 if (c == '+')
336 c = *s++;
337 }
338 if ((base == 0 || base == 16) &&
339 c == '0' && (*s == 'x' || *s == 'X')) {
340 c = s[1];
341 s += 2;
342 base = 16;
343 }
344 if (base == 0)
345 base = c == '0' ? 8 : 10;
346 qbase = unsigned(base);
347 cutoff = qulonglong(ULLONG_MAX) / qbase;
348 cutlim = qulonglong(ULLONG_MAX) % qbase;
349 for (acc = 0, any = 0;; c = *s++) {
350 if (!isascii(c))
351 break;
352 if (isdigit(c))
353 c -= '0';
354 else if (isalpha(c))
355 c -= isupper(c) ? 'A' - 10 : 'a' - 10;
356 else
357 break;
358 if (c >= base)
359 break;
360 if (any < 0 || acc > cutoff || (acc == cutoff && c > cutlim))
361 any = -1;
362 else {
363 any = 1;
364 acc *= qbase;
365 acc += c;
366 }
367 }
368 if (any == 0) {
369 if (ok != 0)
370 *ok = false;
371 } else if (any < 0) {
372 acc = ULLONG_MAX;
373 if (ok != 0)
374 *ok = false;
375 }
376 if (endptr != 0)
377 *endptr = (any ? s - 1 : nptr);
378 return acc;
379}
380
381
382// "$FreeBSD: src/lib/libc/stdlib/strtoll.c,v 1.5.2.1 2001/03/02 09:45:20 obrien Exp $";
383
384
385/*
386 * Convert a string to a long long integer.
387 *
388 * Ignores `locale' stuff. Assumes that the upper and lower case
389 * alphabets and digits are each contiguous.
390 */
391qlonglong qstrtoll(const char *nptr, const char **endptr, register int base, bool *ok)
392{
393 register const char *s;
394 register qulonglong acc;
395 register unsigned char c;
396 register qulonglong qbase, cutoff;
397 register int neg, any, cutlim;
398
399 /*
400 * Skip white space and pick up leading +/- sign if any.
401 * If base is 0, allow 0x for hex and 0 for octal, else
402 * assume decimal; if base is already 16, allow 0x.
403 */
404 s = nptr;
405 do {
406 c = *s++;
407 } while (isspace(c));
408 if (c == '-') {
409 neg = 1;
410 c = *s++;
411 } else {
412 neg = 0;
413 if (c == '+')
414 c = *s++;
415 }
416 if ((base == 0 || base == 16) &&
417 c == '0' && (*s == 'x' || *s == 'X')) {
418 c = s[1];
419 s += 2;
420 base = 16;
421 }
422 if (base == 0)
423 base = c == '0' ? 8 : 10;
424
425 /*
426 * Compute the cutoff value between legal numbers and illegal
427 * numbers. That is the largest legal value, divided by the
428 * base. An input number that is greater than this value, if
429 * followed by a legal input character, is too big. One that
430 * is equal to this value may be valid or not; the limit
431 * between valid and invalid numbers is then based on the last
432 * digit. For instance, if the range for quads is
433 * [-9223372036854775808..9223372036854775807] and the input base
434 * is 10, cutoff will be set to 922337203685477580 and cutlim to
435 * either 7 (neg==0) or 8 (neg==1), meaning that if we have
436 * accumulated a value > 922337203685477580, or equal but the
437 * next digit is > 7 (or 8), the number is too big, and we will
438 * return a range error.
439 *
440 * Set any if any `digits' consumed; make it negative to indicate
441 * overflow.
442 */
443 qbase = unsigned(base);
444 cutoff = neg ? qulonglong(0-(LLONG_MIN + LLONG_MAX)) + LLONG_MAX : LLONG_MAX;
445 cutlim = cutoff % qbase;
446 cutoff /= qbase;
447 for (acc = 0, any = 0;; c = *s++) {
448 if (!isascii(c))
449 break;
450 if (isdigit(c))
451 c -= '0';
452 else if (isalpha(c))
453 c -= isupper(c) ? 'A' - 10 : 'a' - 10;
454 else
455 break;
456 if (c >= base)
457 break;
458 if (any < 0 || acc > cutoff || (acc == cutoff && c > cutlim))
459 any = -1;
460 else {
461 any = 1;
462 acc *= qbase;
463 acc += c;
464 }
465 }
466 if (any < 0) {
467 acc = neg ? LLONG_MIN : LLONG_MAX;
468 if (ok != 0)
469 *ok = false;
470 } else if (neg) {
471 acc = (~acc) + 1;
472 }
473 if (endptr != 0)
474 *endptr = (any >= 0 ? s - 1 : nptr);
475
476 if (ok != 0)
477 *ok = any > 0;
478
479 return acc;
480}
481
482#ifndef QT_QLOCALE_USES_FCVT
483
484/* From: NetBSD: strtod.c,v 1.26 1998/02/03 18:44:21 perry Exp */
485/* $FreeBSD: src/lib/libc/stdlib/netbsd_strtod.c,v 1.2.2.2 2001/03/02 17:14:15 tegge Exp $ */
486
487/* Please send bug reports to
488 David M. Gay
489 AT&T Bell Laboratories, Room 2C-463
490 600 Mountain Avenue
491 Murray Hill, NJ 07974-2070
492 U.S.A.
493 dmg@research.att.com or research!dmg
494 */
495
496/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
497 *
498 * This strtod returns a nearest machine number to the input decimal
499 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
500 * broken by the IEEE round-even rule. Otherwise ties are broken by
501 * biased rounding (add half and chop).
502 *
503 * Inspired loosely by William D. Clinger's paper "How to Read Floating
504 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
505 *
506 * Modifications:
507 *
508 * 1. We only require IEEE, IBM, or VAX double-precision
509 * arithmetic (not IEEE double-extended).
510 * 2. We get by with floating-point arithmetic in a case that
511 * Clinger missed -- when we're computing d * 10^n
512 * for a small integer d and the integer n is not too
513 * much larger than 22 (the maximum integer k for which
514 * we can represent 10^k exactly), we may be able to
515 * compute (d*10^k) * 10^(e-k) with just one roundoff.
516 * 3. Rather than a bit-at-a-time adjustment of the binary
517 * result in the hard case, we use floating-point
518 * arithmetic to determine the adjustment to within
519 * one bit; only in really hard cases do we need to
520 * compute a second residual.
521 * 4. Because of 3., we don't need a large table of powers of 10
522 * for ten-to-e (just some small tables, e.g. of 10^k
523 * for 0 <= k <= 22).
524 */
525
526/*
527 * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
528 * significant byte has the lowest address.
529 * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
530 * significant byte has the lowest address.
531 * #define Long int on machines with 32-bit ints and 64-bit longs.
532 * #define Sudden_Underflow for IEEE-format machines without gradual
533 * underflow (i.e., that flush to zero on underflow).
534 * #define IBM for IBM mainframe-style floating-point arithmetic.
535 * #define VAX for VAX-style floating-point arithmetic.
536 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
537 * #define No_leftright to omit left-right logic in fast floating-point
538 * computation of dtoa.
539 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
540 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
541 * that use extended-precision instructions to compute rounded
542 * products and quotients) with IBM.
543 * #define ROUND_BIASED for IEEE-format with biased rounding.
544 * #define Inaccurate_Divide for IEEE-format with correctly rounded
545 * products but inaccurate quotients, e.g., for Intel i860.
546 * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision
547 * integer arithmetic. Whether this speeds things up or slows things
548 * down depends on the machine and the number being converted.
549 * #define KR_headers for old-style C function headers.
550 * #define Bad_float_h if your system lacks a float.h or if it does not
551 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
552 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
553 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
554 * if memory is available and otherwise does something you deem
555 * appropriate. If MALLOC is undefined, malloc will be invoked
556 * directly -- and assumed always to succeed.
557 */
558
559#if defined(LIBC_SCCS) && !defined(lint)
560__RCSID("$NetBSD: strtod.c,v 1.26 1998/02/03 18:44:21 perry Exp $");
561#endif /* LIBC_SCCS and not lint */
562
563/*
564#if defined(__m68k__) || defined(__sparc__) || defined(__i386__) || \
565 defined(__mips__) || defined(__ns32k__) || defined(__alpha__) || \
566 defined(__powerpc__) || defined(Q_OS_WIN) || defined(Q_OS_DARWIN) || defined(Q_OS_MAC) || \
567 defined(mips) || defined(Q_OS_AIX) || defined(Q_OS_SOLARIS)
568# define IEEE_BIG_OR_LITTLE_ENDIAN 1
569#endif
570*/
571
572// *All* of our architectures have IEEE arithmetic, don't they?
573#define IEEE_BIG_OR_LITTLE_ENDIAN 1
574
575#ifdef __arm32__
576/*
577 * Although the CPU is little endian the FP has different
578 * byte and word endianness. The byte order is still little endian
579 * but the word order is big endian.
580 */
581#define IEEE_BIG_OR_LITTLE_ENDIAN
582#endif
583
584#ifdef vax
585#define VAX
586#endif
587
588#define Long qint32
589#define ULong quint32
590
591#define MALLOC malloc
592
593#ifdef BSD_QDTOA_DEBUG
594QT_BEGIN_INCLUDE_NAMESPACE
595#include <stdio.h>
596QT_END_INCLUDE_NAMESPACE
597
598#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
599#endif
600
601#ifdef Unsigned_Shifts
602#define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
603#else
604#define Sign_Extend(a,b) /*no-op*/
605#endif
606
607#if (defined(IEEE_BIG_OR_LITTLE_ENDIAN) + defined(VAX) + defined(IBM)) != 1
608#error Exactly one of IEEE_BIG_OR_LITTLE_ENDIAN, VAX, or IBM should be defined.
609#endif
610
611static inline ULong _getWord0(const NEEDS_VOLATILE double x)
612{
613 const NEEDS_VOLATILE uchar *ptr = reinterpret_cast<const NEEDS_VOLATILE uchar *>(&x);
614 if (QSysInfo::ByteOrder == QSysInfo::BigEndian) {
615 return (ptr[0]<<24) + (ptr[1]<<16) + (ptr[2]<<8) + ptr[3];
616 } else {
617 return (ptr[7]<<24) + (ptr[6]<<16) + (ptr[5]<<8) + ptr[4];
618 }
619}
620
621static inline void _setWord0(NEEDS_VOLATILE double *x, ULong l)
622{
623 NEEDS_VOLATILE uchar *ptr = reinterpret_cast<NEEDS_VOLATILE uchar *>(x);
624 if (QSysInfo::ByteOrder == QSysInfo::BigEndian) {
625 ptr[0] = uchar(l>>24);
626 ptr[1] = uchar(l>>16);
627 ptr[2] = uchar(l>>8);
628 ptr[3] = uchar(l);
629 } else {
630 ptr[7] = uchar(l>>24);
631 ptr[6] = uchar(l>>16);
632 ptr[5] = uchar(l>>8);
633 ptr[4] = uchar(l);
634 }
635}
636
637static inline ULong _getWord1(const NEEDS_VOLATILE double x)
638{
639 const NEEDS_VOLATILE uchar *ptr = reinterpret_cast<const NEEDS_VOLATILE uchar *>(&x);
640 if (QSysInfo::ByteOrder == QSysInfo::BigEndian) {
641 return (ptr[4]<<24) + (ptr[5]<<16) + (ptr[6]<<8) + ptr[7];
642 } else {
643 return (ptr[3]<<24) + (ptr[2]<<16) + (ptr[1]<<8) + ptr[0];
644 }
645}
646static inline void _setWord1(NEEDS_VOLATILE double *x, ULong l)
647{
648 NEEDS_VOLATILE uchar *ptr = reinterpret_cast<uchar NEEDS_VOLATILE *>(x);
649 if (QSysInfo::ByteOrder == QSysInfo::BigEndian) {
650 ptr[4] = uchar(l>>24);
651 ptr[5] = uchar(l>>16);
652 ptr[6] = uchar(l>>8);
653 ptr[7] = uchar(l);
654 } else {
655 ptr[3] = uchar(l>>24);
656 ptr[2] = uchar(l>>16);
657 ptr[1] = uchar(l>>8);
658 ptr[0] = uchar(l);
659 }
660}
661
662static inline ULong getWord0(const NEEDS_VOLATILE double x)
663{
664#ifdef QT_ARMFPA
665 return _getWord1(x);
666#else
667 return _getWord0(x);
668#endif
669}
670
671static inline void setWord0(NEEDS_VOLATILE double *x, ULong l)
672{
673#ifdef QT_ARMFPA
674 _setWord1(x, l);
675#else
676 _setWord0(x, l);
677#endif
678}
679
680static inline ULong getWord1(const NEEDS_VOLATILE double x)
681{
682#ifdef QT_ARMFPA
683 return _getWord0(x);
684#else
685 return _getWord1(x);
686#endif
687}
688
689static inline void setWord1(NEEDS_VOLATILE double *x, ULong l)
690{
691#ifdef QT_ARMFPA
692 _setWord0(x, l);
693#else
694 _setWord1(x, l);
695#endif
696}
697
698static inline void Storeinc(ULong *&a, const ULong &b, const ULong &c)
699{
700
701 *a = (ushort(b) << 16) | ushort(c);
702 ++a;
703}
704
705/* #define P DBL_MANT_DIG */
706/* Ten_pmax = floor(P*log(2)/log(5)) */
707/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
708/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
709/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
710
711#if defined(IEEE_BIG_OR_LITTLE_ENDIAN)
712#define Exp_shift 20
713#define Exp_shift1 20
714#define Exp_msk1 0x100000
715#define Exp_msk11 0x100000
716#define Exp_mask 0x7ff00000
717#define P 53
718#define Bias 1023
719#define IEEE_Arith
720#define Emin (-1022)
721#define Exp_1 0x3ff00000
722#define Exp_11 0x3ff00000
723#define Ebits 11
724#define Frac_mask 0xfffff
725#define Frac_mask1 0xfffff
726#define Ten_pmax 22
727#define Bletch 0x10
728#define Bndry_mask 0xfffff
729#define Bndry_mask1 0xfffff
730#if defined(LSB) && defined(Q_OS_VXWORKS)
731#undef LSB
732#endif
733#define LSB 1
734#define Sign_bit 0x80000000
735#define Log2P 1
736#define Tiny0 0
737#define Tiny1 1
738#define Quick_max 14
739#define Int_max 14
740#define Infinite(x) (getWord0(x) == 0x7ff00000) /* sufficient test for here */
741#else
742#undef Sudden_Underflow
743#define Sudden_Underflow
744#ifdef IBM
745#define Exp_shift 24
746#define Exp_shift1 24
747#define Exp_msk1 0x1000000
748#define Exp_msk11 0x1000000
749#define Exp_mask 0x7f000000
750#define P 14
751#define Bias 65
752#define Exp_1 0x41000000
753#define Exp_11 0x41000000
754#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
755#define Frac_mask 0xffffff
756#define Frac_mask1 0xffffff
757#define Bletch 4
758#define Ten_pmax 22
759#define Bndry_mask 0xefffff
760#define Bndry_mask1 0xffffff
761#define LSB 1
762#define Sign_bit 0x80000000
763#define Log2P 4
764#define Tiny0 0x100000
765#define Tiny1 0
766#define Quick_max 14
767#define Int_max 15
768#else /* VAX */
769#define Exp_shift 23
770#define Exp_shift1 7
771#define Exp_msk1 0x80
772#define Exp_msk11 0x800000
773#define Exp_mask 0x7f80
774#define P 56
775#define Bias 129
776#define Exp_1 0x40800000
777#define Exp_11 0x4080
778#define Ebits 8
779#define Frac_mask 0x7fffff
780#define Frac_mask1 0xffff007f
781#define Ten_pmax 24
782#define Bletch 2
783#define Bndry_mask 0xffff007f
784#define Bndry_mask1 0xffff007f
785#define LSB 0x10000
786#define Sign_bit 0x8000
787#define Log2P 1
788#define Tiny0 0x80
789#define Tiny1 0
790#define Quick_max 15
791#define Int_max 15
792#endif
793#endif
794
795#ifndef IEEE_Arith
796#define ROUND_BIASED
797#endif
798
799#ifdef RND_PRODQUOT
800#define rounded_product(a,b) a = rnd_prod(a, b)
801#define rounded_quotient(a,b) a = rnd_quot(a, b)
802extern double rnd_prod(double, double), rnd_quot(double, double);
803#else
804#define rounded_product(a,b) a *= b
805#define rounded_quotient(a,b) a /= b
806#endif
807
808#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
809#define Big1 0xffffffff
810
811#ifndef Just_16
812/* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
813 * This makes some inner loops simpler and sometimes saves work
814 * during multiplications, but it often seems to make things slightly
815 * slower. Hence the default is now to store 32 bits per Long.
816 */
817#ifndef Pack_32
818#define Pack_32
819#endif
820#endif
821
822#define Kmax 15
823
824struct
825Bigint {
826 struct Bigint *next;
827 int k, maxwds, sign, wds;
828 ULong x[1];
829};
830
831 typedef struct Bigint Bigint;
832
833static Bigint *Balloc(int k)
834{
835 int x;
836 Bigint *rv;
837
838 x = 1 << k;
839 rv = static_cast<Bigint *>(MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long)));
840 Q_CHECK_PTR(rv);
841 rv->k = k;
842 rv->maxwds = x;
843 rv->sign = rv->wds = 0;
844 return rv;
845}
846
847static void Bfree(Bigint *v)
848{
849 free(v);
850}
851
852#define Bcopy(x,y) memcpy(reinterpret_cast<char *>(&x->sign), reinterpret_cast<char *>(&y->sign), \
853y->wds*sizeof(Long) + 2*sizeof(int))
854
855/* multiply by m and add a */
856static Bigint *multadd(Bigint *b, int m, int a)
857{
858 int i, wds;
859 ULong *x, y;
860#ifdef Pack_32
861 ULong xi, z;
862#endif
863 Bigint *b1;
864
865 wds = b->wds;
866 x = b->x;
867 i = 0;
868 do {
869#ifdef Pack_32
870 xi = *x;
871 y = (xi & 0xffff) * m + a;
872 z = (xi >> 16) * m + (y >> 16);
873 a = (z >> 16);
874 *x++ = (z << 16) + (y & 0xffff);
875#else
876 y = *x * m + a;
877 a = (y >> 16);
878 *x++ = y & 0xffff;
879#endif
880 }
881 while(++i < wds);
882 if (a) {
883 if (wds >= b->maxwds) {
884 b1 = Balloc(b->k+1);
885 Bcopy(b1, b);
886 Bfree(b);
887 b = b1;
888 }
889 b->x[wds++] = a;
890 b->wds = wds;
891 }
892 return b;
893}
894
895static Bigint *s2b(const char *s, int nd0, int nd, ULong y9)
896{
897 Bigint *b;
898 int i, k;
899 Long x, y;
900
901 x = (nd + 8) / 9;
902 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
903#ifdef Pack_32
904 b = Balloc(k);
905 b->x[0] = y9;
906 b->wds = 1;
907#else
908 b = Balloc(k+1);
909 b->x[0] = y9 & 0xffff;
910 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
911#endif
912
913 i = 9;
914 if (9 < nd0) {
915 s += 9;
916 do b = multadd(b, 10, *s++ - '0');
917 while(++i < nd0);
918 s++;
919 }
920 else
921 s += 10;
922 for(; i < nd; i++)
923 b = multadd(b, 10, *s++ - '0');
924 return b;
925}
926
927static int hi0bits(ULong x)
928{
929 int k = 0;
930
931 if (!(x & 0xffff0000)) {
932 k = 16;
933 x <<= 16;
934 }
935 if (!(x & 0xff000000)) {
936 k += 8;
937 x <<= 8;
938 }
939 if (!(x & 0xf0000000)) {
940 k += 4;
941 x <<= 4;
942 }
943 if (!(x & 0xc0000000)) {
944 k += 2;
945 x <<= 2;
946 }
947 if (!(x & 0x80000000)) {
948 k++;
949 if (!(x & 0x40000000))
950 return 32;
951 }
952 return k;
953}
954
955static int lo0bits(ULong *y)
956{
957 int k;
958 ULong x = *y;
959
960 if (x & 7) {
961 if (x & 1)
962 return 0;
963 if (x & 2) {
964 *y = x >> 1;
965 return 1;
966 }
967 *y = x >> 2;
968 return 2;
969 }
970 k = 0;
971 if (!(x & 0xffff)) {
972 k = 16;
973 x >>= 16;
974 }
975 if (!(x & 0xff)) {
976 k += 8;
977 x >>= 8;
978 }
979 if (!(x & 0xf)) {
980 k += 4;
981 x >>= 4;
982 }
983 if (!(x & 0x3)) {
984 k += 2;
985 x >>= 2;
986 }
987 if (!(x & 1)) {
988 k++;
989 x >>= 1;
990 if (!x & 1)
991 return 32;
992 }
993 *y = x;
994 return k;
995}
996
997static Bigint *i2b(int i)
998{
999 Bigint *b;
1000
1001 b = Balloc(1);
1002 b->x[0] = i;
1003 b->wds = 1;
1004 return b;
1005}
1006
1007static Bigint *mult(Bigint *a, Bigint *b)
1008{
1009 Bigint *c;
1010 int k, wa, wb, wc;
1011 ULong carry, y, z;
1012 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
1013#ifdef Pack_32
1014 ULong z2;
1015#endif
1016
1017 if (a->wds < b->wds) {
1018 c = a;
1019 a = b;
1020 b = c;
1021 }
1022 k = a->k;
1023 wa = a->wds;
1024 wb = b->wds;
1025 wc = wa + wb;
1026 if (wc > a->maxwds)
1027 k++;
1028 c = Balloc(k);
1029 for(x = c->x, xa = x + wc; x < xa; x++)
1030 *x = 0;
1031 xa = a->x;
1032 xae = xa + wa;
1033 xb = b->x;
1034 xbe = xb + wb;
1035 xc0 = c->x;
1036#ifdef Pack_32
1037 for(; xb < xbe; xb++, xc0++) {
1038 if ((y = *xb & 0xffff) != 0) {
1039 x = xa;
1040 xc = xc0;
1041 carry = 0;
1042 do {
1043 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
1044 carry = z >> 16;
1045 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
1046 carry = z2 >> 16;
1047 Storeinc(xc, z2, z);
1048 }
1049 while(x < xae);
1050 *xc = carry;
1051 }
1052 if ((y = *xb >> 16) != 0) {
1053 x = xa;
1054 xc = xc0;
1055 carry = 0;
1056 z2 = *xc;
1057 do {
1058 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
1059 carry = z >> 16;
1060 Storeinc(xc, z, z2);
1061 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
1062 carry = z2 >> 16;
1063 }
1064 while(x < xae);
1065 *xc = z2;
1066 }
1067 }
1068#else
1069 for(; xb < xbe; xc0++) {
1070 if (y = *xb++) {
1071 x = xa;
1072 xc = xc0;
1073 carry = 0;
1074 do {
1075 z = *x++ * y + *xc + carry;
1076 carry = z >> 16;
1077 *xc++ = z & 0xffff;
1078 }
1079 while(x < xae);
1080 *xc = carry;
1081 }
1082 }
1083#endif
1084 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
1085 c->wds = wc;
1086 return c;
1087}
1088
1089static Bigint *p5s;
1090
1091struct p5s_deleter
1092{
1093 ~p5s_deleter()
1094 {
1095 while (p5s) {
1096 Bigint *next = p5s->next;
1097 Bfree(p5s);
1098 p5s = next;
1099 }
1100 }
1101};
1102
1103static Bigint *pow5mult(Bigint *b, int k)
1104{
1105 Bigint *b1, *p5, *p51;
1106 int i;
1107 static const int p05[3] = { 5, 25, 125 };
1108
1109 if ((i = k & 3) != 0)
1110#if defined(Q_OS_IRIX) && defined(Q_CC_GNU)
1111 {
1112 // work around a bug on 64 bit IRIX gcc
1113 int *p = (int *) p05;
1114 b = multadd(b, p[i-1], 0);
1115 }
1116#else
1117 b = multadd(b, p05[i-1], 0);
1118#endif
1119
1120 if (!(k >>= 2))
1121 return b;
1122 if (!(p5 = p5s)) {
1123 /* first time */
1124 static p5s_deleter deleter;
1125 p5 = p5s = i2b(625);
1126 p5->next = 0;
1127 }
1128 for(;;) {
1129 if (k & 1) {
1130 b1 = mult(b, p5);
1131 Bfree(b);
1132 b = b1;
1133 }
1134 if (!(k >>= 1))
1135 break;
1136 if (!(p51 = p5->next)) {
1137 p51 = p5->next = mult(p5,p5);
1138 p51->next = 0;
1139 }
1140 p5 = p51;
1141 }
1142 return b;
1143}
1144
1145static Bigint *lshift(Bigint *b, int k)
1146{
1147 int i, k1, n, n1;
1148 Bigint *b1;
1149 ULong *x, *x1, *xe, z;
1150
1151#ifdef Pack_32
1152 n = k >> 5;
1153#else
1154 n = k >> 4;
1155#endif
1156 k1 = b->k;
1157 n1 = n + b->wds + 1;
1158 for(i = b->maxwds; n1 > i; i <<= 1)
1159 k1++;
1160 b1 = Balloc(k1);
1161 x1 = b1->x;
1162 for(i = 0; i < n; i++)
1163 *x1++ = 0;
1164 x = b->x;
1165 xe = x + b->wds;
1166#ifdef Pack_32
1167 if (k &= 0x1f) {
1168 k1 = 32 - k;
1169 z = 0;
1170 do {
1171 *x1++ = *x << k | z;
1172 z = *x++ >> k1;
1173 }
1174 while(x < xe);
1175 if ((*x1 = z) != 0)
1176 ++n1;
1177 }
1178#else
1179 if (k &= 0xf) {
1180 k1 = 16 - k;
1181 z = 0;
1182 do {
1183 *x1++ = *x << k & 0xffff | z;
1184 z = *x++ >> k1;
1185 }
1186 while(x < xe);
1187 if (*x1 = z)
1188 ++n1;
1189 }
1190#endif
1191 else do
1192 *x1++ = *x++;
1193 while(x < xe);
1194 b1->wds = n1 - 1;
1195 Bfree(b);
1196 return b1;
1197}
1198
1199static int cmp(Bigint *a, Bigint *b)
1200{
1201 ULong *xa, *xa0, *xb, *xb0;
1202 int i, j;
1203
1204 i = a->wds;
1205 j = b->wds;
1206#ifdef BSD_QDTOA_DEBUG
1207 if (i > 1 && !a->x[i-1])
1208 Bug("cmp called with a->x[a->wds-1] == 0");
1209 if (j > 1 && !b->x[j-1])
1210 Bug("cmp called with b->x[b->wds-1] == 0");
1211#endif
1212 if (i -= j)
1213 return i;
1214 xa0 = a->x;
1215 xa = xa0 + j;
1216 xb0 = b->x;
1217 xb = xb0 + j;
1218 for(;;) {
1219 if (*--xa != *--xb)
1220 return *xa < *xb ? -1 : 1;
1221 if (xa <= xa0)
1222 break;
1223 }
1224 return 0;
1225}
1226
1227static Bigint *diff(Bigint *a, Bigint *b)
1228{
1229 Bigint *c;
1230 int i, wa, wb;
1231 Long borrow, y; /* We need signed shifts here. */
1232 ULong *xa, *xae, *xb, *xbe, *xc;
1233#ifdef Pack_32
1234 Long z;
1235#endif
1236
1237 i = cmp(a,b);
1238 if (!i) {
1239 c = Balloc(0);
1240 c->wds = 1;
1241 c->x[0] = 0;
1242 return c;
1243 }
1244 if (i < 0) {
1245 c = a;
1246 a = b;
1247 b = c;
1248 i = 1;
1249 }
1250 else
1251 i = 0;
1252 c = Balloc(a->k);
1253 c->sign = i;
1254 wa = a->wds;
1255 xa = a->x;
1256 xae = xa + wa;
1257 wb = b->wds;
1258 xb = b->x;
1259 xbe = xb + wb;
1260 xc = c->x;
1261 borrow = 0;
1262#ifdef Pack_32
1263 do {
1264 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
1265 borrow = y >> 16;
1266 Sign_Extend(borrow, y);
1267 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
1268 borrow = z >> 16;
1269 Sign_Extend(borrow, z);
1270 Storeinc(xc, z, y);
1271 }
1272 while(xb < xbe);
1273 while(xa < xae) {
1274 y = (*xa & 0xffff) + borrow;
1275 borrow = y >> 16;
1276 Sign_Extend(borrow, y);
1277 z = (*xa++ >> 16) + borrow;
1278 borrow = z >> 16;
1279 Sign_Extend(borrow, z);
1280 Storeinc(xc, z, y);
1281 }
1282#else
1283 do {
1284 y = *xa++ - *xb++ + borrow;
1285 borrow = y >> 16;
1286 Sign_Extend(borrow, y);
1287 *xc++ = y & 0xffff;
1288 }
1289 while(xb < xbe);
1290 while(xa < xae) {
1291 y = *xa++ + borrow;
1292 borrow = y >> 16;
1293 Sign_Extend(borrow, y);
1294 *xc++ = y & 0xffff;
1295 }
1296#endif
1297 while(!*--xc)
1298 wa--;
1299 c->wds = wa;
1300 return c;
1301}
1302
1303static double ulp(double x)
1304{
1305 Long L;
1306 double a;
1307
1308 L = (getWord0(x) & Exp_mask) - (P-1)*Exp_msk1;
1309#ifndef Sudden_Underflow
1310 if (L > 0) {
1311#endif
1312#ifdef IBM
1313 L |= Exp_msk1 >> 4;
1314#endif
1315 setWord0(&a, L);
1316 setWord1(&a, 0);
1317#ifndef Sudden_Underflow
1318 }
1319 else {
1320 L = -L >> Exp_shift;
1321 if (L < Exp_shift) {
1322 setWord0(&a, 0x80000 >> L);
1323 setWord1(&a, 0);
1324 }
1325 else {
1326 setWord0(&a, 0);
1327 L -= Exp_shift;
1328 setWord1(&a, L >= 31 ? 1U : 1U << (31 - L));
1329 }
1330 }
1331#endif
1332 return a;
1333}
1334
1335static double b2d(Bigint *a, int *e)
1336{
1337 ULong *xa, *xa0, w, y, z;
1338 int k;
1339 double d;
1340
1341 xa0 = a->x;
1342 xa = xa0 + a->wds;
1343 y = *--xa;
1344#ifdef BSD_QDTOA_DEBUG
1345 if (!y) Bug("zero y in b2d");
1346#endif
1347 k = hi0bits(y);
1348 *e = 32 - k;
1349#ifdef Pack_32
1350 if (k < Ebits) {
1351 setWord0(&d, Exp_1 | y >> (Ebits - k));
1352 w = xa > xa0 ? *--xa : 0;
1353 setWord1(&d, y << ((32-Ebits) + k) | w >> (Ebits - k));
1354 goto ret_d;
1355 }
1356 z = xa > xa0 ? *--xa : 0;
1357 if (k -= Ebits) {
1358 setWord0(&d, Exp_1 | y << k | z >> (32 - k));
1359 y = xa > xa0 ? *--xa : 0;
1360 setWord1(&d, z << k | y >> (32 - k));
1361 }
1362 else {
1363 setWord0(&d, Exp_1 | y);
1364 setWord1(&d, z);
1365 }
1366#else
1367 if (k < Ebits + 16) {
1368 z = xa > xa0 ? *--xa : 0;
1369 setWord0(&d, Exp_1 | y << k - Ebits | z >> Ebits + 16 - k);
1370 w = xa > xa0 ? *--xa : 0;
1371 y = xa > xa0 ? *--xa : 0;
1372 setWord1(&d, z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k);
1373 goto ret_d;
1374 }
1375 z = xa > xa0 ? *--xa : 0;
1376 w = xa > xa0 ? *--xa : 0;
1377 k -= Ebits + 16;
1378 setWord0(&d, Exp_1 | y << k + 16 | z << k | w >> 16 - k);
1379 y = xa > xa0 ? *--xa : 0;
1380 setWord1(&d, w << k + 16 | y << k);
1381#endif
1382 ret_d:
1383 return d;
1384}
1385
1386static Bigint *d2b(double d, int *e, int *bits)
1387{
1388 Bigint *b;
1389 int de, i, k;
1390 ULong *x, y, z;
1391
1392#ifdef Pack_32
1393 b = Balloc(1);
1394#else
1395 b = Balloc(2);
1396#endif
1397 x = b->x;
1398
1399 z = getWord0(d) & Frac_mask;
1400 setWord0(&d, getWord0(d) & 0x7fffffff); /* clear sign bit, which we ignore */
1401#ifdef Sudden_Underflow
1402 de = (int)(getWord0(d) >> Exp_shift);
1403#ifndef IBM
1404 z |= Exp_msk11;
1405#endif
1406#else
1407 if ((de = int(getWord0(d) >> Exp_shift)) != 0)
1408 z |= Exp_msk1;
1409#endif
1410#ifdef Pack_32
1411 if ((y = getWord1(d)) != 0) {
1412 if ((k = lo0bits(&y)) != 0) {
1413 x[0] = y | z << (32 - k);
1414 z >>= k;
1415 }
1416 else
1417 x[0] = y;
1418 i = b->wds = (x[1] = z) ? 2 : 1;
1419 }
1420 else {
1421#ifdef BSD_QDTOA_DEBUG
1422 if (!z)
1423 Bug("Zero passed to d2b");
1424#endif
1425 k = lo0bits(&z);
1426 x[0] = z;
1427 i = b->wds = 1;
1428 k += 32;
1429 }
1430#else
1431 if (y = getWord1(d)) {
1432 if (k = lo0bits(&y))
1433 if (k >= 16) {
1434 x[0] = y | z << 32 - k & 0xffff;
1435 x[1] = z >> k - 16 & 0xffff;
1436 x[2] = z >> k;
1437 i = 2;
1438 }
1439 else {
1440 x[0] = y & 0xffff;
1441 x[1] = y >> 16 | z << 16 - k & 0xffff;
1442 x[2] = z >> k & 0xffff;
1443 x[3] = z >> k+16;
1444 i = 3;
1445 }
1446 else {
1447 x[0] = y & 0xffff;
1448 x[1] = y >> 16;
1449 x[2] = z & 0xffff;
1450 x[3] = z >> 16;
1451 i = 3;
1452 }
1453 }
1454 else {
1455#ifdef BSD_QDTOA_DEBUG
1456 if (!z)
1457 Bug("Zero passed to d2b");
1458#endif
1459 k = lo0bits(&z);
1460 if (k >= 16) {
1461 x[0] = z;
1462 i = 0;
1463 }
1464 else {
1465 x[0] = z & 0xffff;
1466 x[1] = z >> 16;
1467 i = 1;
1468 }
1469 k += 32;
1470 }
1471 while(!x[i])
1472 --i;
1473 b->wds = i + 1;
1474#endif
1475#ifndef Sudden_Underflow
1476 if (de) {
1477#endif
1478#ifdef IBM
1479 *e = (de - Bias - (P-1) << 2) + k;
1480 *bits = 4*P + 8 - k - hi0bits(getWord0(d) & Frac_mask);
1481#else
1482 *e = de - Bias - (P-1) + k;
1483 *bits = P - k;
1484#endif
1485#ifndef Sudden_Underflow
1486 }
1487 else {
1488 *e = de - Bias - (P-1) + 1 + k;
1489#ifdef Pack_32
1490 *bits = 32*i - hi0bits(x[i-1]);
1491#else
1492 *bits = (i+2)*16 - hi0bits(x[i]);
1493#endif
1494 }
1495#endif
1496 return b;
1497}
1498
1499static double ratio(Bigint *a, Bigint *b)
1500{
1501 double da, db;
1502 int k, ka, kb;
1503
1504 da = b2d(a, &ka);
1505 db = b2d(b, &kb);
1506#ifdef Pack_32
1507 k = ka - kb + 32*(a->wds - b->wds);
1508#else
1509 k = ka - kb + 16*(a->wds - b->wds);
1510#endif
1511#ifdef IBM
1512 if (k > 0) {
1513 setWord0(&da, getWord0(da) + (k >> 2)*Exp_msk1);
1514 if (k &= 3)
1515 da *= 1 << k;
1516 }
1517 else {
1518 k = -k;
1519 setWord0(&db, getWord0(db) + (k >> 2)*Exp_msk1);
1520 if (k &= 3)
1521 db *= 1 << k;
1522 }
1523#else
1524 if (k > 0)
1525 setWord0(&da, getWord0(da) + k*Exp_msk1);
1526 else {
1527 k = -k;
1528 setWord0(&db, getWord0(db) + k*Exp_msk1);
1529 }
1530#endif
1531 return da / db;
1532}
1533
1534static const double tens[] = {
1535 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1536 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1537 1e20, 1e21, 1e22
1538#ifdef VAX
1539 , 1e23, 1e24
1540#endif
1541};
1542
1543#ifdef IEEE_Arith
1544static const double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1545static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1546#define n_bigtens 5
1547#else
1548#ifdef IBM
1549static const double bigtens[] = { 1e16, 1e32, 1e64 };
1550static const double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1551#define n_bigtens 3
1552#else
1553static const double bigtens[] = { 1e16, 1e32 };
1554static const double tinytens[] = { 1e-16, 1e-32 };
1555#define n_bigtens 2
1556#endif
1557#endif
1558
1559/*
1560 The pre-release gcc3.3 shipped with SuSE 8.2 has a bug which causes
1561 the comparison 1e-100 == 0.0 to return true. As a workaround, we
1562 compare it to a global variable containing 0.0, which produces
1563 correct assembler output.
1564
1565 ### consider detecting the broken compilers and using the static
1566 ### double for these, and use a #define for all working compilers
1567*/
1568static double g_double_zero = 0.0;
1569
1570Q_CORE_EXPORT double qstrtod(const char *s00, const char **se, bool *ok)
1571{
1572 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1573 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1574 const char *s, *s0, *s1;
1575 double aadj, aadj1, adj, rv, rv0;
1576 Long L;
1577 ULong y, z;
1578 Bigint *bb1, *bd0;
1579 Bigint *bb = NULL, *bd = NULL, *bs = NULL, *delta = NULL;/* pacify gcc */
1580
1581 /*
1582 #ifndef KR_headers
1583 const char decimal_point = localeconv()->decimal_point[0];
1584 #else
1585 const char decimal_point = '.';
1586 #endif */
1587 if (ok != 0)
1588 *ok = true;
1589
1590 const char decimal_point = '.';
1591
1592 sign = nz0 = nz = 0;
1593 rv = 0.;
1594
1595
1596 for(s = s00; isspace(uchar(*s)); s++)
1597 ;
1598
1599 if (*s == '-') {
1600 sign = 1;
1601 s++;
1602 } else if (*s == '+') {
1603 s++;
1604 }
1605
1606 if (*s == '\0') {
1607 s = s00;
1608 goto ret;
1609 }
1610
1611 if (*s == '0') {
1612 nz0 = 1;
1613 while(*++s == '0') ;
1614 if (!*s)
1615 goto ret;
1616 }
1617 s0 = s;
1618 y = z = 0;
1619 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1620 if (nd < 9)
1621 y = 10*y + c - '0';
1622 else if (nd < 16)
1623 z = 10*z + c - '0';
1624 nd0 = nd;
1625 if (c == decimal_point) {
1626 c = *++s;
1627 if (!nd) {
1628 for(; c == '0'; c = *++s)
1629 nz++;
1630 if (c > '0' && c <= '9') {
1631 s0 = s;
1632 nf += nz;
1633 nz = 0;
1634 goto have_dig;
1635 }
1636 goto dig_done;
1637 }
1638 for(; c >= '0' && c <= '9'; c = *++s) {
1639 have_dig:
1640 nz++;
1641 if (c -= '0') {
1642 nf += nz;
1643 for(i = 1; i < nz; i++)
1644 if (nd++ < 9)
1645 y *= 10;
1646 else if (nd <= DBL_DIG + 1)
1647 z *= 10;
1648 if (nd++ < 9)
1649 y = 10*y + c;
1650 else if (nd <= DBL_DIG + 1)
1651 z = 10*z + c;
1652 nz = 0;
1653 }
1654 }
1655 }
1656 dig_done:
1657 e = 0;
1658 if (c == 'e' || c == 'E') {
1659 if (!nd && !nz && !nz0) {
1660 s = s00;
1661 goto ret;
1662 }
1663 s00 = s;
1664 esign = 0;
1665 switch(c = *++s) {
1666 case '-':
1667 esign = 1;
1668 case '+':
1669 c = *++s;
1670 }
1671 if (c >= '0' && c <= '9') {
1672 while(c == '0')
1673 c = *++s;
1674 if (c > '0' && c <= '9') {
1675 L = c - '0';
1676 s1 = s;
1677 while((c = *++s) >= '0' && c <= '9')
1678 L = 10*L + c - '0';
1679 if (s - s1 > 8 || L > 19999)
1680 /* Avoid confusion from exponents
1681 * so large that e might overflow.
1682 */
1683 e = 19999; /* safe for 16 bit ints */
1684 else
1685 e = int(L);
1686 if (esign)
1687 e = -e;
1688 }
1689 else
1690 e = 0;
1691 }
1692 else
1693 s = s00;
1694 }
1695 if (!nd) {
1696 if (!nz && !nz0)
1697 s = s00;
1698 goto ret;
1699 }
1700 e1 = e -= nf;
1701
1702 /* Now we have nd0 digits, starting at s0, followed by a
1703 * decimal point, followed by nd-nd0 digits. The number we're
1704 * after is the integer represented by those digits times
1705 * 10**e */
1706
1707 if (!nd0)
1708 nd0 = nd;
1709 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1710 rv = y;
1711 if (k > 9)
1712#if defined(Q_OS_IRIX) && defined(Q_CC_GNU)
1713 {
1714 // work around a bug on 64 bit IRIX gcc
1715 double *t = (double *) tens;
1716 rv = t[k - 9] * rv + z;
1717 }
1718#else
1719 rv = tens[k - 9] * rv + z;
1720#endif
1721
1722 bd0 = 0;
1723 if (nd <= DBL_DIG
1724#ifndef RND_PRODQUOT
1725 && FLT_ROUNDS == 1
1726#endif
1727 ) {
1728 if (!e)
1729 goto ret;
1730 if (e > 0) {
1731 if (e <= Ten_pmax) {
1732#ifdef VAX
1733 goto vax_ovfl_check;
1734#else
1735 /* rv = */ rounded_product(rv, tens[e]);
1736 goto ret;
1737#endif
1738 }
1739 i = DBL_DIG - nd;
1740 if (e <= Ten_pmax + i) {
1741 /* A fancier test would sometimes let us do
1742 * this for larger i values.
1743 */
1744 e -= i;
1745 rv *= tens[i];
1746#ifdef VAX
1747 /* VAX exponent range is so narrow we must
1748 * worry about overflow here...
1749 */
1750 vax_ovfl_check:
1751 setWord0(&rv, getWord0(rv) - P*Exp_msk1);
1752 /* rv = */ rounded_product(rv, tens[e]);
1753 if ((getWord0(rv) & Exp_mask)
1754 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1755 goto ovfl;
1756 setWord0(&rv, getWord0(rv) + P*Exp_msk1);
1757#else
1758 /* rv = */ rounded_product(rv, tens[e]);
1759#endif
1760 goto ret;
1761 }
1762 }
1763#ifndef Inaccurate_Divide
1764 else if (e >= -Ten_pmax) {
1765 /* rv = */ rounded_quotient(rv, tens[-e]);
1766 goto ret;
1767 }
1768#endif
1769 }
1770 e1 += nd - k;
1771
1772 /* Get starting approximation = rv * 10**e1 */
1773
1774 if (e1 > 0) {
1775 if ((i = e1 & 15) != 0)
1776 rv *= tens[i];
1777 if (e1 &= ~15) {
1778 if (e1 > DBL_MAX_10_EXP) {
1779 ovfl:
1780 // errno = ERANGE;
1781 if (ok != 0)
1782 *ok = false;
1783#ifdef __STDC__
1784 rv = HUGE_VAL;
1785#else
1786 /* Can't trust HUGE_VAL */
1787#ifdef IEEE_Arith
1788 setWord0(&rv, Exp_mask);
1789 setWord1(&rv, 0);
1790#else
1791 setWord0(&rv, Big0);
1792 setWord1(&rv, Big1);
1793#endif
1794#endif
1795 if (bd0)
1796 goto retfree;
1797 goto ret;
1798 }
1799 if (e1 >>= 4) {
1800 for(j = 0; e1 > 1; j++, e1 >>= 1)
1801 if (e1 & 1)
1802 rv *= bigtens[j];
1803 /* The last multiplication could overflow. */
1804 setWord0(&rv, getWord0(rv) - P*Exp_msk1);
1805 rv *= bigtens[j];
1806 if ((z = getWord0(rv) & Exp_mask)
1807 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1808 goto ovfl;
1809 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1810 /* set to largest number */
1811 /* (Can't trust DBL_MAX) */
1812 setWord0(&rv, Big0);
1813 setWord1(&rv, Big1);
1814 }
1815 else
1816 setWord0(&rv, getWord0(rv) + P*Exp_msk1);
1817 }
1818
1819 }
1820 }
1821 else if (e1 < 0) {
1822 e1 = -e1;
1823 if ((i = e1 & 15) != 0)
1824 rv /= tens[i];
1825 if (e1 &= ~15) {
1826 e1 >>= 4;
1827 if (e1 >= 1 << n_bigtens)
1828 goto undfl;
1829 for(j = 0; e1 > 1; j++, e1 >>= 1)
1830 if (e1 & 1)
1831 rv *= tinytens[j];
1832 /* The last multiplication could underflow. */
1833 rv0 = rv;
1834 rv *= tinytens[j];
1835 if (rv == g_double_zero)
1836 {
1837 rv = 2.*rv0;
1838 rv *= tinytens[j];
1839 if (rv == g_double_zero)
1840 {
1841 undfl:
1842 rv = 0.;
1843 // errno = ERANGE;
1844 if (ok != 0)
1845 *ok = false;
1846 if (bd0)
1847 goto retfree;
1848 goto ret;
1849 }
1850 setWord0(&rv, Tiny0);
1851 setWord1(&rv, Tiny1);
1852 /* The refinement below will clean
1853 * this approximation up.
1854 */
1855 }
1856 }
1857 }
1858
1859 /* Now the hard part -- adjusting rv to the correct value.*/
1860
1861 /* Put digits into bd: true value = bd * 10^e */
1862
1863 bd0 = s2b(s0, nd0, nd, y);
1864
1865 for(;;) {
1866 bd = Balloc(bd0->k);
1867 Bcopy(bd, bd0);
1868 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1869 bs = i2b(1);
1870
1871 if (e >= 0) {
1872 bb2 = bb5 = 0;
1873 bd2 = bd5 = e;
1874 }
1875 else {
1876 bb2 = bb5 = -e;
1877 bd2 = bd5 = 0;
1878 }
1879 if (bbe >= 0)
1880 bb2 += bbe;
1881 else
1882 bd2 -= bbe;
1883 bs2 = bb2;
1884#ifdef Sudden_Underflow
1885#ifdef IBM
1886 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1887#else
1888 j = P + 1 - bbbits;
1889#endif
1890#else
1891 i = bbe + bbbits - 1; /* logb(rv) */
1892 if (i < Emin) /* denormal */
1893 j = bbe + (P-Emin);
1894 else
1895 j = P + 1 - bbbits;
1896#endif
1897 bb2 += j;
1898 bd2 += j;
1899 i = bb2 < bd2 ? bb2 : bd2;
1900 if (i > bs2)
1901 i = bs2;
1902 if (i > 0) {
1903 bb2 -= i;
1904 bd2 -= i;
1905 bs2 -= i;
1906 }
1907 if (bb5 > 0) {
1908 bs = pow5mult(bs, bb5);
1909 bb1 = mult(bs, bb);
1910 Bfree(bb);
1911 bb = bb1;
1912 }
1913 if (bb2 > 0)
1914 bb = lshift(bb, bb2);
1915 if (bd5 > 0)
1916 bd = pow5mult(bd, bd5);
1917 if (bd2 > 0)
1918 bd = lshift(bd, bd2);
1919 if (bs2 > 0)
1920 bs = lshift(bs, bs2);
1921 delta = diff(bb, bd);
1922 dsign = delta->sign;
1923 delta->sign = 0;
1924 i = cmp(delta, bs);
1925 if (i < 0) {
1926 /* Error is less than half an ulp -- check for
1927 * special case of mantissa a power of two.
1928 */
1929 if (dsign || getWord1(rv) || getWord0(rv) & Bndry_mask)
1930 break;
1931 delta = lshift(delta,Log2P);
1932 if (cmp(delta, bs) > 0)
1933 goto drop_down;
1934 break;
1935 }
1936 if (i == 0) {
1937 /* exactly half-way between */
1938 if (dsign) {
1939 if ((getWord0(rv) & Bndry_mask1) == Bndry_mask1
1940 && getWord1(rv) == 0xffffffff) {
1941 /*boundary case -- increment exponent*/
1942 setWord0(&rv, (getWord0(rv) & Exp_mask)
1943 + Exp_msk1
1944#ifdef IBM
1945 | Exp_msk1 >> 4
1946#endif
1947 );
1948 setWord1(&rv, 0);
1949 break;
1950 }
1951 }
1952 else if (!(getWord0(rv) & Bndry_mask) && !getWord1(rv)) {
1953 drop_down:
1954 /* boundary case -- decrement exponent */
1955#ifdef Sudden_Underflow
1956 L = getWord0(rv) & Exp_mask;
1957#ifdef IBM
1958 if (L < Exp_msk1)
1959#else
1960 if (L <= Exp_msk1)
1961#endif
1962 goto undfl;
1963 L -= Exp_msk1;
1964#else
1965 L = (getWord0(rv) & Exp_mask) - Exp_msk1;
1966#endif
1967 setWord0(&rv, L | Bndry_mask1);
1968 setWord1(&rv, 0xffffffff);
1969#ifdef IBM
1970 goto cont;
1971#else
1972 break;
1973#endif
1974 }
1975#ifndef ROUND_BIASED
1976 if (!(getWord1(rv) & LSB))
1977 break;
1978#endif
1979 if (dsign)
1980 rv += ulp(rv);
1981#ifndef ROUND_BIASED
1982 else {
1983 rv -= ulp(rv);
1984#ifndef Sudden_Underflow
1985 if (rv == g_double_zero)
1986 goto undfl;
1987#endif
1988 }
1989#endif
1990 break;
1991 }
1992 if ((aadj = ratio(delta, bs)) <= 2.) {
1993 if (dsign)
1994 aadj = aadj1 = 1.;
1995 else if (getWord1(rv) || getWord0(rv) & Bndry_mask) {
1996#ifndef Sudden_Underflow
1997 if (getWord1(rv) == Tiny1 && !getWord0(rv))
1998 goto undfl;
1999#endif
2000 aadj = 1.;
2001 aadj1 = -1.;
2002 }
2003 else {
2004 /* special case -- power of FLT_RADIX to be */
2005 /* rounded down... */
2006
2007 if (aadj < 2./FLT_RADIX)
2008 aadj = 1./FLT_RADIX;
2009 else
2010 aadj *= 0.5;
2011 aadj1 = -aadj;
2012 }
2013 }
2014 else {
2015 aadj *= 0.5;
2016 aadj1 = dsign ? aadj : -aadj;
2017#ifdef Check_FLT_ROUNDS
2018 switch(FLT_ROUNDS) {
2019 case 2: /* towards +infinity */
2020 aadj1 -= 0.5;
2021 break;
2022 case 0: /* towards 0 */
2023 case 3: /* towards -infinity */
2024 aadj1 += 0.5;
2025 }
2026#else
2027 if (FLT_ROUNDS == 0)
2028 aadj1 += 0.5;
2029#endif
2030 }
2031 y = getWord0(rv) & Exp_mask;
2032
2033 /* Check for overflow */
2034
2035 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2036 rv0 = rv;
2037 setWord0(&rv, getWord0(rv) - P*Exp_msk1);
2038 adj = aadj1 * ulp(rv);
2039 rv += adj;
2040 if ((getWord0(rv) & Exp_mask) >=
2041 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2042 if (getWord0(rv0) == Big0 && getWord1(rv0) == Big1)
2043 goto ovfl;
2044 setWord0(&rv, Big0);
2045 setWord1(&rv, Big1);
2046 goto cont;
2047 }
2048 else
2049 setWord0(&rv, getWord0(rv) + P*Exp_msk1);
2050 }
2051 else {
2052#ifdef Sudden_Underflow
2053 if ((getWord0(rv) & Exp_mask) <= P*Exp_msk1) {
2054 rv0 = rv;
2055 setWord0(&rv, getWord0(rv) + P*Exp_msk1);
2056 adj = aadj1 * ulp(rv);
2057 rv += adj;
2058#ifdef IBM
2059 if ((getWord0(rv) & Exp_mask) < P*Exp_msk1)
2060#else
2061 if ((getWord0(rv) & Exp_mask) <= P*Exp_msk1)
2062#endif
2063 {
2064 if (getWord0(rv0) == Tiny0
2065 && getWord1(rv0) == Tiny1)
2066 goto undfl;
2067 setWord0(&rv, Tiny0);
2068 setWord1(&rv, Tiny1);
2069 goto cont;
2070 }
2071 else
2072 setWord0(&rv, getWord0(rv) - P*Exp_msk1);
2073 }
2074 else {
2075 adj = aadj1 * ulp(rv);
2076 rv += adj;
2077 }
2078#else
2079 /* Compute adj so that the IEEE rounding rules will
2080 * correctly round rv + adj in some half-way cases.
2081 * If rv * ulp(rv) is denormalized (i.e.,
2082 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2083 * trouble from bits lost to denormalization;
2084 * example: 1.2e-307 .
2085 */
2086 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
2087 aadj1 = int(aadj + 0.5);
2088 if (!dsign)
2089 aadj1 = -aadj1;
2090 }
2091 adj = aadj1 * ulp(rv);
2092 rv += adj;
2093#endif
2094 }
2095 z = getWord0(rv) & Exp_mask;
2096 if (y == z) {
2097 /* Can we stop now? */
2098 L = Long(aadj);
2099 aadj -= L;
2100 /* The tolerances below are conservative. */
2101 if (dsign || getWord1(rv) || getWord0(rv) & Bndry_mask) {
2102 if (aadj < .4999999 || aadj > .5000001)
2103 break;
2104 }
2105 else if (aadj < .4999999/FLT_RADIX)
2106 break;
2107 }
2108 cont:
2109 Bfree(bb);
2110 Bfree(bd);
2111 Bfree(bs);
2112 Bfree(delta);
2113 }
2114 retfree:
2115 Bfree(bb);
2116 Bfree(bd);
2117 Bfree(bs);
2118 Bfree(bd0);
2119 Bfree(delta);
2120 ret:
2121 if (se)
2122 *se = s;
2123 return sign ? -rv : rv;
2124}
2125
2126static int quorem(Bigint *b, Bigint *S)
2127{
2128 int n;
2129 Long borrow, y;
2130 ULong carry, q, ys;
2131 ULong *bx, *bxe, *sx, *sxe;
2132#ifdef Pack_32
2133 Long z;
2134 ULong si, zs;
2135#endif
2136
2137 n = S->wds;
2138#ifdef BSD_QDTOA_DEBUG
2139 /*debug*/ if (b->wds > n)
2140 /*debug*/ Bug("oversize b in quorem");
2141#endif
2142 if (b->wds < n)
2143 return 0;
2144 sx = S->x;
2145 sxe = sx + --n;
2146 bx = b->x;
2147 bxe = bx + n;
2148 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2149#ifdef BSD_QDTOA_DEBUG
2150 /*debug*/ if (q > 9)
2151 /*debug*/ Bug("oversized quotient in quorem");
2152#endif
2153 if (q) {
2154 borrow = 0;
2155 carry = 0;
2156 do {
2157#ifdef Pack_32
2158 si = *sx++;
2159 ys = (si & 0xffff) * q + carry;
2160 zs = (si >> 16) * q + (ys >> 16);
2161 carry = zs >> 16;
2162 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
2163 borrow = y >> 16;
2164 Sign_Extend(borrow, y);
2165 z = (*bx >> 16) - (zs & 0xffff) + borrow;
2166 borrow = z >> 16;
2167 Sign_Extend(borrow, z);
2168 Storeinc(bx, z, y);
2169#else
2170 ys = *sx++ * q + carry;
2171 carry = ys >> 16;
2172 y = *bx - (ys & 0xffff) + borrow;
2173 borrow = y >> 16;
2174 Sign_Extend(borrow, y);
2175 *bx++ = y & 0xffff;
2176#endif
2177 }
2178 while(sx <= sxe);
2179 if (!*bxe) {
2180 bx = b->x;
2181 while(--bxe > bx && !*bxe)
2182 --n;
2183 b->wds = n;
2184 }
2185 }
2186 if (cmp(b, S) >= 0) {
2187 q++;
2188 borrow = 0;
2189 carry = 0;
2190 bx = b->x;
2191 sx = S->x;
2192 do {
2193#ifdef Pack_32
2194 si = *sx++;
2195 ys = (si & 0xffff) + carry;
2196 zs = (si >> 16) + (ys >> 16);
2197 carry = zs >> 16;
2198 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
2199 borrow = y >> 16;
2200 Sign_Extend(borrow, y);
2201 z = (*bx >> 16) - (zs & 0xffff) + borrow;
2202 borrow = z >> 16;
2203 Sign_Extend(borrow, z);
2204 Storeinc(bx, z, y);
2205#else
2206 ys = *sx++ + carry;
2207 carry = ys >> 16;
2208 y = *bx - (ys & 0xffff) + borrow;
2209 borrow = y >> 16;
2210 Sign_Extend(borrow, y);
2211 *bx++ = y & 0xffff;
2212#endif
2213 }
2214 while(sx <= sxe);
2215 bx = b->x;
2216 bxe = bx + n;
2217 if (!*bxe) {
2218 while(--bxe > bx && !*bxe)
2219 --n;
2220 b->wds = n;
2221 }
2222 }
2223 return q;
2224}
2225
2226/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2227 *
2228 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2229 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
2230 *
2231 * Modifications:
2232 * 1. Rather than iterating, we use a simple numeric overestimate
2233 * to determine k = floor(log10(d)). We scale relevant
2234 * quantities using O(log2(k)) rather than O(k) multiplications.
2235 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2236 * try to generate digits strictly left to right. Instead, we
2237 * compute with fewer bits and propagate the carry if necessary
2238 * when rounding the final digit up. This is often faster.
2239 * 3. Under the assumption that input will be rounded nearest,
2240 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2241 * That is, we allow equality in stopping tests when the
2242 * round-nearest rule will give the same floating-point value
2243 * as would satisfaction of the stopping test with strict
2244 * inequality.
2245 * 4. We remove common factors of powers of 2 from relevant
2246 * quantities.
2247 * 5. When converting floating-point integers less than 1e16,
2248 * we use floating-point arithmetic rather than resorting
2249 * to multiple-precision integers.
2250 * 6. When asked to produce fewer than 15 digits, we first try
2251 * to get by with floating-point arithmetic; we resort to
2252 * multiple-precision integer arithmetic only if we cannot
2253 * guarantee that the floating-point calculation has given
2254 * the correctly rounded result. For k requested digits and
2255 * "uniformly" distributed input, the probability is
2256 * something like 10^(k-15) that we must resort to the Long
2257 * calculation.
2258 */
2259
2260#if defined(Q_OS_WIN) && defined (Q_CC_GNU) && !defined(_clear87) // See QTBUG-7576
2261extern "C" {
2262__attribute__ ((dllimport)) unsigned int __cdecl __MINGW_NOTHROW _control87 (unsigned int unNew, unsigned int unMask);
2263__attribute__ ((dllimport)) unsigned int __cdecl __MINGW_NOTHROW _clearfp (void); /* Clear the FPU status word */
2264}
2265# define _clear87 _clearfp
2266#endif
2267
2268/* This actually sometimes returns a pointer to a string literal
2269 cast to a char*. Do NOT try to modify the return value. */
2270
2271Q_CORE_EXPORT char *qdtoa ( double d, int mode, int ndigits, int *decpt, int *sign, char **rve, char **resultp)
2272{
2273 // Some values of the floating-point control word can cause _qdtoa to crash with an underflow.
2274 // We set a safe value here.
2275#ifdef Q_OS_WIN
2276 _clear87();
2277 unsigned int oldbits = _control87(0, 0);
2278#ifndef MCW_EM
2279# ifdef _MCW_EM
2280# define MCW_EM _MCW_EM
2281# else
2282# define MCW_EM 0x0008001F
2283# endif
2284#endif
2285 _control87(MCW_EM, MCW_EM);
2286#endif
2287
2288#if defined(Q_OS_LINUX) && !defined(__UCLIBC__)
2289 fenv_t envp;
2290 feholdexcept(&envp);
2291#endif
2292
2293 char *s = _qdtoa(d, mode, ndigits, decpt, sign, rve, resultp);
2294
2295#ifdef Q_OS_WIN
2296 _clear87();
2297#ifndef _M_X64
2298 _control87(oldbits, 0xFFFFF);
2299#else
2300 _control87(oldbits, _MCW_EM|_MCW_DN|_MCW_RC);
2301#endif //_M_X64
2302#endif //Q_OS_WIN
2303
2304#if defined(Q_OS_LINUX) && !defined(__UCLIBC__)
2305 fesetenv(&envp);
2306#endif
2307
2308 return s;
2309}
2310
2311static char *_qdtoa( NEEDS_VOLATILE double d, int mode, int ndigits, int *decpt, int *sign, char **rve, char **resultp)
2312{
2313 /*
2314 Arguments ndigits, decpt, sign are similar to those
2315 of ecvt and fcvt; trailing zeros are suppressed from
2316 the returned string. If not null, *rve is set to point
2317 to the end of the return value. If d is +-Infinity or NaN,
2318 then *decpt is set to 9999.
2319
2320 mode:
2321 0 ==> shortest string that yields d when read in
2322 and rounded to nearest.
2323 1 ==> like 0, but with Steele & White stopping rule;
2324 e.g. with IEEE P754 arithmetic , mode 0 gives
2325 1e23 whereas mode 1 gives 9.999999999999999e22.
2326 2 ==> max(1,ndigits) significant digits. This gives a
2327 return value similar to that of ecvt, except
2328 that trailing zeros are suppressed.
2329 3 ==> through ndigits past the decimal point. This
2330 gives a return value similar to that from fcvt,
2331 except that trailing zeros are suppressed, and
2332 ndigits can be negative.
2333 4-9 should give the same return values as 2-3, i.e.,
2334 4 <= mode <= 9 ==> same return as mode
2335 2 + (mode & 1). These modes are mainly for
2336 debugging; often they run slower but sometimes
2337 faster than modes 2-3.
2338 4,5,8,9 ==> left-to-right digit generation.
2339 6-9 ==> don't try fast floating-point estimate
2340 (if applicable).
2341
2342 Values of mode other than 0-9 are treated as mode 0.
2343
2344 Sufficient space is allocated to the return value
2345 to hold the suppressed trailing zeros.
2346 */
2347
2348 int bbits, b2, b5, be, dig, i, ieps, ilim0,
2349 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2350 try_quick;
2351 int ilim = 0, ilim1 = 0, spec_case = 0; /* pacify gcc */
2352 Long L;
2353#ifndef Sudden_Underflow
2354 int denorm;
2355 ULong x;
2356#endif
2357 Bigint *b, *b1, *delta, *mhi, *S;
2358 Bigint *mlo = NULL; /* pacify gcc */
2359 double d2;
2360 double ds, eps;
2361 char *s, *s0;
2362
2363 if (getWord0(d) & Sign_bit) {
2364 /* set sign for everything, including 0's and NaNs */
2365 *sign = 1;
2366 setWord0(&d, getWord0(d) & ~Sign_bit); /* clear sign bit */
2367 }
2368 else
2369 *sign = 0;
2370
2371#if defined(IEEE_Arith) + defined(VAX)
2372#ifdef IEEE_Arith
2373 if ((getWord0(d) & Exp_mask) == Exp_mask)
2374#else
2375 if (getWord0(d) == 0x8000)
2376#endif
2377 {
2378 /* Infinity or NaN */
2379 *decpt = 9999;
2380 s =
2381#ifdef IEEE_Arith
2382 !getWord1(d) && !(getWord0(d) & 0xfffff) ? const_cast<char*>("Infinity") :
2383#endif
2384 const_cast<char*>("NaN");
2385 if (rve)
2386 *rve =
2387#ifdef IEEE_Arith
2388 s[3] ? s + 8 :
2389#endif
2390 s + 3;
2391 return s;
2392 }
2393#endif
2394#ifdef IBM
2395 d += 0; /* normalize */
2396#endif
2397 if (d == g_double_zero)
2398 {
2399 *decpt = 1;
2400 s = const_cast<char*>("0");
2401 if (rve)
2402 *rve = s + 1;
2403 return s;
2404 }
2405
2406 b = d2b(d, &be, &bbits);
2407#ifdef Sudden_Underflow
2408 i = (int)(getWord0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2409#else
2410 if ((i = int(getWord0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) {
2411#endif
2412 d2 = d;
2413 setWord0(&d2, getWord0(d2) & Frac_mask1);
2414 setWord0(&d2, getWord0(d2) | Exp_11);
2415#ifdef IBM
2416 if (j = 11 - hi0bits(getWord0(d2) & Frac_mask))
2417 d2 /= 1 << j;
2418#endif
2419
2420 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2421 * log10(x) = log(x) / log(10)
2422 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2423 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2424 *
2425 * This suggests computing an approximation k to log10(d) by
2426 *
2427 * k = (i - Bias)*0.301029995663981
2428 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2429 *
2430 * We want k to be too large rather than too small.
2431 * The error in the first-order Taylor series approximation
2432 * is in our favor, so we just round up the constant enough
2433 * to compensate for any error in the multiplication of
2434 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2435 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2436 * adding 1e-13 to the constant term more than suffices.
2437 * Hence we adjust the constant term to 0.1760912590558.
2438 * (We could get a more accurate k by invoking log10,
2439 * but this is probably not worthwhile.)
2440 */
2441
2442 i -= Bias;
2443#ifdef IBM
2444 i <<= 2;
2445 i += j;
2446#endif
2447#ifndef Sudden_Underflow
2448 denorm = 0;
2449 }
2450 else {
2451 /* d is denormalized */
2452
2453 i = bbits + be + (Bias + (P-1) - 1);
2454 x = i > 32 ? getWord0(d) << (64 - i) | getWord1(d) >> (i - 32)
2455 : getWord1(d) << (32 - i);
2456 d2 = x;
2457 setWord0(&d2, getWord0(d2) - 31*Exp_msk1); /* adjust exponent */
2458 i -= (Bias + (P-1) - 1) + 1;
2459 denorm = 1;
2460 }
2461#endif
2462 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2463 k = int(ds);
2464 if (ds < 0. && ds != k)
2465 k--; /* want k = floor(ds) */
2466 k_check = 1;
2467 if (k >= 0 && k <= Ten_pmax) {
2468 if (d < tens[k])
2469 k--;
2470 k_check = 0;
2471 }
2472 j = bbits - i - 1;
2473 if (j >= 0) {
2474 b2 = 0;
2475 s2 = j;
2476 }
2477 else {
2478 b2 = -j;
2479 s2 = 0;
2480 }
2481 if (k >= 0) {
2482 b5 = 0;
2483 s5 = k;
2484 s2 += k;
2485 }
2486 else {
2487 b2 -= k;
2488 b5 = -k;
2489 s5 = 0;
2490 }
2491 if (mode < 0 || mode > 9)
2492 mode = 0;
2493 try_quick = 1;
2494 if (mode > 5) {
2495 mode -= 4;
2496 try_quick = 0;
2497 }
2498 leftright = 1;
2499 switch(mode) {
2500 case 0:
2501 case 1:
2502 ilim = ilim1 = -1;
2503 i = 18;
2504 ndigits = 0;
2505 break;
2506 case 2:
2507 leftright = 0;
2508 /* no break */
2509 case 4:
2510 if (ndigits <= 0)
2511 ndigits = 1;
2512 ilim = ilim1 = i = ndigits;
2513 break;
2514 case 3:
2515 leftright = 0;
2516 /* no break */
2517 case 5:
2518 i = ndigits + k + 1;
2519 ilim = i;
2520 ilim1 = i - 1;
2521 if (i <= 0)
2522 i = 1;
2523 }
2524 QT_TRY {
2525 *resultp = static_cast<char *>(malloc(i + 1));
2526 Q_CHECK_PTR(*resultp);
2527 } QT_CATCH(...) {
2528 Bfree(b);
2529 QT_RETHROW;
2530 }
2531 s = s0 = *resultp;
2532
2533 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2534
2535 /* Try to get by with floating-point arithmetic. */
2536
2537 i = 0;
2538 d2 = d;
2539 k0 = k;
2540 ilim0 = ilim;
2541 ieps = 2; /* conservative */
2542 if (k > 0) {
2543 ds = tens[k&0xf];
2544 j = k >> 4;
2545 if (j & Bletch) {
2546 /* prevent overflows */
2547 j &= Bletch - 1;
2548 d /= bigtens[n_bigtens-1];
2549 ieps++;
2550 }
2551 for(; j; j >>= 1, i++)
2552 if (j & 1) {
2553 ieps++;
2554 ds *= bigtens[i];
2555 }
2556 d /= ds;
2557 }
2558 else if ((j1 = -k) != 0) {
2559 d *= tens[j1 & 0xf];
2560 for(j = j1 >> 4; j; j >>= 1, i++)
2561 if (j & 1) {
2562 ieps++;
2563 d *= bigtens[i];
2564 }
2565 }
2566 if (k_check && d < 1. && ilim > 0) {
2567 if (ilim1 <= 0)
2568 goto fast_failed;
2569 ilim = ilim1;
2570 k--;
2571 d *= 10.;
2572 ieps++;
2573 }
2574 eps = ieps*d + 7.;
2575 setWord0(&eps, getWord0(eps) - (P-1)*Exp_msk1);
2576 if (ilim == 0) {
2577 S = mhi = 0;
2578 d -= 5.;
2579 if (d > eps)
2580 goto one_digit;
2581 if (d < -eps)
2582 goto no_digits;
2583 goto fast_failed;
2584 }
2585#ifndef No_leftright
2586 if (leftright) {
2587 /* Use Steele & White method of only
2588 * generating digits needed.
2589 */
2590 eps = 0.5/tens[ilim-1] - eps;
2591 for(i = 0;;) {
2592 L = Long(d);
2593 d -= L;
2594 *s++ = '0' + int(L);
2595 if (d < eps)
2596 goto ret1;
2597 if (1. - d < eps)
2598 goto bump_up;
2599 if (++i >= ilim)
2600 break;
2601 eps *= 10.;
2602 d *= 10.;
2603 }
2604 }
2605 else {
2606#endif
2607 /* Generate ilim digits, then fix them up. */
2608#if defined(Q_OS_IRIX) && defined(Q_CC_GNU)
2609 // work around a bug on 64 bit IRIX gcc
2610 double *t = (double *) tens;
2611 eps *= t[ilim-1];
2612#else
2613 eps *= tens[ilim-1];
2614#endif
2615 for(i = 1;; i++, d *= 10.) {
2616 L = Long(d);
2617 d -= L;
2618 *s++ = '0' + int(L);
2619 if (i == ilim) {
2620 if (d > 0.5 + eps)
2621 goto bump_up;
2622 else if (d < 0.5 - eps) {
2623 while(*--s == '0') {}
2624 s++;
2625 goto ret1;
2626 }
2627 break;
2628 }
2629 }
2630#ifndef No_leftright
2631 }
2632#endif
2633 fast_failed:
2634 s = s0;
2635 d = d2;
2636 k = k0;
2637 ilim = ilim0;
2638 }
2639
2640 /* Do we have a "small" integer? */
2641
2642 if (be >= 0 && k <= Int_max) {
2643 /* Yes. */
2644 ds = tens[k];
2645 if (ndigits < 0 && ilim <= 0) {
2646 S = mhi = 0;
2647 if (ilim < 0 || d <= 5*ds)
2648 goto no_digits;
2649 goto one_digit;
2650 }
2651 for(i = 1;; i++) {
2652 L = Long(d / ds);
2653 d -= L*ds;
2654#ifdef Check_FLT_ROUNDS
2655 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2656 if (d < 0) {
2657 L--;
2658 d += ds;
2659 }
2660#endif
2661 *s++ = '0' + int(L);
2662 if (i == ilim) {
2663 d += d;
2664 if (d > ds || (d == ds && L & 1)) {
2665 bump_up:
2666 while(*--s == '9')
2667 if (s == s0) {
2668 k++;
2669 *s = '0';
2670 break;
2671 }
2672 ++*s++;
2673 }
2674 break;
2675 }
2676 if ((d *= 10.) == g_double_zero)
2677 break;
2678 }
2679 goto ret1;
2680 }
2681
2682 m2 = b2;
2683 m5 = b5;
2684 mhi = mlo = 0;
2685 if (leftright) {
2686 if (mode < 2) {
2687 i =
2688#ifndef Sudden_Underflow
2689 denorm ? be + (Bias + (P-1) - 1 + 1) :
2690#endif
2691#ifdef IBM
2692 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2693#else
2694 1 + P - bbits;
2695#endif
2696 }
2697 else {
2698 j = ilim - 1;
2699 if (m5 >= j)
2700 m5 -= j;
2701 else {
2702 s5 += j -= m5;
2703 b5 += j;
2704 m5 = 0;
2705 }
2706 if ((i = ilim) < 0) {
2707 m2 -= i;
2708 i = 0;
2709 }
2710 }
2711 b2 += i;
2712 s2 += i;
2713 mhi = i2b(1);
2714 }
2715 if (m2 > 0 && s2 > 0) {
2716 i = m2 < s2 ? m2 : s2;
2717 b2 -= i;
2718 m2 -= i;
2719 s2 -= i;
2720 }
2721 if (b5 > 0) {
2722 if (leftright) {
2723 if (m5 > 0) {
2724 mhi = pow5mult(mhi, m5);
2725 b1 = mult(mhi, b);
2726 Bfree(b);
2727 b = b1;
2728 }
2729 if ((j = b5 - m5) != 0)
2730 b = pow5mult(b, j);
2731 }
2732 else
2733 b = pow5mult(b, b5);
2734 }
2735 S = i2b(1);
2736 if (s5 > 0)
2737 S = pow5mult(S, s5);
2738
2739 /* Check for special case that d is a normalized power of 2. */
2740
2741 if (mode < 2) {
2742 if (!getWord1(d) && !(getWord0(d) & Bndry_mask)
2743#ifndef Sudden_Underflow
2744 && getWord0(d) & Exp_mask
2745#endif
2746 ) {
2747 /* The special case */
2748 b2 += Log2P;
2749 s2 += Log2P;
2750 spec_case = 1;
2751 }
2752 else
2753 spec_case = 0;
2754 }
2755
2756 /* Arrange for convenient computation of quotients:
2757 * shift left if necessary so divisor has 4 leading 0 bits.
2758 *
2759 * Perhaps we should just compute leading 28 bits of S once
2760 * and for all and pass them and a shift to quorem, so it
2761 * can do shifts and ors to compute the numerator for q.
2762 */
2763#ifdef Pack_32
2764 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0)
2765 i = 32 - i;
2766#else
2767 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
2768 i = 16 - i;
2769#endif
2770 if (i > 4) {
2771 i -= 4;
2772 b2 += i;
2773 m2 += i;
2774 s2 += i;
2775 }
2776 else if (i < 4) {
2777 i += 28;
2778 b2 += i;
2779 m2 += i;
2780 s2 += i;
2781 }
2782 if (b2 > 0)
2783 b = lshift(b, b2);
2784 if (s2 > 0)
2785 S = lshift(S, s2);
2786 if (k_check) {
2787 if (cmp(b,S) < 0) {
2788 k--;
2789 b = multadd(b, 10, 0); /* we botched the k estimate */
2790 if (leftright)
2791 mhi = multadd(mhi, 10, 0);
2792 ilim = ilim1;
2793 }
2794 }
2795 if (ilim <= 0 && mode > 2) {
2796 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2797 /* no digits, fcvt style */
2798 no_digits:
2799 k = -1 - ndigits;
2800 goto ret;
2801 }
2802 one_digit:
2803 *s++ = '1';
2804 k++;
2805 goto ret;
2806 }
2807 if (leftright) {
2808 if (m2 > 0)
2809 mhi = lshift(mhi, m2);
2810
2811 /* Compute mlo -- check for special case
2812 * that d is a normalized power of 2.
2813 */
2814
2815 mlo = mhi;
2816 if (spec_case) {
2817 mhi = Balloc(mhi->k);
2818 Bcopy(mhi, mlo);
2819 mhi = lshift(mhi, Log2P);
2820 }
2821
2822 for(i = 1;;i++) {
2823 dig = quorem(b,S) + '0';
2824 /* Do we yet have the shortest decimal string
2825 * that will round to d?
2826 */
2827 j = cmp(b, mlo);
2828 delta = diff(S, mhi);
2829 j1 = delta->sign ? 1 : cmp(b, delta);
2830 Bfree(delta);
2831#ifndef ROUND_BIASED
2832 if (j1 == 0 && !mode && !(getWord1(d) & 1)) {
2833 if (dig == '9')
2834 goto round_9_up;
2835 if (j > 0)
2836 dig++;
2837 *s++ = dig;
2838 goto ret;
2839 }
2840#endif
2841 if (j < 0 || (j == 0 && !mode
2842#ifndef ROUND_BIASED
2843 && !(getWord1(d) & 1)
2844#endif
2845 )) {
2846 if (j1 > 0) {
2847 b = lshift(b, 1);
2848 j1 = cmp(b, S);
2849 if ((j1 > 0 || (j1 == 0 && dig & 1))
2850 && dig++ == '9')
2851 goto round_9_up;
2852 }
2853 *s++ = dig;
2854 goto ret;
2855 }
2856 if (j1 > 0) {
2857 if (dig == '9') { /* possible if i == 1 */
2858 round_9_up:
2859 *s++ = '9';
2860 goto roundoff;
2861 }
2862 *s++ = dig + 1;
2863 goto ret;
2864 }
2865 *s++ = dig;
2866 if (i == ilim)
2867 break;
2868 b = multadd(b, 10, 0);
2869 if (mlo == mhi)
2870 mlo = mhi = multadd(mhi, 10, 0);
2871 else {
2872 mlo = multadd(mlo, 10, 0);
2873 mhi = multadd(mhi, 10, 0);
2874 }
2875 }
2876 }
2877 else
2878 for(i = 1;; i++) {
2879 *s++ = dig = quorem(b,S) + '0';
2880 if (i >= ilim)
2881 break;
2882 b = multadd(b, 10, 0);
2883 }
2884
2885 /* Round off last digit */
2886
2887 b = lshift(b, 1);
2888 j = cmp(b, S);
2889 if (j > 0 || (j == 0 && dig & 1)) {
2890 roundoff:
2891 while(*--s == '9')
2892 if (s == s0) {
2893 k++;
2894 *s++ = '1';
2895 goto ret;
2896 }
2897 ++*s++;
2898 }
2899 else {
2900 while(*--s == '0') {}
2901 s++;
2902 }
2903 ret:
2904 Bfree(S);
2905 if (mhi) {
2906 if (mlo && mlo != mhi)
2907 Bfree(mlo);
2908 Bfree(mhi);
2909 }
2910 ret1:
2911 Bfree(b);
2912 if (s == s0) { /* don't return empty string */
2913 *s++ = '0';
2914 k = 0;
2915 }
2916 *s = 0;
2917 *decpt = k + 1;
2918 if (rve)
2919 *rve = s;
2920 return s0;
2921}
2922#else
2923// NOT thread safe!
2924
2925#include <errno.h>
2926
2927Q_CORE_EXPORT char *qdtoa( double d, int mode, int ndigits, int *decpt, int *sign, char **rve, char **resultp)
2928{
2929 if(rve)
2930 *rve = 0;
2931
2932 char *res;
2933 if (mode == 0)
2934 ndigits = 80;
2935
2936 if (mode == 3)
2937 res = fcvt(d, ndigits, decpt, sign);
2938 else
2939 res = ecvt(d, ndigits, decpt, sign);
2940
2941 int n = qstrlen(res);
2942 if (mode == 0) { // remove trailing 0's
2943 const int stop = qMax(1, *decpt);
2944 int i;
2945 for (i = n-1; i >= stop; --i) {
2946 if (res[i] != '0')
2947 break;
2948 }
2949 n = i + 1;
2950 }
2951 *resultp = static_cast<char*>(malloc(n + 1));
2952 Q_CHECK_PTR(resultp);
2953 qstrncpy(*resultp, res, n + 1);
2954 return *resultp;
2955}
2956
2957Q_CORE_EXPORT double qstrtod(const char *s00, const char **se, bool *ok)
2958{
2959 double ret = strtod((char*)s00, (char**)se);
2960 if (ok) {
2961 if((ret == 0.0l && errno == ERANGE)
2962 || ret == HUGE_VAL || ret == -HUGE_VAL)
2963 *ok = false;
2964 else
2965 *ok = true; // the result will be that we don't report underflow in this case
2966 }
2967 return ret;
2968}
2969
2970#endif // QT_QLOCALE_USES_FCVT
2971
2972QT_END_NAMESPACE
2973