mysql中 decimal 的实现_mysql decimal实现

mysql中定义deciaml

/**
   intg is the number of *decimal* digits (NOT number of decimal_digit_t's !)
   before the point
   frac is the number of decimal digits after the point
   len  is the length of buf (length of allocated space) in decimal_digit_t's,
   not in bytes
   sign false means positive, true means negative
   buf  is an array of decimal_digit_t's
*/
typedef struct st_decimal_t {
 int intg;
 int frac;
 int len;
 my_bool sign;
 decimal_digit_t *buf;
} decimal_t;

具体实现感觉是是一个1000000000为基数的定点浮点数。

 

从中扣出来的代码实现如下：

 

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <fcntl.h>
#include <malloc.h>
#include <assert.h>

#ifndef rdtsc
//#define rdtscfreq 2.1280493e9
#define rdtscfreq 1
#define rdtsc(low,high)          \
 __asm__ __volatile__("rdtsc" : "=a" (low), "=d" (high))

static inline  double rdtscdiff(unsigned int low1,unsigned int high1,unsigned int low2,unsigned int high2)
{
 const  double MAXintplus1 = 4294967296.0;
 unsigned int low,high;
 if (low2<low1) {
  low = 0xffffffff - low1 + low2;
  high = high2 - high1 -1;
 }
 else {
  high = high2 - high1;
  low = low2 - low1;
 }
// printf("%d,",low);
 return (double)high * MAXintplus1 + (double)low ;
}
static inline double mdtime(int id)
{
 static unsigned int high0,low0,high1,low1;
 if(id){
  rdtsc(low1,high1);
  return rdtscdiff(low0,high0,low1,high1)/rdtscfreq; 
 }
 else{
  rdtsc(low0,high0);
  return 0.0;
 }
}

#endif

#define E_DEC_OK                0
#define E_DEC_TRUNCATED         1
#define E_DEC_OVERFLOW          2
#define E_DEC_DIV_ZERO          4
#define E_DEC_BAD_NUM           8
#define E_DEC_OOM              16

#define E_DEC_ERROR            31
#define E_DEC_FATAL_ERROR      30

typedef int int32;
typedef enum{
 TRUNCATE=0,
 HALF_EVEN,
 HALF_UP,
 CEILING,
 FLOOR
}decimal_round_mode;
typedef int  my_bool; /* Small bool */
typedef int32 decimal_digit_t;

typedef unsigned long long int ulonglong; /* ulong or unsigned long long */
typedef long long int longlong;

/**
   intg is the number of *decimal* digits (NOT number of decimal_digit_t's !)
   before the point
   frac is the number of decimal digits after the point
   len  is the length of buf (length of allocated space) in decimal_digit_t's,
   not in bytes
   sign false means positive, true means negative
   buf  is an array of decimal_digit_t's
*/
typedef struct st_decimal_t {
 int intg;
 int frac;
 int len;
 my_bool sign;
 decimal_digit_t *buf;
} decimal_t;

typedef decimal_digit_t dec1;
typedef longlong      dec2;

#define DIG_PER_DEC1 9
#define DIG_MASK     100000000
#define DIG_BASE     1000000000
#define DIG_MAX      (DIG_BASE-1)
#define DIG_BASE2    ((dec2)DIG_BASE * (dec2)DIG_BASE)
#define ROUND_UP(X)  (((X)+DIG_PER_DEC1-1)/DIG_PER_DEC1)

static const dec1 powers10[DIG_PER_DEC1+1]={
 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000};
#define DBUG_ASSERT(A) assert(A)
static const int dig2bytes[DIG_PER_DEC1+1]={0, 1, 1, 2, 2, 3, 3, 4, 4, 4};
static const dec1 frac_max[DIG_PER_DEC1-1]={
 900000000, 990000000, 999000000,
 999900000, 999990000, 999999000,
 999999900, 999999990 };
static double scaler10[]= {
 1.0, 1e10, 1e20, 1e30, 1e40, 1e50, 1e60, 1e70, 1e80, 1e90
};
static double scaler1[]= {
 1.0, 10.0, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9
};

#define sanity(d) DBUG_ASSERT((d)->len >0 && ((d)->buf[0] |    \
             (d)->buf[(d)->len-1] | 1))

#define FIX_INTG_FRAC_ERROR(len, intg1, frac1, error) \
 do             \
 {             \
  if (unlikely(intg1+frac1 > (len)))    \
  {            \
            if (unlikely(intg1 > (len)))    \
            {           \
    intg1=(len);       \
    frac1=0;        \
    error=E_DEC_OVERFLOW;     \
            }           \
            else          \
            {           \
    frac1=(len)-intg1;      \
    error=E_DEC_TRUNCATED;     \
            }           \
  }            \
  else           \
            error=E_DEC_OK;        \
 } while(0)

#define ADD(to, from1, from2, carry)  /* assume carry <= 1 */ \
 do               \
 {               \
  dec1 a=(from1)+(from2)+(carry);       \
  DBUG_ASSERT((carry) <= 1);        \
  if (((carry)= a >= DIG_BASE)) /* no division here! */ \
            a-=DIG_BASE;          \
  (to)=a;             \
 } while(0)

#define ADD2(to, from1, from2, carry)   \
 do           \
 {           \
  dec2 a=((dec2)(from1))+(from2)+(carry); \
  if (((carry)= a >= DIG_BASE))   \
            a-=DIG_BASE;      \
  if (unlikely(a >= DIG_BASE))   \
  {          \
            a-=DIG_BASE;      \
            carry++;       \
  }          \
  (to)=(dec1) a;       \
 } while(0)

#define SUB(to, from1, from2, carry) /* to=from1-from2 */ \
 do              \
 {              \
  dec1 a=(from1)-(from2)-(carry);      \
  if (((carry)= a < 0))        \
            a+=DIG_BASE;         \
  (to)=a;            \
 } while(0)

#define SUB2(to, from1, from2, carry) /* to=from1-from2 */ \
 do              \
 {              \
  dec1 a=(from1)-(from2)-(carry);      \
  if (((carry)= a < 0))        \
            a+=DIG_BASE;         \
  if (unlikely(a < 0))        \
  {             \
            a+=DIG_BASE;         \
            carry++;          \
  }             \
  (to)=a;            \
 } while(0)

#define decimal_make_zero(dec)        do {  \
  (dec)->buf[0]=0;      \
  (dec)->intg=1;       \
  (dec)->frac=0;       \
  (dec)->sign=0;       \
 } while(0)

#define swap_variables(t, a, b) { t dummy; dummy= a; a= b; b= dummy; }
#define set_if_bigger(a,b)  do { if ((a) < (b)) (a)=(b); } while(0)
#define set_if_smaller(a,b) do { if ((a) > (b)) (a)=(b); } while(0)

static dec1 *remove_leading_zeroes(decimal_t *from, int *intg_result)
{
 int intg= from->intg, i;
 dec1 *buf0= from->buf;
 i= ((intg - 1) % DIG_PER_DEC1) + 1;
 while (intg > 0 && *buf0 == 0)
 {
  intg-= i;
  i= DIG_PER_DEC1;
  buf0++;
 }
 if (intg > 0)
 {
  for (i= (intg - 1) % DIG_PER_DEC1; *buf0 < powers10[i--]; intg--) ;
  DBUG_ASSERT(intg > 0);
 }
 else
  intg=0;
 *intg_result= intg;
 return buf0;
}

int decimal_intg(decimal_t *from)
{
 int res;
 remove_leading_zeroes(from, &res);
 return res;
}

#define __builtin_expect(x, expected_value) (x)

#define likely(x) __builtin_expect((x),1)
#define unlikely(x) __builtin_expect((x),0)

 

#define max(a, b) ((a) > (b) ? (a) : (b))
#define min(a, b) ((a) < (b) ? (a) : (b))

/*
  Get maximum value for given precision and scale

  SYNOPSIS
  max_decimal()
  precision/scale - see decimal_bin_size() below
  to              - decimal where where the result will be stored
  to->buf and to->len must be set.
*/

void max_decimal(int precision, int frac, decimal_t *to)
{
 int intpart;
 dec1 *buf= to->buf;
 DBUG_ASSERT(precision && precision >= frac);

 to->sign= 0;
 if ((intpart= to->intg= (precision - frac)))
 {
  int firstdigits= intpart % DIG_PER_DEC1;
  if (firstdigits)
   *buf++= powers10[firstdigits] - 1; /* get 9 99 999 ... */
  for(intpart/= DIG_PER_DEC1; intpart; intpart--)
   *buf++= DIG_MAX;
 }

 if ((to->frac= frac))
 {
  int lastdigits= frac % DIG_PER_DEC1;
  for(frac/= DIG_PER_DEC1; frac; frac--)
   *buf++= DIG_MAX;
  if (lastdigits)
   *buf= frac_max[lastdigits - 1];
 }
}

/* to=from1-from2.
   if to==0, return -1/0/+1 - the result of the comparison */
static int do_sub(decimal_t *from1, decimal_t *from2, decimal_t *to)
{
 int intg1=ROUND_UP(from1->intg), intg2=ROUND_UP(from2->intg),
  frac1=ROUND_UP(from1->frac), frac2=ROUND_UP(from2->frac);
 int frac0=max(frac1, frac2), error;
 dec1 *buf1, *buf2, *buf0, *stop1, *stop2, *start1, *start2, carry=0;

 /* let carry:=1 if from2 > from1 */
 start1=buf1=from1->buf; stop1=buf1+intg1;
 start2=buf2=from2->buf; stop2=buf2+intg2;
 if (unlikely(*buf1 == 0))
 {
  while (buf1 < stop1 && *buf1 == 0)
   buf1++;
  start1=buf1;
  intg1= (int) (stop1-buf1);
 }
 if (unlikely(*buf2 == 0))
 {
  while (buf2 < stop2 && *buf2 == 0)
   buf2++;
  start2=buf2;
  intg2= (int) (stop2-buf2);
 }
 if (intg2 > intg1)
  carry=1;
 else if (intg2 == intg1)
 {
  dec1 *end1= stop1 + (frac1 - 1);
  dec1 *end2= stop2 + (frac2 - 1);
  while (unlikely((buf1 <= end1) && (*end1 == 0)))
   end1--;
  while (unlikely((buf2 <= end2) && (*end2 == 0)))
   end2--;
  frac1= (int) (end1 - stop1) + 1;
  frac2= (int) (end2 - stop2) + 1;
  while (buf1 <=end1 && buf2 <= end2 && *buf1 == *buf2)
   buf1++, buf2++;
  if (buf1 <= end1)
  {
   if (buf2 <= end2)
    carry= *buf2 > *buf1;
   else
    carry= 0;
  }
  else
  {
   if (buf2 <= end2)
    carry=1;
   else /* short-circuit everything: from1 == from2 */
   {
    if (to == 0) /* decimal_cmp() */
     return 0;
    decimal_make_zero(to);
    return E_DEC_OK;
   }
  }
 }

 if (to == 0) /* decimal_cmp() */
  return carry == from1->sign ? 1 : -1;

 sanity(to);

 to->sign=from1->sign;

 /* ensure that always from1 > from2 (and intg1 >= intg2) */
 if (carry)
 {
  swap_variables(decimal_t *,from1,from1);
  swap_variables(dec1 *,start1, start2);
  swap_variables(int,intg1,intg2);
  swap_variables(int,frac1,frac2);
  to->sign= 1 - to->sign;
 }

 FIX_INTG_FRAC_ERROR(to->len, intg1, frac0, error);
 buf0=to->buf+intg1+frac0;

 to->frac=max(from1->frac, from2->frac);
 to->intg=intg1*DIG_PER_DEC1;
 if (unlikely(error))
 {
  set_if_smaller(to->frac, frac0*DIG_PER_DEC1);
  set_if_smaller(frac1, frac0);
  set_if_smaller(frac2, frac0);
  set_if_smaller(intg2, intg1);
 }
 carry=0;

 /* part 1 - max(frac) ... min (frac) */
 if (frac1 > frac2)
 {
  buf1=start1+intg1+frac1;
  stop1=start1+intg1+frac2;
  buf2=start2+intg2+frac2;
  while (frac0-- > frac1)
   *--buf0=0;
  while (buf1 > stop1)
   *--buf0=*--buf1;
 }
 else
 {
  buf1=start1+intg1+frac1;
  buf2=start2+intg2+frac2;
  stop2=start2+intg2+frac1;
  while (frac0-- > frac2)
   *--buf0=0;
  while (buf2 > stop2)
  {
   SUB(*--buf0, 0, *--buf2, carry);
  }
 }

 /* part 2 - min(frac) ... intg2 */
 while (buf2 > start2)
 {
  SUB(*--buf0, *--buf1, *--buf2, carry);
 }

 /* part 3 - intg2 ... intg1 */
 while (carry && buf1 > start1)
 {
  SUB(*--buf0, *--buf1, 0, carry);
 }

 while (buf1 > start1)
  *--buf0=*--buf1;

 while (buf0 > to->buf)
  *--buf0=0;

 return error;
}

static int do_add(decimal_t *from1, decimal_t *from2, decimal_t *to)
{
 int intg1=ROUND_UP(from1->intg), intg2=ROUND_UP(from2->intg),
  frac1=ROUND_UP(from1->frac), frac2=ROUND_UP(from2->frac),
  frac0=max(frac1, frac2), intg0=max(intg1, intg2), error;
 dec1 *buf1, *buf2, *buf0, *stop, *stop2, x, carry;

 sanity(to);

 /* is there a need for extra word because of carry ? */
 x=intg1 > intg2 ? from1->buf[0] :
  intg2 > intg1 ? from2->buf[0] :
  from1->buf[0] + from2->buf[0] ;
 if (unlikely(x > DIG_MAX-1)) /* yes, there is */
 {
  intg0++;
  to->buf[0]=0; /* safety */
 }

 FIX_INTG_FRAC_ERROR(to->len, intg0, frac0, error);
 if (unlikely(error == E_DEC_OVERFLOW))
 {
  max_decimal(to->len * DIG_PER_DEC1, 0, to);
  return error;
 }

 buf0=to->buf+intg0+frac0;

 to->sign=from1->sign;
 to->frac=max(from1->frac, from2->frac);
 to->intg=intg0*DIG_PER_DEC1;
 if (unlikely(error))
 {
  set_if_smaller(to->frac, frac0*DIG_PER_DEC1);
  set_if_smaller(frac1, frac0);
  set_if_smaller(frac2, frac0);
  set_if_smaller(intg1, intg0);
  set_if_smaller(intg2, intg0);
 }

 /* part 1 - max(frac) ... min (frac) */
 if (frac1 > frac2)
 {
  buf1=from1->buf+intg1+frac1;
  stop=from1->buf+intg1+frac2;
  buf2=from2->buf+intg2+frac2;
  stop2=from1->buf+(intg1 > intg2 ? intg1-intg2 : 0);
 }
 else
 {
  buf1=from2->buf+intg2+frac2;
  stop=from2->buf+intg2+frac1;
  buf2=from1->buf+intg1+frac1;
  stop2=from2->buf+(intg2 > intg1 ? intg2-intg1 : 0);
 }
 while (buf1 > stop)
  *--buf0=*--buf1;

 /* part 2 - min(frac) ... min(intg) */
 carry=0;
 while (buf1 > stop2)
 {
  ADD(*--buf0, *--buf1, *--buf2, carry);
 }

 /* part 3 - min(intg) ... max(intg) */
 buf1= intg1 > intg2 ? ((stop=from1->buf)+intg1-intg2) :
  ((stop=from2->buf)+intg2-intg1) ;
 while (buf1 > stop)
 {
  ADD(*--buf0, *--buf1, 0, carry);
 }

 if (unlikely(carry))
  *--buf0=1;
 DBUG_ASSERT(buf0 == to->buf || buf0 == to->buf+1);

 return error;
}

int decimal_add(decimal_t *from1, decimal_t *from2, decimal_t *to)
{
 if (likely(from1->sign == from2->sign))
  return do_add(from1, from2, to);
 return do_sub(from1, from2, to);
}
int decimal_sub(decimal_t *from1, decimal_t *from2, decimal_t *to)
{
 if (likely(from1->sign == from2->sign))
  return do_sub(from1, from2, to);
 return do_add(from1, from2, to);
}
/*
  multiply two decimals

  SYNOPSIS
  decimal_mul()
  from1, from2 - factors
  to      - product

  RETURN VALUE
  E_DEC_OK/E_DEC_TRUNCATED/E_DEC_OVERFLOW;

  NOTES
  in this implementation, with sizeof(dec1)=4 we have DIG_PER_DEC1=9,
  and 63-digit number will take only 7 dec1 words (basically a 7-digit
  "base 999999999" number).  Thus there's no need in fast multiplication
  algorithms, 7-digit numbers can be multiplied with a naive O(n*n)
  method.

  XXX if this library is to be used with huge numbers of thousands of
  digits, fast multiplication must be implemented.
*/
int decimal_mul(decimal_t *from1, decimal_t *from2, decimal_t *to)
{
 int intg1=ROUND_UP(from1->intg), intg2=ROUND_UP(from2->intg),
  frac1=ROUND_UP(from1->frac), frac2=ROUND_UP(from2->frac),
  intg0=ROUND_UP(from1->intg+from2->intg),
  frac0=frac1+frac2, error, i, j, d_to_move;
 dec1 *buf1=from1->buf+intg1, *buf2=from2->buf+intg2, *buf0,
  *start2, *stop2, *stop1, *start0, carry;

 sanity(to);

 i=intg0;                                       /* save 'ideal' values */
 j=frac0;
 FIX_INTG_FRAC_ERROR(to->len, intg0, frac0, error);  /* bound size */
 to->sign=from1->sign != from2->sign;
 to->frac=from1->frac+from2->frac;              /* store size in digits */
 to->intg=intg0*DIG_PER_DEC1;

 if (unlikely(error))
 {
  set_if_smaller(to->frac, frac0*DIG_PER_DEC1);
  set_if_smaller(to->intg, intg0*DIG_PER_DEC1);
  if (unlikely(i > intg0))                     /* bounded integer-part */
  {
   i-=intg0;
   j=i >> 1;
   intg1-= j;
   intg2-=i-j;
   frac1=frac2=0; /* frac0 is already 0 here */
  }
  else                                         /* bounded fract part */
  {
   j-=frac0;
   i=j >> 1;
   if (frac1 <= frac2)
   {
    frac1-= i;
    frac2-=j-i;
   }
   else
   {
    frac2-= i;
    frac1-=j-i;
   }
  }
 }
 start0=to->buf+intg0+frac0-1;
 start2=buf2+frac2-1;
 stop1=buf1-intg1;
 stop2=buf2-intg2;

 bzero(to->buf, (intg0+frac0)*sizeof(dec1));

 for (buf1+=frac1-1; buf1 >= stop1; buf1--, start0--)
 {
  carry=0;
  for (buf0=start0, buf2=start2; buf2 >= stop2; buf2--, buf0--)
  {
   dec1 hi, lo;
   dec2 p= ((dec2)*buf1) * ((dec2)*buf2);
   hi=(dec1)(p/DIG_BASE);
   lo=(dec1)(p-((dec2)hi)*DIG_BASE);
   ADD2(*buf0, *buf0, lo, carry);
   carry+=hi;
  }
  if (carry)
  {
   if (buf0 < to->buf)
    return E_DEC_OVERFLOW;
   ADD2(*buf0, *buf0, 0, carry);
  }
  for (buf0--; carry; buf0--)
  {
   if (buf0 < to->buf)
    return E_DEC_OVERFLOW;
   ADD(*buf0, *buf0, 0, carry);
  }
 }

 /* Now we have to check for -0.000 case */
 if (to->sign)
 {
  dec1 *buf= to->buf;
  dec1 *end= to->buf + intg0 + frac0;
  DBUG_ASSERT(buf != end);
  for (;;)
  {
   if (*buf)
    break;
   if (++buf == end)
   {
    /* We got decimal zero */
    decimal_make_zero(to);
    break;
   }
  }
 }
 buf1= to->buf;
 d_to_move= intg0 + ROUND_UP(to->frac);
 while (!*buf1 && (to->intg > DIG_PER_DEC1))
 {
  buf1++;
  to->intg-= DIG_PER_DEC1;
  d_to_move--;
 }
 if (to->buf < buf1)
 {
  dec1 *cur_d= to->buf;
  for (; d_to_move--; cur_d++, buf1++)
   *cur_d= *buf1;
 }
 return error;
}

/*
  naive division algorithm (Knuth's Algorithm D in 4.3.1) -
  it's ok for short numbers
  also we're using alloca() to allocate a temporary buffer

  XXX if this library is to be used with huge numbers of thousands of
  digits, fast division must be implemented and alloca should be
  changed to malloc (or at least fallback to malloc if alloca() fails)
  but then, decimal_mul() should be rewritten too :(
*/
#define UNINIT_VAR(x) x= x

#define my_alloca(SZ) alloca((size_t) (SZ))
#define my_afree(PTR) {}

static int do_div_mod(decimal_t *from1, decimal_t *from2,
       decimal_t *to, decimal_t *mod, int scale_incr)
{
 int frac1=ROUND_UP(from1->frac)*DIG_PER_DEC1, prec1=from1->intg+frac1,
  frac2=ROUND_UP(from2->frac)*DIG_PER_DEC1, prec2=from2->intg+frac2,
  UNINIT_VAR(error), i, intg0, frac0, len1, len2, dintg, div_mod=(!mod);
 dec1 *buf0, *buf1=from1->buf, *buf2=from2->buf, *tmp1,
  *start2, *stop2, *stop1, *stop0, norm2, carry, *start1, dcarry;
 dec2 norm_factor, x, guess, y;

 if (mod)
  to=mod;

 sanity(to);

 /* removing all the leading zeroes */
 i= ((prec2 - 1) % DIG_PER_DEC1) + 1;
 while (prec2 > 0 && *buf2 == 0)
 {
  prec2-= i;
  i= DIG_PER_DEC1;
  buf2++;
 }
 if (prec2 <= 0) /* short-circuit everything: from2 == 0 */
  return E_DEC_DIV_ZERO;
 for (i= (prec2 - 1) % DIG_PER_DEC1; *buf2 < powers10[i--]; prec2--) ;
 DBUG_ASSERT(prec2 > 0);

 i=((prec1-1) % DIG_PER_DEC1)+1;
 while (prec1 > 0 && *buf1 == 0)
 {
  prec1-=i;
  i=DIG_PER_DEC1;
  buf1++;
 }
 if (prec1 <= 0)
 { /* short-circuit everything: from1 == 0 */
  decimal_make_zero(to);
  return E_DEC_OK;
 }
 for (i=(prec1-1) % DIG_PER_DEC1; *buf1 < powers10[i--]; prec1--) ;
 DBUG_ASSERT(prec1 > 0);

 /* let's fix scale_incr, taking into account frac1,frac2 increase */
 if ((scale_incr-= frac1 - from1->frac + frac2 - from2->frac) < 0)
  scale_incr=0;

 dintg=(prec1-frac1)-(prec2-frac2)+(*buf1 >= *buf2);
 if (dintg < 0)
 {
  dintg/=DIG_PER_DEC1;
  intg0=0;
 }
 else
  intg0=ROUND_UP(dintg);
 if (mod)
 {
  /* we're calculating N1 % N2.
     The result will have
     frac=max(frac1, frac2), as for subtraction
     intg=intg2
  */
  to->sign=from1->sign;
  to->frac=max(from1->frac, from2->frac);
  frac0=0;
 }
 else
 {
  /*
    we're calculating N1/N2. N1 is in the buf1, has prec1 digits
    N2 is in the buf2, has prec2 digits. Scales are frac1 and
    frac2 accordingly.
    Thus, the result will have
    frac = ROUND_UP(frac1+frac2+scale_incr)
    and
    intg = (prec1-frac1) - (prec2-frac2) + 1
    prec = intg+frac
  */
  frac0=ROUND_UP(frac1+frac2+scale_incr);
  FIX_INTG_FRAC_ERROR(to->len, intg0, frac0, error);
  to->sign=from1->sign != from2->sign;
  to->intg=intg0*DIG_PER_DEC1;
  to->frac=frac0*DIG_PER_DEC1;
 }
 buf0=to->buf;
 stop0=buf0+intg0+frac0;
 if (likely(div_mod))
  while (dintg++ < 0)
   *buf0++=0;

 len1=(i=ROUND_UP(prec1))+ROUND_UP(2*frac2+scale_incr+1) + 1;
 set_if_bigger(len1, 3);
 if (!(tmp1=(dec1 *)my_alloca(len1*sizeof(dec1))))
  return E_DEC_OOM;
 memcpy(tmp1, buf1, i*sizeof(dec1));
 bzero(tmp1+i, (len1-i)*sizeof(dec1));

 start1=tmp1;
 stop1=start1+len1;
 start2=buf2;
 stop2=buf2+ROUND_UP(prec2)-1;

 /* removing end zeroes */
 while (*stop2 == 0 && stop2 >= start2)
  stop2--;
 len2= (int) (stop2++ - start2);

 /*
   calculating norm2 (normalized *start2) - we need *start2 to be large
   (at least > DIG_BASE/2), but unlike Knuth's Alg. D we don't want to
   normalize input numbers (as we don't make a copy of the divisor).
   Thus we normalize first dec1 of buf2 only, and we'll normalize *start1
   on the fly for the purpose of guesstimation only.
   It's also faster, as we're saving on normalization of buf2
 */
 norm_factor=DIG_BASE/(*start2+1);
 norm2=(dec1)(norm_factor*start2[0]);
 if (likely(len2>0))
  norm2+=(dec1)(norm_factor*start2[1]/DIG_BASE);

 if (*start1 < *start2)
  dcarry=*start1++;
 else
  dcarry=0;

 /* main loop */
 for (; buf0 < stop0; buf0++)
 {
  /* short-circuit, if possible */
  if (unlikely(dcarry == 0 && *start1 < *start2))
   guess=0;
  else
  {
   /* D3: make a guess */
   x=start1[0]+((dec2)dcarry)*DIG_BASE;
   y=start1[1];
   guess=(norm_factor*x+norm_factor*y/DIG_BASE)/norm2;
   if (unlikely(guess >= DIG_BASE))
    guess=DIG_BASE-1;
   if (likely(len2>0))
   {
    /* hmm, this is a suspicious trick - I removed normalization here */
    if (start2[1]*guess > (x-guess*start2[0])*DIG_BASE+y)
     guess--;
    if (unlikely(start2[1]*guess > (x-guess*start2[0])*DIG_BASE+y))
     guess--;
    DBUG_ASSERT(start2[1]*guess <= (x-guess*start2[0])*DIG_BASE+y);
   }

   /* D4: multiply and subtract */
   buf2=stop2;
   buf1=start1+len2;
   DBUG_ASSERT(buf1 < stop1);
   for (carry=0; buf2 > start2; buf1--)
   {
    dec1 hi, lo;
    x=guess * (*--buf2);
    hi=(dec1)(x/DIG_BASE);
    lo=(dec1)(x-((dec2)hi)*DIG_BASE);
    SUB2(*buf1, *buf1, lo, carry);
    carry+=hi;
   }
   carry= dcarry < carry;

   /* D5: check the remainder */
   if (unlikely(carry))
   {
    /* D6: correct the guess */
    guess--;
    buf2=stop2;
    buf1=start1+len2;
    for (carry=0; buf2 > start2; buf1--)
    {
     ADD(*buf1, *buf1, *--buf2, carry);
    }
   }
  }
  if (likely(div_mod))
   *buf0=(dec1)guess;
  dcarry= *start1;
  start1++;
 }
 if (mod)
 {
  /*
    now the result is in tmp1, it has
    intg=prec1-frac1
    frac=max(frac1, frac2)=to->frac
  */
  if (dcarry)
   *--start1=dcarry;
  buf0=to->buf;
  intg0=(int) (ROUND_UP(prec1-frac1)-(start1-tmp1));
  frac0=ROUND_UP(to->frac);
  error=E_DEC_OK;
  if (unlikely(frac0==0 && intg0==0))
  {
   decimal_make_zero(to);
   goto done;
  }
  if (intg0<=0)
  {
   if (unlikely(-intg0 >= to->len))
   {
    decimal_make_zero(to);
    error=E_DEC_TRUNCATED;
    goto done;
   }
   stop1=start1+frac0;
   frac0+=intg0;
   to->intg=0;
   while (intg0++ < 0)
    *buf0++=0;
  }
  else
  {
   if (unlikely(intg0 > to->len))
   {
    frac0=0;
    intg0=to->len;
    error=E_DEC_OVERFLOW;
    goto done;
   }
   DBUG_ASSERT(intg0 <= ROUND_UP(from2->intg));
   stop1=start1+frac0+intg0;
   to->intg=min(intg0*DIG_PER_DEC1, from2->intg);
  }
  if (unlikely(intg0+frac0 > to->len))
  {
   stop1-=frac0+intg0-to->len;
   frac0=to->len-intg0;
   to->frac=frac0*DIG_PER_DEC1;
   error=E_DEC_TRUNCATED;
  }
  DBUG_ASSERT(buf0 + (stop1 - start1) <= to->buf + to->len);
  while (start1 < stop1)
   *buf0++=*start1++;
 }
done:
 my_afree(tmp1);
 return error;
}

/*
  division of two decimals

  SYNOPSIS
  decimal_div()
  from1   - dividend
  from2   - divisor
  to      - quotient

  RETURN VALUE
  E_DEC_OK/E_DEC_TRUNCATED/E_DEC_OVERFLOW/E_DEC_DIV_ZERO;

  NOTES
  see do_div_mod()
*/

int
decimal_div(decimal_t *from1, decimal_t *from2, decimal_t *to, int scale_incr)
{
 return do_div_mod(from1, from2, to, 0, scale_incr);
}

/*
  Convert string to decimal

  SYNOPSIS
  internal_str2decl()
  from    - value to convert. Doesn't have to be \0 terminated!
  to      - decimal where where the result will be stored
  to->buf and to->len must be set.
  end     - Pointer to pointer to end of string. Will on return be
  set to the char after the last used character
  fixed   - use to->intg, to->frac as limits for input number

  NOTE
  to->intg and to->frac can be modified even when fixed=1
  (but only decreased, in this case)

  RETURN VALUE
  E_DEC_OK/E_DEC_TRUNCATED/E_DEC_OVERFLOW/E_DEC_BAD_NUM/E_DEC_OOM
  In case of E_DEC_FATAL_ERROR *to is set to decimal zero
  (to make error handling easier)
*/

int internal_str2dec(const char *from, decimal_t *to, char **end, my_bool fixed)
{
 const char *s= from, *s1, *endp, *end_of_string= *end;
 int i, intg, frac, error, intg1, frac1;
 dec1 x,*buf;
 sanity(to);

 error= E_DEC_BAD_NUM;                         /* In case of bad number */
 // while (s < end_of_string && my_isspace(&my_charset_latin1, *s))
 //  s++;
 if (s == end_of_string)
  goto fatal_error;

 if ((to->sign= (*s == '-')))
  s++;
 else if (*s == '+')
  s++;

 s1=s;
 while (s < end_of_string && (*s>='0'&&*s<='9'))
  s++;
 intg= (int) (s-s1);
 if (s < end_of_string && *s=='.')
 {
  endp= s+1;
  while (endp < end_of_string && (*endp>='0'&&*endp<='9'))
   endp++;
  frac= (int) (endp - s - 1);
 }
 else
 {
  frac= 0;
  endp= s;
 }

 *end= (char*) endp;

 if (frac+intg == 0)
  goto fatal_error;

 error= 0;
 if (fixed)
 {
  if (frac > to->frac)
  {
   error=E_DEC_TRUNCATED;
   frac=to->frac;
  }
  if (intg > to->intg)
  {
   error=E_DEC_OVERFLOW;
   intg=to->intg;
  }
  intg1=ROUND_UP(intg);
  frac1=ROUND_UP(frac);
  if (intg1+frac1 > to->len)
  {
   error= E_DEC_OOM;
   goto fatal_error;
  }
 }
 else
 {
  intg1=ROUND_UP(intg);
  frac1=ROUND_UP(frac);
  FIX_INTG_FRAC_ERROR(to->len, intg1, frac1, error);
  if (unlikely(error))
  {
   frac=frac1*DIG_PER_DEC1;
   if (error == E_DEC_OVERFLOW)
    intg=intg1*DIG_PER_DEC1;
  }
 }
 /* Error is guranteed to be set here */
 to->intg=intg;
 to->frac=frac;

 buf=to->buf+intg1;
 s1=s;

 for (x=0, i=0; intg; intg--)
 {
  x+= (*--s - '0')*powers10[i];

  if (unlikely(++i == DIG_PER_DEC1))
  {
   *--buf=x;
   x=0;
   i=0;
  }
 }
 if (i)
  *--buf=x;

 buf=to->buf+intg1;
 for (x=0, i=0; frac; frac--)
 {
  x= (*++s1 - '0') + x*10;

  if (unlikely(++i == DIG_PER_DEC1))
  {
   *buf++=x;
   x=0;
   i=0;
  }
 }
 if (i)
  *buf=x*powers10[DIG_PER_DEC1-i];

 /* Handle exponent */
 // if (endp+1 < end_of_string && (*endp == 'e' || *endp == 'E'))
 // {
 //  int str_error;
 //  longlong exponent= my_strtoll10(endp+1, (char**) &end_of_string,
 //          &str_error);

 //  if (end_of_string != endp +1)               /* If at least one digit */
 //  {
 //   *end= (char*) end_of_string;
 //   if (str_error > 0)
 //   {
 //    error= E_DEC_BAD_NUM;
 //    goto fatal_error;
 //   }
 //   if (exponent > INT_MAX/2 || (str_error == 0 && exponent < 0))
 //   {
 //    error= E_DEC_OVERFLOW;
 //    goto fatal_error;
 //   }
 //   if (exponent < INT_MIN/2 && error != E_DEC_OVERFLOW)
 //   {
 //    error= E_DEC_TRUNCATED;
 //    goto fatal_error;
 //   }
 //   if (error != E_DEC_OVERFLOW)
 //    error= decimal_shift(to, (int) exponent);
 //  }
 // }
 return error;

fatal_error:
 decimal_make_zero(to);
 return error;
}

int full= 0;
decimal_t a, b, c;
char buf1[100], buf2[100], buf3[100];
#if VaxAsm

char *strend(s)
const char *s;
{
  asm("locc $0,$65535,*4(ap)");
  asm("movl r1,r0");
}

#else /* ~VaxAsm */

char *strend(register const char *s)
{
  while (*s++);
  return (char*) (s-1);
}

#endif /* VaxAsm */

#define string2decimal(A,B,C) internal_str2dec((A), (B), (C), 0)
#define string2decimal_fixed(A,B,C) internal_str2dec((A), (B), (C), 1)
void dump_decimal(decimal_t *d)
{
  int i;
  printf("/* intg=%d, frac=%d, sign=%d, buf[]={", d->intg, d->frac, d->sign);
  for (i=0; i < ROUND_UP(d->frac)+ROUND_UP(d->intg)-1; i++)
    printf("%09d, ", d->buf[i]);
  printf("%09d} */ ", d->buf[i]);
}

void check_result_code(int actual, int want)
{
  if (actual != want)
  {
    printf("\n^^^^^^^^^^^^^ must return %d\n", want);
    exit(1);
  }
}

int double2decimal(double from, decimal_t *to)
{
  /* TODO: fix it, when we'll have dtoa */
  char buff[400], *end;
  int length, res;

  length= sprintf(buff, "%.16G", from);
  end= buff+length;
  res= string2decimal(buff, to, &end);
  return(res);
}

void test_s2d(const char *s, const char *orig, int ex)
{
  char s1[100], *end;
  int res;
  end= strend(s);
  res= string2decimal(s, &a, &end);
}

void test_f2d(double from, int ex)
{
  int res;

  res=double2decimal(from, &a);
   printf("\n");
}

void mysqltest ()
{
 double tt;
 
 /* a.buf=(decimal_digit_t *)buf1; */
 /* a.len=sizeof(buf1)/sizeof(dec1); */
 /* b.buf=(decimal_digit_t *)buf2; */
 /* b.len=sizeof(buf2)/sizeof(dec1); */
 /* c.buf=(decimal_digit_t *)buf3; */
 /* c.len=sizeof(buf3)/sizeof(dec1); */
 printf("====================\n");
 
  
 test_s2d("12345", "12345", 0);
 test_s2d("12345.", "12345", 0);
 test_s2d("123.45", "123.45", 0);
 test_s2d("-123.45", "-123.45", 0);
 test_s2d("1.00012345000098765", "1.00012345000098765", 0);
 test_s2d("-1.000000012345000098765", "-1.000000012345000098765", 0);
 test_s2d("1234500009876.5", "1234500009876.5", 0);
 a.len=1;
 test_s2d("123450000098765", "98765", 2);
 test_s2d("123450.000098765", "123450", 1);

 test_f2d(12345, 0);
 test_f2d(1.0/3, 0);
 test_f2d(-123.45, 0);
 test_f2d(0.00012345000098765, 0);
 test_f2d(1234500009876.5, 0); 
 
  
 // mdtime(0);
 // f += e;
 // tt= mdtime(1);
 // printf(" += cycle:%15f \n", tt);

 // LongDecimal a(123456789,5);
 // LongDecimal b(987654321,4);
 // mdtime(0);
 // b *= a;
 // tt= mdtime(1);
 // printf(" *= cycle:%15f \n", tt);

 // LongDecimal c(123456789,5);
 // LongDecimal d(987654321,4);
 // mdtime(0);
 // c /= d;
 // tt= mdtime(1);
 // printf(" /= cycle:%15f \n", tt);
 
}

int main()
{
 printf("%s", "=================\n");
 
 mysqltest();
 
}