40 #include <drizzled/internal/m_string.h>
56 #define DTOA_OVERFLOW 9999
58 static double my_strtod_int(
const char *,
char **,
int *);
59 static char *dtoa(
double,
int,
int,
int *,
int *,
char **);
92 size_t my_fcvt(
double x,
int precision,
char *to,
bool *error)
94 int decpt, sign, len, i;
95 char *res, *src, *end, *dst= to;
96 assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
98 res= dtoa(x, 5, precision, &decpt, &sign, &end);
100 if (decpt == DTOA_OVERFLOW)
120 for (i= decpt; i < 0; i++)
124 for (i= 1; i <= len; i++)
127 if (i == decpt && i < len)
138 for (i= precision - max(0, (len - decpt)); i > 0; i--)
214 size_t my_gcvt(
double x, my_gcvt_arg_type type,
int width,
char *to,
217 int decpt, sign, len, exp_len;
218 char *res, *src, *end, *dst= to, *dend= dst + width;
219 bool have_space, force_e_format;
220 assert(width > 0 && to != NULL);
226 res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? min(width, DBL_DIG) :
227 min(width, FLT_DIG), &decpt, &sign, &end);
229 if (decpt == DTOA_OVERFLOW)
250 exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
261 have_space= (decpt <= 0 ? len - decpt + 2 :
262 decpt > 0 && decpt < len ? len + 1 :
269 force_e_format= (decpt <= 0 && width <= 2 - decpt && width >= 3 + exp_len);
301 ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
302 (len > 1 || !force_e_format)))) &&
309 (!have_space || (decpt >= -MAX_DECPT_FOR_F_FORMAT + 1 &&
310 (decpt <= MAX_DECPT_FOR_F_FORMAT || len > decpt))))
315 width-= (decpt < len) + (decpt <= 0 ? 1 - decpt : 0);
333 res= dtoa(x, 5, width - decpt, &decpt, &sign, &end);
349 if (sign && dst < dend)
355 if (len > 0 && dst < dend)
357 for (; decpt < 0 && dst < dend; decpt++)
361 for (i= 1; i <= len && dst < dend; i++)
364 if (i == decpt && i < len && dst < dend)
367 while (i++ <= decpt && dst < dend)
399 res= dtoa(x, 4, width, &decpt, &sign, &end);
409 if (sign && dst < dend)
413 if (len > 1 && dst < dend)
416 while (src < end && dst < dend)
421 if (decpt_sign && dst < dend)
424 if (decpt >= 100 && dst < dend)
426 *dst++= decpt / 100 +
'0';
429 *dst++= decpt / 10 +
'0';
431 else if (decpt >= 10 && dst < dend)
432 *dst++= decpt / 10 +
'0';
434 *dst++= decpt % 10 +
'0';
463 double my_strtod(
const char *str,
char **end,
int *error)
466 assert(str != NULL && end != NULL && *end != NULL && error != NULL);
468 res= my_strtod_int(str, end, error);
469 return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX);
473 double my_atof(
const char *nptr)
476 const char *end= nptr+65535;
477 return (my_strtod(nptr, (
char**) &end, &error));
535 typedef int32_t Long;
536 typedef uint32_t ULong;
537 typedef int64_t LLong;
538 typedef uint64_t ULLong;
540 typedef union {
double d; ULong L[2]; }
U;
542 #if defined(WORDS_BIGENDIAN) || (defined(__FLOAT_WORD_ORDER) && \
543 (__FLOAT_WORD_ORDER == __BIG_ENDIAN))
544 #define word0(x) ((U*)&x)->L[0]
545 #define word1(x) ((U*)&x)->L[1]
547 #define word0(x) ((U*)&x)->L[1]
548 #define word1(x) ((U*)&x)->L[0]
551 #define dval(x) ((U*)&x)->d
560 #define Exp_shift1 20
561 #define Exp_msk1 0x100000
562 #define Exp_mask 0x7ff00000
566 #define Exp_1 0x3ff00000
567 #define Exp_11 0x3ff00000
569 #define Frac_mask 0xfffff
570 #define Frac_mask1 0xfffff
573 #define Bndry_mask 0xfffff
574 #define Bndry_mask1 0xfffff
576 #define Sign_bit 0x80000000
584 #define Flt_Rounds FLT_ROUNDS
590 #define rounded_product(a,b) a*= b
591 #define rounded_quotient(a,b) a/= b
593 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
594 #define Big1 0xffffffff
595 #define FFFFFFFF 0xffffffffUL
612 dst->sign= src->sign;
615 assert(dst->maxwds >= src->wds);
617 memcpy(dst->p.x, src->p.x, src->wds*
sizeof(ULong));
622 static Bigint *Balloc(
int k)
625 Bigint* rv= (Bigint*) malloc(
sizeof(Bigint));
627 rv->p.x= (ULong*)malloc(x *
sizeof(ULong));
632 rv->sign= rv->wds= 0;
643 static void Bfree(Bigint *v)
655 static Bigint *multadd(Bigint *b,
int m,
int a)
668 y= *x * (ULLong)m + carry;
670 *x++= (ULong)(y & FFFFFFFF);
675 if (wds >= b->maxwds)
682 b->p.x[wds++]= (ULong) carry;
689 static Bigint *s2b(
const char *s,
int nd0,
int nd, ULong y9)
696 for (k= 0, y= 1; x > y; y <<= 1, k++) ;
706 b= multadd(b, 10, *s++ -
'0');
713 b= multadd(b, 10, *s++ -
'0');
718 static int hi0bits(ULong x)
722 if (!(x & 0xffff0000))
727 if (!(x & 0xff000000))
732 if (!(x & 0xf0000000))
737 if (!(x & 0xc0000000))
742 if (!(x & 0x80000000))
745 if (!(x & 0x40000000))
752 static int lo0bits(ULong *y)
804 static Bigint *i2b(
int i)
817 static Bigint *mult(Bigint *a, Bigint *b)
821 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
838 for (x= c->p.x, xa= x + wc; x < xa; x++)
845 for (; xb < xbe; xc0++)
854 z= *x++ * (ULLong)y + *xc + carry;
856 *xc++= (ULong) (z & FFFFFFFF);
862 for (xc0= c->p.x, xc= xc0 + wc; wc > 0 && !*--xc; --wc) ;
873 static ULong powers5[]=
881 2242703233UL, 762134875UL, 1262UL,
883 3211403009UL, 1849224548UL, 3668416493UL, 3913284084UL, 1593091UL,
885 781532673UL, 64985353UL, 253049085UL, 594863151UL, 3553621484UL,
886 3288652808UL, 3167596762UL, 2788392729UL, 3911132675UL, 590UL,
888 2553183233UL, 3201533787UL, 3638140786UL, 303378311UL, 1809731782UL,
889 3477761648UL, 3583367183UL, 649228654UL, 2915460784UL, 487929380UL,
890 1011012442UL, 1677677582UL, 3428152256UL, 1710878487UL, 1438394610UL,
891 2161952759UL, 4100910556UL, 1608314830UL, 349175UL
895 static Bigint p5_a[]=
898 { { powers5 }, 1, 1, 0, 1 },
899 { { powers5 + 1 }, 1, 1, 0, 1 },
900 { { powers5 + 2 }, 1, 2, 0, 2 },
901 { { powers5 + 4 }, 2, 3, 0, 3 },
902 { { powers5 + 7 }, 3, 5, 0, 5 },
903 { { powers5 + 12 }, 4, 10, 0, 10 },
904 { { powers5 + 22 }, 5, 19, 0, 19 }
907 #define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
909 static Bigint *pow5mult(Bigint *b,
int k)
911 Bigint *b1, *p5, *p51;
913 static int p05[3]= { 5, 25, 125 };
916 b= multadd(b, p05[i-1], 0);
932 if (p5 < p5_a + P5A_MAX)
934 else if (p5 == p5_a + P5A_MAX)
947 static Bigint *lshift(Bigint *b,
int k)
951 ULong *x, *x1, *xe, z;
956 for (i= b->maxwds; n1 > i; i<<= 1)
960 for (i= 0; i < n; i++)
987 static int cmp(Bigint *a, Bigint *b)
989 ULong *xa, *xa0, *xb, *xb0;
1003 return *xa < *xb ? -1 : 1;
1011 static Bigint *diff(Bigint *a, Bigint *b)
1015 ULong *xa, *xae, *xb, *xbe, *xc;
1047 y= (ULLong)*xa++ - *xb++ - borrow;
1048 borrow= y >> 32 & (ULong)1;
1049 *xc++= (ULong) (y & FFFFFFFF);
1055 borrow= y >> 32 & (ULong)1;
1056 *xc++= (ULong) (y & FFFFFFFF);
1065 static double ulp(
double x)
1070 L= (word0(x) & Exp_mask) - (P - 1)*Exp_msk1;
1077 static double b2d(Bigint *a,
int *e)
1079 ULong *xa, *xa0, w, y, z;
1092 d0= Exp_1 | y >> (Ebits - k);
1093 w= xa > xa0 ? *--xa : 0;
1094 d1= y << ((32-Ebits) + k) | w >> (Ebits - k);
1097 z= xa > xa0 ? *--xa : 0;
1100 d0= Exp_1 | y << k | z >> (32 - k);
1101 y= xa > xa0 ? *--xa : 0;
1102 d1= z << k | y >> (32 - k);
1116 static Bigint *d2b(
double d,
int *e,
int *bits)
1130 if ((de= (
int)(d0 >> Exp_shift)))
1134 if ((k= lo0bits(&y)))
1136 x[0]= y | z << (32 - k);
1141 i= b->wds= (x[1]= z) ? 2 : 1;
1152 *e= de - Bias - (P-1) + k;
1157 *e= de - Bias - (P-1) + 1 + k;
1158 *bits= 32*i - hi0bits(x[i-1]);
1166 static double ratio(Bigint *a, Bigint *b)
1171 dval(da)= b2d(a, &ka);
1172 dval(db)= b2d(b, &kb);
1173 k= ka - kb + 32*(a->wds - b->wds);
1175 word0(da)+= k*Exp_msk1;
1179 word0(db)+= k*Exp_msk1;
1181 return dval(da) / dval(db);
1184 static const double tens[] =
1186 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1187 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1191 static const double bigtens[]= { 1e16, 1e32, 1e64, 1e128, 1e256 };
1192 static const double tinytens[]=
1193 { 1e-16, 1e-32, 1e-64, 1e-128,
1194 9007199254740992.*9007199254740992.e-256
1200 #define Scale_Bit 0x10
1232 static double my_strtod_int(
const char *s00,
char **se,
int *error)
1235 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1236 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1237 const char *s, *s0, *s1, *end = *se;
1238 double aadj, aadj1, adj, rv, rv0;
1241 Bigint *bb= NULL, *bb1= NULL, *bd= NULL, *bd0= NULL, *bs= NULL, *delta= NULL;
1243 int inexact, oldinexact;
1251 for (s= s00; s < end; s++)
1276 while (++s < end && *s ==
'0') ;
1282 for (nd= nf= 0; s < end && (c= *s) >=
'0' && c <=
'9'; nd++, s++)
1288 if (s < end - 1 && c ==
'.')
1293 for (; s < end && c ==
'0'; c= *++s)
1295 if (s < end && c >
'0' && c <=
'9')
1304 for (; s < end && c >=
'0' && c <=
'9'; c = *++s)
1311 for (i= 1; i < nz; i++)
1314 else if (nd <= DBL_DIG + 1)
1318 else if (nd <= DBL_DIG + 1)
1326 if (s < end && (c ==
'e' || c ==
'E'))
1328 if (!nd && !nz && !nz0)
1339 if (s < end && c >=
'0' && c <=
'9')
1341 while (s < end && c ==
'0')
1343 if (s < end && c >
'0' && c <=
'9') {
1346 while (++s < end && (c= *s) >=
'0' && c <=
'9')
1348 if (s - s1 > 8 || L > 19999)
1385 k= nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1391 oldinexact = get_inexact();
1393 dval(rv)= tens[k - 9] * dval(rv) + z;
1405 rounded_product(dval(rv), tens[e]);
1409 if (e <= Ten_pmax + i)
1417 rounded_product(dval(rv), tens[e]);
1421 #ifndef Inaccurate_Divide
1422 else if (e >= -Ten_pmax)
1424 rounded_quotient(dval(rv), tens[-e]);
1434 oldinexact= get_inexact();
1446 if (e1 > DBL_MAX_10_EXP)
1451 word0(rv)= Exp_mask;
1456 dval(rv0)*= dval(rv0);
1463 for(j= 0; e1 > 1; j++, e1>>= 1)
1465 dval(rv)*= bigtens[j];
1467 word0(rv)-= P*Exp_msk1;
1468 dval(rv)*= bigtens[j];
1469 if ((z= word0(rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - P))
1471 if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P))
1478 word0(rv)+= P*Exp_msk1;
1488 if (e1 >= 1 << n_bigtens)
1492 for(j= 0; e1 > 0; j++, e1>>= 1)
1494 dval(rv)*= tinytens[j];
1495 if (scale && (j = 2 * P + 1 - ((word0(rv) & Exp_mask) >> Exp_shift)) > 0)
1502 word0(rv)= (P + 2) * Exp_msk1;
1504 word0(rv)&= 0xffffffff << (j - 32);
1507 word1(rv)&= 0xffffffff << j;
1524 bd0= s2b(s0, nd0, nd, y);
1530 bb= d2b(dval(rv), &bbe, &bbbits);
1557 i= bb2 < bd2 ? bb2 : bd2;
1568 bs= pow5mult(bs, bb5);
1574 bb= lshift(bb, bb2);
1576 bd= pow5mult(bd, bd5);
1578 bd= lshift(bd, bd2);
1580 bs= lshift(bs, bs2);
1581 delta= diff(bb, bd);
1592 if (dsign || word1(rv) || word0(rv) & Bndry_mask ||
1593 (word0(rv) & Exp_mask) <= (2 * P + 1) * Exp_msk1)
1596 if (!delta->x[0] && delta->wds <= 1)
1601 if (!delta->p.x[0] && delta->wds <= 1)
1609 delta= lshift(delta, Log2P);
1610 if (cmp(delta, bs) > 0)
1619 if ((word0(rv) & Bndry_mask1) == Bndry_mask1 &&
1621 ((scale && (y = word0(rv) & Exp_mask) <= 2 * P * Exp_msk1) ?
1622 (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
1626 word0(rv)= (word0(rv) & Exp_mask) + Exp_msk1;
1632 else if (!(word0(rv) & Bndry_mask) && !word1(rv))
1638 L= word0(rv) & Exp_mask;
1639 if (L <= (2 *P + 1) * Exp_msk1)
1641 if (L > (P + 2) * Exp_msk1)
1648 L= (word0(rv) & Exp_mask) - Exp_msk1;
1649 word0(rv)= L | Bndry_mask1;
1650 word1(rv)= 0xffffffff;
1653 if (!(word1(rv) & LSB))
1656 dval(rv)+= ulp(dval(rv));
1659 dval(rv)-= ulp(dval(rv));
1666 if ((aadj= ratio(delta, bs)) <= 2.)
1670 else if (word1(rv) || word0(rv) & Bndry_mask)
1672 if (word1(rv) == Tiny1 && !word0(rv))
1680 if (aadj < 2. / FLT_RADIX)
1681 aadj= 1. / FLT_RADIX;
1690 aadj1= dsign ? aadj : -aadj;
1691 if (Flt_Rounds == 0)
1694 y= word0(rv) & Exp_mask;
1698 if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
1700 dval(rv0)= dval(rv);
1701 word0(rv)-= P * Exp_msk1;
1702 adj= aadj1 * ulp(dval(rv));
1704 if ((word0(rv) & Exp_mask) >= Exp_msk1 * (DBL_MAX_EXP + Bias - P))
1706 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1713 word0(rv)+= P * Exp_msk1;
1717 if (scale && y <= 2 * P * Exp_msk1)
1719 if (aadj <= 0x7fffffff)
1721 if ((z= (ULong) aadj) <= 0)
1724 aadj1= dsign ? aadj : -aadj;
1726 word0(aadj1)+= (2 * P + 1) * Exp_msk1 - y;
1728 adj = aadj1 * ulp(dval(rv));
1731 z= word0(rv) & Exp_mask;
1740 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1742 if (aadj < .4999999 || aadj > .5000001)
1745 else if (aadj < .4999999 / FLT_RADIX)
1760 word0(rv0)= Exp_1 + (70 << Exp_shift);
1765 else if (!oldinexact)
1770 word0(rv0)= Exp_1 - 2 * P * Exp_msk1;
1772 dval(rv)*= dval(rv0);
1775 if (inexact && !(word0(rv) & Exp_mask))
1779 dval(rv0)*= dval(rv0);
1790 return sign ? -dval(rv) : dval(rv);
1794 static int quorem(Bigint *b, Bigint *S)
1797 ULong *bx, *bxe, q, *sx, *sxe;
1798 ULLong borrow, carry, y, ys;
1807 q= *bxe / (*sxe + 1);
1814 ys= *sx++ * (ULLong)q + carry;
1816 y= *bx - (ys & FFFFFFFF) - borrow;
1817 borrow= y >> 32 & (ULong)1;
1818 *bx++= (ULong) (y & FFFFFFFF);
1824 while (--bxe > bx && !*bxe)
1840 y= *bx - (ys & FFFFFFFF) - borrow;
1841 borrow= y >> 32 & (ULong)1;
1842 *bx++= (ULong) (y & FFFFFFFF);
1849 while (--bxe > bx && !*bxe)
1893 static char *dtoa(
double d,
int mode,
int ndigits,
int *decpt,
int *sign,
1929 int bbits, b2, b5, be, dig, i, ieps, ilim=0, ilim0, ilim1= 0,
1930 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1931 spec_case, try_quick;
1935 Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
1939 if (word0(d) & Sign_bit)
1943 word0(d) &= ~Sign_bit;
1949 if ((((word0(d) & Exp_mask) == Exp_mask) && ((*decpt= DTOA_OVERFLOW) != 0)) ||
1950 (!dval(d) && ((*decpt= 1) != 0)))
1953 char *res= (
char*) malloc(2);
1962 b= d2b(dval(d), &be, &bbits);
1963 if ((i= (
int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
1966 word0(d2) &= Frac_mask1;
1967 word0(d2) |= Exp_11;
1999 i= bbits + be + (Bias + (P-1) - 1);
2000 x= i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
2001 : word1(d) << (32 - i);
2003 word0(d2)-= 31*Exp_msk1;
2004 i-= (Bias + (P-1) - 1) + 1;
2007 ds= (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2009 if (ds < 0. && ds != k)
2012 if (k >= 0 && k <= Ten_pmax)
2014 if (dval(d) < tens[k])
2041 if (mode < 0 || mode > 9)
2065 ilim= ilim1= i= ndigits;
2077 s= s0= (
char*)malloc(i+1);
2080 if (ilim >= 0 && ilim <= Quick_max && try_quick)
2096 dval(d)/= bigtens[n_bigtens-1];
2099 for (; j; j>>= 1, i++)
2111 dval(d)*= tens[j1 & 0xf];
2112 for (j= j1 >> 4; j; j>>= 1, i++)
2117 dval(d)*= bigtens[i];
2121 if (k_check && dval(d) < 1. && ilim > 0)
2130 dval(eps)= ieps*dval(d) + 7.;
2131 word0(eps)-= (P-1)*Exp_msk1;
2136 if (dval(d) > dval(eps))
2138 if (dval(d) < -dval(eps))
2145 dval(eps)= 0.5/tens[ilim-1] - dval(eps);
2151 if (dval(d) < dval(eps))
2153 if (1. - dval(d) < dval(eps))
2164 dval(eps)*= tens[ilim-1];
2165 for (i= 1;; i++, dval(d)*= 10.)
2173 if (dval(d) > 0.5 + dval(eps))
2175 else if (dval(d) < 0.5 - dval(eps))
2177 while (*--s ==
'0') {}
2194 if (be >= 0 && k <= Int_max)
2198 if (ndigits < 0 && ilim <= 0)
2201 if (ilim < 0 || dval(d) <= 5*ds)
2205 for (i= 1;; i++, dval(d)*= 10.)
2207 L= (Long)(dval(d) / ds);
2217 if (dval(d) > ds || (dval(d) == ds && L & 1))
2240 i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
2245 if (m2 > 0 && s2 > 0)
2247 i= m2 < s2 ? m2 : s2;
2258 mhi= pow5mult(mhi, m5);
2276 if ((mode < 2 || leftright)
2279 if (!word1(d) && !(word0(d) & Bndry_mask) &&
2280 word0(d) & (Exp_mask & ~Exp_msk1)
2298 if ((i= ((s5 ? 32 - hi0bits(S->p.x[S->wds-1]) : 1) + s2) & 0x1f))
2324 b= multadd(b, 10, 0);
2326 mhi= multadd(mhi, 10, 0);
2330 if (ilim <= 0 && (mode == 3 || mode == 5))
2332 if (ilim < 0 || cmp(b,S= multadd(S,5,0)) <= 0)
2347 mhi= lshift(mhi, m2);
2356 mhi= Balloc(mhi->k);
2358 mhi= lshift(mhi, Log2P);
2363 dig= quorem(b,S) +
'0';
2366 delta= diff(S, mhi);
2367 j1= delta->sign ? 1 : cmp(b, delta);
2369 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
2379 if (j < 0 || (j == 0 && mode != 1 && !(word1(d) & 1)))
2381 if (!b->p.x[0] && b->wds <= 1)
2389 if ((j1 > 0 || (j1 == 0 && dig & 1))
2411 b= multadd(b, 10, 0);
2413 mlo= mhi= multadd(mhi, 10, 0);
2416 mlo= multadd(mlo, 10, 0);
2417 mhi= multadd(mhi, 10, 0);
2424 *s++= dig= quorem(b,S) +
'0';
2425 if (!b->p.x[0] && b->wds <= 1)
2431 b= multadd(b, 10, 0);
2438 if (j > 0 || (j == 0 && dig & 1))
2452 while (*--s ==
'0') {}
2459 if (mlo && mlo != mhi)