vasnprintf.c

gnu tar 源码包。 tar 软件是 Unix 系统下的一个打包软件

	    s -= 2;
	  }
	if (msd >= 0x2)
	  {
	    msd = msd >> 1;
	    s -= 1;
	  }
      }
      /* 0 <= s < GMP_LIMB_BITS.
	 Copy b, shifting it left by s bits.  */
      if (s > 0)
	{
	  tmp_roomptr = (mp_limb_t *) malloc (b_len * sizeof (mp_limb_t));
	  if (tmp_roomptr == NULL)
	    {
	      free (roomptr);
	      return NULL;
	    }
	  {
	    const mp_limb_t *sourceptr = b_ptr;
	    mp_limb_t *destptr = tmp_roomptr;
	    mp_twolimb_t accu = 0;
	    size_t count;
	    for (count = b_len; count > 0; count--)
	      {
		accu += (mp_twolimb_t) *sourceptr++ << s;
		*destptr++ = (mp_limb_t) accu;
		accu = accu >> GMP_LIMB_BITS;
	      }
	    /* accu must be zero, since that was how s was determined.  */
	    if (accu != 0)
	      abort ();
	  }
	  b_ptr = tmp_roomptr;
	}
      /* Copy a, shifting it left by s bits, yields r.
	 Memory layout:
	 At the beginning: r = roomptr[0..a_len],
	 at the end: r = roomptr[0..b_len-1], q = roomptr[b_len..a_len]  */
      r_ptr = roomptr;
      if (s == 0)
	{
	  memcpy (r_ptr, a_ptr, a_len * sizeof (mp_limb_t));
	  r_ptr[a_len] = 0;
	}
      else
	{
	  const mp_limb_t *sourceptr = a_ptr;
	  mp_limb_t *destptr = r_ptr;
	  mp_twolimb_t accu = 0;
	  size_t count;
	  for (count = a_len; count > 0; count--)
	    {
	      accu += (mp_twolimb_t) *sourceptr++ << s;
	      *destptr++ = (mp_limb_t) accu;
	      accu = accu >> GMP_LIMB_BITS;
	    }
	  *destptr++ = (mp_limb_t) accu;
	}
      q_ptr = roomptr + b_len;
      q_len = a_len - b_len + 1; /* q will have m-n+1 limbs */
      {
	size_t j = a_len - b_len; /* m-n */
	mp_limb_t b_msd = b_ptr[b_len - 1]; /* b[n-1] */
	mp_limb_t b_2msd = b_ptr[b_len - 2]; /* b[n-2] */
	mp_twolimb_t b_msdd = /* b[n-1]*beta+b[n-2] */
	  ((mp_twolimb_t) b_msd << GMP_LIMB_BITS) | b_2msd;
	/* Division loop, traversed m-n+1 times.
	   j counts down, b is unchanged, beta/2 <= b[n-1] < beta.  */
	for (;;)
	  {
	    mp_limb_t q_star;
	    mp_limb_t c1;
	    if (r_ptr[j + b_len] < b_msd) /* r[j+n] < b[n-1] ? */
	      {
		/* Divide r[j+n]*beta+r[j+n-1] by b[n-1], no overflow.  */
		mp_twolimb_t num =
		  ((mp_twolimb_t) r_ptr[j + b_len] << GMP_LIMB_BITS)
		  | r_ptr[j + b_len - 1];
		q_star = num / b_msd;
		c1 = num % b_msd;
	      }
	    else
	      {
		/* Overflow, hence r[j+n]*beta+r[j+n-1] >= beta*b[n-1].  */
		q_star = (mp_limb_t)~(mp_limb_t)0; /* q* = beta-1 */
		/* Test whether r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] >= beta
		   <==> r[j+n]*beta+r[j+n-1] + b[n-1] >= beta*b[n-1]+beta
		   <==> b[n-1] < floor((r[j+n]*beta+r[j+n-1]+b[n-1])/beta)
		        {<= beta !}.
		   If yes, jump directly to the subtraction loop.
		   (Otherwise, r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] < beta
		    <==> floor((r[j+n]*beta+r[j+n-1]+b[n-1])/beta) = b[n-1] ) */
		if (r_ptr[j + b_len] > b_msd
		    || (c1 = r_ptr[j + b_len - 1] + b_msd) < b_msd)
		  /* r[j+n] >= b[n-1]+1 or
		     r[j+n] = b[n-1] and the addition r[j+n-1]+b[n-1] gives a
		     carry.  */
		  goto subtract;
	      }
	    /* q_star = q*,
	       c1 = (r[j+n]*beta+r[j+n-1]) - q* * b[n-1] (>=0, <beta).  */
	    {
	      mp_twolimb_t c2 = /* c1*beta+r[j+n-2] */
		((mp_twolimb_t) c1 << GMP_LIMB_BITS) | r_ptr[j + b_len - 2];
	      mp_twolimb_t c3 = /* b[n-2] * q* */
		(mp_twolimb_t) b_2msd * (mp_twolimb_t) q_star;
	      /* While c2 < c3, increase c2 and decrease c3.
		 Consider c3-c2.  While it is > 0, decrease it by
		 b[n-1]*beta+b[n-2].  Because of b[n-1]*beta+b[n-2] >= beta^2/2
		 this can happen only twice.  */
	      if (c3 > c2)
		{
		  q_star = q_star - 1; /* q* := q* - 1 */
		  if (c3 - c2 > b_msdd)
		    q_star = q_star - 1; /* q* := q* - 1 */
		}
	    }
	    if (q_star > 0)
	      subtract:
	      {
		/* Subtract r := r - b * q* * beta^j.  */
		mp_limb_t cr;
		{
		  const mp_limb_t *sourceptr = b_ptr;
		  mp_limb_t *destptr = r_ptr + j;
		  mp_twolimb_t carry = 0;
		  size_t count;
		  for (count = b_len; count > 0; count--)
		    {
		      /* Here 0 <= carry <= q*.  */
		      carry =
			carry
			+ (mp_twolimb_t) q_star * (mp_twolimb_t) *sourceptr++
			+ (mp_limb_t) ~(*destptr);
		      /* Here 0 <= carry <= beta*q* + beta-1.  */
		      *destptr++ = ~(mp_limb_t) carry;
		      carry = carry >> GMP_LIMB_BITS; /* <= q* */
		    }
		  cr = (mp_limb_t) carry;
		}
		/* Subtract cr from r_ptr[j + b_len], then forget about
		   r_ptr[j + b_len].  */
		if (cr > r_ptr[j + b_len])
		  {
		    /* Subtraction gave a carry.  */
		    q_star = q_star - 1; /* q* := q* - 1 */
		    /* Add b back.  */
		    {
		      const mp_limb_t *sourceptr = b_ptr;
		      mp_limb_t *destptr = r_ptr + j;
		      mp_limb_t carry = 0;
		      size_t count;
		      for (count = b_len; count > 0; count--)
			{
			  mp_limb_t source1 = *sourceptr++;
			  mp_limb_t source2 = *destptr;
			  *destptr++ = source1 + source2 + carry;
			  carry =
			    (carry
			     ? source1 >= (mp_limb_t) ~source2
			     : source1 > (mp_limb_t) ~source2);
			}
		    }
		    /* Forget about the carry and about r[j+n].  */
		  }
	      }
	    /* q* is determined.  Store it as q[j].  */
	    q_ptr[j] = q_star;
	    if (j == 0)
	      break;
	    j--;
	  }
      }
      r_len = b_len;
      /* Normalise q.  */
      if (q_ptr[q_len - 1] == 0)
	q_len--;
# if 0 /* Not needed here, since we need r only to compare it with b/2, and
	  b is shifted left by s bits.  */
      /* Shift r right by s bits.  */
      if (s > 0)
	{
	  mp_limb_t ptr = r_ptr + r_len;
	  mp_twolimb_t accu = 0;
	  size_t count;
	  for (count = r_len; count > 0; count--)
	    {
	      accu = (mp_twolimb_t) (mp_limb_t) accu << GMP_LIMB_BITS;
	      accu += (mp_twolimb_t) *--ptr << (GMP_LIMB_BITS - s);
	      *ptr = (mp_limb_t) (accu >> GMP_LIMB_BITS);
	    }
	}
# endif
      /* Normalise r.  */
      while (r_len > 0 && r_ptr[r_len - 1] == 0)
	r_len--;
    }
  /* Compare r << 1 with b.  */
  if (r_len > b_len)
    goto increment_q;
  {
    size_t i;
    for (i = b_len;;)
      {
	mp_limb_t r_i =
	  (i <= r_len && i > 0 ? r_ptr[i - 1] >> (GMP_LIMB_BITS - 1) : 0)
	  | (i < r_len ? r_ptr[i] << 1 : 0);
	mp_limb_t b_i = (i < b_len ? b_ptr[i] : 0);
	if (r_i > b_i)
	  goto increment_q;
	if (r_i < b_i)
	  goto keep_q;
	if (i == 0)
	  break;
	i--;
      }
  }
  if (q_len > 0 && ((q_ptr[0] & 1) != 0))
    /* q is odd.  */
    increment_q:
    {
      size_t i;
      for (i = 0; i < q_len; i++)
	if (++(q_ptr[i]) != 0)
	  goto keep_q;
      q_ptr[q_len++] = 1;
    }
  keep_q:
  if (tmp_roomptr != NULL)
    free (tmp_roomptr);
  q->limbs = q_ptr;
  q->nlimbs = q_len;
  return roomptr;
}

/* Convert a bignum a >= 0, multiplied with 10^extra_zeroes, to decimal
   representation.
   Destroys the contents of a.
   Return the allocated memory - containing the decimal digits in low-to-high
   order, terminated with a NUL character - in case of success, NULL in case
   of memory allocation failure.  */
static char *
convert_to_decimal (mpn_t a, size_t extra_zeroes)
{
  mp_limb_t *a_ptr = a.limbs;
  size_t a_len = a.nlimbs;
  /* 0.03345 is slightly larger than log(2)/(9*log(10)).  */
  size_t c_len = 9 * ((size_t)(a_len * (GMP_LIMB_BITS * 0.03345f)) + 1);
  char *c_ptr = (char *) malloc (xsum (c_len, extra_zeroes));
  if (c_ptr != NULL)
    {
      char *d_ptr = c_ptr;
      for (; extra_zeroes > 0; extra_zeroes--)
	*d_ptr++ = '0';
      while (a_len > 0)
	{
	  /* Divide a by 10^9, in-place.  */
	  mp_limb_t remainder = 0;
	  mp_limb_t *ptr = a_ptr + a_len;
	  size_t count;
	  for (count = a_len; count > 0; count--)
	    {
	      mp_twolimb_t num =
		((mp_twolimb_t) remainder << GMP_LIMB_BITS) | *--ptr;
	      *ptr = num / 1000000000;
	      remainder = num % 1000000000;
	    }
	  /* Store the remainder as 9 decimal digits.  */
	  for (count = 9; count > 0; count--)
	    {
	      *d_ptr++ = '0' + (remainder % 10);
	      remainder = remainder / 10;
	    }
	  /* Normalize a.  */
	  if (a_ptr[a_len - 1] == 0)
	    a_len--;
	}
      /* Remove leading zeroes.  */
      while (d_ptr > c_ptr && d_ptr[-1] == '0')
	d_ptr--;
      /* But keep at least one zero.  */
      if (d_ptr == c_ptr)
	*d_ptr++ = '0';
      /* Terminate the string.  */
      *d_ptr = '\0';
    }
  return c_ptr;
}

# if NEED_PRINTF_LONG_DOUBLE

/* Assuming x is finite and >= 0:
   write x as x = 2^e * m, where m is a bignum.
   Return the allocated memory in case of success, NULL in case of memory
   allocation failure.  */
static void *
decode_long_double (long double x, int *ep, mpn_t *mp)
{
  mpn_t m;
  int exp;
  long double y;
  size_t i;

  /* Allocate memory for result.  */
  m.nlimbs = (LDBL_MANT_BIT + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
  m.limbs = (mp_limb_t *) malloc (m.nlimbs * sizeof (mp_limb_t));
  if (m.limbs == NULL)
    return NULL;
  /* Split into exponential part and mantissa.  */
  y = frexpl (x, &exp);
  if (!(y >= 0.0L && y < 1.0L))
    abort ();
  /* x = 2^exp * y = 2^(exp - LDBL_MANT_BIT) * (y * LDBL_MANT_BIT), and the
     latter is an integer.  */
  /* Convert the mantissa (y * LDBL_MANT_BIT) to a sequence of limbs.
     I'm not sure whether it's safe to cast a 'long double' value between
     2^31 and 2^32 to 'unsigned int', therefore play safe and cast only
     'long double' values between 0 and 2^16 (to 'unsigned int' or 'int',
     doesn't matter).  */
#  if (LDBL_MANT_BIT % GMP_LIMB_BITS) != 0
#   if (LDBL_MANT_BIT % GMP_LIMB_BITS) > GMP_LIMB_BITS / 2
    {
      mp_limb_t hi, lo;
      y *= (mp_limb_t) 1 << (LDBL_MANT_BIT % (GMP_LIMB_BITS / 2));
      hi = (int) y;
      y -= hi;
      if (!(y >= 0.0L && y < 1.0L))
	abort ();
      y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
      lo = (int) y;
      y -= lo;
      if (!(y >= 0.0L && y < 1.0L))
	abort ();
      m.limbs[LDBL_MANT_BIT / GMP_LIMB_BITS] = (hi << (GMP_LIMB_BITS / 2)) | lo;
    }
#   else
    {
      mp_limb_t d;
      y *= (mp_limb_t) 1 << (LDBL_MANT_BIT % GMP_LIMB_BITS);
      d = (int) y;
      y -= d;
      if (!(y >= 0.0L && y < 1.0L))
	abort ();
      m.limbs[LDBL_MANT_BIT / GMP_LIMB_BITS] = d;
    }
#   endif
#  endif
  for (i = LDBL_MANT_BIT / GMP_LIMB_BITS; i > 0; )
    {
      mp_limb_t hi, lo;
      y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
      hi = (int) y;
      y -= hi;
      if (!(y >= 0.0L && y < 1.0L))
	abort ();
      y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
      lo = (int) y;
      y -= lo;
      if (!(y >= 0.0L && y < 1.0L))
	abort ();
      m.limbs[--i] = (hi << (GMP_LIMB_BITS / 2)) | lo;
    }
#if 0 /* On FreeBSD 6.1/x86, 'long double' numbers sometimes have excess
         precision.  */
  if (!(y == 0.0L))
    abort ();
#endif
  /* Normalise.  */
  while (m.nlimbs > 0 && m.limbs[m.nlimbs - 1] == 0)
    m.nlimbs--;
  *mp = m;
  *ep = exp - LDBL_MANT_BIT;
  return m.limbs;
}

# endif

# if NEED_PRINTF_DOUBLE

/* Assuming x is finite and >= 0:
   write x as x = 2^e * m, where m is a bignum.
   Return the allocated memory in case of success, NULL in case of memory
   allocation failure.  */
static void *
decode_double (double x, int *ep, mpn_t *mp)
{
  mpn_t m;
  int exp;
  double y;
  size_t i;

  /* Allocate memory for result.  */
  m.nlimbs = (DBL_MANT_BIT + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
  m.limbs = (mp_limb_t *) malloc (m.nlimbs * sizeof (mp_limb_t));
  if (m.limbs == NULL)
    return NULL;
  /* Split into exponential part and mantissa.  */
  y = frexp (x, &exp);
  if (!(y >= 0.0 && y < 1.0))
    abort ();
  /* x = 2^exp * y = 2^(exp - DBL_MANT_BIT) * (y * DBL_MANT_BIT), and the
     latter is an integer.  */
  /* Convert the mantissa (y * DBL_MANT_BIT) to a sequence of limbs.
     I'm not sure whether it's safe to cast a 'double' value between
     2^31 and 2^32 to 'unsigned int', therefore play safe and cast only
     'double' values between 0 and 2^16 (to 'unsigned int' or 'int',
     doesn't matter).  */
#  if (DBL_MANT_BIT % GMP_LIMB_BITS) != 0
#   if (DBL_MANT_BIT % GMP_LIMB_BITS) > GMP_LIMB_BITS / 2
    {
      mp_limb_t hi, lo;
      y *= (mp_limb_t) 1 << (DBL_MANT_BIT % (GMP_LIMB_BITS / 2));
      hi = (int) y;
      y -= hi;
      if (!(y >= 0.0 && y < 1.0))
	abort ();
      y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
      lo = (int) y;
      y -= lo;
      if (!(y >= 0.0 && y < 1.0))
	abort ();
      m.limbs[DBL_MANT_BIT / GMP_LIMB_BITS] = (hi << (GMP_LIMB_BITS / 2)) | lo;
    }
#   else
    {
      mp_limb_t d;
      y *= (mp_limb_t) 1 << (DBL_MANT_BIT % GMP_LIMB_BITS);
      d = (int) y;
      y -= d;
      if (!(y >= 0.0 && y < 1.0))
	abort ();
      m.limbs[DBL_MANT_BIT / GMP_LIMB_BITS] = d;
    }
#   endif
#  endif
  for (i = DBL_MANT_BIT / GMP_LIMB_BITS; i > 0; )
    {
      mp_limb_t hi, lo;
      y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
      hi = (int) y;
      y -= hi;
      if (!(y >= 0.0 && y < 1.0))
	abort ();
      y *= (mp_limb_t) 1 << (GMP_LIMB_BITS / 2);
      lo = (int) y;
      y -= lo;
      if (!(y >= 0.0 && y < 1.0))
	abort ();
      m.limbs[--i] = (hi << (GMP_LIMB_BITS / 2)) | lo;
    }
  if (!(y == 0.0))
    abort ();
  /* Normalise.  */
  while (m.nlimbs > 0 && m.limbs[m.nlimbs - 1] == 0)
    m.nlimbs--;
  *mp = m;
  *ep = exp - DBL_MANT_BIT;
  return m.limbs;
}

# endif

/* Assuming x = 2^e * m is finite and >= 0, and n is an integer:
   Returns the decimal representation of round (x * 10^n).
   Return the allocated memory - containing the decimal digits in low-to-high
   order, terminated with a NUL character - in case of success, NULL in case
   of memory allocation failure.  */
static char *
scale10_round_decimal_decoded (int e, mpn_t m, void *memory, int n)
{
  int s;
  size_t extra_zeroes;
  unsigned int abs_n;
  unsigned int abs_s;
  mp_limb_t *pow5_ptr;
  size_t pow5_len;
  unsigned int s_limbs;
  unsigned int s_bits;
  mpn_t pow5;
  mpn_t z;
  void *z_memory;
  char *digits;

  if (memory == NULL)
    return NULL;
  /* x = 2^e * m, hence
     y = round (2^e * 10^n * m) = round (2^(e+n) * 5^n * m)
       = round (2^s * 5^n * m).  */
  s = e + n;
  extra_zeroes = 0;
  /* Factor out a common power of 10 if possible.  */
  if (s > 0 && n > 0)
    {
      extra_zeroes = (s < n ? s : n);
      s -= extra_zeroes;
      n -= extra_zeroes;
    }
  /* Here y = round (2^s * 5^n * m) * 10^extra_zeroes.
     Before converting to decimal, we need to compute
     z = round (2^s * 5^n * m).  */
  /* Compute 5^|n|, possibly shifted by |s| bits if n and s have the same
     sign.  2.322 is slightly larger than log(5)/log(2).  */
  abs_n = (n >= 0 ? n : -n);
  abs_s = (s >= 0 ? s : -s);
  pow5_ptr = (mp_limb_t *) malloc (((int)(abs_n * (2.322f / GMP_LIMB_BITS)) + 1
				    + abs_s / GMP_LIMB_BITS + 1)