00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 #include <config.h>
00039
00040 #include <drizzled/internal/m_string.h>
00041
00042 #include <float.h>
00043
00044 #include <cstdlib>
00045 #include <cerrno>
00046 #include <algorithm>
00047
00048 using namespace std;
00049
00050 namespace drizzled
00051 {
00052 namespace internal
00053 {
00054
00055
00056 #define DTOA_OVERFLOW 9999
00057
00058 static double my_strtod_int(const char *, char **, int *);
00059 static char *dtoa(double, int, int, int *, int *, char **);
00060
00092 size_t my_fcvt(double x, int precision, char *to, bool *error)
00093 {
00094 int decpt, sign, len, i;
00095 char *res, *src, *end, *dst= to;
00096 assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
00097
00098 res= dtoa(x, 5, precision, &decpt, &sign, &end);
00099
00100 if (decpt == DTOA_OVERFLOW)
00101 {
00102 free(res);
00103 *to++= '0';
00104 *to= '\0';
00105 if (error != NULL)
00106 *error= true;
00107 return 1;
00108 }
00109
00110 src= res;
00111 len= end - src;
00112
00113 if (sign)
00114 *dst++= '-';
00115
00116 if (decpt <= 0)
00117 {
00118 *dst++= '0';
00119 *dst++= '.';
00120 for (i= decpt; i < 0; i++)
00121 *dst++= '0';
00122 }
00123
00124 for (i= 1; i <= len; i++)
00125 {
00126 *dst++= *src++;
00127 if (i == decpt && i < len)
00128 *dst++= '.';
00129 }
00130 while (i++ <= decpt)
00131 *dst++= '0';
00132
00133 if (precision > 0)
00134 {
00135 if (len <= decpt)
00136 *dst++= '.';
00137
00138 for (i= precision - max(0, (len - decpt)); i > 0; i--)
00139 *dst++= '0';
00140 }
00141
00142 *dst= '\0';
00143 if (error != NULL)
00144 *error= false;
00145
00146 free(res);
00147
00148 return dst - to;
00149 }
00150
00214 size_t my_gcvt(double x, my_gcvt_arg_type type, int width, char *to,
00215 bool *error)
00216 {
00217 int decpt, sign, len, exp_len;
00218 char *res, *src, *end, *dst= to, *dend= dst + width;
00219 bool have_space, force_e_format;
00220 assert(width > 0 && to != NULL);
00221
00222
00223 if (x < 0.)
00224 width--;
00225
00226 res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? min(width, DBL_DIG) :
00227 min(width, FLT_DIG), &decpt, &sign, &end);
00228
00229 if (decpt == DTOA_OVERFLOW)
00230 {
00231 free(res);
00232 *to++= '0';
00233 *to= '\0';
00234 if (error != NULL)
00235 *error= true;
00236 return 1;
00237 }
00238
00239 if (error != NULL)
00240 *error= false;
00241
00242 src= res;
00243 len= end - res;
00244
00245
00246
00247
00248
00249
00250 exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261 have_space= (decpt <= 0 ? len - decpt + 2 :
00262 decpt > 0 && decpt < len ? len + 1 :
00263 decpt) <= width;
00264
00265
00266
00267
00268
00269 force_e_format= (decpt <= 0 && width <= 2 - decpt && width >= 3 + exp_len);
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296 if ((have_space ||
00297
00298
00299
00300
00301 ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
00302 (len > 1 || !force_e_format)))) &&
00303 !force_e_format)) &&
00304
00305
00306
00307
00308
00309 (!have_space || (decpt >= -MAX_DECPT_FOR_F_FORMAT + 1 &&
00310 (decpt <= MAX_DECPT_FOR_F_FORMAT || len > decpt))))
00311 {
00312
00313 int i;
00314
00315 width-= (decpt < len) + (decpt <= 0 ? 1 - decpt : 0);
00316
00317
00318 if (width < len)
00319 {
00320 if (width < decpt)
00321 {
00322 if (error != NULL)
00323 *error= true;
00324 width= decpt;
00325 }
00326
00327
00328
00329
00330
00331
00332 free(res);
00333 res= dtoa(x, 5, width - decpt, &decpt, &sign, &end);
00334 src= res;
00335 len= end - res;
00336 }
00337
00338 if (len == 0)
00339 {
00340
00341 *dst++= '0';
00342 goto end;
00343 }
00344
00345
00346
00347
00348
00349 if (sign && dst < dend)
00350 *dst++= '-';
00351 if (decpt <= 0)
00352 {
00353 if (dst < dend)
00354 *dst++= '0';
00355 if (len > 0 && dst < dend)
00356 *dst++= '.';
00357 for (; decpt < 0 && dst < dend; decpt++)
00358 *dst++= '0';
00359 }
00360
00361 for (i= 1; i <= len && dst < dend; i++)
00362 {
00363 *dst++= *src++;
00364 if (i == decpt && i < len && dst < dend)
00365 *dst++= '.';
00366 }
00367 while (i++ <= decpt && dst < dend)
00368 *dst++= '0';
00369 }
00370 else
00371 {
00372
00373 int decpt_sign= 0;
00374
00375 if (--decpt < 0)
00376 {
00377 decpt= -decpt;
00378 width--;
00379 decpt_sign= 1;
00380 }
00381 width-= 1 + exp_len;
00382
00383 if (len > 1)
00384 width--;
00385
00386 if (width <= 0)
00387 {
00388
00389 if (error != NULL)
00390 *error= true;
00391 width= 0;
00392 }
00393
00394
00395 if (width < len)
00396 {
00397
00398 free(res);
00399 res= dtoa(x, 4, width, &decpt, &sign, &end);
00400 src= res;
00401 len= end - res;
00402 if (--decpt < 0)
00403 decpt= -decpt;
00404 }
00405
00406
00407
00408
00409 if (sign && dst < dend)
00410 *dst++= '-';
00411 if (dst < dend)
00412 *dst++= *src++;
00413 if (len > 1 && dst < dend)
00414 {
00415 *dst++= '.';
00416 while (src < end && dst < dend)
00417 *dst++= *src++;
00418 }
00419 if (dst < dend)
00420 *dst++= 'e';
00421 if (decpt_sign && dst < dend)
00422 *dst++= '-';
00423
00424 if (decpt >= 100 && dst < dend)
00425 {
00426 *dst++= decpt / 100 + '0';
00427 decpt%= 100;
00428 if (dst < dend)
00429 *dst++= decpt / 10 + '0';
00430 }
00431 else if (decpt >= 10 && dst < dend)
00432 *dst++= decpt / 10 + '0';
00433 if (dst < dend)
00434 *dst++= decpt % 10 + '0';
00435
00436 }
00437
00438 end:
00439 free(res);
00440 *dst= '\0';
00441
00442 return dst - to;
00443 }
00444
00463 double my_strtod(const char *str, char **end, int *error)
00464 {
00465 double res;
00466 assert(str != NULL && end != NULL && *end != NULL && error != NULL);
00467
00468 res= my_strtod_int(str, end, error);
00469 return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX);
00470 }
00471
00472
00473 double my_atof(const char *nptr)
00474 {
00475 int error;
00476 const char *end= nptr+65535;
00477 return (my_strtod(nptr, (char**) &end, &error));
00478 }
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535 typedef int32_t Long;
00536 typedef uint32_t ULong;
00537 typedef int64_t LLong;
00538 typedef uint64_t ULLong;
00539
00540 typedef union { double d; ULong L[2]; } U;
00541
00542 #if defined(WORDS_BIGENDIAN) || (defined(__FLOAT_WORD_ORDER) && \
00543 (__FLOAT_WORD_ORDER == __BIG_ENDIAN))
00544 #define word0(x) ((U*)&x)->L[0]
00545 #define word1(x) ((U*)&x)->L[1]
00546 #else
00547 #define word0(x) ((U*)&x)->L[1]
00548 #define word1(x) ((U*)&x)->L[0]
00549 #endif
00550
00551 #define dval(x) ((U*)&x)->d
00552
00553
00554
00555
00556
00557
00558
00559 #define Exp_shift 20
00560 #define Exp_shift1 20
00561 #define Exp_msk1 0x100000
00562 #define Exp_mask 0x7ff00000
00563 #define P 53
00564 #define Bias 1023
00565 #define Emin (-1022)
00566 #define Exp_1 0x3ff00000
00567 #define Exp_11 0x3ff00000
00568 #define Ebits 11
00569 #define Frac_mask 0xfffff
00570 #define Frac_mask1 0xfffff
00571 #define Ten_pmax 22
00572 #define Bletch 0x10
00573 #define Bndry_mask 0xfffff
00574 #define Bndry_mask1 0xfffff
00575 #define LSB 1
00576 #define Sign_bit 0x80000000
00577 #define Log2P 1
00578 #define Tiny1 1
00579 #define Quick_max 14
00580 #define Int_max 14
00581
00582 #ifndef Flt_Rounds
00583 #ifdef FLT_ROUNDS
00584 #define Flt_Rounds FLT_ROUNDS
00585 #else
00586 #define Flt_Rounds 1
00587 #endif
00588 #endif
00589
00590 #define rounded_product(a,b) a*= b
00591 #define rounded_quotient(a,b) a/= b
00592
00593 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
00594 #define Big1 0xffffffff
00595 #define FFFFFFFF 0xffffffffUL
00596
00597
00598
00599 typedef struct Bigint
00600 {
00601 union {
00602 ULong *x;
00603 struct Bigint *next;
00604 } p;
00605 int k;
00606 int maxwds;
00607 int sign;
00608 int wds;
00609 } Bigint;
00610
00611 static Bigint *Bcopy(Bigint* dst, Bigint* src)
00612 {
00613 dst->sign= src->sign;
00614 dst->wds= src->wds;
00615
00616 assert(dst->maxwds >= src->wds);
00617
00618 memcpy(dst->p.x, src->p.x, src->wds*sizeof(ULong));
00619
00620 return dst;
00621 }
00622
00623 static Bigint *Balloc(int k)
00624 {
00625 Bigint *rv;
00626
00627
00628
00629 int x= 1 << k;
00630 rv= (Bigint*) malloc(sizeof(Bigint));
00631
00632 rv->p.x= (ULong*)malloc(x * sizeof(ULong));
00633
00634 rv->k= k;
00635 rv->maxwds= x;
00636
00637 rv->sign= rv->wds= 0;
00638
00639 return rv;
00640 }
00641
00642
00643
00644
00645
00646
00647
00648 static void Bfree(Bigint *v)
00649 {
00650 if(!v)
00651 return;
00652 free(v->p.x);
00653 free(v);
00654 }
00655
00656
00657
00658
00659
00660 static Bigint *multadd(Bigint *b, int m, int a)
00661 {
00662 int i, wds;
00663 ULong *x;
00664 ULLong carry, y;
00665 Bigint *b1;
00666
00667 wds= b->wds;
00668 x= b->p.x;
00669 i= 0;
00670 carry= a;
00671 do
00672 {
00673 y= *x * (ULLong)m + carry;
00674 carry= y >> 32;
00675 *x++= (ULong)(y & FFFFFFFF);
00676 }
00677 while (++i < wds);
00678 if (carry)
00679 {
00680 if (wds >= b->maxwds)
00681 {
00682 b1= Balloc(b->k+1);
00683 Bcopy(b1, b);
00684 Bfree(b);
00685 b= b1;
00686 }
00687 b->p.x[wds++]= (ULong) carry;
00688 b->wds= wds;
00689 }
00690 return b;
00691 }
00692
00693
00694 static Bigint *s2b(const char *s, int nd0, int nd, ULong y9)
00695 {
00696 Bigint *b;
00697 int i, k;
00698 Long x, y;
00699
00700 x= (nd + 8) / 9;
00701 for (k= 0, y= 1; x > y; y <<= 1, k++) ;
00702 b= Balloc(k);
00703 b->p.x[0]= y9;
00704 b->wds= 1;
00705
00706 i= 9;
00707 if (9 < nd0)
00708 {
00709 s+= 9;
00710 do
00711 b= multadd(b, 10, *s++ - '0');
00712 while (++i < nd0);
00713 s++;
00714 }
00715 else
00716 s+= 10;
00717 for(; i < nd; i++)
00718 b= multadd(b, 10, *s++ - '0');
00719 return b;
00720 }
00721
00722
00723 static int hi0bits(ULong x)
00724 {
00725 int k= 0;
00726
00727 if (!(x & 0xffff0000))
00728 {
00729 k= 16;
00730 x<<= 16;
00731 }
00732 if (!(x & 0xff000000))
00733 {
00734 k+= 8;
00735 x<<= 8;
00736 }
00737 if (!(x & 0xf0000000))
00738 {
00739 k+= 4;
00740 x<<= 4;
00741 }
00742 if (!(x & 0xc0000000))
00743 {
00744 k+= 2;
00745 x<<= 2;
00746 }
00747 if (!(x & 0x80000000))
00748 {
00749 k++;
00750 if (!(x & 0x40000000))
00751 return 32;
00752 }
00753 return k;
00754 }
00755
00756
00757 static int lo0bits(ULong *y)
00758 {
00759 int k;
00760 ULong x= *y;
00761
00762 if (x & 7)
00763 {
00764 if (x & 1)
00765 return 0;
00766 if (x & 2)
00767 {
00768 *y= x >> 1;
00769 return 1;
00770 }
00771 *y= x >> 2;
00772 return 2;
00773 }
00774 k= 0;
00775 if (!(x & 0xffff))
00776 {
00777 k= 16;
00778 x>>= 16;
00779 }
00780 if (!(x & 0xff))
00781 {
00782 k+= 8;
00783 x>>= 8;
00784 }
00785 if (!(x & 0xf))
00786 {
00787 k+= 4;
00788 x>>= 4;
00789 }
00790 if (!(x & 0x3))
00791 {
00792 k+= 2;
00793 x>>= 2;
00794 }
00795 if (!(x & 1))
00796 {
00797 k++;
00798 x>>= 1;
00799 if (!x)
00800 return 32;
00801 }
00802 *y= x;
00803 return k;
00804 }
00805
00806
00807
00808
00809 static Bigint *i2b(int i)
00810 {
00811 Bigint *b;
00812
00813 b= Balloc(1);
00814 b->p.x[0]= i;
00815 b->wds= 1;
00816 return b;
00817 }
00818
00819
00820
00821
00822 static Bigint *mult(Bigint *a, Bigint *b)
00823 {
00824 Bigint *c;
00825 int k, wa, wb, wc;
00826 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
00827 ULong y;
00828 ULLong carry, z;
00829
00830 if (a->wds < b->wds)
00831 {
00832 c= a;
00833 a= b;
00834 b= c;
00835 }
00836 k= a->k;
00837 wa= a->wds;
00838 wb= b->wds;
00839 wc= wa + wb;
00840 if (wc > a->maxwds)
00841 k++;
00842 c= Balloc(k);
00843 for (x= c->p.x, xa= x + wc; x < xa; x++)
00844 *x= 0;
00845 xa= a->p.x;
00846 xae= xa + wa;
00847 xb= b->p.x;
00848 xbe= xb + wb;
00849 xc0= c->p.x;
00850 for (; xb < xbe; xc0++)
00851 {
00852 if ((y= *xb++))
00853 {
00854 x= xa;
00855 xc= xc0;
00856 carry= 0;
00857 do
00858 {
00859 z= *x++ * (ULLong)y + *xc + carry;
00860 carry= z >> 32;
00861 *xc++= (ULong) (z & FFFFFFFF);
00862 }
00863 while (x < xae);
00864 *xc= (ULong) carry;
00865 }
00866 }
00867 for (xc0= c->p.x, xc= xc0 + wc; wc > 0 && !*--xc; --wc) ;
00868 c->wds= wc;
00869 return c;
00870 }
00871
00872
00873
00874
00875
00876
00877
00878 static ULong powers5[]=
00879 {
00880 625UL,
00881
00882 390625UL,
00883
00884 2264035265UL, 35UL,
00885
00886 2242703233UL, 762134875UL, 1262UL,
00887
00888 3211403009UL, 1849224548UL, 3668416493UL, 3913284084UL, 1593091UL,
00889
00890 781532673UL, 64985353UL, 253049085UL, 594863151UL, 3553621484UL,
00891 3288652808UL, 3167596762UL, 2788392729UL, 3911132675UL, 590UL,
00892
00893 2553183233UL, 3201533787UL, 3638140786UL, 303378311UL, 1809731782UL,
00894 3477761648UL, 3583367183UL, 649228654UL, 2915460784UL, 487929380UL,
00895 1011012442UL, 1677677582UL, 3428152256UL, 1710878487UL, 1438394610UL,
00896 2161952759UL, 4100910556UL, 1608314830UL, 349175UL
00897 };
00898
00899
00900 static Bigint p5_a[]=
00901 {
00902
00903 { { powers5 }, 1, 1, 0, 1 },
00904 { { powers5 + 1 }, 1, 1, 0, 1 },
00905 { { powers5 + 2 }, 1, 2, 0, 2 },
00906 { { powers5 + 4 }, 2, 3, 0, 3 },
00907 { { powers5 + 7 }, 3, 5, 0, 5 },
00908 { { powers5 + 12 }, 4, 10, 0, 10 },
00909 { { powers5 + 22 }, 5, 19, 0, 19 }
00910 };
00911
00912 #define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
00913
00914 static Bigint *pow5mult(Bigint *b, int k)
00915 {
00916 Bigint *b1, *p5, *p51;
00917 int i;
00918 static int p05[3]= { 5, 25, 125 };
00919
00920 if ((i= k & 3))
00921 b= multadd(b, p05[i-1], 0);
00922
00923 if (!(k>>= 2))
00924 return b;
00925 p5= p5_a;
00926 for (;;)
00927 {
00928 if (k & 1)
00929 {
00930 b1= mult(b, p5);
00931 Bfree(b);
00932 b= b1;
00933 }
00934 if (!(k>>= 1))
00935 break;
00936
00937 if (p5 < p5_a + P5A_MAX)
00938 ++p5;
00939 else if (p5 == p5_a + P5A_MAX)
00940 p5= mult(p5, p5);
00941 else
00942 {
00943 p51= mult(p5, p5);
00944 Bfree(p5);
00945 p5= p51;
00946 }
00947 }
00948 return b;
00949 }
00950
00951
00952 static Bigint *lshift(Bigint *b, int k)
00953 {
00954 int i, k1, n, n1;
00955 Bigint *b1;
00956 ULong *x, *x1, *xe, z;
00957
00958 n= k >> 5;
00959 k1= b->k;
00960 n1= n + b->wds + 1;
00961 for (i= b->maxwds; n1 > i; i<<= 1)
00962 k1++;
00963 b1= Balloc(k1);
00964 x1= b1->p.x;
00965 for (i= 0; i < n; i++)
00966 *x1++= 0;
00967 x= b->p.x;
00968 xe= x + b->wds;
00969 if (k&= 0x1f)
00970 {
00971 k1= 32 - k;
00972 z= 0;
00973 do
00974 {
00975 *x1++= *x << k | z;
00976 z= *x++ >> k1;
00977 }
00978 while (x < xe);
00979 if ((*x1= z))
00980 ++n1;
00981 }
00982 else
00983 do
00984 *x1++= *x++;
00985 while (x < xe);
00986 b1->wds= n1 - 1;
00987 Bfree(b);
00988 return b1;
00989 }
00990
00991
00992 static int cmp(Bigint *a, Bigint *b)
00993 {
00994 ULong *xa, *xa0, *xb, *xb0;
00995 int i, j;
00996
00997 i= a->wds;
00998 j= b->wds;
00999 if (i-= j)
01000 return i;
01001 xa0= a->p.x;
01002 xa= xa0 + j;
01003 xb0= b->p.x;
01004 xb= xb0 + j;
01005 for (;;)
01006 {
01007 if (*--xa != *--xb)
01008 return *xa < *xb ? -1 : 1;
01009 if (xa <= xa0)
01010 break;
01011 }
01012 return 0;
01013 }
01014
01015
01016 static Bigint *diff(Bigint *a, Bigint *b)
01017 {
01018 Bigint *c;
01019 int i, wa, wb;
01020 ULong *xa, *xae, *xb, *xbe, *xc;
01021 ULLong borrow, y;
01022
01023 i= cmp(a,b);
01024 if (!i)
01025 {
01026 c= Balloc(0);
01027 c->wds= 1;
01028 c->p.x[0]= 0;
01029 return c;
01030 }
01031 if (i < 0)
01032 {
01033 c= a;
01034 a= b;
01035 b= c;
01036 i= 1;
01037 }
01038 else
01039 i= 0;
01040 c= Balloc(a->k);
01041 c->sign= i;
01042 wa= a->wds;
01043 xa= a->p.x;
01044 xae= xa + wa;
01045 wb= b->wds;
01046 xb= b->p.x;
01047 xbe= xb + wb;
01048 xc= c->p.x;
01049 borrow= 0;
01050 do
01051 {
01052 y= (ULLong)*xa++ - *xb++ - borrow;
01053 borrow= y >> 32 & (ULong)1;
01054 *xc++= (ULong) (y & FFFFFFFF);
01055 }
01056 while (xb < xbe);
01057 while (xa < xae)
01058 {
01059 y= *xa++ - borrow;
01060 borrow= y >> 32 & (ULong)1;
01061 *xc++= (ULong) (y & FFFFFFFF);
01062 }
01063 while (!*--xc)
01064 wa--;
01065 c->wds= wa;
01066 return c;
01067 }
01068
01069
01070 static double ulp(double x)
01071 {
01072 Long L;
01073 double a;
01074
01075 L= (word0(x) & Exp_mask) - (P - 1)*Exp_msk1;
01076 word0(a) = L;
01077 word1(a) = 0;
01078 return dval(a);
01079 }
01080
01081
01082 static double b2d(Bigint *a, int *e)
01083 {
01084 ULong *xa, *xa0, w, y, z;
01085 int k;
01086 double d;
01087 #define d0 word0(d)
01088 #define d1 word1(d)
01089
01090 xa0= a->p.x;
01091 xa= xa0 + a->wds;
01092 y= *--xa;
01093 k= hi0bits(y);
01094 *e= 32 - k;
01095 if (k < Ebits)
01096 {
01097 d0= Exp_1 | y >> (Ebits - k);
01098 w= xa > xa0 ? *--xa : 0;
01099 d1= y << ((32-Ebits) + k) | w >> (Ebits - k);
01100 goto ret_d;
01101 }
01102 z= xa > xa0 ? *--xa : 0;
01103 if (k-= Ebits)
01104 {
01105 d0= Exp_1 | y << k | z >> (32 - k);
01106 y= xa > xa0 ? *--xa : 0;
01107 d1= z << k | y >> (32 - k);
01108 }
01109 else
01110 {
01111 d0= Exp_1 | y;
01112 d1= z;
01113 }
01114 ret_d:
01115 #undef d0
01116 #undef d1
01117 return dval(d);
01118 }
01119
01120
01121 static Bigint *d2b(double d, int *e, int *bits)
01122 {
01123 Bigint *b;
01124 int de, k;
01125 ULong *x, y, z;
01126 int i;
01127 #define d0 word0(d)
01128 #define d1 word1(d)
01129
01130 b= Balloc(1);
01131 x= b->p.x;
01132
01133 z= d0 & Frac_mask;
01134 d0 &= 0x7fffffff;
01135 if ((de= (int)(d0 >> Exp_shift)))
01136 z|= Exp_msk1;
01137 if ((y= d1))
01138 {
01139 if ((k= lo0bits(&y)))
01140 {
01141 x[0]= y | z << (32 - k);
01142 z>>= k;
01143 }
01144 else
01145 x[0]= y;
01146 i= b->wds= (x[1]= z) ? 2 : 1;
01147 }
01148 else
01149 {
01150 k= lo0bits(&z);
01151 x[0]= z;
01152 i= b->wds= 1;
01153 k+= 32;
01154 }
01155 if (de)
01156 {
01157 *e= de - Bias - (P-1) + k;
01158 *bits= P - k;
01159 }
01160 else
01161 {
01162 *e= de - Bias - (P-1) + 1 + k;
01163 *bits= 32*i - hi0bits(x[i-1]);
01164 }
01165 return b;
01166 #undef d0
01167 #undef d1
01168 }
01169
01170
01171 static double ratio(Bigint *a, Bigint *b)
01172 {
01173 double da, db;
01174 int k, ka, kb;
01175
01176 dval(da)= b2d(a, &ka);
01177 dval(db)= b2d(b, &kb);
01178 k= ka - kb + 32*(a->wds - b->wds);
01179 if (k > 0)
01180 word0(da)+= k*Exp_msk1;
01181 else
01182 {
01183 k= -k;
01184 word0(db)+= k*Exp_msk1;
01185 }
01186 return dval(da) / dval(db);
01187 }
01188
01189 static const double tens[] =
01190 {
01191 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
01192 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
01193 1e20, 1e21, 1e22
01194 };
01195
01196 static const double bigtens[]= { 1e16, 1e32, 1e64, 1e128, 1e256 };
01197 static const double tinytens[]=
01198 { 1e-16, 1e-32, 1e-64, 1e-128,
01199 9007199254740992.*9007199254740992.e-256
01200 };
01201
01202
01203
01204
01205 #define Scale_Bit 0x10
01206 #define n_bigtens 5
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237 static double my_strtod_int(const char *s00, char **se, int *error)
01238 {
01239 int scale;
01240 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
01241 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
01242 const char *s, *s0, *s1, *end = *se;
01243 double aadj, aadj1, adj, rv, rv0;
01244 Long L;
01245 ULong y, z;
01246 Bigint *bb= NULL, *bb1= NULL, *bd= NULL, *bd0= NULL, *bs= NULL, *delta= NULL;
01247 #ifdef SET_INEXACT
01248 int inexact, oldinexact;
01249 #endif
01250
01251 c= 0;
01252 *error= 0;
01253
01254 sign= nz0= nz= 0;
01255 dval(rv)= 0.;
01256 for (s= s00; s < end; s++)
01257 switch (*s) {
01258 case '-':
01259 sign= 1;
01260
01261 case '+':
01262 s++;
01263 goto break2;
01264 case '\t':
01265 case '\n':
01266 case '\v':
01267 case '\f':
01268 case '\r':
01269 case ' ':
01270 continue;
01271 default:
01272 goto break2;
01273 }
01274 break2:
01275 if (s >= end)
01276 goto ret0;
01277
01278 if (*s == '0')
01279 {
01280 nz0= 1;
01281 while (++s < end && *s == '0') ;
01282 if (s >= end)
01283 goto ret;
01284 }
01285 s0= s;
01286 y= z= 0;
01287 for (nd= nf= 0; s < end && (c= *s) >= '0' && c <= '9'; nd++, s++)
01288 if (nd < 9)
01289 y= 10*y + c - '0';
01290 else if (nd < 16)
01291 z= 10*z + c - '0';
01292 nd0= nd;
01293 if (s < end - 1 && c == '.')
01294 {
01295 c= *++s;
01296 if (!nd)
01297 {
01298 for (; s < end && c == '0'; c= *++s)
01299 nz++;
01300 if (s < end && c > '0' && c <= '9')
01301 {
01302 s0= s;
01303 nf+= nz;
01304 nz= 0;
01305 goto have_dig;
01306 }
01307 goto dig_done;
01308 }
01309 for (; s < end && c >= '0' && c <= '9'; c = *++s)
01310 {
01311 have_dig:
01312 nz++;
01313 if (c-= '0')
01314 {
01315 nf+= nz;
01316 for (i= 1; i < nz; i++)
01317 if (nd++ < 9)
01318 y*= 10;
01319 else if (nd <= DBL_DIG + 1)
01320 z*= 10;
01321 if (nd++ < 9)
01322 y= 10*y + c;
01323 else if (nd <= DBL_DIG + 1)
01324 z= 10*z + c;
01325 nz= 0;
01326 }
01327 }
01328 }
01329 dig_done:
01330 e= 0;
01331 if (s < end && (c == 'e' || c == 'E'))
01332 {
01333 if (!nd && !nz && !nz0)
01334 goto ret0;
01335 s00= s;
01336 esign= 0;
01337 if (++s < end)
01338 switch (c= *s) {
01339 case '-':
01340 esign= 1;
01341 case '+':
01342 c= *++s;
01343 }
01344 if (s < end && c >= '0' && c <= '9')
01345 {
01346 while (s < end && c == '0')
01347 c= *++s;
01348 if (s < end && c > '0' && c <= '9') {
01349 L= c - '0';
01350 s1= s;
01351 while (++s < end && (c= *s) >= '0' && c <= '9')
01352 L= 10*L + c - '0';
01353 if (s - s1 > 8 || L > 19999)
01354
01355
01356
01357 e= 19999;
01358 else
01359 e= (int)L;
01360 if (esign)
01361 e= -e;
01362 }
01363 else
01364 e= 0;
01365 }
01366 else
01367 s= s00;
01368 }
01369 if (!nd)
01370 {
01371 if (!nz && !nz0)
01372 {
01373 ret0:
01374 s= s00;
01375 sign= 0;
01376 }
01377 goto ret;
01378 }
01379 e1= e -= nf;
01380
01381
01382
01383
01384
01385
01386
01387
01388 if (!nd0)
01389 nd0= nd;
01390 k= nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
01391 dval(rv)= y;
01392 if (k > 9)
01393 {
01394 #ifdef SET_INEXACT
01395 if (k > DBL_DIG)
01396 oldinexact = get_inexact();
01397 #endif
01398 dval(rv)= tens[k - 9] * dval(rv) + z;
01399 }
01400 bd0= 0;
01401 if (nd <= DBL_DIG
01402 )
01403 {
01404 if (!e)
01405 goto ret;
01406 if (e > 0)
01407 {
01408 if (e <= Ten_pmax)
01409 {
01410 rounded_product(dval(rv), tens[e]);
01411 goto ret;
01412 }
01413 i= DBL_DIG - nd;
01414 if (e <= Ten_pmax + i)
01415 {
01416
01417
01418
01419
01420 e-= i;
01421 dval(rv)*= tens[i];
01422 rounded_product(dval(rv), tens[e]);
01423 goto ret;
01424 }
01425 }
01426 #ifndef Inaccurate_Divide
01427 else if (e >= -Ten_pmax)
01428 {
01429 rounded_quotient(dval(rv), tens[-e]);
01430 goto ret;
01431 }
01432 #endif
01433 }
01434 e1+= nd - k;
01435
01436 #ifdef SET_INEXACT
01437 inexact= 1;
01438 if (k <= DBL_DIG)
01439 oldinexact= get_inexact();
01440 #endif
01441 scale= 0;
01442
01443
01444
01445 if (e1 > 0)
01446 {
01447 if ((i= e1 & 15))
01448 dval(rv)*= tens[i];
01449 if (e1&= ~15)
01450 {
01451 if (e1 > DBL_MAX_10_EXP)
01452 {
01453 ovfl:
01454 *error= EOVERFLOW;
01455
01456 word0(rv)= Exp_mask;
01457 word1(rv)= 0;
01458 #ifdef SET_INEXACT
01459
01460 dval(rv0)= 1e300;
01461 dval(rv0)*= dval(rv0);
01462 #endif
01463 if (bd0)
01464 goto retfree;
01465 goto ret;
01466 }
01467 e1>>= 4;
01468 for(j= 0; e1 > 1; j++, e1>>= 1)
01469 if (e1 & 1)
01470 dval(rv)*= bigtens[j];
01471
01472 word0(rv)-= P*Exp_msk1;
01473 dval(rv)*= bigtens[j];
01474 if ((z= word0(rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - P))
01475 goto ovfl;
01476 if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P))
01477 {
01478
01479 word0(rv)= Big0;
01480 word1(rv)= Big1;
01481 }
01482 else
01483 word0(rv)+= P*Exp_msk1;
01484 }
01485 }
01486 else if (e1 < 0)
01487 {
01488 e1= -e1;
01489 if ((i= e1 & 15))
01490 dval(rv)/= tens[i];
01491 if ((e1>>= 4))
01492 {
01493 if (e1 >= 1 << n_bigtens)
01494 goto undfl;
01495 if (e1 & Scale_Bit)
01496 scale= 2 * P;
01497 for(j= 0; e1 > 0; j++, e1>>= 1)
01498 if (e1 & 1)
01499 dval(rv)*= tinytens[j];
01500 if (scale && (j = 2 * P + 1 - ((word0(rv) & Exp_mask) >> Exp_shift)) > 0)
01501 {
01502
01503 if (j >= 32)
01504 {
01505 word1(rv)= 0;
01506 if (j >= 53)
01507 word0(rv)= (P + 2) * Exp_msk1;
01508 else
01509 word0(rv)&= 0xffffffff << (j - 32);
01510 }
01511 else
01512 word1(rv)&= 0xffffffff << j;
01513 }
01514 if (!dval(rv))
01515 {
01516 undfl:
01517 dval(rv)= 0.;
01518 if (bd0)
01519 goto retfree;
01520 goto ret;
01521 }
01522 }
01523 }
01524
01525
01526
01527
01528
01529 bd0= s2b(s0, nd0, nd, y);
01530
01531 for(;;)
01532 {
01533 bd= Balloc(bd0->k);
01534 Bcopy(bd, bd0);
01535 bb= d2b(dval(rv), &bbe, &bbbits);
01536 bs= i2b(1);
01537
01538 if (e >= 0)
01539 {
01540 bb2= bb5= 0;
01541 bd2= bd5= e;
01542 }
01543 else
01544 {
01545 bb2= bb5= -e;
01546 bd2= bd5= 0;
01547 }
01548 if (bbe >= 0)
01549 bb2+= bbe;
01550 else
01551 bd2-= bbe;
01552 bs2= bb2;
01553 j= bbe - scale;
01554 i= j + bbbits - 1;
01555 if (i < Emin)
01556 j+= P - Emin;
01557 else
01558 j= P + 1 - bbbits;
01559 bb2+= j;
01560 bd2+= j;
01561 bd2+= scale;
01562 i= bb2 < bd2 ? bb2 : bd2;
01563 if (i > bs2)
01564 i= bs2;
01565 if (i > 0)
01566 {
01567 bb2-= i;
01568 bd2-= i;
01569 bs2-= i;
01570 }
01571 if (bb5 > 0)
01572 {
01573 bs= pow5mult(bs, bb5);
01574 bb1= mult(bs, bb);
01575 Bfree(bb);
01576 bb= bb1;
01577 }
01578 if (bb2 > 0)
01579 bb= lshift(bb, bb2);
01580 if (bd5 > 0)
01581 bd= pow5mult(bd, bd5);
01582 if (bd2 > 0)
01583 bd= lshift(bd, bd2);
01584 if (bs2 > 0)
01585 bs= lshift(bs, bs2);
01586 delta= diff(bb, bd);
01587 dsign= delta->sign;
01588 delta->sign= 0;
01589 i= cmp(delta, bs);
01590
01591 if (i < 0)
01592 {
01593
01594
01595
01596
01597 if (dsign || word1(rv) || word0(rv) & Bndry_mask ||
01598 (word0(rv) & Exp_mask) <= (2 * P + 1) * Exp_msk1)
01599 {
01600 #ifdef SET_INEXACT
01601 if (!delta->x[0] && delta->wds <= 1)
01602 inexact= 0;
01603 #endif
01604 break;
01605 }
01606 if (!delta->p.x[0] && delta->wds <= 1)
01607 {
01608
01609 #ifdef SET_INEXACT
01610 inexact= 0;
01611 #endif
01612 break;
01613 }
01614 delta= lshift(delta, Log2P);
01615 if (cmp(delta, bs) > 0)
01616 goto drop_down;
01617 break;
01618 }
01619 if (i == 0)
01620 {
01621
01622 if (dsign)
01623 {
01624 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 &&
01625 word1(rv) ==
01626 ((scale && (y = word0(rv) & Exp_mask) <= 2 * P * Exp_msk1) ?
01627 (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
01628 0xffffffff))
01629 {
01630
01631 word0(rv)= (word0(rv) & Exp_mask) + Exp_msk1;
01632 word1(rv) = 0;
01633 dsign = 0;
01634 break;
01635 }
01636 }
01637 else if (!(word0(rv) & Bndry_mask) && !word1(rv))
01638 {
01639 drop_down:
01640
01641 if (scale)
01642 {
01643 L= word0(rv) & Exp_mask;
01644 if (L <= (2 *P + 1) * Exp_msk1)
01645 {
01646 if (L > (P + 2) * Exp_msk1)
01647
01648 break;
01649
01650 goto undfl;
01651 }
01652 }
01653 L= (word0(rv) & Exp_mask) - Exp_msk1;
01654 word0(rv)= L | Bndry_mask1;
01655 word1(rv)= 0xffffffff;
01656 break;
01657 }
01658 if (!(word1(rv) & LSB))
01659 break;
01660 if (dsign)
01661 dval(rv)+= ulp(dval(rv));
01662 else
01663 {
01664 dval(rv)-= ulp(dval(rv));
01665 if (!dval(rv))
01666 goto undfl;
01667 }
01668 dsign= 1 - dsign;
01669 break;
01670 }
01671 if ((aadj= ratio(delta, bs)) <= 2.)
01672 {
01673 if (dsign)
01674 aadj= aadj1= 1.;
01675 else if (word1(rv) || word0(rv) & Bndry_mask)
01676 {
01677 if (word1(rv) == Tiny1 && !word0(rv))
01678 goto undfl;
01679 aadj= 1.;
01680 aadj1= -1.;
01681 }
01682 else
01683 {
01684
01685 if (aadj < 2. / FLT_RADIX)
01686 aadj= 1. / FLT_RADIX;
01687 else
01688 aadj*= 0.5;
01689 aadj1= -aadj;
01690 }
01691 }
01692 else
01693 {
01694 aadj*= 0.5;
01695 aadj1= dsign ? aadj : -aadj;
01696 if (Flt_Rounds == 0)
01697 aadj1+= 0.5;
01698 }
01699 y= word0(rv) & Exp_mask;
01700
01701
01702
01703 if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
01704 {
01705 dval(rv0)= dval(rv);
01706 word0(rv)-= P * Exp_msk1;
01707 adj= aadj1 * ulp(dval(rv));
01708 dval(rv)+= adj;
01709 if ((word0(rv) & Exp_mask) >= Exp_msk1 * (DBL_MAX_EXP + Bias - P))
01710 {
01711 if (word0(rv0) == Big0 && word1(rv0) == Big1)
01712 goto ovfl;
01713 word0(rv)= Big0;
01714 word1(rv)= Big1;
01715 goto cont;
01716 }
01717 else
01718 word0(rv)+= P * Exp_msk1;
01719 }
01720 else
01721 {
01722 if (scale && y <= 2 * P * Exp_msk1)
01723 {
01724 if (aadj <= 0x7fffffff)
01725 {
01726 if ((z= (ULong) aadj) <= 0)
01727 z= 1;
01728 aadj= z;
01729 aadj1= dsign ? aadj : -aadj;
01730 }
01731 word0(aadj1)+= (2 * P + 1) * Exp_msk1 - y;
01732 }
01733 adj = aadj1 * ulp(dval(rv));
01734 dval(rv) += adj;
01735 }
01736 z= word0(rv) & Exp_mask;
01737 #ifndef SET_INEXACT
01738 if (!scale)
01739 if (y == z)
01740 {
01741
01742 L= (Long)aadj;
01743 aadj-= L;
01744
01745 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
01746 {
01747 if (aadj < .4999999 || aadj > .5000001)
01748 break;
01749 }
01750 else if (aadj < .4999999 / FLT_RADIX)
01751 break;
01752 }
01753 #endif
01754 cont:
01755 Bfree(bb);
01756 Bfree(bd);
01757 Bfree(bs);
01758 Bfree(delta);
01759 }
01760 #ifdef SET_INEXACT
01761 if (inexact)
01762 {
01763 if (!oldinexact)
01764 {
01765 word0(rv0)= Exp_1 + (70 << Exp_shift);
01766 word1(rv0)= 0;
01767 dval(rv0)+= 1.;
01768 }
01769 }
01770 else if (!oldinexact)
01771 clear_inexact();
01772 #endif
01773 if (scale)
01774 {
01775 word0(rv0)= Exp_1 - 2 * P * Exp_msk1;
01776 word1(rv0)= 0;
01777 dval(rv)*= dval(rv0);
01778 }
01779 #ifdef SET_INEXACT
01780 if (inexact && !(word0(rv) & Exp_mask))
01781 {
01782
01783 dval(rv0)= 1e-300;
01784 dval(rv0)*= dval(rv0);
01785 }
01786 #endif
01787 retfree:
01788 Bfree(bb);
01789 Bfree(bd);
01790 Bfree(bs);
01791 Bfree(bd0);
01792 Bfree(delta);
01793 ret:
01794 *se= (char *)s;
01795 return sign ? -dval(rv) : dval(rv);
01796 }
01797
01798
01799 static int quorem(Bigint *b, Bigint *S)
01800 {
01801 int n;
01802 ULong *bx, *bxe, q, *sx, *sxe;
01803 ULLong borrow, carry, y, ys;
01804
01805 n= S->wds;
01806 if (b->wds < n)
01807 return 0;
01808 sx= S->p.x;
01809 sxe= sx + --n;
01810 bx= b->p.x;
01811 bxe= bx + n;
01812 q= *bxe / (*sxe + 1);
01813 if (q)
01814 {
01815 borrow= 0;
01816 carry= 0;
01817 do
01818 {
01819 ys= *sx++ * (ULLong)q + carry;
01820 carry= ys >> 32;
01821 y= *bx - (ys & FFFFFFFF) - borrow;
01822 borrow= y >> 32 & (ULong)1;
01823 *bx++= (ULong) (y & FFFFFFFF);
01824 }
01825 while (sx <= sxe);
01826 if (!*bxe)
01827 {
01828 bx= b->p.x;
01829 while (--bxe > bx && !*bxe)
01830 --n;
01831 b->wds= n;
01832 }
01833 }
01834 if (cmp(b, S) >= 0)
01835 {
01836 q++;
01837 borrow= 0;
01838 carry= 0;
01839 bx= b->p.x;
01840 sx= S->p.x;
01841 do
01842 {
01843 ys= *sx++ + carry;
01844 carry= ys >> 32;
01845 y= *bx - (ys & FFFFFFFF) - borrow;
01846 borrow= y >> 32 & (ULong)1;
01847 *bx++= (ULong) (y & FFFFFFFF);
01848 }
01849 while (sx <= sxe);
01850 bx= b->p.x;
01851 bxe= bx + n;
01852 if (!*bxe)
01853 {
01854 while (--bxe > bx && !*bxe)
01855 --n;
01856 b->wds= n;
01857 }
01858 }
01859 return q;
01860 }
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898 static char *dtoa(double d, int mode, int ndigits, int *decpt, int *sign,
01899 char **rve)
01900 {
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934 int bbits, b2, b5, be, dig, i, ieps, ilim=0, ilim0, ilim1= 0,
01935 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
01936 spec_case, try_quick;
01937 Long L;
01938 int denorm;
01939 ULong x;
01940 Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
01941 double d2, ds, eps;
01942 char *s, *s0;
01943
01944 if (word0(d) & Sign_bit)
01945 {
01946
01947 *sign= 1;
01948 word0(d) &= ~Sign_bit;
01949 }
01950 else
01951 *sign= 0;
01952
01953
01954 if ((((word0(d) & Exp_mask) == Exp_mask) && ((*decpt= DTOA_OVERFLOW) != 0)) ||
01955 (!dval(d) && ((*decpt= 1) != 0)))
01956 {
01957
01958 char *res= (char*) malloc(2);
01959 res[0]= '0';
01960 res[1]= '\0';
01961 if (rve)
01962 *rve= res + 1;
01963 return res;
01964 }
01965
01966
01967 b= d2b(dval(d), &be, &bbits);
01968 if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
01969 {
01970 dval(d2)= dval(d);
01971 word0(d2) &= Frac_mask1;
01972 word0(d2) |= Exp_11;
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997 i-= Bias;
01998 denorm= 0;
01999 }
02000 else
02001 {
02002
02003
02004 i= bbits + be + (Bias + (P-1) - 1);
02005 x= i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
02006 : word1(d) << (32 - i);
02007 dval(d2)= x;
02008 word0(d2)-= 31*Exp_msk1;
02009 i-= (Bias + (P-1) - 1) + 1;
02010 denorm= 1;
02011 }
02012 ds= (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
02013 k= (int)ds;
02014 if (ds < 0. && ds != k)
02015 k--;
02016 k_check= 1;
02017 if (k >= 0 && k <= Ten_pmax)
02018 {
02019 if (dval(d) < tens[k])
02020 k--;
02021 k_check= 0;
02022 }
02023 j= bbits - i - 1;
02024 if (j >= 0)
02025 {
02026 b2= 0;
02027 s2= j;
02028 }
02029 else
02030 {
02031 b2= -j;
02032 s2= 0;
02033 }
02034 if (k >= 0)
02035 {
02036 b5= 0;
02037 s5= k;
02038 s2+= k;
02039 }
02040 else
02041 {
02042 b2-= k;
02043 b5= -k;
02044 s5= 0;
02045 }
02046 if (mode < 0 || mode > 9)
02047 mode= 0;
02048
02049 try_quick= 1;
02050
02051 if (mode > 5)
02052 {
02053 mode-= 4;
02054 try_quick= 0;
02055 }
02056 leftright= 1;
02057 switch (mode) {
02058 case 0:
02059 case 1:
02060 ilim= ilim1= -1;
02061 i= 18;
02062 ndigits= 0;
02063 break;
02064 case 2:
02065 leftright= 0;
02066
02067 case 4:
02068 if (ndigits <= 0)
02069 ndigits= 1;
02070 ilim= ilim1= i= ndigits;
02071 break;
02072 case 3:
02073 leftright= 0;
02074
02075 case 5:
02076 i= ndigits + k + 1;
02077 ilim= i;
02078 ilim1= i - 1;
02079 if (i <= 0)
02080 i= 1;
02081 }
02082 s= s0= (char*)malloc(i+1);
02083
02084
02085 if (ilim >= 0 && ilim <= Quick_max && try_quick)
02086 {
02087
02088 i= 0;
02089 dval(d2)= dval(d);
02090 k0= k;
02091 ilim0= ilim;
02092 ieps= 2;
02093 if (k > 0)
02094 {
02095 ds= tens[k&0xf];
02096 j= k >> 4;
02097 if (j & Bletch)
02098 {
02099
02100 j&= Bletch - 1;
02101 dval(d)/= bigtens[n_bigtens-1];
02102 ieps++;
02103 }
02104 for (; j; j>>= 1, i++)
02105 {
02106 if (j & 1)
02107 {
02108 ieps++;
02109 ds*= bigtens[i];
02110 }
02111 }
02112 dval(d)/= ds;
02113 }
02114 else if ((j1= -k))
02115 {
02116 dval(d)*= tens[j1 & 0xf];
02117 for (j= j1 >> 4; j; j>>= 1, i++)
02118 {
02119 if (j & 1)
02120 {
02121 ieps++;
02122 dval(d)*= bigtens[i];
02123 }
02124 }
02125 }
02126 if (k_check && dval(d) < 1. && ilim > 0)
02127 {
02128 if (ilim1 <= 0)
02129 goto fast_failed;
02130 ilim= ilim1;
02131 k--;
02132 dval(d)*= 10.;
02133 ieps++;
02134 }
02135 dval(eps)= ieps*dval(d) + 7.;
02136 word0(eps)-= (P-1)*Exp_msk1;
02137 if (ilim == 0)
02138 {
02139 S= mhi= 0;
02140 dval(d)-= 5.;
02141 if (dval(d) > dval(eps))
02142 goto one_digit;
02143 if (dval(d) < -dval(eps))
02144 goto no_digits;
02145 goto fast_failed;
02146 }
02147 if (leftright)
02148 {
02149
02150 dval(eps)= 0.5/tens[ilim-1] - dval(eps);
02151 for (i= 0;;)
02152 {
02153 L= (Long) dval(d);
02154 dval(d)-= L;
02155 *s++= '0' + (int)L;
02156 if (dval(d) < dval(eps))
02157 goto ret1;
02158 if (1. - dval(d) < dval(eps))
02159 goto bump_up;
02160 if (++i >= ilim)
02161 break;
02162 dval(eps)*= 10.;
02163 dval(d)*= 10.;
02164 }
02165 }
02166 else
02167 {
02168
02169 dval(eps)*= tens[ilim-1];
02170 for (i= 1;; i++, dval(d)*= 10.)
02171 {
02172 L= (Long)(dval(d));
02173 if (!(dval(d)-= L))
02174 ilim= i;
02175 *s++= '0' + (int)L;
02176 if (i == ilim)
02177 {
02178 if (dval(d) > 0.5 + dval(eps))
02179 goto bump_up;
02180 else if (dval(d) < 0.5 - dval(eps))
02181 {
02182 while (*--s == '0') {}
02183 s++;
02184 goto ret1;
02185 }
02186 break;
02187 }
02188 }
02189 }
02190 fast_failed:
02191 s= s0;
02192 dval(d)= dval(d2);
02193 k= k0;
02194 ilim= ilim0;
02195 }
02196
02197
02198
02199 if (be >= 0 && k <= Int_max)
02200 {
02201
02202 ds= tens[k];
02203 if (ndigits < 0 && ilim <= 0)
02204 {
02205 S= mhi= 0;
02206 if (ilim < 0 || dval(d) <= 5*ds)
02207 goto no_digits;
02208 goto one_digit;
02209 }
02210 for (i= 1;; i++, dval(d)*= 10.)
02211 {
02212 L= (Long)(dval(d) / ds);
02213 dval(d)-= L*ds;
02214 *s++= '0' + (int)L;
02215 if (!dval(d))
02216 {
02217 break;
02218 }
02219 if (i == ilim)
02220 {
02221 dval(d)+= dval(d);
02222 if (dval(d) > ds || (dval(d) == ds && L & 1))
02223 {
02224 bump_up:
02225 while (*--s == '9')
02226 if (s == s0)
02227 {
02228 k++;
02229 *s= '0';
02230 break;
02231 }
02232 ++*s++;
02233 }
02234 break;
02235 }
02236 }
02237 goto ret1;
02238 }
02239
02240 m2= b2;
02241 m5= b5;
02242 mhi= mlo= 0;
02243 if (leftright)
02244 {
02245 i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
02246 b2+= i;
02247 s2+= i;
02248 mhi= i2b(1);
02249 }
02250 if (m2 > 0 && s2 > 0)
02251 {
02252 i= m2 < s2 ? m2 : s2;
02253 b2-= i;
02254 m2-= i;
02255 s2-= i;
02256 }
02257 if (b5 > 0)
02258 {
02259 if (leftright)
02260 {
02261 if (m5 > 0)
02262 {
02263 mhi= pow5mult(mhi, m5);
02264 b1= mult(mhi, b);
02265 Bfree(b);
02266 b= b1;
02267 }
02268 if ((j= b5 - m5))
02269 b= pow5mult(b, j);
02270 }
02271 else
02272 b= pow5mult(b, b5);
02273 }
02274 S= i2b(1);
02275 if (s5 > 0)
02276 S= pow5mult(S, s5);
02277
02278
02279
02280 spec_case= 0;
02281 if ((mode < 2 || leftright)
02282 )
02283 {
02284 if (!word1(d) && !(word0(d) & Bndry_mask) &&
02285 word0(d) & (Exp_mask & ~Exp_msk1)
02286 )
02287 {
02288
02289 b2+= Log2P;
02290 s2+= Log2P;
02291 spec_case= 1;
02292 }
02293 }
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303 if ((i= ((s5 ? 32 - hi0bits(S->p.x[S->wds-1]) : 1) + s2) & 0x1f))
02304 i= 32 - i;
02305 if (i > 4)
02306 {
02307 i-= 4;
02308 b2+= i;
02309 m2+= i;
02310 s2+= i;
02311 }
02312 else if (i < 4)
02313 {
02314 i+= 28;
02315 b2+= i;
02316 m2+= i;
02317 s2+= i;
02318 }
02319 if (b2 > 0)
02320 b= lshift(b, b2);
02321 if (s2 > 0)
02322 S= lshift(S, s2);
02323 if (k_check)
02324 {
02325 if (cmp(b,S) < 0)
02326 {
02327 k--;
02328
02329 b= multadd(b, 10, 0);
02330 if (leftright)
02331 mhi= multadd(mhi, 10, 0);
02332 ilim= ilim1;
02333 }
02334 }
02335 if (ilim <= 0 && (mode == 3 || mode == 5))
02336 {
02337 if (ilim < 0 || cmp(b,S= multadd(S,5,0)) <= 0)
02338 {
02339
02340 no_digits:
02341 k= -1 - ndigits;
02342 goto ret;
02343 }
02344 one_digit:
02345 *s++= '1';
02346 k++;
02347 goto ret;
02348 }
02349 if (leftright)
02350 {
02351 if (m2 > 0)
02352 mhi= lshift(mhi, m2);
02353
02354
02355
02356
02357
02358 mlo= mhi;
02359 if (spec_case)
02360 {
02361 mhi= Balloc(mhi->k);
02362 Bcopy(mhi, mlo);
02363 mhi= lshift(mhi, Log2P);
02364 }
02365
02366 for (i= 1;;i++)
02367 {
02368 dig= quorem(b,S) + '0';
02369
02370 j= cmp(b, mlo);
02371 delta= diff(S, mhi);
02372 j1= delta->sign ? 1 : cmp(b, delta);
02373 Bfree(delta);
02374 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
02375 )
02376 {
02377 if (dig == '9')
02378 goto round_9_up;
02379 if (j > 0)
02380 dig++;
02381 *s++= dig;
02382 goto ret;
02383 }
02384 if (j < 0 || (j == 0 && mode != 1 && !(word1(d) & 1)))
02385 {
02386 if (!b->p.x[0] && b->wds <= 1)
02387 {
02388 goto accept_dig;
02389 }
02390 if (j1 > 0)
02391 {
02392 b= lshift(b, 1);
02393 j1= cmp(b, S);
02394 if ((j1 > 0 || (j1 == 0 && dig & 1))
02395 && dig++ == '9')
02396 goto round_9_up;
02397 }
02398 accept_dig:
02399 *s++= dig;
02400 goto ret;
02401 }
02402 if (j1 > 0)
02403 {
02404 if (dig == '9')
02405 {
02406 round_9_up:
02407 *s++= '9';
02408 goto roundoff;
02409 }
02410 *s++= dig + 1;
02411 goto ret;
02412 }
02413 *s++= dig;
02414 if (i == ilim)
02415 break;
02416 b= multadd(b, 10, 0);
02417 if (mlo == mhi)
02418 mlo= mhi= multadd(mhi, 10, 0);
02419 else
02420 {
02421 mlo= multadd(mlo, 10, 0);
02422 mhi= multadd(mhi, 10, 0);
02423 }
02424 }
02425 }
02426 else
02427 for (i= 1;; i++)
02428 {
02429 *s++= dig= quorem(b,S) + '0';
02430 if (!b->p.x[0] && b->wds <= 1)
02431 {
02432 goto ret;
02433 }
02434 if (i >= ilim)
02435 break;
02436 b= multadd(b, 10, 0);
02437 }
02438
02439
02440
02441 b= lshift(b, 1);
02442 j= cmp(b, S);
02443 if (j > 0 || (j == 0 && dig & 1))
02444 {
02445 roundoff:
02446 while (*--s == '9')
02447 if (s == s0)
02448 {
02449 k++;
02450 *s++= '1';
02451 goto ret;
02452 }
02453 ++*s++;
02454 }
02455 else
02456 {
02457 while (*--s == '0') {}
02458 s++;
02459 }
02460 ret:
02461 Bfree(S);
02462 if (mhi)
02463 {
02464 if (mlo && mlo != mhi)
02465 Bfree(mlo);
02466 Bfree(mhi);
02467 }
02468 ret1:
02469 Bfree(b);
02470 *s= 0;
02471 *decpt= k + 1;
02472 if (rve)
02473 *rve= s;
02474 return s0;
02475 }
02476
02477 }
02478 }