1 : /****************************************************************
2 : *
3 : * The author of this software is David M. Gay.
4 : *
5 : * Copyright (c) 1991 by AT&T.
6 : *
7 : * Permission to use, copy, modify, and distribute this software for any
8 : * purpose without fee is hereby granted, provided that this entire notice
9 : * is included in all copies of any software which is or includes a copy
10 : * or modification of this software and in all copies of the supporting
11 : * documentation for such software.
12 : *
13 : * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14 : * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
15 : * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16 : * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
17 : *
18 : ***************************************************************/
19 :
20 : /* Please send bug reports to
21 : David M. Gay
22 : AT&T Bell Laboratories, Room 2C-463
23 : 600 Mountain Avenue
24 : Murray Hill, NJ 07974-2070
25 : U.S.A.
26 : dmg@research.att.com or research!dmg
27 : */
28 :
29 : /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
30 : *
31 : * This strtod returns a nearest machine number to the input decimal
32 : * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
33 : * broken by the IEEE round-even rule. Otherwise ties are broken by
34 : * biased rounding (add half and chop).
35 : *
36 : * Inspired loosely by William D. Clinger's paper "How to Read Floating
37 : * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
38 : *
39 : * Modifications:
40 : *
41 : * 1. We only require IEEE, IBM, or VAX double-precision
42 : * arithmetic (not IEEE double-extended).
43 : * 2. We get by with floating-point arithmetic in a case that
44 : * Clinger missed -- when we're computing d * 10^n
45 : * for a small integer d and the integer n is not too
46 : * much larger than 22 (the maximum integer k for which
47 : * we can represent 10^k exactly), we may be able to
48 : * compute (d*10^k) * 10^(e-k) with just one roundoff.
49 : * 3. Rather than a bit-at-a-time adjustment of the binary
50 : * result in the hard case, we use floating-point
51 : * arithmetic to determine the adjustment to within
52 : * one bit; only in really hard cases do we need to
53 : * compute a second residual.
54 : * 4. Because of 3., we don't need a large table of powers of 10
55 : * for ten-to-e (just some small tables, e.g. of 10^k
56 : * for 0 <= k <= 22).
57 : */
58 :
59 : /*
60 : * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
61 : * significant byte has the lowest address.
62 : * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
63 : * significant byte has the lowest address.
64 : * #define Long int on machines with 32-bit ints and 64-bit longs.
65 : * #define Sudden_Underflow for IEEE-format machines without gradual
66 : * underflow (i.e., that flush to zero on underflow).
67 : * #define IBM for IBM mainframe-style floating-point arithmetic.
68 : * #define VAX for VAX-style floating-point arithmetic.
69 : * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
70 : * #define No_leftright to omit left-right logic in fast floating-point
71 : * computation of dtoa.
72 : * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
73 : * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
74 : * that use extended-precision instructions to compute rounded
75 : * products and quotients) with IBM.
76 : * #define ROUND_BIASED for IEEE-format with biased rounding.
77 : * #define Inaccurate_Divide for IEEE-format with correctly rounded
78 : * products but inaccurate quotients, e.g., for Intel i860.
79 : * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision
80 : * integer arithmetic. Whether this speeds things up or slows things
81 : * down depends on the machine and the number being converted.
82 : * #define KR_headers for old-style C function headers.
83 : * #define Bad_float_h if your system lacks a float.h or if it does not
84 : * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
85 : * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
86 : * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
87 : * if memory is available and otherwise does something you deem
88 : * appropriate. If MALLOC is undefined, malloc will be invoked
89 : * directly -- and assumed always to succeed.
90 : */
91 :
92 : /* $Id: zend_strtod.c 277398 2009-03-18 10:18:10Z dmitry $ */
93 :
94 : #include <zend_operators.h>
95 : #include <zend_strtod.h>
96 :
97 : #ifdef ZTS
98 : #include <TSRM.h>
99 : #endif
100 :
101 : #include <stddef.h>
102 : #include <stdio.h>
103 : #include <ctype.h>
104 : #include <stdarg.h>
105 : #include <string.h>
106 : #include <stdlib.h>
107 : #include <math.h>
108 :
109 : #ifdef HAVE_LOCALE_H
110 : #include <locale.h>
111 : #endif
112 :
113 : #ifdef HAVE_SYS_TYPES_H
114 : #include <sys/types.h>
115 : #endif
116 :
117 : #if defined(HAVE_INTTYPES_H)
118 : #include <inttypes.h>
119 : #elif defined(HAVE_STDINT_H)
120 : #include <stdint.h>
121 : #endif
122 :
123 : #ifndef HAVE_INT32_T
124 : # if SIZEOF_INT == 4
125 : typedef int int32_t;
126 : # elif SIZEOF_LONG == 4
127 : typedef long int int32_t;
128 : # endif
129 : #endif
130 :
131 : #ifndef HAVE_UINT32_T
132 : # if SIZEOF_INT == 4
133 : typedef unsigned int uint32_t;
134 : # elif SIZEOF_LONG == 4
135 : typedef unsigned long int uint32_t;
136 : # endif
137 : #endif
138 :
139 : #if (defined(__APPLE__) || defined(__APPLE_CC__)) && (defined(__BIG_ENDIAN__) || defined(__LITTLE_ENDIAN__))
140 : # if defined(__LITTLE_ENDIAN__)
141 : # undef WORDS_BIGENDIAN
142 : # else
143 : # if defined(__BIG_ENDIAN__)
144 : # define WORDS_BIGENDIAN
145 : # endif
146 : # endif
147 : #endif
148 :
149 : #ifdef WORDS_BIGENDIAN
150 : #define IEEE_BIG_ENDIAN
151 : #else
152 : #define IEEE_LITTLE_ENDIAN
153 : #endif
154 :
155 : #if defined(__arm__) && !defined(__VFP_FP__)
156 : /*
157 : * * Although the CPU is little endian the FP has different
158 : * * byte and word endianness. The byte order is still little endian
159 : * * but the word order is big endian.
160 : * */
161 : #define IEEE_BIG_ENDIAN
162 : #undef IEEE_LITTLE_ENDIAN
163 : #endif
164 :
165 : #ifdef __vax__
166 : #define VAX
167 : #endif
168 :
169 : #if defined(_MSC_VER)
170 : #define int32_t __int32
171 : #define uint32_t unsigned __int32
172 : #define IEEE_LITTLE_ENDIAN
173 : #endif
174 :
175 : #define Long int32_t
176 : #define ULong uint32_t
177 :
178 : #ifdef __cplusplus
179 : #include "malloc.h"
180 : #include "memory.h"
181 : #else
182 : #ifndef KR_headers
183 : #include "stdlib.h"
184 : #include "string.h"
185 : #include "locale.h"
186 : #else
187 : #include "malloc.h"
188 : #include "memory.h"
189 : #endif
190 : #endif
191 :
192 : #ifdef MALLOC
193 : #ifdef KR_headers
194 : extern char *MALLOC();
195 : #else
196 : extern void *MALLOC(size_t);
197 : #endif
198 : #else
199 : #define MALLOC malloc
200 : #endif
201 :
202 : #include "ctype.h"
203 : #include "errno.h"
204 :
205 : #ifdef Bad_float_h
206 : #ifdef IEEE_BIG_ENDIAN
207 : #define IEEE_ARITHMETIC
208 : #endif
209 : #ifdef IEEE_LITTLE_ENDIAN
210 : #define IEEE_ARITHMETIC
211 : #endif
212 :
213 : #ifdef IEEE_ARITHMETIC
214 : #define DBL_DIG 15
215 : #define DBL_MAX_10_EXP 308
216 : #define DBL_MAX_EXP 1024
217 : #define FLT_RADIX 2
218 : #define FLT_ROUNDS 1
219 : #define DBL_MAX 1.7976931348623157e+308
220 : #endif
221 :
222 : #ifdef IBM
223 : #define DBL_DIG 16
224 : #define DBL_MAX_10_EXP 75
225 : #define DBL_MAX_EXP 63
226 : #define FLT_RADIX 16
227 : #define FLT_ROUNDS 0
228 : #define DBL_MAX 7.2370055773322621e+75
229 : #endif
230 :
231 : #ifdef VAX
232 : #define DBL_DIG 16
233 : #define DBL_MAX_10_EXP 38
234 : #define DBL_MAX_EXP 127
235 : #define FLT_RADIX 2
236 : #define FLT_ROUNDS 1
237 : #define DBL_MAX 1.7014118346046923e+38
238 : #endif
239 :
240 :
241 : #ifndef LONG_MAX
242 : #define LONG_MAX 2147483647
243 : #endif
244 : #else
245 : #include "float.h"
246 : #endif
247 : #ifndef __MATH_H__
248 : #include "math.h"
249 : #endif
250 :
251 : BEGIN_EXTERN_C()
252 :
253 : #ifndef CONST
254 : #ifdef KR_headers
255 : #define CONST /* blank */
256 : #else
257 : #define CONST const
258 : #endif
259 : #endif
260 :
261 : #ifdef Unsigned_Shifts
262 : #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
263 : #else
264 : #define Sign_Extend(a,b) /*no-op*/
265 : #endif
266 :
267 : #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + \
268 : defined(IBM) != 1
269 : Exactly one of IEEE_LITTLE_ENDIAN IEEE_BIG_ENDIAN, VAX, or
270 : IBM should be defined.
271 : #endif
272 :
273 : typedef union {
274 : double d;
275 : ULong ul[2];
276 : } _double;
277 : #define value(x) ((x).d)
278 : #ifdef IEEE_LITTLE_ENDIAN
279 : #define word0(x) ((x).ul[1])
280 : #define word1(x) ((x).ul[0])
281 : #else
282 : #define word0(x) ((x).ul[0])
283 : #define word1(x) ((x).ul[1])
284 : #endif
285 :
286 : /* The following definition of Storeinc is appropriate for MIPS processors.
287 : * An alternative that might be better on some machines is
288 : * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
289 : */
290 : #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__)
291 : #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
292 : ((unsigned short *)a)[0] = (unsigned short)c, a++)
293 : #else
294 : #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
295 : ((unsigned short *)a)[1] = (unsigned short)c, a++)
296 : #endif
297 :
298 : /* #define P DBL_MANT_DIG */
299 : /* Ten_pmax = floor(P*log(2)/log(5)) */
300 : /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
301 : /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
302 : /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
303 :
304 : #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN)
305 : #define Exp_shift 20
306 : #define Exp_shift1 20
307 : #define Exp_msk1 0x100000
308 : #define Exp_msk11 0x100000
309 : #define Exp_mask 0x7ff00000
310 : #define P 53
311 : #define Bias 1023
312 : #define IEEE_Arith
313 : #define Emin (-1022)
314 : #define Exp_1 0x3ff00000
315 : #define Exp_11 0x3ff00000
316 : #define Ebits 11
317 : #define Frac_mask 0xfffff
318 : #define Frac_mask1 0xfffff
319 : #define Ten_pmax 22
320 : #define Bletch 0x10
321 : #define Bndry_mask 0xfffff
322 : #define Bndry_mask1 0xfffff
323 : #define LSB 1
324 : #define Sign_bit 0x80000000
325 : #define Log2P 1
326 : #define Tiny0 0
327 : #define Tiny1 1
328 : #define Quick_max 14
329 : #define Int_max 14
330 : #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
331 : #else
332 : #undef Sudden_Underflow
333 : #define Sudden_Underflow
334 : #ifdef IBM
335 : #define Exp_shift 24
336 : #define Exp_shift1 24
337 : #define Exp_msk1 0x1000000
338 : #define Exp_msk11 0x1000000
339 : #define Exp_mask 0x7f000000
340 : #define P 14
341 : #define Bias 65
342 : #define Exp_1 0x41000000
343 : #define Exp_11 0x41000000
344 : #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
345 : #define Frac_mask 0xffffff
346 : #define Frac_mask1 0xffffff
347 : #define Bletch 4
348 : #define Ten_pmax 22
349 : #define Bndry_mask 0xefffff
350 : #define Bndry_mask1 0xffffff
351 : #define LSB 1
352 : #define Sign_bit 0x80000000
353 : #define Log2P 4
354 : #define Tiny0 0x100000
355 : #define Tiny1 0
356 : #define Quick_max 14
357 : #define Int_max 15
358 : #else /* VAX */
359 : #define Exp_shift 23
360 : #define Exp_shift1 7
361 : #define Exp_msk1 0x80
362 : #define Exp_msk11 0x800000
363 : #define Exp_mask 0x7f80
364 : #define P 56
365 : #define Bias 129
366 : #define Exp_1 0x40800000
367 : #define Exp_11 0x4080
368 : #define Ebits 8
369 : #define Frac_mask 0x7fffff
370 : #define Frac_mask1 0xffff007f
371 : #define Ten_pmax 24
372 : #define Bletch 2
373 : #define Bndry_mask 0xffff007f
374 : #define Bndry_mask1 0xffff007f
375 : #define LSB 0x10000
376 : #define Sign_bit 0x8000
377 : #define Log2P 1
378 : #define Tiny0 0x80
379 : #define Tiny1 0
380 : #define Quick_max 15
381 : #define Int_max 15
382 : #endif
383 : #endif
384 :
385 : #ifndef IEEE_Arith
386 : #define ROUND_BIASED
387 : #endif
388 :
389 : #ifdef RND_PRODQUOT
390 : #define rounded_product(a,b) a = rnd_prod(a, b)
391 : #define rounded_quotient(a,b) a = rnd_quot(a, b)
392 : #ifdef KR_headers
393 : extern double rnd_prod(), rnd_quot();
394 : #else
395 : extern double rnd_prod(double, double), rnd_quot(double, double);
396 : #endif
397 : #else
398 : #define rounded_product(a,b) a *= b
399 : #define rounded_quotient(a,b) a /= b
400 : #endif
401 :
402 : #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
403 : #define Big1 0xffffffff
404 :
405 : #ifndef Just_16
406 : /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
407 : * * This makes some inner loops simpler and sometimes saves work
408 : * * during multiplications, but it often seems to make things slightly
409 : * * slower. Hence the default is now to store 32 bits per Long.
410 : * */
411 : #ifndef Pack_32
412 : #define Pack_32
413 : #endif
414 : #endif
415 :
416 : #define Kmax 15
417 :
418 : struct Bigint {
419 : struct Bigint *next;
420 : int k, maxwds, sign, wds;
421 : ULong x[1];
422 : };
423 :
424 : typedef struct Bigint Bigint;
425 :
426 : /* static variables, multithreading fun! */
427 : static Bigint *freelist[Kmax+1];
428 : static Bigint *p5s;
429 :
430 : static void destroy_freelist(void);
431 :
432 : #ifdef ZTS
433 :
434 : static MUTEX_T dtoa_mutex;
435 : static MUTEX_T pow5mult_mutex;
436 :
437 : #define _THREAD_PRIVATE_MUTEX_LOCK(x) tsrm_mutex_lock(x);
438 : #define _THREAD_PRIVATE_MUTEX_UNLOCK(x) tsrm_mutex_unlock(x);
439 :
440 : #else
441 :
442 : #define _THREAD_PRIVATE_MUTEX_LOCK(x)
443 : #define _THREAD_PRIVATE_MUTEX_UNLOCK(x)
444 :
445 : #endif /* ZTS */
446 :
447 : ZEND_API int zend_startup_strtod(void) /* {{{ */
448 17633 : {
449 : #ifdef ZTS
450 : dtoa_mutex = tsrm_mutex_alloc();
451 : pow5mult_mutex = tsrm_mutex_alloc();
452 : #endif
453 17633 : return 1;
454 : }
455 : /* }}} */
456 : ZEND_API int zend_shutdown_strtod(void) /* {{{ */
457 17665 : {
458 17665 : destroy_freelist();
459 : #ifdef ZTS
460 : tsrm_mutex_free(dtoa_mutex);
461 : dtoa_mutex = NULL;
462 :
463 : tsrm_mutex_free(pow5mult_mutex);
464 : pow5mult_mutex = NULL;
465 : #endif
466 17665 : return 1;
467 : }
468 : /* }}} */
469 :
470 : static Bigint * Balloc(int k)
471 274338 : {
472 : int x;
473 : Bigint *rv;
474 :
475 274338 : if (k > Kmax) {
476 0 : zend_error(E_ERROR, "Balloc() allocation exceeds list boundary");
477 : }
478 :
479 : _THREAD_PRIVATE_MUTEX_LOCK(dtoa_mutex);
480 274338 : if ((rv = freelist[k])) {
481 270089 : freelist[k] = rv->next;
482 : } else {
483 4249 : x = 1 << k;
484 4249 : rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long));
485 4249 : if (!rv) {
486 : _THREAD_PRIVATE_MUTEX_UNLOCK(dtoa_mutex);
487 0 : zend_error(E_ERROR, "Balloc() failed to allocate memory");
488 : }
489 4249 : rv->k = k;
490 4249 : rv->maxwds = x;
491 : }
492 : _THREAD_PRIVATE_MUTEX_UNLOCK(dtoa_mutex);
493 274338 : rv->sign = rv->wds = 0;
494 274338 : return rv;
495 : }
496 :
497 : static void Bfree(Bigint *v)
498 273250 : {
499 273250 : if (v) {
500 : _THREAD_PRIVATE_MUTEX_LOCK(dtoa_mutex);
501 273242 : v->next = freelist[v->k];
502 273242 : freelist[v->k] = v;
503 : _THREAD_PRIVATE_MUTEX_UNLOCK(dtoa_mutex);
504 : }
505 273250 : }
506 :
507 : #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
508 : y->wds*sizeof(Long) + 2*sizeof(int))
509 :
510 : /* return value is only used as a simple string, so mis-aligned parts
511 : * inside the Bigint are not at risk on strict align architectures
512 : */
513 128895 : static char * rv_alloc(int i) {
514 : int j, k, *r;
515 :
516 128895 : j = sizeof(ULong);
517 128895 : for(k = 0;
518 258289 : sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
519 499 : j <<= 1) {
520 499 : k++;
521 : }
522 128895 : r = (int*)Balloc(k);
523 128895 : *r = k;
524 128895 : return (char *)(r+1);
525 : }
526 :
527 :
528 : static char * nrv_alloc(char *s, char **rve, int n)
529 621 : {
530 : char *rv, *t;
531 :
532 621 : t = rv = rv_alloc(n);
533 1927 : while((*t = *s++) !=0) {
534 685 : t++;
535 : }
536 621 : if (rve) {
537 0 : *rve = t;
538 : }
539 621 : return rv;
540 : }
541 :
542 : static Bigint * multadd(Bigint *b, int m, int a) /* multiply by m and add a */
543 50926 : {
544 : int i, wds;
545 : ULong *x, y;
546 : #ifdef Pack_32
547 : ULong xi, z;
548 : #endif
549 : Bigint *b1;
550 :
551 50926 : wds = b->wds;
552 50926 : x = b->x;
553 50926 : i = 0;
554 : do {
555 : #ifdef Pack_32
556 139530 : xi = *x;
557 139530 : y = (xi & 0xffff) * m + a;
558 139530 : z = (xi >> 16) * m + (y >> 16);
559 139530 : a = (int)(z >> 16);
560 139530 : *x++ = (z << 16) + (y & 0xffff);
561 : #else
562 : y = *x * m + a;
563 : a = (int)(y >> 16);
564 : *x++ = y & 0xffff;
565 : #endif
566 : }
567 139530 : while(++i < wds);
568 50926 : if (a) {
569 4603 : if (wds >= b->maxwds) {
570 0 : b1 = Balloc(b->k+1);
571 0 : Bcopy(b1, b);
572 0 : Bfree(b);
573 0 : b = b1;
574 : }
575 4603 : b->x[wds++] = a;
576 4603 : b->wds = wds;
577 : }
578 50926 : return b;
579 : }
580 :
581 : static int hi0bits(ULong x)
582 1568 : {
583 1568 : int k = 0;
584 :
585 1568 : if (!(x & 0xffff0000)) {
586 817 : k = 16;
587 817 : x <<= 16;
588 : }
589 1568 : if (!(x & 0xff000000)) {
590 750 : k += 8;
591 750 : x <<= 8;
592 : }
593 1568 : if (!(x & 0xf0000000)) {
594 864 : k += 4;
595 864 : x <<= 4;
596 : }
597 1568 : if (!(x & 0xc0000000)) {
598 776 : k += 2;
599 776 : x <<= 2;
600 : }
601 1568 : if (!(x & 0x80000000)) {
602 916 : k++;
603 916 : if (!(x & 0x40000000)) {
604 0 : return 32;
605 : }
606 : }
607 1568 : return k;
608 : }
609 :
610 : static int lo0bits(ULong *y)
611 130223 : {
612 : int k;
613 130223 : ULong x = *y;
614 :
615 130223 : if (x & 7) {
616 104502 : if (x & 1) {
617 59866 : return 0;
618 : }
619 44636 : if (x & 2) {
620 30686 : *y = x >> 1;
621 30686 : return 1;
622 : }
623 13950 : *y = x >> 2;
624 13950 : return 2;
625 : }
626 25721 : k = 0;
627 25721 : if (!(x & 0xffff)) {
628 6533 : k = 16;
629 6533 : x >>= 16;
630 : }
631 25721 : if (!(x & 0xff)) {
632 4323 : k += 8;
633 4323 : x >>= 8;
634 : }
635 25721 : if (!(x & 0xf)) {
636 11687 : k += 4;
637 11687 : x >>= 4;
638 : }
639 25721 : if (!(x & 0x3)) {
640 13740 : k += 2;
641 13740 : x >>= 2;
642 : }
643 25721 : if (!(x & 1)) {
644 14546 : k++;
645 14546 : x >>= 1;
646 14546 : if (!(x & 1)) {
647 0 : return 32;
648 : }
649 : }
650 25721 : *y = x;
651 25721 : return k;
652 : }
653 :
654 : static Bigint * i2b(int i)
655 2697 : {
656 : Bigint *b;
657 :
658 2697 : b = Balloc(1);
659 2697 : b->x[0] = i;
660 2697 : b->wds = 1;
661 2697 : return b;
662 : }
663 :
664 : static Bigint * mult(Bigint *a, Bigint *b)
665 1485 : {
666 : Bigint *c;
667 : int k, wa, wb, wc;
668 : ULong carry, y, z;
669 : ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
670 : #ifdef Pack_32
671 : ULong z2;
672 : #endif
673 :
674 1485 : if (a->wds < b->wds) {
675 325 : c = a;
676 325 : a = b;
677 325 : b = c;
678 : }
679 1485 : k = a->k;
680 1485 : wa = a->wds;
681 1485 : wb = b->wds;
682 1485 : wc = wa + wb;
683 1485 : if (wc > a->maxwds) {
684 549 : k++;
685 : }
686 1485 : c = Balloc(k);
687 6997 : for(x = c->x, xa = x + wc; x < xa; x++) {
688 5512 : *x = 0;
689 : }
690 1485 : xa = a->x;
691 1485 : xae = xa + wa;
692 1485 : xb = b->x;
693 1485 : xbe = xb + wb;
694 1485 : xc0 = c->x;
695 : #ifdef Pack_32
696 3713 : for(; xb < xbe; xb++, xc0++) {
697 2228 : if ((y = *xb & 0xffff)) {
698 2225 : x = xa;
699 2225 : xc = xc0;
700 2225 : carry = 0;
701 : do {
702 7140 : z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
703 7140 : carry = z >> 16;
704 7140 : z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
705 7140 : carry = z2 >> 16;
706 7140 : Storeinc(xc, z2, z);
707 : }
708 7140 : while(x < xae);
709 2225 : *xc = carry;
710 : }
711 2228 : if ((y = *xb >> 16)) {
712 1374 : x = xa;
713 1374 : xc = xc0;
714 1374 : carry = 0;
715 1374 : z2 = *xc;
716 : do {
717 5259 : z = (*x & 0xffff) * y + (*xc >> 16) + carry;
718 5259 : carry = z >> 16;
719 5259 : Storeinc(xc, z, z2);
720 5259 : z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
721 5259 : carry = z2 >> 16;
722 : }
723 5259 : while(x < xae);
724 1374 : *xc = z2;
725 : }
726 : }
727 : #else
728 : for(; xb < xbe; xc0++) {
729 : if (y = *xb++) {
730 : x = xa;
731 : xc = xc0;
732 : carry = 0;
733 : do {
734 : z = *x++ * y + *xc + carry;
735 : carry = z >> 16;
736 : *xc++ = z & 0xffff;
737 : }
738 : while(x < xae);
739 : *xc = carry;
740 : }
741 : }
742 : #endif
743 1485 : for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
744 1485 : c->wds = wc;
745 1485 : return c;
746 : }
747 :
748 : static Bigint * s2b (CONST char *s, int nd0, int nd, ULong y9)
749 1927 : {
750 : Bigint *b;
751 : int i, k;
752 : Long x, y;
753 :
754 1927 : x = (nd + 8) / 9;
755 1927 : for(k = 0, y = 1; x > y; y <<= 1, k++) ;
756 : #ifdef Pack_32
757 1927 : b = Balloc(k);
758 1927 : b->x[0] = y9;
759 1927 : b->wds = 1;
760 : #else
761 : b = Balloc(k+1);
762 : b->x[0] = y9 & 0xffff;
763 : b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
764 : #endif
765 :
766 1927 : i = 9;
767 1927 : if (9 < nd0) {
768 1745 : s += 9;
769 35967 : do b = multadd(b, 10, *s++ - '0');
770 35967 : while(++i < nd0);
771 1745 : s++;
772 : } else {
773 182 : s += 10;
774 : }
775 4045 : for(; i < nd; i++) {
776 2118 : b = multadd(b, 10, *s++ - '0');
777 : }
778 1927 : return b;
779 : }
780 :
781 : static Bigint * pow5mult(Bigint *b, int k)
782 632 : {
783 : Bigint *b1, *p5, *p51;
784 : int i;
785 : static int p05[3] = { 5, 25, 125 };
786 :
787 : _THREAD_PRIVATE_MUTEX_LOCK(pow5mult_mutex);
788 632 : if ((i = k & 3)) {
789 583 : b = multadd(b, p05[i-1], 0);
790 : }
791 :
792 632 : if (!(k >>= 2)) {
793 : _THREAD_PRIVATE_MUTEX_UNLOCK(pow5mult_mutex);
794 334 : return b;
795 : }
796 298 : if (!(p5 = p5s)) {
797 : /* first time */
798 298 : p5 = p5s = i2b(625);
799 298 : p5->next = 0;
800 : }
801 : for(;;) {
802 1096 : if (k & 1) {
803 532 : b1 = mult(b, p5);
804 532 : Bfree(b);
805 532 : b = b1;
806 : }
807 1096 : if (!(k >>= 1)) {
808 298 : break;
809 : }
810 798 : if (!(p51 = p5->next)) {
811 798 : if (!(p51 = p5->next)) {
812 798 : p51 = p5->next = mult(p5,p5);
813 798 : p51->next = 0;
814 : }
815 : }
816 798 : p5 = p51;
817 798 : }
818 : _THREAD_PRIVATE_MUTEX_UNLOCK(pow5mult_mutex);
819 298 : return b;
820 : }
821 :
822 :
823 : static Bigint *lshift(Bigint *b, int k)
824 5213 : {
825 : int i, k1, n, n1;
826 : Bigint *b1;
827 : ULong *x, *x1, *xe, z;
828 :
829 : #ifdef Pack_32
830 5213 : n = k >> 5;
831 : #else
832 : n = k >> 4;
833 : #endif
834 5213 : k1 = b->k;
835 5213 : n1 = n + b->wds + 1;
836 9921 : for(i = b->maxwds; n1 > i; i <<= 1) {
837 4708 : k1++;
838 : }
839 5213 : b1 = Balloc(k1);
840 5213 : x1 = b1->x;
841 9793 : for(i = 0; i < n; i++) {
842 4580 : *x1++ = 0;
843 : }
844 5213 : x = b->x;
845 5213 : xe = x + b->wds;
846 : #ifdef Pack_32
847 5213 : if (k &= 0x1f) {
848 5144 : k1 = 32 - k;
849 5144 : z = 0;
850 : do {
851 8879 : *x1++ = *x << k | z;
852 8879 : z = *x++ >> k1;
853 : }
854 8879 : while(x < xe);
855 5144 : if ((*x1 = z)) {
856 1129 : ++n1;
857 : }
858 : }
859 : #else
860 : if (k &= 0xf) {
861 : k1 = 16 - k;
862 : z = 0;
863 : do {
864 : *x1++ = *x << k & 0xffff | z;
865 : z = *x++ >> k1;
866 : }
867 : while(x < xe);
868 : if (*x1 = z) {
869 : ++n1;
870 : }
871 : }
872 : #endif
873 : else do
874 109 : *x1++ = *x++;
875 109 : while(x < xe);
876 5213 : b1->wds = n1 - 1;
877 5213 : Bfree(b);
878 5213 : return b1;
879 : }
880 :
881 : static int cmp(Bigint *a, Bigint *b)
882 14295 : {
883 : ULong *xa, *xa0, *xb, *xb0;
884 : int i, j;
885 :
886 14295 : i = a->wds;
887 14295 : j = b->wds;
888 : #ifdef DEBUG
889 : if (i > 1 && !a->x[i-1])
890 : Bug("cmp called with a->x[a->wds-1] == 0");
891 : if (j > 1 && !b->x[j-1])
892 : Bug("cmp called with b->x[b->wds-1] == 0");
893 : #endif
894 14295 : if (i -= j)
895 164 : return i;
896 14131 : xa0 = a->x;
897 14131 : xa = xa0 + j;
898 14131 : xb0 = b->x;
899 14131 : xb = xb0 + j;
900 : for(;;) {
901 17662 : if (*--xa != *--xb)
902 13873 : return *xa < *xb ? -1 : 1;
903 3789 : if (xa <= xa0)
904 258 : break;
905 3531 : }
906 258 : return 0;
907 : }
908 :
909 :
910 : static Bigint * diff(Bigint *a, Bigint *b)
911 1949 : {
912 : Bigint *c;
913 : int i, wa, wb;
914 : Long borrow, y; /* We need signed shifts here. */
915 : ULong *xa, *xae, *xb, *xbe, *xc;
916 : #ifdef Pack_32
917 : Long z;
918 : #endif
919 :
920 1949 : i = cmp(a,b);
921 1949 : if (!i) {
922 134 : c = Balloc(0);
923 134 : c->wds = 1;
924 134 : c->x[0] = 0;
925 134 : return c;
926 : }
927 1815 : if (i < 0) {
928 1358 : c = a;
929 1358 : a = b;
930 1358 : b = c;
931 1358 : i = 1;
932 : } else {
933 457 : i = 0;
934 : }
935 1815 : c = Balloc(a->k);
936 1815 : c->sign = i;
937 1815 : wa = a->wds;
938 1815 : xa = a->x;
939 1815 : xae = xa + wa;
940 1815 : wb = b->wds;
941 1815 : xb = b->x;
942 1815 : xbe = xb + wb;
943 1815 : xc = c->x;
944 1815 : borrow = 0;
945 : #ifdef Pack_32
946 : do {
947 6875 : y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
948 6875 : borrow = y >> 16;
949 : Sign_Extend(borrow, y);
950 6875 : z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
951 6875 : borrow = z >> 16;
952 : Sign_Extend(borrow, z);
953 6875 : Storeinc(xc, z, y);
954 6875 : } while(xb < xbe);
955 3652 : while(xa < xae) {
956 22 : y = (*xa & 0xffff) + borrow;
957 22 : borrow = y >> 16;
958 : Sign_Extend(borrow, y);
959 22 : z = (*xa++ >> 16) + borrow;
960 22 : borrow = z >> 16;
961 : Sign_Extend(borrow, z);
962 22 : Storeinc(xc, z, y);
963 : }
964 : #else
965 : do {
966 : y = *xa++ - *xb++ + borrow;
967 : borrow = y >> 16;
968 : Sign_Extend(borrow, y);
969 : *xc++ = y & 0xffff;
970 : } while(xb < xbe);
971 : while(xa < xae) {
972 : y = *xa++ + borrow;
973 : borrow = y >> 16;
974 : Sign_Extend(borrow, y);
975 : *xc++ = y & 0xffff;
976 : }
977 : #endif
978 6683 : while(!*--xc) {
979 3053 : wa--;
980 : }
981 1815 : c->wds = wa;
982 1815 : return c;
983 : }
984 :
985 : static double ulp (double _x)
986 746 : {
987 : volatile _double x;
988 : register Long L;
989 : volatile _double a;
990 :
991 746 : value(x) = _x;
992 746 : L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
993 : #ifndef Sudden_Underflow
994 746 : if (L > 0) {
995 : #endif
996 : #ifdef IBM
997 : L |= Exp_msk1 >> 4;
998 : #endif
999 746 : word0(a) = L;
1000 746 : word1(a) = 0;
1001 : #ifndef Sudden_Underflow
1002 : }
1003 : else {
1004 0 : L = -L >> Exp_shift;
1005 0 : if (L < Exp_shift) {
1006 0 : word0(a) = 0x80000 >> L;
1007 0 : word1(a) = 0;
1008 : }
1009 : else {
1010 0 : word0(a) = 0;
1011 0 : L -= Exp_shift;
1012 0 : word1(a) = L >= 31 ? 1 : 1 << (31 - L);
1013 : }
1014 : }
1015 : #endif
1016 746 : return value(a);
1017 : }
1018 :
1019 : static double
1020 : b2d
1021 : #ifdef KR_headers
1022 : (a, e) Bigint *a; int *e;
1023 : #else
1024 : (Bigint *a, int *e)
1025 : #endif
1026 1462 : {
1027 : ULong *xa, *xa0, w, y, z;
1028 : int k;
1029 : volatile _double d;
1030 : #ifdef VAX
1031 : ULong d0, d1;
1032 : #else
1033 : #define d0 word0(d)
1034 : #define d1 word1(d)
1035 : #endif
1036 :
1037 1462 : xa0 = a->x;
1038 1462 : xa = xa0 + a->wds;
1039 1462 : y = *--xa;
1040 : #ifdef DEBUG
1041 : if (!y) Bug("zero y in b2d");
1042 : #endif
1043 1462 : k = hi0bits(y);
1044 1462 : *e = 32 - k;
1045 : #ifdef Pack_32
1046 1462 : if (k < Ebits) {
1047 500 : d0 = Exp_1 | y >> (Ebits - k);
1048 500 : w = xa > xa0 ? *--xa : 0;
1049 500 : d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1050 500 : goto ret_d;
1051 : }
1052 962 : z = xa > xa0 ? *--xa : 0;
1053 962 : if (k -= Ebits) {
1054 912 : d0 = Exp_1 | y << k | z >> (32 - k);
1055 912 : y = xa > xa0 ? *--xa : 0;
1056 912 : d1 = z << k | y >> (32 - k);
1057 : }
1058 : else {
1059 50 : d0 = Exp_1 | y;
1060 50 : d1 = z;
1061 : }
1062 : #else
1063 : if (k < Ebits + 16) {
1064 : z = xa > xa0 ? *--xa : 0;
1065 : d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1066 : w = xa > xa0 ? *--xa : 0;
1067 : y = xa > xa0 ? *--xa : 0;
1068 : d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1069 : goto ret_d;
1070 : }
1071 : z = xa > xa0 ? *--xa : 0;
1072 : w = xa > xa0 ? *--xa : 0;
1073 : k -= Ebits + 16;
1074 : d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1075 : y = xa > xa0 ? *--xa : 0;
1076 : d1 = w << k + 16 | y << k;
1077 : #endif
1078 1462 : ret_d:
1079 : #ifdef VAX
1080 : word0(d) = d0 >> 16 | d0 << 16;
1081 : word1(d) = d1 >> 16 | d1 << 16;
1082 : #else
1083 : #undef d0
1084 : #undef d1
1085 : #endif
1086 1462 : return value(d);
1087 : }
1088 :
1089 :
1090 : static Bigint * d2b(double _d, int *e, int *bits)
1091 130223 : {
1092 : Bigint *b;
1093 : int de, i, k;
1094 : ULong *x, y, z;
1095 : volatile _double d;
1096 : #ifdef VAX
1097 : ULong d0, d1;
1098 : #endif
1099 :
1100 130223 : value(d) = _d;
1101 : #ifdef VAX
1102 : d0 = word0(d) >> 16 | word0(d) << 16;
1103 : d1 = word1(d) >> 16 | word1(d) << 16;
1104 : #else
1105 : #define d0 word0(d)
1106 : #define d1 word1(d)
1107 : #endif
1108 :
1109 : #ifdef Pack_32
1110 130223 : b = Balloc(1);
1111 : #else
1112 : b = Balloc(2);
1113 : #endif
1114 130223 : x = b->x;
1115 :
1116 130223 : z = d0 & Frac_mask;
1117 130223 : d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1118 : #ifdef Sudden_Underflow
1119 : de = (int)(d0 >> Exp_shift);
1120 : #ifndef IBM
1121 : z |= Exp_msk11;
1122 : #endif
1123 : #else
1124 130223 : if ((de = (int)(d0 >> Exp_shift)))
1125 130205 : z |= Exp_msk1;
1126 : #endif
1127 : #ifdef Pack_32
1128 130223 : if ((y = d1)) {
1129 121750 : if ((k = lo0bits(&y))) {
1130 61956 : x[0] = y | (z << (32 - k));
1131 61956 : z >>= k;
1132 : } else {
1133 59794 : x[0] = y;
1134 : }
1135 121750 : i = b->wds = (x[1] = z) ? 2 : 1;
1136 : } else {
1137 : #ifdef DEBUG
1138 : if (!z)
1139 : Bug("Zero passed to d2b");
1140 : #endif
1141 8473 : k = lo0bits(&z);
1142 8473 : x[0] = z;
1143 8473 : i = b->wds = 1;
1144 8473 : k += 32;
1145 : }
1146 : #else
1147 : if (y = d1) {
1148 : if (k = lo0bits(&y)) {
1149 : if (k >= 16) {
1150 : x[0] = y | z << 32 - k & 0xffff;
1151 : x[1] = z >> k - 16 & 0xffff;
1152 : x[2] = z >> k;
1153 : i = 2;
1154 : } else {
1155 : x[0] = y & 0xffff;
1156 : x[1] = y >> 16 | z << 16 - k & 0xffff;
1157 : x[2] = z >> k & 0xffff;
1158 : x[3] = z >> k+16;
1159 : i = 3;
1160 : }
1161 : } else {
1162 : x[0] = y & 0xffff;
1163 : x[1] = y >> 16;
1164 : x[2] = z & 0xffff;
1165 : x[3] = z >> 16;
1166 : i = 3;
1167 : }
1168 : } else {
1169 : #ifdef DEBUG
1170 : if (!z)
1171 : Bug("Zero passed to d2b");
1172 : #endif
1173 : k = lo0bits(&z);
1174 : if (k >= 16) {
1175 : x[0] = z;
1176 : i = 0;
1177 : } else {
1178 : x[0] = z & 0xffff;
1179 : x[1] = z >> 16;
1180 : i = 1;
1181 : }
1182 : k += 32;
1183 : }
1184 : while(!x[i])
1185 : --i;
1186 : b->wds = i + 1;
1187 : #endif
1188 : #ifndef Sudden_Underflow
1189 130223 : if (de) {
1190 : #endif
1191 : #ifdef IBM
1192 : *e = (de - Bias - (P-1) << 2) + k;
1193 : *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1194 : #else
1195 130205 : *e = de - Bias - (P-1) + k;
1196 130205 : *bits = P - k;
1197 : #endif
1198 : #ifndef Sudden_Underflow
1199 : } else {
1200 18 : *e = de - Bias - (P-1) + 1 + k;
1201 : #ifdef Pack_32
1202 18 : *bits = 32*i - hi0bits(x[i-1]);
1203 : #else
1204 : *bits = (i+2)*16 - hi0bits(x[i]);
1205 : #endif
1206 : }
1207 : #endif
1208 130223 : return b;
1209 : }
1210 : #undef d0
1211 : #undef d1
1212 :
1213 :
1214 : static double ratio (Bigint *a, Bigint *b)
1215 731 : {
1216 : volatile _double da, db;
1217 : int k, ka, kb;
1218 :
1219 731 : value(da) = b2d(a, &ka);
1220 731 : value(db) = b2d(b, &kb);
1221 : #ifdef Pack_32
1222 731 : k = ka - kb + 32*(a->wds - b->wds);
1223 : #else
1224 : k = ka - kb + 16*(a->wds - b->wds);
1225 : #endif
1226 : #ifdef IBM
1227 : if (k > 0) {
1228 : word0(da) += (k >> 2)*Exp_msk1;
1229 : if (k &= 3) {
1230 : da *= 1 << k;
1231 : }
1232 : } else {
1233 : k = -k;
1234 : word0(db) += (k >> 2)*Exp_msk1;
1235 : if (k &= 3)
1236 : db *= 1 << k;
1237 : }
1238 : #else
1239 731 : if (k > 0) {
1240 225 : word0(da) += k*Exp_msk1;
1241 : } else {
1242 506 : k = -k;
1243 506 : word0(db) += k*Exp_msk1;
1244 : }
1245 : #endif
1246 731 : return value(da) / value(db);
1247 : }
1248 :
1249 : static CONST double
1250 : tens[] = {
1251 : 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1252 : 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1253 : 1e20, 1e21, 1e22
1254 : #ifdef VAX
1255 : , 1e23, 1e24
1256 : #endif
1257 : };
1258 :
1259 : #ifdef IEEE_Arith
1260 : static CONST double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1261 : static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1262 : #define n_bigtens 5
1263 : #else
1264 : #ifdef IBM
1265 : static CONST double bigtens[] = { 1e16, 1e32, 1e64 };
1266 : static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1267 : #define n_bigtens 3
1268 : #else
1269 : static CONST double bigtens[] = { 1e16, 1e32 };
1270 : static CONST double tinytens[] = { 1e-16, 1e-32 };
1271 : #define n_bigtens 2
1272 : #endif
1273 : #endif
1274 :
1275 :
1276 : static int quorem(Bigint *b, Bigint *S)
1277 12672 : {
1278 : int n;
1279 : Long borrow, y;
1280 : ULong carry, q, ys;
1281 : ULong *bx, *bxe, *sx, *sxe;
1282 : #ifdef Pack_32
1283 : Long z;
1284 : ULong si, zs;
1285 : #endif
1286 :
1287 12672 : n = S->wds;
1288 : #ifdef DEBUG
1289 : /*debug*/ if (b->wds > n)
1290 : /*debug*/ Bug("oversize b in quorem");
1291 : #endif
1292 12672 : if (b->wds < n)
1293 3133 : return 0;
1294 9539 : sx = S->x;
1295 9539 : sxe = sx + --n;
1296 9539 : bx = b->x;
1297 9539 : bxe = bx + n;
1298 9539 : q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1299 : #ifdef DEBUG
1300 : /*debug*/ if (q > 9)
1301 : /*debug*/ Bug("oversized quotient in quorem");
1302 : #endif
1303 9539 : if (q) {
1304 8360 : borrow = 0;
1305 8360 : carry = 0;
1306 : do {
1307 : #ifdef Pack_32
1308 32892 : si = *sx++;
1309 32892 : ys = (si & 0xffff) * q + carry;
1310 32892 : zs = (si >> 16) * q + (ys >> 16);
1311 32892 : carry = zs >> 16;
1312 32892 : y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1313 32892 : borrow = y >> 16;
1314 : Sign_Extend(borrow, y);
1315 32892 : z = (*bx >> 16) - (zs & 0xffff) + borrow;
1316 32892 : borrow = z >> 16;
1317 : Sign_Extend(borrow, z);
1318 32892 : Storeinc(bx, z, y);
1319 : #else
1320 : ys = *sx++ * q + carry;
1321 : carry = ys >> 16;
1322 : y = *bx - (ys & 0xffff) + borrow;
1323 : borrow = y >> 16;
1324 : Sign_Extend(borrow, y);
1325 : *bx++ = y & 0xffff;
1326 : #endif
1327 : }
1328 32892 : while(sx <= sxe);
1329 8360 : if (!*bxe) {
1330 0 : bx = b->x;
1331 0 : while(--bxe > bx && !*bxe)
1332 0 : --n;
1333 0 : b->wds = n;
1334 : }
1335 : }
1336 9539 : if (cmp(b, S) >= 0) {
1337 141 : q++;
1338 141 : borrow = 0;
1339 141 : carry = 0;
1340 141 : bx = b->x;
1341 141 : sx = S->x;
1342 : do {
1343 : #ifdef Pack_32
1344 508 : si = *sx++;
1345 508 : ys = (si & 0xffff) + carry;
1346 508 : zs = (si >> 16) + (ys >> 16);
1347 508 : carry = zs >> 16;
1348 508 : y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1349 508 : borrow = y >> 16;
1350 : Sign_Extend(borrow, y);
1351 508 : z = (*bx >> 16) - (zs & 0xffff) + borrow;
1352 508 : borrow = z >> 16;
1353 : Sign_Extend(borrow, z);
1354 508 : Storeinc(bx, z, y);
1355 : #else
1356 : ys = *sx++ + carry;
1357 : carry = ys >> 16;
1358 : y = *bx - (ys & 0xffff) + borrow;
1359 : borrow = y >> 16;
1360 : Sign_Extend(borrow, y);
1361 : *bx++ = y & 0xffff;
1362 : #endif
1363 : }
1364 508 : while(sx <= sxe);
1365 141 : bx = b->x;
1366 141 : bxe = bx + n;
1367 141 : if (!*bxe) {
1368 449 : while(--bxe > bx && !*bxe)
1369 171 : --n;
1370 139 : b->wds = n;
1371 : }
1372 : }
1373 9539 : return q;
1374 : }
1375 :
1376 : static void destroy_freelist(void)
1377 17665 : {
1378 : int i;
1379 : Bigint *tmp;
1380 :
1381 : _THREAD_PRIVATE_MUTEX_LOCK(dtoa_mutex);
1382 300305 : for (i = 0; i <= Kmax; i++) {
1383 282640 : Bigint **listp = &freelist[i];
1384 568433 : while ((tmp = *listp) != NULL) {
1385 3153 : *listp = tmp->next;
1386 3153 : free(tmp);
1387 : }
1388 282640 : freelist[i] = NULL;
1389 : }
1390 : _THREAD_PRIVATE_MUTEX_UNLOCK(dtoa_mutex);
1391 :
1392 17665 : }
1393 :
1394 :
1395 : ZEND_API void zend_freedtoa(char *s)
1396 128895 : {
1397 128895 : Bigint *b = (Bigint *)((int *)s - 1);
1398 128895 : b->maxwds = 1 << (b->k = *(int*)b);
1399 128895 : Bfree(b);
1400 128895 : }
1401 :
1402 : /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1403 : *
1404 : * Inspired by "How to Print Floating-Point Numbers Accurately" by
1405 : * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1406 : *
1407 : * Modifications:
1408 : * 1. Rather than iterating, we use a simple numeric overestimate
1409 : * to determine k = floor(log10(d)). We scale relevant
1410 : * quantities using O(log2(k)) rather than O(k) multiplications.
1411 : * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1412 : * try to generate digits strictly left to right. Instead, we
1413 : * compute with fewer bits and propagate the carry if necessary
1414 : * when rounding the final digit up. This is often faster.
1415 : * 3. Under the assumption that input will be rounded nearest,
1416 : * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1417 : * That is, we allow equality in stopping tests when the
1418 : * round-nearest rule will give the same floating-point value
1419 : * as would satisfaction of the stopping test with strict
1420 : * inequality.
1421 : * 4. We remove common factors of powers of 2 from relevant
1422 : * quantities.
1423 : * 5. When converting floating-point integers less than 1e16,
1424 : * we use floating-point arithmetic rather than resorting
1425 : * to multiple-precision integers.
1426 : * 6. When asked to produce fewer than 15 digits, we first try
1427 : * to get by with floating-point arithmetic; we resort to
1428 : * multiple-precision integer arithmetic only if we cannot
1429 : * guarantee that the floating-point calculation has given
1430 : * the correctly rounded result. For k requested digits and
1431 : * "uniformly" distributed input, the probability is
1432 : * something like 10^(k-15) that we must resort to the Long
1433 : * calculation.
1434 : */
1435 :
1436 : ZEND_API char * zend_dtoa(double _d, int mode, int ndigits, int *decpt, int *sign, char **rve)
1437 128895 : {
1438 : /* Arguments ndigits, decpt, sign are similar to those
1439 : of ecvt and fcvt; trailing zeros are suppressed from
1440 : the returned string. If not null, *rve is set to point
1441 : to the end of the return value. If d is +-Infinity or NaN,
1442 : then *decpt is set to 9999.
1443 :
1444 : mode:
1445 : 0 ==> shortest string that yields d when read in
1446 : and rounded to nearest.
1447 : 1 ==> like 0, but with Steele & White stopping rule;
1448 : e.g. with IEEE P754 arithmetic , mode 0 gives
1449 : 1e23 whereas mode 1 gives 9.999999999999999e22.
1450 : 2 ==> max(1,ndigits) significant digits. This gives a
1451 : return value similar to that of ecvt, except
1452 : that trailing zeros are suppressed.
1453 : 3 ==> through ndigits past the decimal point. This
1454 : gives a return value similar to that from fcvt,
1455 : except that trailing zeros are suppressed, and
1456 : ndigits can be negative.
1457 : 4-9 should give the same return values as 2-3, i.e.,
1458 : 4 <= mode <= 9 ==> same return as mode
1459 : 2 + (mode & 1). These modes are mainly for
1460 : debugging; often they run slower but sometimes
1461 : faster than modes 2-3.
1462 : 4,5,8,9 ==> left-to-right digit generation.
1463 : 6-9 ==> don't try fast floating-point estimate
1464 : (if applicable).
1465 :
1466 : Values of mode other than 0-9 are treated as mode 0.
1467 :
1468 : Sufficient space is allocated to the return value
1469 : to hold the suppressed trailing zeros.
1470 : */
1471 :
1472 128895 : int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1,
1473 : j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1474 128895 : spec_case = 0, try_quick;
1475 : Long L;
1476 : #ifndef Sudden_Underflow
1477 : int denorm;
1478 : ULong x;
1479 : #endif
1480 : Bigint *b, *b1, *delta, *mlo, *mhi, *S, *tmp;
1481 : double ds;
1482 : char *s, *s0;
1483 : volatile _double d, d2, eps;
1484 :
1485 128895 : value(d) = _d;
1486 :
1487 128895 : if (word0(d) & Sign_bit) {
1488 : /* set sign for everything, including 0's and NaNs */
1489 3285 : *sign = 1;
1490 3285 : word0(d) &= ~Sign_bit; /* clear sign bit */
1491 : }
1492 : else
1493 125610 : *sign = 0;
1494 :
1495 : #if defined(IEEE_Arith) + defined(VAX)
1496 : #ifdef IEEE_Arith
1497 128895 : if ((word0(d) & Exp_mask) == Exp_mask)
1498 : #else
1499 : if (word0(d) == 0x8000)
1500 : #endif
1501 : {
1502 : /* Infinity or NaN */
1503 12 : *decpt = 9999;
1504 : #ifdef IEEE_Arith
1505 12 : if (!word1(d) && !(word0(d) & 0xfffff))
1506 8 : return nrv_alloc("Infinity", rve, 8);
1507 : #endif
1508 4 : return nrv_alloc("NaN", rve, 3);
1509 : }
1510 : #endif
1511 : #ifdef IBM
1512 : value(d) += 0; /* normalize */
1513 : #endif
1514 128883 : if (!value(d)) {
1515 609 : *decpt = 1;
1516 609 : return nrv_alloc("0", rve, 1);
1517 : }
1518 :
1519 128274 : b = d2b(value(d), &be, &bbits);
1520 : #ifdef Sudden_Underflow
1521 : i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1522 : #else
1523 128274 : if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
1524 : #endif
1525 128262 : value(d2) = value(d);
1526 128262 : word0(d2) &= Frac_mask1;
1527 128262 : word0(d2) |= Exp_11;
1528 : #ifdef IBM
1529 : if (j = 11 - hi0bits(word0(d2) & Frac_mask))
1530 : value(d2) /= 1 << j;
1531 : #endif
1532 :
1533 : /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1534 : * log10(x) = log(x) / log(10)
1535 : * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1536 : * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1537 : *
1538 : * This suggests computing an approximation k to log10(d) by
1539 : *
1540 : * k = (i - Bias)*0.301029995663981
1541 : * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1542 : *
1543 : * We want k to be too large rather than too small.
1544 : * The error in the first-order Taylor series approximation
1545 : * is in our favor, so we just round up the constant enough
1546 : * to compensate for any error in the multiplication of
1547 : * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1548 : * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1549 : * adding 1e-13 to the constant term more than suffices.
1550 : * Hence we adjust the constant term to 0.1760912590558.
1551 : * (We could get a more accurate k by invoking log10,
1552 : * but this is probably not worthwhile.)
1553 : */
1554 :
1555 128262 : i -= Bias;
1556 : #ifdef IBM
1557 : i <<= 2;
1558 : i += j;
1559 : #endif
1560 : #ifndef Sudden_Underflow
1561 128262 : denorm = 0;
1562 : }
1563 : else {
1564 : /* d is denormalized */
1565 :
1566 12 : i = bbits + be + (Bias + (P-1) - 1);
1567 12 : x = i > 32 ? (word0(d) << (64 - i)) | (word1(d) >> (i - 32))
1568 : : (word1(d) << (32 - i));
1569 12 : value(d2) = x;
1570 12 : word0(d2) -= 31*Exp_msk1; /* adjust exponent */
1571 12 : i -= (Bias + (P-1) - 1) + 1;
1572 12 : denorm = 1;
1573 : }
1574 : #endif
1575 128274 : ds = (value(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1576 128274 : k = (int)ds;
1577 128274 : if (ds < 0. && ds != k)
1578 100635 : k--; /* want k = floor(ds) */
1579 128274 : k_check = 1;
1580 128274 : if (k >= 0 && k <= Ten_pmax) {
1581 27521 : if (value(d) < tens[k])
1582 3554 : k--;
1583 27521 : k_check = 0;
1584 : }
1585 128274 : j = bbits - i - 1;
1586 128274 : if (j >= 0) {
1587 124042 : b2 = 0;
1588 124042 : s2 = j;
1589 : }
1590 : else {
1591 4232 : b2 = -j;
1592 4232 : s2 = 0;
1593 : }
1594 128274 : if (k >= 0) {
1595 24198 : b5 = 0;
1596 24198 : s5 = k;
1597 24198 : s2 += k;
1598 : }
1599 : else {
1600 104076 : b2 -= k;
1601 104076 : b5 = -k;
1602 104076 : s5 = 0;
1603 : }
1604 128274 : if (mode < 0 || mode > 9)
1605 0 : mode = 0;
1606 128274 : try_quick = 1;
1607 128274 : if (mode > 5) {
1608 0 : mode -= 4;
1609 0 : try_quick = 0;
1610 : }
1611 128274 : leftright = 1;
1612 128274 : switch(mode) {
1613 : case 0:
1614 : case 1:
1615 0 : ilim = ilim1 = -1;
1616 0 : i = 18;
1617 0 : ndigits = 0;
1618 0 : break;
1619 : case 2:
1620 24099 : leftright = 0;
1621 : /* no break */
1622 : case 4:
1623 24099 : if (ndigits <= 0)
1624 0 : ndigits = 1;
1625 24099 : ilim = ilim1 = i = ndigits;
1626 24099 : break;
1627 : case 3:
1628 104175 : leftright = 0;
1629 : /* no break */
1630 : case 5:
1631 104175 : i = ndigits + k + 1;
1632 104175 : ilim = i;
1633 104175 : ilim1 = i - 1;
1634 104175 : if (i <= 0)
1635 43 : i = 1;
1636 : }
1637 128274 : s = s0 = rv_alloc(i);
1638 :
1639 128274 : if (ilim >= 0 && ilim <= Quick_max && try_quick) {
1640 :
1641 : /* Try to get by with floating-point arithmetic. */
1642 :
1643 128002 : i = 0;
1644 128002 : value(d2) = value(d);
1645 128002 : k0 = k;
1646 128002 : ilim0 = ilim;
1647 128002 : ieps = 2; /* conservative */
1648 128002 : if (k > 0) {
1649 10064 : ds = tens[k&0xf];
1650 10064 : j = k >> 4;
1651 10064 : if (j & Bletch) {
1652 : /* prevent overflows */
1653 4 : j &= Bletch - 1;
1654 4 : value(d) /= bigtens[n_bigtens-1];
1655 4 : ieps++;
1656 : }
1657 10255 : for(; j; j >>= 1, i++)
1658 191 : if (j & 1) {
1659 126 : ieps++;
1660 126 : ds *= bigtens[i];
1661 : }
1662 10064 : value(d) /= ds;
1663 : }
1664 117938 : else if ((j1 = -k)) {
1665 104022 : value(d) *= tens[j1 & 0xf];
1666 104222 : for(j = j1 >> 4; j; j >>= 1, i++)
1667 200 : if (j & 1) {
1668 116 : ieps++;
1669 116 : value(d) *= bigtens[i];
1670 : }
1671 : }
1672 128002 : if (k_check && value(d) < 1. && ilim > 0) {
1673 37 : if (ilim1 <= 0)
1674 0 : goto fast_failed;
1675 37 : ilim = ilim1;
1676 37 : k--;
1677 37 : value(d) *= 10.;
1678 37 : ieps++;
1679 : }
1680 128002 : value(eps) = ieps*value(d) + 7.;
1681 128002 : word0(eps) -= (P-1)*Exp_msk1;
1682 128002 : if (ilim == 0) {
1683 8 : S = mhi = 0;
1684 8 : value(d) -= 5.;
1685 8 : if (value(d) > value(eps))
1686 2 : goto one_digit;
1687 6 : if (value(d) < -value(eps))
1688 6 : goto no_digits;
1689 0 : goto fast_failed;
1690 : }
1691 : #ifndef No_leftright
1692 127994 : if (leftright) {
1693 : /* Use Steele & White method of only
1694 : * generating digits needed.
1695 : */
1696 0 : value(eps) = 0.5/tens[ilim-1] - value(eps);
1697 0 : for(i = 0;;) {
1698 0 : L = value(d);
1699 0 : value(d) -= L;
1700 0 : *s++ = '0' + (int)L;
1701 0 : if (value(d) < value(eps))
1702 0 : goto ret1;
1703 0 : if (1. - value(d) < value(eps))
1704 0 : goto bump_up;
1705 0 : if (++i >= ilim)
1706 0 : break;
1707 0 : value(eps) *= 10.;
1708 0 : value(d) *= 10.;
1709 0 : }
1710 : }
1711 : else {
1712 : #endif
1713 : /* Generate ilim digits, then fix them up. */
1714 127994 : value(eps) *= tens[ilim-1];
1715 1145669 : for(i = 1;; i++, value(d) *= 10.) {
1716 1145669 : L = value(d);
1717 1145669 : value(d) -= L;
1718 1145669 : *s++ = '0' + (int)L;
1719 1145669 : if (i == ilim) {
1720 127994 : if (value(d) > 0.5 + value(eps))
1721 59907 : goto bump_up;
1722 68087 : else if (value(d) < 0.5 - value(eps)) {
1723 276969 : while(*--s == '0');
1724 67764 : s++;
1725 67764 : goto ret1;
1726 : }
1727 323 : break;
1728 : }
1729 1017675 : }
1730 : #ifndef No_leftright
1731 : }
1732 : #endif
1733 323 : fast_failed:
1734 323 : s = s0;
1735 323 : value(d) = value(d2);
1736 323 : k = k0;
1737 323 : ilim = ilim0;
1738 : }
1739 :
1740 : /* Do we have a "small" integer? */
1741 :
1742 595 : if (be >= 0 && k <= Int_max) {
1743 : /* Yes. */
1744 145 : ds = tens[k];
1745 145 : if (ndigits < 0 && ilim <= 0) {
1746 0 : S = mhi = 0;
1747 0 : if (ilim < 0 || value(d) <= 5*ds)
1748 : goto no_digits;
1749 0 : goto one_digit;
1750 : }
1751 1398 : for(i = 1;; i++) {
1752 1398 : L = value(d) / ds;
1753 1398 : value(d) -= L*ds;
1754 : #ifdef Check_FLT_ROUNDS
1755 : /* If FLT_ROUNDS == 2, L will usually be high by 1 */
1756 : if (value(d) < 0) {
1757 : L--;
1758 : value(d) += ds;
1759 : }
1760 : #endif
1761 1398 : *s++ = '0' + (int)L;
1762 1398 : if (i == ilim) {
1763 0 : value(d) += value(d);
1764 0 : if (value(d) > ds || (value(d) == ds && (L & 1))) {
1765 59907 : bump_up:
1766 256667 : while(*--s == '9')
1767 136862 : if (s == s0) {
1768 9 : k++;
1769 9 : *s = '0';
1770 9 : break;
1771 : }
1772 59907 : ++*s++;
1773 : }
1774 59907 : break;
1775 : }
1776 1398 : if (!(value(d) *= 10.))
1777 145 : break;
1778 1253 : }
1779 60052 : goto ret1;
1780 : }
1781 :
1782 450 : m2 = b2;
1783 450 : m5 = b5;
1784 450 : mhi = mlo = 0;
1785 450 : if (leftright) {
1786 0 : if (mode < 2) {
1787 0 : i =
1788 : #ifndef Sudden_Underflow
1789 : denorm ? be + (Bias + (P-1) - 1 + 1) :
1790 : #endif
1791 : #ifdef IBM
1792 : 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
1793 : #else
1794 : 1 + P - bbits;
1795 : #endif
1796 : }
1797 : else {
1798 0 : j = ilim - 1;
1799 0 : if (m5 >= j)
1800 0 : m5 -= j;
1801 : else {
1802 0 : s5 += j -= m5;
1803 0 : b5 += j;
1804 0 : m5 = 0;
1805 : }
1806 0 : if ((i = ilim) < 0) {
1807 0 : m2 -= i;
1808 0 : i = 0;
1809 : }
1810 : }
1811 0 : b2 += i;
1812 0 : s2 += i;
1813 0 : mhi = i2b(1);
1814 : }
1815 450 : if (m2 > 0 && s2 > 0) {
1816 382 : i = m2 < s2 ? m2 : s2;
1817 382 : b2 -= i;
1818 382 : m2 -= i;
1819 382 : s2 -= i;
1820 : }
1821 450 : if (b5 > 0) {
1822 339 : if (leftright) {
1823 0 : if (m5 > 0) {
1824 0 : mhi = pow5mult(mhi, m5);
1825 0 : b1 = mult(mhi, b);
1826 0 : Bfree(b);
1827 0 : b = b1;
1828 : }
1829 0 : if ((j = b5 - m5)) {
1830 0 : b = pow5mult(b, j);
1831 : }
1832 : } else {
1833 339 : b = pow5mult(b, b5);
1834 : }
1835 : }
1836 450 : S = i2b(1);
1837 450 : if (s5 > 0)
1838 88 : S = pow5mult(S, s5);
1839 : /* Check for special case that d is a normalized power of 2. */
1840 :
1841 450 : if (mode < 2) {
1842 0 : if (!word1(d) && !(word0(d) & Bndry_mask)
1843 : #ifndef Sudden_Underflow
1844 : && word0(d) & Exp_mask
1845 : #endif
1846 : ) {
1847 : /* The special case */
1848 0 : b2 += Log2P;
1849 0 : s2 += Log2P;
1850 0 : spec_case = 1;
1851 : } else {
1852 0 : spec_case = 0;
1853 : }
1854 : }
1855 :
1856 : /* Arrange for convenient computation of quotients:
1857 : * shift left if necessary so divisor has 4 leading 0 bits.
1858 : *
1859 : * Perhaps we should just compute leading 28 bits of S once
1860 : * and for all and pass them and a shift to quorem, so it
1861 : * can do shifts and ors to compute the numerator for q.
1862 : */
1863 : #ifdef Pack_32
1864 450 : if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
1865 450 : i = 32 - i;
1866 : #else
1867 : if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf))
1868 : i = 16 - i;
1869 : #endif
1870 450 : if (i > 4) {
1871 433 : i -= 4;
1872 433 : b2 += i;
1873 433 : m2 += i;
1874 433 : s2 += i;
1875 : }
1876 17 : else if (i < 4) {
1877 17 : i += 28;
1878 17 : b2 += i;
1879 17 : m2 += i;
1880 17 : s2 += i;
1881 : }
1882 450 : if (b2 > 0)
1883 450 : b = lshift(b, b2);
1884 450 : if (s2 > 0)
1885 450 : S = lshift(S, s2);
1886 450 : if (k_check) {
1887 379 : if (cmp(b,S) < 0) {
1888 1 : k--;
1889 1 : b = multadd(b, 10, 0); /* we botched the k estimate */
1890 1 : if (leftright)
1891 0 : mhi = multadd(mhi, 10, 0);
1892 1 : ilim = ilim1;
1893 : }
1894 : }
1895 450 : if (ilim <= 0 && mode > 2) {
1896 35 : if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
1897 : /* no digits, fcvt style */
1898 41 : no_digits:
1899 41 : k = -1 - ndigits;
1900 41 : goto ret;
1901 : }
1902 2 : one_digit:
1903 2 : *s++ = '1';
1904 2 : k++;
1905 2 : goto ret;
1906 : }
1907 415 : if (leftright) {
1908 0 : if (m2 > 0)
1909 0 : mhi = lshift(mhi, m2);
1910 :
1911 : /* Compute mlo -- check for special case
1912 : * that d is a normalized power of 2.
1913 : */
1914 :
1915 0 : mlo = mhi;
1916 0 : if (spec_case) {
1917 0 : mhi = Balloc(mhi->k);
1918 0 : Bcopy(mhi, mlo);
1919 0 : mhi = lshift(mhi, Log2P);
1920 : }
1921 :
1922 0 : for(i = 1;;i++) {
1923 0 : dig = quorem(b,S) + '0';
1924 : /* Do we yet have the shortest decimal string
1925 : * that will round to d?
1926 : */
1927 0 : j = cmp(b, mlo);
1928 0 : delta = diff(S, mhi);
1929 0 : j1 = delta->sign ? 1 : cmp(b, delta);
1930 0 : Bfree(delta);
1931 : #ifndef ROUND_BIASED
1932 0 : if (j1 == 0 && !mode && !(word1(d) & 1)) {
1933 0 : if (dig == '9')
1934 0 : goto round_9_up;
1935 0 : if (j > 0)
1936 0 : dig++;
1937 0 : *s++ = dig;
1938 0 : goto ret;
1939 : }
1940 : #endif
1941 0 : if (j < 0 || (j == 0 && !mode
1942 : #ifndef ROUND_BIASED
1943 : && !(word1(d) & 1)
1944 : #endif
1945 : )) {
1946 0 : if (j1 > 0) {
1947 0 : b = lshift(b, 1);
1948 0 : j1 = cmp(b, S);
1949 0 : if ((j1 > 0 || (j1 == 0 && (dig & 1)))
1950 : && dig++ == '9')
1951 0 : goto round_9_up;
1952 : }
1953 0 : *s++ = dig;
1954 0 : goto ret;
1955 : }
1956 0 : if (j1 > 0) {
1957 0 : if (dig == '9') { /* possible if i == 1 */
1958 0 : round_9_up:
1959 0 : *s++ = '9';
1960 0 : goto roundoff;
1961 : }
1962 0 : *s++ = dig + 1;
1963 0 : goto ret;
1964 : }
1965 0 : *s++ = dig;
1966 0 : if (i == ilim)
1967 0 : break;
1968 0 : b = multadd(b, 10, 0);
1969 0 : if (mlo == mhi)
1970 0 : mlo = mhi = multadd(mhi, 10, 0);
1971 : else {
1972 0 : mlo = multadd(mlo, 10, 0);
1973 0 : mhi = multadd(mhi, 10, 0);
1974 : }
1975 0 : }
1976 : }
1977 : else
1978 12672 : for(i = 1;; i++) {
1979 12672 : *s++ = dig = quorem(b,S) + '0';
1980 12672 : if (i >= ilim)
1981 415 : break;
1982 12257 : b = multadd(b, 10, 0);
1983 12257 : }
1984 :
1985 : /* Round off last digit */
1986 :
1987 415 : b = lshift(b, 1);
1988 415 : j = cmp(b, S);
1989 728 : if (j > 0 || (j == 0 && (dig & 1))) {
1990 313 : roundoff:
1991 627 : while(*--s == '9')
1992 1 : if (s == s0) {
1993 0 : k++;
1994 0 : *s++ = '1';
1995 0 : goto ret;
1996 : }
1997 313 : ++*s++;
1998 : }
1999 : else {
2000 2948 : while(*--s == '0');
2001 102 : s++;
2002 : }
2003 458 : ret:
2004 458 : Bfree(S);
2005 458 : if (mhi) {
2006 0 : if (mlo && mlo != mhi)
2007 0 : Bfree(mlo);
2008 0 : Bfree(mhi);
2009 : }
2010 128274 : ret1:
2011 :
2012 : _THREAD_PRIVATE_MUTEX_LOCK(pow5mult_mutex);
2013 256881 : while (p5s) {
2014 333 : tmp = p5s;
2015 333 : p5s = p5s->next;
2016 333 : free(tmp);
2017 : }
2018 : _THREAD_PRIVATE_MUTEX_UNLOCK(pow5mult_mutex);
2019 :
2020 128274 : Bfree(b);
2021 :
2022 128274 : if (s == s0) { /* don't return empty string */
2023 41 : *s++ = '0';
2024 41 : k = 0;
2025 : }
2026 128274 : *s = 0;
2027 128274 : *decpt = k + 1;
2028 128274 : if (rve)
2029 104915 : *rve = s;
2030 128274 : return s0;
2031 : }
2032 :
2033 : ZEND_API double zend_strtod (CONST char *s00, char **se)
2034 214894 : {
2035 : int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
2036 : e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
2037 : CONST char *s, *s0, *s1;
2038 : double aadj, aadj1, adj;
2039 : volatile _double rv, rv0;
2040 : Long L;
2041 : ULong y, z;
2042 : Bigint *bb, *bb1, *bd, *bd0, *bs, *delta, *tmp;
2043 : double result;
2044 :
2045 214894 : CONST char decimal_point = '.';
2046 :
2047 214894 : sign = nz0 = nz = 0;
2048 214894 : value(rv) = 0.;
2049 :
2050 :
2051 214894 : for(s = s00; isspace((unsigned char) *s); s++)
2052 : ;
2053 :
2054 214894 : if (*s == '-') {
2055 981 : sign = 1;
2056 981 : s++;
2057 213913 : } else if (*s == '+') {
2058 20 : s++;
2059 : }
2060 :
2061 214894 : if (*s == '\0') {
2062 82 : s = s00;
2063 82 : goto ret;
2064 : }
2065 :
2066 214812 : if (*s == '0') {
2067 201192 : nz0 = 1;
2068 201869 : while(*++s == '0') ;
2069 201192 : if (!*s)
2070 112 : goto ret;
2071 : }
2072 214700 : s0 = s;
2073 214700 : y = z = 0;
2074 303533 : for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
2075 88833 : if (nd < 9)
2076 49112 : y = 10*y + c - '0';
2077 39721 : else if (nd < 16)
2078 16120 : z = 10*z + c - '0';
2079 214700 : nd0 = nd;
2080 214700 : if (c == decimal_point) {
2081 209610 : c = *++s;
2082 209610 : if (!nd) {
2083 224132 : for(; c == '0'; c = *++s)
2084 22339 : nz++;
2085 201793 : if (c > '0' && c <= '9') {
2086 201634 : s0 = s;
2087 201634 : nf += nz;
2088 201634 : nz = 0;
2089 201634 : goto have_dig;
2090 : }
2091 159 : goto dig_done;
2092 : }
2093 1613574 : for(; c >= '0' && c <= '9'; c = *++s) {
2094 1605757 : have_dig:
2095 1605757 : nz++;
2096 1605757 : if (c -= '0') {
2097 1103294 : nf += nz;
2098 1180299 : for(i = 1; i < nz; i++)
2099 77005 : if (nd++ < 9)
2100 76246 : y *= 10;
2101 759 : else if (nd <= DBL_DIG + 1)
2102 511 : z *= 10;
2103 1103294 : if (nd++ < 9)
2104 1100620 : y = 10*y + c;
2105 2674 : else if (nd <= DBL_DIG + 1)
2106 1192 : z = 10*z + c;
2107 1103294 : nz = 0;
2108 : }
2109 : }
2110 : }
2111 214700 : dig_done:
2112 214700 : e = 0;
2113 214700 : if (c == 'e' || c == 'E') {
2114 2737 : if (!nd && !nz && !nz0) {
2115 0 : s = s00;
2116 0 : goto ret;
2117 : }
2118 2737 : s00 = s;
2119 2737 : esign = 0;
2120 2737 : switch(c = *++s) {
2121 : case '-':
2122 724 : esign = 1;
2123 : case '+':
2124 942 : c = *++s;
2125 : }
2126 5473 : if (c >= '0' && c <= '9') {
2127 5488 : while(c == '0')
2128 16 : c = *++s;
2129 5462 : if (c > '0' && c <= '9') {
2130 2726 : L = c - '0';
2131 2726 : s1 = s;
2132 6894 : while((c = *++s) >= '0' && c <= '9')
2133 1442 : L = 10*L + c - '0';
2134 2727 : if (s - s1 > 8 || L > 19999)
2135 : /* Avoid confusion from exponents
2136 : * so large that e might overflow.
2137 : */
2138 1 : e = 19999; /* safe for 16 bit ints */
2139 : else
2140 2725 : e = (int)L;
2141 2726 : if (esign)
2142 724 : e = -e;
2143 : }
2144 : else
2145 10 : e = 0;
2146 : }
2147 : else
2148 1 : s = s00;
2149 : }
2150 214700 : if (!nd) {
2151 554 : if (!nz && !nz0)
2152 375 : s = s00;
2153 554 : goto ret;
2154 : }
2155 214146 : e1 = e -= nf;
2156 :
2157 : /* Now we have nd0 digits, starting at s0, followed by a
2158 : * decimal point, followed by nd-nd0 digits. The number we're
2159 : * after is the integer represented by those digits times
2160 : * 10**e */
2161 :
2162 214146 : if (!nd0)
2163 201634 : nd0 = nd;
2164 214146 : k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
2165 214146 : value(rv) = y;
2166 214146 : if (k > 9)
2167 3961 : value(rv) = tens[k - 9] * value(rv) + z;
2168 214146 : bd0 = 0;
2169 214146 : if (nd <= DBL_DIG
2170 : #ifndef RND_PRODQUOT
2171 : && FLT_ROUNDS == 1
2172 : #endif
2173 : ) {
2174 212324 : if (!e)
2175 3337 : goto ret;
2176 208987 : if (e > 0) {
2177 1655 : if (e <= Ten_pmax) {
2178 : #ifdef VAX
2179 : goto vax_ovfl_check;
2180 : #else
2181 1566 : /* value(rv) = */ rounded_product(value(rv),
2182 : tens[e]);
2183 1566 : goto ret;
2184 : #endif
2185 : }
2186 89 : i = DBL_DIG - nd;
2187 89 : if (e <= Ten_pmax + i) {
2188 : /* A fancier test would sometimes let us do
2189 : * this for larger i values.
2190 : */
2191 25 : e -= i;
2192 25 : value(rv) *= tens[i];
2193 : #ifdef VAX
2194 : /* VAX exponent range is so narrow we must
2195 : * worry about overflow here...
2196 : */
2197 : vax_ovfl_check:
2198 : word0(rv) -= P*Exp_msk1;
2199 : /* value(rv) = */ rounded_product(value(rv),
2200 : tens[e]);
2201 : if ((word0(rv) & Exp_mask)
2202 : > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
2203 : goto ovfl;
2204 : word0(rv) += P*Exp_msk1;
2205 : #else
2206 25 : /* value(rv) = */ rounded_product(value(rv),
2207 : tens[e]);
2208 : #endif
2209 25 : goto ret;
2210 : }
2211 : }
2212 : #ifndef Inaccurate_Divide
2213 207332 : else if (e >= -Ten_pmax) {
2214 207275 : /* value(rv) = */ rounded_quotient(value(rv),
2215 : tens[-e]);
2216 207275 : goto ret;
2217 : }
2218 : #endif
2219 : }
2220 1943 : e1 += nd - k;
2221 :
2222 : /* Get starting approximation = rv * 10**e1 */
2223 :
2224 1943 : if (e1 > 0) {
2225 1729 : if ((i = e1 & 15))
2226 1676 : value(rv) *= tens[i];
2227 1729 : if (e1 &= ~15) {
2228 794 : if (e1 > DBL_MAX_10_EXP) {
2229 14 : ovfl:
2230 14 : errno = ERANGE;
2231 : #ifndef Bad_float_h
2232 14 : value(rv) = HUGE_VAL;
2233 : #else
2234 : /* Can't trust HUGE_VAL */
2235 : #ifdef IEEE_Arith
2236 : word0(rv) = Exp_mask;
2237 : word1(rv) = 0;
2238 : #else
2239 : word0(rv) = Big0;
2240 : word1(rv) = Big1;
2241 : #endif
2242 : #endif
2243 14 : if (bd0)
2244 0 : goto retfree;
2245 14 : goto ret;
2246 : }
2247 783 : if (e1 >>= 4) {
2248 873 : for(j = 0; e1 > 1; j++, e1 >>= 1)
2249 90 : if (e1 & 1)
2250 23 : value(rv) *= bigtens[j];
2251 : /* The last multiplication could overflow. */
2252 783 : word0(rv) -= P*Exp_msk1;
2253 783 : value(rv) *= bigtens[j];
2254 783 : if ((z = word0(rv) & Exp_mask)
2255 : > Exp_msk1*(DBL_MAX_EXP+Bias-P))
2256 3 : goto ovfl;
2257 780 : if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
2258 : /* set to largest number */
2259 : /* (Can't trust DBL_MAX) */
2260 0 : word0(rv) = Big0;
2261 0 : word1(rv) = Big1;
2262 : }
2263 : else
2264 780 : word0(rv) += P*Exp_msk1;
2265 : }
2266 :
2267 : }
2268 : }
2269 214 : else if (e1 < 0) {
2270 154 : e1 = -e1;
2271 154 : if ((i = e1 & 15))
2272 149 : value(rv) /= tens[i];
2273 154 : if (e1 &= ~15) {
2274 71 : e1 >>= 4;
2275 71 : if (e1 >= 1 << n_bigtens)
2276 1 : goto undfl;
2277 142 : for(j = 0; e1 > 1; j++, e1 >>= 1)
2278 72 : if (e1 & 1)
2279 13 : value(rv) *= tinytens[j];
2280 : /* The last multiplication could underflow. */
2281 70 : value(rv0) = value(rv);
2282 70 : value(rv) *= tinytens[j];
2283 70 : if (!value(rv)) {
2284 2 : value(rv) = 2.*value(rv0);
2285 2 : value(rv) *= tinytens[j];
2286 2 : if (!value(rv)) {
2287 3 : undfl:
2288 3 : value(rv) = 0.;
2289 3 : errno = ERANGE;
2290 3 : if (bd0)
2291 1 : goto retfree;
2292 2 : goto ret;
2293 : }
2294 1 : word0(rv) = Tiny0;
2295 1 : word1(rv) = Tiny1;
2296 : /* The refinement below will clean
2297 : * this approximation up.
2298 : */
2299 : }
2300 : }
2301 : }
2302 :
2303 : /* Now the hard part -- adjusting rv to the correct value.*/
2304 :
2305 : /* Put digits into bd: true value = bd * 10^e */
2306 :
2307 1927 : bd0 = s2b(s0, nd0, nd, y);
2308 :
2309 : for(;;) {
2310 1949 : bd = Balloc(bd0->k);
2311 1949 : Bcopy(bd, bd0);
2312 1949 : bb = d2b(value(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
2313 1949 : bs = i2b(1);
2314 :
2315 1949 : if (e >= 0) {
2316 1794 : bb2 = bb5 = 0;
2317 1794 : bd2 = bd5 = e;
2318 : }
2319 : else {
2320 155 : bb2 = bb5 = -e;
2321 155 : bd2 = bd5 = 0;
2322 : }
2323 1949 : if (bbe >= 0)
2324 1825 : bb2 += bbe;
2325 : else
2326 124 : bd2 -= bbe;
2327 1949 : bs2 = bb2;
2328 : #ifdef Sudden_Underflow
2329 : #ifdef IBM
2330 : j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
2331 : #else
2332 : j = P + 1 - bbbits;
2333 : #endif
2334 : #else
2335 1949 : i = bbe + bbbits - 1; /* logb(rv) */
2336 1949 : if (i < Emin) /* denormal */
2337 6 : j = bbe + (P-Emin);
2338 : else
2339 1943 : j = P + 1 - bbbits;
2340 : #endif
2341 1949 : bb2 += j;
2342 1949 : bd2 += j;
2343 1949 : i = bb2 < bd2 ? bb2 : bd2;
2344 1949 : if (i > bs2)
2345 210 : i = bs2;
2346 1949 : if (i > 0) {
2347 1921 : bb2 -= i;
2348 1921 : bd2 -= i;
2349 1921 : bs2 -= i;
2350 : }
2351 1949 : if (bb5 > 0) {
2352 155 : bs = pow5mult(bs, bb5);
2353 155 : bb1 = mult(bs, bb);
2354 155 : Bfree(bb);
2355 155 : bb = bb1;
2356 : }
2357 1949 : if (bb2 > 0)
2358 1949 : bb = lshift(bb, bb2);
2359 1949 : if (bd5 > 0)
2360 50 : bd = pow5mult(bd, bd5);
2361 1949 : if (bd2 > 0)
2362 210 : bd = lshift(bd, bd2);
2363 1949 : if (bs2 > 0)
2364 1675 : bs = lshift(bs, bs2);
2365 1949 : delta = diff(bb, bd);
2366 1949 : dsign = delta->sign;
2367 1949 : delta->sign = 0;
2368 1949 : i = cmp(delta, bs);
2369 1949 : if (i < 0) {
2370 : /* Error is less than half an ulp -- check for
2371 : * special case of mantissa a power of two.
2372 : */
2373 1185 : if (dsign || word1(rv) || word0(rv) & Bndry_mask)
2374 : break;
2375 64 : delta = lshift(delta,Log2P);
2376 64 : if (cmp(delta, bs) > 0)
2377 0 : goto drop_down;
2378 64 : break;
2379 : }
2380 764 : if (i == 0) {
2381 : /* exactly half-way between */
2382 33 : if (dsign) {
2383 30 : if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2384 : && word1(rv) == 0xffffffff) {
2385 : /*boundary case -- increment exponent*/
2386 0 : word0(rv) = (word0(rv) & Exp_mask)
2387 : + Exp_msk1
2388 : #ifdef IBM
2389 : | Exp_msk1 >> 4
2390 : #endif
2391 : ;
2392 0 : word1(rv) = 0;
2393 0 : break;
2394 : }
2395 : }
2396 3 : else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2397 0 : drop_down:
2398 : /* boundary case -- decrement exponent */
2399 : #ifdef Sudden_Underflow
2400 : L = word0(rv) & Exp_mask;
2401 : #ifdef IBM
2402 : if (L < Exp_msk1)
2403 : #else
2404 : if (L <= Exp_msk1)
2405 : #endif
2406 : goto undfl;
2407 : L -= Exp_msk1;
2408 : #else
2409 0 : L = (word0(rv) & Exp_mask) - Exp_msk1;
2410 : #endif
2411 0 : word0(rv) = L | Bndry_mask1;
2412 0 : word1(rv) = 0xffffffff;
2413 : #ifdef IBM
2414 : goto cont;
2415 : #else
2416 0 : break;
2417 : #endif
2418 : }
2419 : #ifndef ROUND_BIASED
2420 33 : if (!(word1(rv) & LSB))
2421 17 : break;
2422 : #endif
2423 16 : if (dsign)
2424 16 : value(rv) += ulp(value(rv));
2425 : #ifndef ROUND_BIASED
2426 : else {
2427 0 : value(rv) -= ulp(value(rv));
2428 : #ifndef Sudden_Underflow
2429 0 : if (!value(rv))
2430 0 : goto undfl;
2431 : #endif
2432 : }
2433 : #endif
2434 16 : break;
2435 : }
2436 731 : if ((aadj = ratio(delta, bs)) <= 2.) {
2437 519 : if (dsign)
2438 477 : aadj = aadj1 = 1.;
2439 83 : else if (word1(rv) || word0(rv) & Bndry_mask) {
2440 : #ifndef Sudden_Underflow
2441 42 : if (word1(rv) == Tiny1 && !word0(rv))
2442 1 : goto undfl;
2443 : #endif
2444 41 : aadj = 1.;
2445 41 : aadj1 = -1.;
2446 : }
2447 : else {
2448 : /* special case -- power of FLT_RADIX to be */
2449 : /* rounded down... */
2450 :
2451 0 : if (aadj < 2./FLT_RADIX)
2452 0 : aadj = 1./FLT_RADIX;
2453 : else
2454 0 : aadj *= 0.5;
2455 0 : aadj1 = -aadj;
2456 : }
2457 : }
2458 : else {
2459 212 : aadj *= 0.5;
2460 212 : aadj1 = dsign ? aadj : -aadj;
2461 : #ifdef Check_FLT_ROUNDS
2462 : switch(FLT_ROUNDS) {
2463 : case 2: /* towards +infinity */
2464 : aadj1 -= 0.5;
2465 : break;
2466 : case 0: /* towards 0 */
2467 : case 3: /* towards -infinity */
2468 : aadj1 += 0.5;
2469 : }
2470 : #else
2471 : if (FLT_ROUNDS == 0)
2472 : aadj1 += 0.5;
2473 : #endif
2474 : }
2475 730 : y = word0(rv) & Exp_mask;
2476 :
2477 : /* Check for overflow */
2478 :
2479 730 : if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2480 0 : value(rv0) = value(rv);
2481 0 : word0(rv) -= P*Exp_msk1;
2482 0 : adj = aadj1 * ulp(value(rv));
2483 0 : value(rv) += adj;
2484 0 : if ((word0(rv) & Exp_mask) >=
2485 : Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2486 0 : if (word0(rv0) == Big0 && word1(rv0) == Big1)
2487 0 : goto ovfl;
2488 0 : word0(rv) = Big0;
2489 0 : word1(rv) = Big1;
2490 0 : goto cont;
2491 : }
2492 : else
2493 0 : word0(rv) += P*Exp_msk1;
2494 : }
2495 : else {
2496 : #ifdef Sudden_Underflow
2497 : if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2498 : value(rv0) = value(rv);
2499 : word0(rv) += P*Exp_msk1;
2500 : adj = aadj1 * ulp(value(rv));
2501 : value(rv) += adj;
2502 : #ifdef IBM
2503 : if ((word0(rv) & Exp_mask) < P*Exp_msk1)
2504 : #else
2505 : if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2506 : #endif
2507 : {
2508 : if (word0(rv0) == Tiny0
2509 : && word1(rv0) == Tiny1)
2510 : goto undfl;
2511 : word0(rv) = Tiny0;
2512 : word1(rv) = Tiny1;
2513 : goto cont;
2514 : }
2515 : else
2516 : word0(rv) -= P*Exp_msk1;
2517 : }
2518 : else {
2519 : adj = aadj1 * ulp(value(rv));
2520 : value(rv) += adj;
2521 : }
2522 : #else
2523 : /* Compute adj so that the IEEE rounding rules will
2524 : * correctly round rv + adj in some half-way cases.
2525 : * If rv * ulp(rv) is denormalized (i.e.,
2526 : * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2527 : * trouble from bits lost to denormalization;
2528 : * example: 1.2e-307 .
2529 : */
2530 730 : if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
2531 0 : aadj1 = (double)(int)(aadj + 0.5);
2532 0 : if (!dsign)
2533 0 : aadj1 = -aadj1;
2534 : }
2535 730 : adj = aadj1 * ulp(value(rv));
2536 730 : value(rv) += adj;
2537 : #endif
2538 : }
2539 730 : z = word0(rv) & Exp_mask;
2540 730 : if (y == z) {
2541 : /* Can we stop now? */
2542 708 : L = aadj;
2543 708 : aadj -= L;
2544 : /* The tolerances below are conservative. */
2545 708 : if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2546 708 : if (aadj < .4999999 || aadj > .5000001)
2547 : break;
2548 : }
2549 0 : else if (aadj < .4999999/FLT_RADIX)
2550 0 : break;
2551 : }
2552 22 : cont:
2553 22 : Bfree(bb);
2554 22 : Bfree(bd);
2555 22 : Bfree(bs);
2556 22 : Bfree(delta);
2557 22 : }
2558 1927 : retfree:
2559 1927 : Bfree(bb);
2560 1927 : Bfree(bd);
2561 1927 : Bfree(bs);
2562 1927 : Bfree(bd0);
2563 1927 : Bfree(delta);
2564 214894 : ret:
2565 214894 : if (se)
2566 205380 : *se = (char *)s;
2567 214894 : result = sign ? -value(rv) : value(rv);
2568 :
2569 : _THREAD_PRIVATE_MUTEX_LOCK(pow5mult_mutex);
2570 430551 : while (p5s) {
2571 763 : tmp = p5s;
2572 763 : p5s = p5s->next;
2573 763 : free(tmp);
2574 : }
2575 : _THREAD_PRIVATE_MUTEX_UNLOCK(pow5mult_mutex);
2576 :
2577 214894 : return result;
2578 : }
2579 :
2580 : ZEND_API double zend_hex_strtod(const char *str, char **endptr)
2581 97 : {
2582 97 : const char *s = str;
2583 : char c;
2584 97 : int any = 0;
2585 97 : double value = 0;
2586 :
2587 97 : if (*s == '0' && (s[1] == 'x' || s[1] == 'X')) {
2588 0 : s += 2;
2589 : }
2590 :
2591 1053 : while ((c = *s++)) {
2592 1683 : if (c >= '0' && c <= '9') {
2593 727 : c -= '0';
2594 288 : } else if (c >= 'A' && c <= 'F') {
2595 59 : c -= 'A' - 10;
2596 170 : } else if (c >= 'a' && c <= 'f') {
2597 73 : c -= 'a' - 10;
2598 : } else {
2599 : break;
2600 : }
2601 :
2602 859 : any = 1;
2603 859 : value = value * 16 + c;
2604 : }
2605 :
2606 97 : if (endptr != NULL) {
2607 0 : *endptr = (char *)(any ? s - 1 : str);
2608 : }
2609 :
2610 97 : return value;
2611 : }
2612 :
2613 : ZEND_API double zend_oct_strtod(const char *str, char **endptr)
2614 83 : {
2615 83 : const char *s = str;
2616 : char c;
2617 83 : double value = 0;
2618 83 : int any = 0;
2619 :
2620 : /* skip leading zero */
2621 83 : s++;
2622 :
2623 1138 : while ((c = *s++)) {
2624 1055 : if (c < '0' || c > '7') {
2625 : /* break and return the current value if the number is not well-formed
2626 : * that's what Linux strtol() does
2627 : */
2628 : break;
2629 : }
2630 972 : value = value * 8 + c - '0';
2631 972 : any = 1;
2632 : }
2633 :
2634 83 : if (endptr != NULL) {
2635 0 : *endptr = (char *)(any ? s - 1 : str);
2636 : }
2637 :
2638 83 : return value;
2639 : }
2640 :
2641 : /*
2642 : * Local variables:
2643 : * tab-width: 4
2644 : * c-basic-offset: 4
2645 : * End:
2646 : * vim600: sw=4 ts=4 fdm=marker
2647 : * vim<600: sw=4 ts=4
2648 : */
|
