obtained gcc-4.6.4.tar.bz2 from upstream website;upstream

verified gcc-4.6.4.tar.bz2.sig;
imported gcc-4.6.4 source tree from verified upstream tarball.

downloading a git-generated archive based on the 'upstream' tag
should provide you with a source tree that is binary identical
to the one extracted from the above tarball.

if you have obtained the source via the command 'git clone',
however, do note that line-endings of files in your working
directory might differ from line-endings of the respective
files in the upstream repository.

22 files changed, 4673 insertions, 0 deletions

diff --git a/libquadmath/printf/_i18n_number.h b/libquadmath/printf/_i18n_number.h
new file mode 100644
index 000000000..80065d88b
--- /dev/null
+++ b/libquadmath/printf/_i18n_number.h
@@ -0,0 +1,139 @@
+/* Copyright (C) 2000, 2004, 2008 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gnu.org>, 2000.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+/* Look up the value of the next multibyte character and return its numerical
+   value if it is one of the digits known in the locale.  If *DECIDED is
+   -1 this means it is not yet decided which form it is and we have to
+   search through all available digits.  Otherwise we know which script
+   the digits are from.  */
+static inline char *
+outdigit_value (char *s, int n)
+{
+  const char *outdigit;
+  size_t dlen;
+
+  assert (0 <= n && n <= 9);
+  outdigit = nl_langinfo (_NL_CTYPE_OUTDIGIT0_MB + n);
+  dlen = strlen (outdigit);
+
+  s -= dlen;
+  while (dlen-- > 0)
+    s[dlen] = outdigit[dlen];
+
+  return s;
+}
+
+/* Look up the value of the next multibyte character and return its numerical
+   value if it is one of the digits known in the locale.  If *DECIDED is
+   -1 this means it is not yet decided which form it is and we have to
+   search through all available digits.  Otherwise we know which script
+   the digits are from.  */
+static inline wchar_t
+outdigitwc_value (int n)
+{
+  assert (0 <= n && n <= 9);
+
+  return nl_langinfo_wc (_NL_CTYPE_OUTDIGIT0_WC + n);
+}
+
+static char *
+_i18n_number_rewrite (char *w, char *rear_ptr, char *end)
+{
+  char decimal[MB_LEN_MAX + 1];
+  char thousands[MB_LEN_MAX + 1];
+
+  /* "to_outpunct" is a map from ASCII decimal point and thousands-sep
+     to their equivalent in locale. This is defined for locales which
+     use extra decimal point and thousands-sep.  */
+  wctrans_t map = wctrans ("to_outpunct");
+  wint_t wdecimal = towctrans (L_('.'), map);
+  wint_t wthousands = towctrans (L_(','), map);
+
+  if (__builtin_expect (map != NULL, 0))
+    {
+      mbstate_t state;
+      memset (&state, '\0', sizeof (state));
+
+      size_t n = wcrtomb (decimal, wdecimal, &state);
+      if (n == (size_t) -1)
+	memcpy (decimal, ".", 2);
+      else
+	decimal[n] = '\0';
+
+      memset (&state, '\0', sizeof (state));
+
+      n = wcrtomb (thousands, wthousands, &state);
+      if (n == (size_t) -1)
+	memcpy (thousands, ",", 2);
+      else
+	thousands[n] = '\0';
+    }
+
+  /* Copy existing string so that nothing gets overwritten.  */
+  char *src;
+  int use_alloca = (rear_ptr - w) < 4096;
+  if (__builtin_expect (use_alloca, 1))
+    src = (char *) alloca (rear_ptr - w);
+  else
+    {
+      src = (char *) malloc (rear_ptr - w);
+      if (src == NULL)
+	/* If we cannot allocate the memory don't rewrite the string.
+	   It is better than nothing.  */
+	return w;
+    }
+
+  memcpy (src, w, rear_ptr - w);
+  char *s = src + (rear_ptr - w);
+
+  w = end;
+
+  /* Process all characters in the string.  */
+  while (--s >= src)
+    {
+      if (*s >= '0' && *s <= '9')
+	{
+	  if (sizeof (char) == 1)
+	    w = (char *) outdigit_value ((char *) w, *s - '0');
+	  else
+	    *--w = (char) outdigitwc_value (*s - '0');
+	}
+      else if (__builtin_expect (map == NULL, 1) || (*s != '.' && *s != ','))
+	*--w = *s;
+      else
+	{
+	  if (sizeof (char) == 1)
+	    {
+	      const char *outpunct = *s == '.' ? decimal : thousands;
+	      size_t dlen = strlen (outpunct);
+
+	      w -= dlen;
+	      while (dlen-- > 0)
+		w[dlen] = outpunct[dlen];
+	    }
+	  else
+	    *--w = *s == '.' ? (char) wdecimal : (char) wthousands;
+	}
+    }
+
+  if (! use_alloca)
+    free (src);
+
+  return w;
+}
diff --git a/libquadmath/printf/_itoa.h b/libquadmath/printf/_itoa.h
new file mode 100644
index 000000000..a0cd2b05c
--- /dev/null
+++ b/libquadmath/printf/_itoa.h
@@ -0,0 +1,77 @@
+/* Internal function for converting integers to ASCII.
+   Copyright (C) 1994-1999,2002,2003,2007 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#ifndef _ITOA_H
+#define _ITOA_H
+
+/* Convert VALUE into ASCII in base BASE (2..16).
+   Write backwards starting the character just before BUFLIM.
+   Return the address of the first (left-to-right) character in the number.
+   Use upper case letters iff UPPER_CASE is nonzero.  */
+
+static const char _itoa_lower_digits[16] = "0123456789abcdef";
+static const char _itoa_upper_digits[16] = "0123456789ABCDEF";
+
+static inline char * __attribute__ ((unused, always_inline))
+_itoa_word (unsigned long value, char *buflim,
+	    unsigned int base, int upper_case)
+{
+  const char *digits = (upper_case ? _itoa_upper_digits : _itoa_lower_digits);
+
+  switch (base)
+    {
+# define SPECIAL(Base)							      \
+    case Base:								      \
+      do								      \
+	*--buflim = digits[value % Base];				      \
+      while ((value /= Base) != 0);					      \
+      break
+
+      SPECIAL (10);
+      SPECIAL (16);
+      SPECIAL (8);
+    default:
+      do
+	*--buflim = digits[value % base];
+      while ((value /= base) != 0);
+    }
+  return buflim;
+}
+
+static inline char * __attribute__ ((unused, always_inline))
+_itoa (uint64_t value, char *buflim,
+       unsigned int base, int upper_case)
+{
+  const char *digits = (upper_case ? _itoa_upper_digits : _itoa_lower_digits);
+
+  switch (base)
+    {
+      SPECIAL (10);
+      SPECIAL (16);
+      SPECIAL (8);
+    default:
+      do
+	*--buflim = digits[value % base];
+      while ((value /= base) != 0);
+    }
+  return buflim;
+}
+# undef SPECIAL
+
+#endif	/* itoa.h */
diff --git a/libquadmath/printf/_itowa.h b/libquadmath/printf/_itowa.h
new file mode 100644
index 000000000..4717b5c65
--- /dev/null
+++ b/libquadmath/printf/_itowa.h
@@ -0,0 +1,83 @@
+/* Internal function for converting integers to ASCII.
+   Copyright (C) 1994,95,96,97,98,99,2002,2003 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#ifndef _ITOWA_H
+#define _ITOWA_H	1
+
+/* Convert VALUE into ASCII in base BASE (2..16).
+   Write backwards starting the character just before BUFLIM.
+   Return the address of the first (left-to-right) character in the number.
+   Use upper case letters iff UPPER_CASE is nonzero.  */
+
+static const wchar_t _itowa_lower_digits[16] = L_("0123456789abcdef");
+static const wchar_t _itowa_upper_digits[16] = L_("0123456789ABCDEF");
+
+static inline wchar_t *
+__attribute__ ((unused, always_inline))
+_itowa_word (unsigned long value, wchar_t *buflim,
+	     unsigned int base, int upper_case)
+{
+  const wchar_t *digits = (upper_case
+			   ? _itowa_upper_digits : _itowa_lower_digits);
+  wchar_t *bp = buflim;
+
+  switch (base)
+    {
+#define SPECIAL(Base)							      \
+    case Base:								      \
+      do								      \
+	*--bp = digits[value % Base];					      \
+      while ((value /= Base) != 0);					      \
+      break
+
+      SPECIAL (10);
+      SPECIAL (16);
+      SPECIAL (8);
+    default:
+      do
+	*--bp = digits[value % base];
+      while ((value /= base) != 0);
+    }
+  return bp;
+}
+
+static inline wchar_t *
+__attribute__ ((unused, always_inline))
+_itowa (uint64_t value, wchar_t *buflim,
+	unsigned int base, int upper_case)
+{
+  const wchar_t *digits = (upper_case
+			   ? _itowa_upper_digits : _itowa_lower_digits);
+  wchar_t *bp = buflim;
+
+  switch (base)
+    {
+      SPECIAL (10);
+      SPECIAL (16);
+      SPECIAL (8);
+    default:
+      do
+	*--bp = digits[value % base];
+      while ((value /= base) != 0);
+    }
+  return bp;
+}
+#undef SPECIAL
+
+#endif	/* itowa.h */
diff --git a/libquadmath/printf/add_n.c b/libquadmath/printf/add_n.c
new file mode 100644
index 000000000..749cf3ee8
--- /dev/null
+++ b/libquadmath/printf/add_n.c
@@ -0,0 +1,62 @@
+/* mpn_add_n -- Add two limb vectors of equal, non-zero length.
+
+Copyright (C) 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+mp_limb_t
+#if __STDC__
+mpn_add_n (mp_ptr res_ptr, mp_srcptr s1_ptr, mp_srcptr s2_ptr, mp_size_t size)
+#else
+mpn_add_n (res_ptr, s1_ptr, s2_ptr, size)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     register mp_srcptr s2_ptr;
+     mp_size_t size;
+#endif
+{
+  register mp_limb_t x, y, cy;
+  register mp_size_t j;
+
+  /* The loop counter and index J goes from -SIZE to -1.  This way
+     the loop becomes faster.  */
+  j = -size;
+
+  /* Offset the base pointers to compensate for the negative indices.  */
+  s1_ptr -= j;
+  s2_ptr -= j;
+  res_ptr -= j;
+
+  cy = 0;
+  do
+    {
+      y = s2_ptr[j];
+      x = s1_ptr[j];
+      y += cy;			/* add previous carry to one addend */
+      cy = (y < cy);		/* get out carry from that addition */
+      y = x + y;		/* add other addend */
+      cy = (y < x) + cy;	/* get out carry from that add, combine */
+      res_ptr[j] = y;
+    }
+  while (++j != 0);
+
+  return cy;
+}
diff --git a/libquadmath/printf/addmul_1.c b/libquadmath/printf/addmul_1.c
new file mode 100644
index 000000000..f527f9848
--- /dev/null
+++ b/libquadmath/printf/addmul_1.c
@@ -0,0 +1,64 @@
+/* mpn_addmul_1 -- multiply the S1_SIZE long limb vector pointed to by S1_PTR
+   by S2_LIMB, add the S1_SIZE least significant limbs of the product to the
+   limb vector pointed to by RES_PTR.  Return the most significant limb of
+   the product, adjusted for carry-out from the addition.
+
+Copyright (C) 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+mp_limb_t
+mpn_addmul_1 (res_ptr, s1_ptr, s1_size, s2_limb)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     mp_size_t s1_size;
+     register mp_limb_t s2_limb;
+{
+  register mp_limb_t cy_limb;
+  register mp_size_t j;
+  register mp_limb_t prod_high, prod_low;
+  register mp_limb_t x;
+
+  /* The loop counter and index J goes from -SIZE to -1.  This way
+     the loop becomes faster.  */
+  j = -s1_size;
+
+  /* Offset the base pointers to compensate for the negative indices.  */
+  res_ptr -= j;
+  s1_ptr -= j;
+
+  cy_limb = 0;
+  do
+    {
+      umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
+
+      prod_low += cy_limb;
+      cy_limb = (prod_low < cy_limb) + prod_high;
+
+      x = res_ptr[j];
+      prod_low = x + prod_low;
+      cy_limb += (prod_low < x);
+      res_ptr[j] = prod_low;
+    }
+  while (++j != 0);
+
+  return cy_limb;
+}
diff --git a/libquadmath/printf/cmp.c b/libquadmath/printf/cmp.c
new file mode 100644
index 000000000..a4be43e2a
--- /dev/null
+++ b/libquadmath/printf/cmp.c
@@ -0,0 +1,56 @@
+/* mpn_cmp -- Compare two low-level natural-number integers.
+
+Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+/* Compare OP1_PTR/OP1_SIZE with OP2_PTR/OP2_SIZE.
+   There are no restrictions on the relative sizes of
+   the two arguments.
+   Return 1 if OP1 > OP2, 0 if they are equal, and -1 if OP1 < OP2.  */
+
+int
+#if __STDC__
+mpn_cmp (mp_srcptr op1_ptr, mp_srcptr op2_ptr, mp_size_t size)
+#else
+mpn_cmp (op1_ptr, op2_ptr, size)
+     mp_srcptr op1_ptr;
+     mp_srcptr op2_ptr;
+     mp_size_t size;
+#endif
+{
+  mp_size_t i;
+  mp_limb_t op1_word, op2_word;
+
+  for (i = size - 1; i >= 0; i--)
+    {
+      op1_word = op1_ptr[i];
+      op2_word = op2_ptr[i];
+      if (op1_word != op2_word)
+	goto diff;
+    }
+  return 0;
+ diff:
+  /* This can *not* be simplified to
+	op2_word - op2_word
+     since that expression might give signed overflow.  */
+  return (op1_word > op2_word) ? 1 : -1;
+}
diff --git a/libquadmath/printf/divrem.c b/libquadmath/printf/divrem.c
new file mode 100644
index 000000000..944d1a0c0
--- /dev/null
+++ b/libquadmath/printf/divrem.c
@@ -0,0 +1,244 @@
+/* mpn_divrem -- Divide natural numbers, producing both remainder and
+   quotient.
+
+Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+/* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
+   the NSIZE-DSIZE least significant quotient limbs at QP
+   and the DSIZE long remainder at NP.  If QEXTRA_LIMBS is
+   non-zero, generate that many fraction bits and append them after the
+   other quotient limbs.
+   Return the most significant limb of the quotient, this is always 0 or 1.
+
+   Preconditions:
+   0. NSIZE >= DSIZE.
+   1. The most significant bit of the divisor must be set.
+   2. QP must either not overlap with the input operands at all, or
+      QP + DSIZE >= NP must hold true.  (This means that it's
+      possible to put the quotient in the high part of NUM, right after the
+      remainder in NUM.
+   3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero.  */
+
+mp_limb_t
+#if __STDC__
+mpn_divrem (mp_ptr qp, mp_size_t qextra_limbs,
+	    mp_ptr np, mp_size_t nsize,
+	    mp_srcptr dp, mp_size_t dsize)
+#else
+mpn_divrem (qp, qextra_limbs, np, nsize, dp, dsize)
+     mp_ptr qp;
+     mp_size_t qextra_limbs;
+     mp_ptr np;
+     mp_size_t nsize;
+     mp_srcptr dp;
+     mp_size_t dsize;
+#endif
+{
+  mp_limb_t most_significant_q_limb = 0;
+
+  switch (dsize)
+    {
+    case 0:
+      /* We are asked to divide by zero, so go ahead and do it!  (To make
+	 the compiler not remove this statement, return the value.)  */
+      return 1 / dsize;
+
+    case 1:
+      {
+	mp_size_t i;
+	mp_limb_t n1;
+	mp_limb_t d;
+
+	d = dp[0];
+	n1 = np[nsize - 1];
+
+	if (n1 >= d)
+	  {
+	    n1 -= d;
+	    most_significant_q_limb = 1;
+	  }
+
+	qp += qextra_limbs;
+	for (i = nsize - 2; i >= 0; i--)
+	  udiv_qrnnd (qp[i], n1, n1, np[i], d);
+	qp -= qextra_limbs;
+
+	for (i = qextra_limbs - 1; i >= 0; i--)
+	  udiv_qrnnd (qp[i], n1, n1, 0, d);
+
+	np[0] = n1;
+      }
+      break;
+
+    case 2:
+      {
+	mp_size_t i;
+	mp_limb_t n1, n0, n2;
+	mp_limb_t d1, d0;
+
+	np += nsize - 2;
+	d1 = dp[1];
+	d0 = dp[0];
+	n1 = np[1];
+	n0 = np[0];
+
+	if (n1 >= d1 && (n1 > d1 || n0 >= d0))
+	  {
+	    sub_ddmmss (n1, n0, n1, n0, d1, d0);
+	    most_significant_q_limb = 1;
+	  }
+
+	for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--)
+	  {
+	    mp_limb_t q;
+	    mp_limb_t r;
+
+	    if (i >= qextra_limbs)
+	      np--;
+	    else
+	      np[0] = 0;
+
+	    if (n1 == d1)
+	      {
+		/* Q should be either 111..111 or 111..110.  Need special
+		   treatment of this rare case as normal division would
+		   give overflow.  */
+		q = ~(mp_limb_t) 0;
+
+		r = n0 + d1;
+		if (r < d1)	/* Carry in the addition? */
+		  {
+		    add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
+		    qp[i] = q;
+		    continue;
+		  }
+		n1 = d0 - (d0 != 0);
+		n0 = -d0;
+	      }
+	    else
+	      {
+		udiv_qrnnd (q, r, n1, n0, d1);
+		umul_ppmm (n1, n0, d0, q);
+	      }
+
+	    n2 = np[0];
+	  q_test:
+	    if (n1 > r || (n1 == r && n0 > n2))
+	      {
+		/* The estimated Q was too large.  */
+		q--;
+
+		sub_ddmmss (n1, n0, n1, n0, 0, d0);
+		r += d1;
+		if (r >= d1)	/* If not carry, test Q again.  */
+		  goto q_test;
+	      }
+
+	    qp[i] = q;
+	    sub_ddmmss (n1, n0, r, n2, n1, n0);
+	  }
+	np[1] = n1;
+	np[0] = n0;
+      }
+      break;
+
+    default:
+      {
+	mp_size_t i;
+	mp_limb_t dX, d1, n0;
+
+	np += nsize - dsize;
+	dX = dp[dsize - 1];
+	d1 = dp[dsize - 2];
+	n0 = np[dsize - 1];
+
+	if (n0 >= dX)
+	  {
+	    if (n0 > dX || mpn_cmp (np, dp, dsize - 1) >= 0)
+	      {
+		mpn_sub_n (np, np, dp, dsize);
+		n0 = np[dsize - 1];
+		most_significant_q_limb = 1;
+	      }
+	  }
+
+	for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--)
+	  {
+	    mp_limb_t q;
+	    mp_limb_t n1, n2;
+	    mp_limb_t cy_limb;
+
+	    if (i >= qextra_limbs)
+	      {
+		np--;
+		n2 = np[dsize];
+	      }
+	    else
+	      {
+		n2 = np[dsize - 1];
+		MPN_COPY_DECR (np + 1, np, dsize);
+		np[0] = 0;
+	      }
+
+	    if (n0 == dX)
+	      /* This might over-estimate q, but it's probably not worth
+		 the extra code here to find out.  */
+	      q = ~(mp_limb_t) 0;
+	    else
+	      {
+		mp_limb_t r;
+
+		udiv_qrnnd (q, r, n0, np[dsize - 1], dX);
+		umul_ppmm (n1, n0, d1, q);
+
+		while (n1 > r || (n1 == r && n0 > np[dsize - 2]))
+		  {
+		    q--;
+		    r += dX;
+		    if (r < dX)	/* I.e. "carry in previous addition?"  */
+		      break;
+		    n1 -= n0 < d1;
+		    n0 -= d1;
+		  }
+	      }
+
+	    /* Possible optimization: We already have (q * n0) and (1 * n1)
+	       after the calculation of q.  Taking advantage of that, we
+	       could make this loop make two iterations less.  */
+
+	    cy_limb = mpn_submul_1 (np, dp, dsize, q);
+
+	    if (n2 != cy_limb)
+	      {
+		mpn_add_n (np, np, dp, dsize);
+		q--;
+	      }
+
+	    qp[i] = q;
+	    n0 = np[dsize - 1];
+	  }
+      }
+    }
+
+  return most_significant_q_limb;
+}
diff --git a/libquadmath/printf/flt1282mpn.c b/libquadmath/printf/flt1282mpn.c
new file mode 100644
index 000000000..0105314ef
--- /dev/null
+++ b/libquadmath/printf/flt1282mpn.c
@@ -0,0 +1,136 @@
+/* Copyright (C) 1995,1996,1997,1998,1999,2002,2003
+	Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <float.h>
+#include <math.h>
+#include <stdlib.h>
+#include "gmp-impl.h"
+
+/* Convert a `__float128' in IEEE854 quad-precision format to a
+   multi-precision integer representing the significand scaled up by its
+   number of bits (113 for long double) and an integral power of two
+   (MPN frexpl). */
+
+mp_size_t
+mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size,
+		      int *expt, int *is_neg,
+		      __float128 value)
+{
+  ieee854_float128 u;
+  u.value = value;
+
+  *is_neg = u.ieee.negative;
+  *expt = (int) u.ieee.exponent - IEEE854_FLOAT128_BIAS;
+
+#if BITS_PER_MP_LIMB == 32
+  res_ptr[0] = u.ieee.mant_low; /* Low-order 32 bits of fraction.  */
+  res_ptr[1] = (u.ieee.mant_low >> 32);
+  res_ptr[2] = u.ieee.mant_high;
+  res_ptr[3] = (u.ieee.mant_high >> 32); /* High-order 32 bits.  */
+  #define N 4
+#elif BITS_PER_MP_LIMB == 64
+  res_ptr[0] = u.ieee.mant_low;
+  res_ptr[1] = u.ieee.mant_high;
+  #define N 2
+#else
+  #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+/* The format does not fill the last limb.  There are some zeros.  */
+#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
+			   - (FLT128_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB)))
+
+  if (u.ieee.exponent == 0)
+    {
+      /* A biased exponent of zero is a special case.
+	 Either it is a zero or it is a denormal number.  */
+      if (res_ptr[0] == 0 && res_ptr[1] == 0
+	  && res_ptr[N - 2] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=4.  */
+	/* It's zero.  */
+	*expt = 0;
+      else
+	{
+	  /* It is a denormal number, meaning it has no implicit leading
+  	     one bit, and its exponent is in fact the format minimum.  */
+	  int cnt;
+
+#if N == 2
+	  if (res_ptr[N - 1] != 0)
+	    {
+	      count_leading_zeros (cnt, res_ptr[N - 1]);
+	      cnt -= NUM_LEADING_ZEROS;
+	      res_ptr[N - 1] = res_ptr[N - 1] << cnt
+			       | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
+	      res_ptr[0] <<= cnt;
+	      *expt = FLT128_MIN_EXP - 1 - cnt;
+	    }
+	  else
+	    {
+	      count_leading_zeros (cnt, res_ptr[0]);
+	      if (cnt >= NUM_LEADING_ZEROS)
+		{
+		  res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS);
+		  res_ptr[0] = 0;
+		}
+	      else
+		{
+		  res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt);
+		  res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt);
+		}
+	      *expt = FLT128_MIN_EXP - 1
+		- (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt;
+	    }
+#else
+	  int j, k, l;
+
+	  for (j = N - 1; j > 0; j--)
+	    if (res_ptr[j] != 0)
+	      break;
+
+	  count_leading_zeros (cnt, res_ptr[j]);
+	  cnt -= NUM_LEADING_ZEROS;
+	  l = N - 1 - j;
+	  if (cnt < 0)
+	    {
+	      cnt += BITS_PER_MP_LIMB;
+	      l--;
+	    }
+	  if (!cnt)
+	    for (k = N - 1; k >= l; k--)
+	      res_ptr[k] = res_ptr[k-l];
+	  else
+	    {
+	      for (k = N - 1; k > l; k--)
+		res_ptr[k] = res_ptr[k-l] << cnt
+			     | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt);
+	      res_ptr[k--] = res_ptr[0] << cnt;
+	    }
+
+	  for (; k >= 0; k--)
+	    res_ptr[k] = 0;
+	  *expt = FLT128_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt;
+#endif
+	}
+    }
+  else
+    /* Add the implicit leading one bit for a normalized number.  */
+    res_ptr[N - 1] |= (mp_limb_t) 1 << (FLT128_MANT_DIG - 1
+					- ((N - 1) * BITS_PER_MP_LIMB));
+
+  return N;
+}
diff --git a/libquadmath/printf/fpioconst.c b/libquadmath/printf/fpioconst.c
new file mode 100644
index 000000000..8c67e6f90
--- /dev/null
+++ b/libquadmath/printf/fpioconst.c
@@ -0,0 +1,463 @@
+/* Table of MP integer constants 10^(2^i), used for floating point <-> decimal.
+   Copyright (C) 1995, 1996, 1997, 1998, 1999, 2002, 2003
+   Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <config.h>
+#include "gmp-impl.h"		/* This defines BITS_PER_MP_LIMB.  */
+#include "fpioconst.h"
+
+/* First page	: 32-bit limbs
+   Second page	: 64-bit limbs
+   Last page	: table of pointers
+ */
+
+#if BITS_PER_MP_LIMB == 32
+
+/* Table with constants of 10^(2^i), i=0..12 for 32-bit limbs.	*/
+
+const mp_limb_t __tens[] =
+{
+#define TENS_P0_IDX	0
+#define TENS_P0_SIZE	3
+  [TENS_P0_IDX] = 0x00000000, 0x00000000, 0x0000000a,
+
+#define TENS_P1_IDX	(TENS_P0_IDX + TENS_P0_SIZE)
+#define TENS_P1_SIZE	3
+  [TENS_P1_IDX] = 0x00000000, 0x00000000, 0x00000064,
+
+#define TENS_P2_IDX	(TENS_P1_IDX + TENS_P1_SIZE)
+#define TENS_P2_SIZE	3
+  [TENS_P2_IDX] = 0x00000000, 0x00000000, 0x00002710,
+
+#define TENS_P3_IDX	(TENS_P2_IDX + TENS_P2_SIZE)
+#define TENS_P3_SIZE	3
+  [TENS_P3_IDX] = 0x00000000, 0x00000000, 0x05f5e100,
+
+#define TENS_P4_IDX	(TENS_P3_IDX + TENS_P3_SIZE)
+#define TENS_P4_SIZE	4
+  [TENS_P4_IDX] = 0x00000000, 0x00000000, 0x6fc10000, 0x002386f2,
+
+#define TENS_P5_IDX	(TENS_P4_IDX + TENS_P4_SIZE)
+#define TENS_P5_SIZE	6
+  [TENS_P5_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x85acef81, 0x2d6d415b,
+  0x000004ee,
+
+#define TENS_P6_IDX	(TENS_P5_IDX + TENS_P5_SIZE)
+#define TENS_P6_SIZE	9
+  [TENS_P6_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0xbf6a1f01,
+  0x6e38ed64, 0xdaa797ed, 0xe93ff9f4, 0x00184f03,
+
+#define TENS_P7_IDX	(TENS_P6_IDX + TENS_P6_SIZE)
+#define TENS_P7_SIZE	16
+  [TENS_P7_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x2e953e01, 0x03df9909, 0x0f1538fd, 0x2374e42f, 0xd3cff5ec,
+  0xc404dc08, 0xbccdb0da, 0xa6337f19, 0xe91f2603, 0x0000024e,
+
+#define TENS_P8_IDX	(TENS_P7_IDX + TENS_P7_SIZE)
+#define TENS_P8_SIZE	29
+  [TENS_P8_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x982e7c01,
+  0xbed3875b, 0xd8d99f72, 0x12152f87, 0x6bde50c6, 0xcf4a6e70, 0xd595d80f,
+  0x26b2716e, 0xadc666b0, 0x1d153624, 0x3c42d35a, 0x63ff540e, 0xcc5573c0,
+  0x65f9ef17, 0x55bc28f2, 0x80dcc7f7, 0xf46eeddc, 0x5fdcefce, 0x000553f7,
+
+#ifndef __NO_LONG_DOUBLE_MATH
+# define TENS_P9_IDX	(TENS_P8_IDX + TENS_P8_SIZE)
+# define TENS_P9_SIZE	56
+  [TENS_P9_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0xfc6cf801, 0x77f27267, 0x8f9546dc, 0x5d96976f, 0xb83a8a97,
+  0xc31e1ad9, 0x46c40513, 0x94e65747, 0xc88976c1, 0x4475b579, 0x28f8733b,
+  0xaa1da1bf, 0x703ed321, 0x1e25cfea, 0xb21a2f22, 0xbc51fb2e, 0x96e14f5d,
+  0xbfa3edac, 0x329c57ae, 0xe7fc7153, 0xc3fc0695, 0x85a91924, 0xf95f635e,
+  0xb2908ee0, 0x93abade4, 0x1366732a, 0x9449775c, 0x69be5b0e, 0x7343afac,
+  0xb099bc81, 0x45a71d46, 0xa2699748, 0x8cb07303, 0x8a0b1f13, 0x8cab8a97,
+  0xc1d238d9, 0x633415d4, 0x0000001c,
+
+# define TENS_P10_IDX	(TENS_P9_IDX + TENS_P9_SIZE)
+# define TENS_P10_SIZE	109
+  [TENS_P10_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x2919f001,
+  0xf55b2b72, 0x6e7c215b, 0x1ec29f86, 0x991c4e87, 0x15c51a88, 0x140ac535,
+  0x4c7d1e1a, 0xcc2cd819, 0x0ed1440e, 0x896634ee, 0x7de16cfb, 0x1e43f61f,
+  0x9fce837d, 0x231d2b9c, 0x233e55c7, 0x65dc60d7, 0xf451218b, 0x1c5cd134,
+  0xc9635986, 0x922bbb9f, 0xa7e89431, 0x9f9f2a07, 0x62be695a, 0x8e1042c4,
+  0x045b7a74, 0x1abe1de3, 0x8ad822a5, 0xba34c411, 0xd814b505, 0xbf3fdeb3,
+  0x8fc51a16, 0xb1b896bc, 0xf56deeec, 0x31fb6bfd, 0xb6f4654b, 0x101a3616,
+  0x6b7595fb, 0xdc1a47fe, 0x80d98089, 0x80bda5a5, 0x9a202882, 0x31eb0f66,
+  0xfc8f1f90, 0x976a3310, 0xe26a7b7e, 0xdf68368a, 0x3ce3a0b8, 0x8e4262ce,
+  0x75a351a2, 0x6cb0b6c9, 0x44597583, 0x31b5653f, 0xc356e38a, 0x35faaba6,
+  0x0190fba0, 0x9fc4ed52, 0x88bc491b, 0x1640114a, 0x005b8041, 0xf4f3235e,
+  0x1e8d4649, 0x36a8de06, 0x73c55349, 0xa7e6bd2a, 0xc1a6970c, 0x47187094,
+  0xd2db49ef, 0x926c3f5b, 0xae6209d4, 0x2d433949, 0x34f4a3c6, 0xd4305d94,
+  0xd9d61a05, 0x00000325,
+
+# define TENS_P11_IDX	(TENS_P10_IDX + TENS_P10_SIZE)
+# define TENS_P11_SIZE	215
+  [TENS_P11_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x1333e001, 0xe3096865, 0xb27d4d3f, 0x49e28dcf, 0xec2e4721,
+  0xee87e354, 0xb6067584, 0x368b8abb, 0xa5e5a191, 0x2ed56d55, 0xfd827773,
+  0xea50d142, 0x51b78db2, 0x98342c9e, 0xc850dabc, 0x866ed6f1, 0x19342c12,
+  0x92794987, 0xd2f869c2, 0x66912e4a, 0x71c7fd8f, 0x57a7842d, 0x235552eb,
+  0xfb7fedcc, 0xf3861ce0, 0x38209ce1, 0x9713b449, 0x34c10134, 0x8c6c54de,
+  0xa7a8289c, 0x2dbb6643, 0xe3cb64f3, 0x8074ff01, 0xe3892ee9, 0x10c17f94,
+  0xa8f16f92, 0xa8281ed6, 0x967abbb3, 0x5a151440, 0x9952fbed, 0x13b41e44,
+  0xafe609c3, 0xa2bca416, 0xf111821f, 0xfb1264b4, 0x91bac974, 0xd6c7d6ab,
+  0x8e48ff35, 0x4419bd43, 0xc4a65665, 0x685e5510, 0x33554c36, 0xab498697,
+  0x0dbd21fe, 0x3cfe491d, 0x982da466, 0xcbea4ca7, 0x9e110c7b, 0x79c56b8a,
+  0x5fc5a047, 0x84d80e2e, 0x1aa9f444, 0x730f203c, 0x6a57b1ab, 0xd752f7a6,
+  0x87a7dc62, 0x944545ff, 0x40660460, 0x77c1a42f, 0xc9ac375d, 0xe866d7ef,
+  0x744695f0, 0x81428c85, 0xa1fc6b96, 0xd7917c7b, 0x7bf03c19, 0x5b33eb41,
+  0x5715f791, 0x8f6cae5f, 0xdb0708fd, 0xb125ac8e, 0x785ce6b7, 0x56c6815b,
+  0x6f46eadb, 0x4eeebeee, 0x195355d8, 0xa244de3c, 0x9d7389c0, 0x53761abd,
+  0xcf99d019, 0xde9ec24b, 0x0d76ce39, 0x70beb181, 0x2e55ecee, 0xd5f86079,
+  0xf56d9d4b, 0xfb8886fb, 0x13ef5a83, 0x408f43c5, 0x3f3389a4, 0xfad37943,
+  0x58ccf45c, 0xf82df846, 0x415c7f3e, 0x2915e818, 0x8b3d5cf4, 0x6a445f27,
+  0xf8dbb57a, 0xca8f0070, 0x8ad803ec, 0xb2e87c34, 0x038f9245, 0xbedd8a6c,
+  0xc7c9dee0, 0x0eac7d56, 0x2ad3fa14, 0xe0de0840, 0xf775677c, 0xf1bd0ad5,
+  0x92be221e, 0x87fa1fb9, 0xce9d04a4, 0xd2c36fa9, 0x3f6f7024, 0xb028af62,
+  0x907855ee, 0xd83e49d6, 0x4efac5dc, 0xe7151aab, 0x77cd8c6b, 0x0a753b7d,
+  0x0af908b4, 0x8c983623, 0xe50f3027, 0x94222771, 0x1d08e2d6, 0xf7e928e6,
+  0xf2ee5ca6, 0x1b61b93c, 0x11eb962b, 0x9648b21c, 0xce2bcba1, 0x34f77154,
+  0x7bbebe30, 0xe526a319, 0x8ce329ac, 0xde4a74d2, 0xb5dc53d5, 0x0009e8b3,
+
+# define TENS_P12_IDX	(TENS_P11_IDX + TENS_P11_SIZE)
+# define TENS_P12_SIZE	428
+  [TENS_P12_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
+  0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x2a67c001,
+  0xd4724e8d, 0x8efe7ae7, 0xf89a1e90, 0xef084117, 0x54e05154, 0x13b1bb51,
+  0x506be829, 0xfb29b172, 0xe599574e, 0xf0da6146, 0x806c0ed3, 0xb86ae5be,
+  0x45155e93, 0xc0591cc2, 0x7e1e7c34, 0x7c4823da, 0x1d1f4cce, 0x9b8ba1e8,
+  0xd6bfdf75, 0xe341be10, 0xc2dfae78, 0x016b67b2, 0x0f237f1a, 0x3dbeabcd,
+  0xaf6a2574, 0xcab3e6d7, 0x142e0e80, 0x61959127, 0x2c234811, 0x87009701,
+  0xcb4bf982, 0xf8169c84, 0x88052f8c, 0x68dde6d4, 0xbc131761, 0xff0b0905,
+  0x54ab9c41, 0x7613b224, 0x1a1c304e, 0x3bfe167b, 0x441c2d47, 0x4f6cea9c,
+  0x78f06181, 0xeb659fb8, 0x30c7ae41, 0x947e0d0e, 0xa1ebcad7, 0xd97d9556,
+  0x2130504d, 0x1a8309cb, 0xf2acd507, 0x3f8ec72a, 0xfd82373a, 0x95a842bc,
+  0x280f4d32, 0xf3618ac0, 0x811a4f04, 0x6dc3a5b4, 0xd3967a1b, 0x15b8c898,
+  0xdcfe388f, 0x454eb2a0, 0x8738b909, 0x10c4e996, 0x2bd9cc11, 0x3297cd0c,
+  0x655fec30, 0xae0725b1, 0xf4090ee8, 0x037d19ee, 0x398c6fed, 0x3b9af26b,
+  0xc994a450, 0xb5341743, 0x75a697b2, 0xac50b9c1, 0x3ccb5b92, 0xffe06205,
+  0xa8329761, 0xdfea5242, 0xeb83cadb, 0xe79dadf7, 0x3c20ee69, 0x1e0a6817,
+  0x7021b97a, 0x743074fa, 0x176ca776, 0x77fb8af6, 0xeca19beb, 0x92baf1de,
+  0xaf63b712, 0xde35c88b, 0xa4eb8f8c, 0xe137d5e9, 0x40b464a0, 0x87d1cde8,
+  0x42923bbd, 0xcd8f62ff, 0x2e2690f3, 0x095edc16, 0x59c89f1b, 0x1fa8fd5d,
+  0x5138753d, 0x390a2b29, 0x80152f18, 0x2dd8d925, 0xf984d83e, 0x7a872e74,
+  0xc19e1faf, 0xed4d542d, 0xecf9b5d0, 0x9462ea75, 0xc53c0adf, 0x0caea134,
+  0x37a2d439, 0xc8fa2e8a, 0x2181327e, 0x6e7bb827, 0x2d240820, 0x50be10e0,
+  0x5893d4b8, 0xab312bb9, 0x1f2b2322, 0x440b3f25, 0xbf627ede, 0x72dac789,
+  0xb608b895, 0x78787e2a, 0x86deb3f0, 0x6fee7aab, 0xbb9373f4, 0x27ecf57b,
+  0xf7d8b57e, 0xfca26a9f, 0x3d04e8d2, 0xc9df13cb, 0x3172826a, 0xcd9e8d7c,
+  0xa8fcd8e0, 0xb2c39497, 0x307641d9, 0x1cc939c1, 0x2608c4cf, 0xb6d1c7bf,
+  0x3d326a7e, 0xeeaf19e6, 0x8e13e25f, 0xee63302b, 0x2dfe6d97, 0x25971d58,
+  0xe41d3cc4, 0x0a80627c, 0xab8db59a, 0x9eea37c8, 0xe90afb77, 0x90ca19cf,
+  0x9ee3352c, 0x3613c850, 0xfe78d682, 0x788f6e50, 0x5b060904, 0xb71bd1a4,
+  0x3fecb534, 0xb32c450c, 0x20c33857, 0xa6e9cfda, 0x0239f4ce, 0x48497187,
+  0xa19adb95, 0xb492ed8a, 0x95aca6a8, 0x4dcd6cd9, 0xcf1b2350, 0xfbe8b12a,
+  0x1a67778c, 0x38eb3acc, 0xc32da383, 0xfb126ab1, 0xa03f40a8, 0xed5bf546,
+  0xe9ce4724, 0x4c4a74fd, 0x73a130d8, 0xd9960e2d, 0xa2ebd6c1, 0x94ab6feb,
+  0x6f233b7c, 0x49126080, 0x8e7b9a73, 0x4b8c9091, 0xd298f999, 0x35e836b5,
+  0xa96ddeff, 0x96119b31, 0x6b0dd9bc, 0xc6cc3f8d, 0x282566fb, 0x72b882e7,
+  0xd6769f3b, 0xa674343d, 0x00fc509b, 0xdcbf7789, 0xd6266a3f, 0xae9641fd,
+  0x4e89541b, 0x11953407, 0x53400d03, 0x8e0dd75a, 0xe5b53345, 0x108f19ad,
+  0x108b89bc, 0x41a4c954, 0xe03b2b63, 0x437b3d7f, 0x97aced8e, 0xcbd66670,
+  0x2c5508c2, 0x650ebc69, 0x5c4f2ef0, 0x904ff6bf, 0x9985a2df, 0x9faddd9e,
+  0x5ed8d239, 0x25585832, 0xe3e51cb9, 0x0ff4f1d4, 0x56c02d9a, 0x8c4ef804,
+  0xc1a08a13, 0x13fd01c8, 0xe6d27671, 0xa7c234f4, 0x9d0176cc, 0xd0d73df2,
+  0x4d8bfa89, 0x544f10cd, 0x2b17e0b2, 0xb70a5c7d, 0xfd86fe49, 0xdf373f41,
+  0x214495bb, 0x84e857fd, 0x00d313d5, 0x0496fcbe, 0xa4ba4744, 0xe8cac982,
+  0xaec29e6e, 0x87ec7038, 0x7000a519, 0xaeee333b, 0xff66e42c, 0x8afd6b25,
+  0x03b4f63b, 0xbd7991dc, 0x5ab8d9c7, 0x2ed4684e, 0x48741a6c, 0xaf06940d,
+  0x2fdc6349, 0xb03d7ecd, 0xe974996f, 0xac7867f9, 0x52ec8721, 0xbcdd9d4a,
+  0x8edd2d00, 0x3557de06, 0x41c759f8, 0x3956d4b9, 0xa75409f2, 0x123cd8a1,
+  0xb6100fab, 0x3e7b21e2, 0x2e8d623b, 0x92959da2, 0xbca35f77, 0x200c03a5,
+  0x35fcb457, 0x1bb6c6e4, 0xf74eb928, 0x3d5d0b54, 0x87cc1d21, 0x4964046f,
+  0x18ae4240, 0xd868b275, 0x8bd2b496, 0x1c5563f4, 0xc234d8f5, 0xf868e970,
+  0xf9151fff, 0xae7be4a2, 0x271133ee, 0xbb0fd922, 0x25254932, 0xa60a9fc0,
+  0x104bcd64, 0x30290145, 0x00000062
+#endif	/* !__NO_LONG_DOUBLE_MATH */
+};
+
+#elif BITS_PER_MP_LIMB == 64
+
+/* Table with constants of 10^(2^i), i=0..12 for 64-bit limbs.	*/
+
+const mp_limb_t __tens[] =
+{
+#define TENS_P0_IDX	0
+#define TENS_P0_SIZE	2
+  [TENS_P0_IDX] = 0x0000000000000000ull, 0x000000000000000aull,
+
+#define TENS_P1_IDX	(TENS_P0_IDX + TENS_P0_SIZE)
+#define TENS_P1_SIZE	2
+  [TENS_P1_IDX] = 0x0000000000000000ull, 0x0000000000000064ull,
+
+#define TENS_P2_IDX	(TENS_P1_IDX + TENS_P1_SIZE)
+#define TENS_P2_SIZE	2
+  [TENS_P2_IDX] = 0x0000000000000000ull, 0x0000000000002710ull,
+
+#define TENS_P3_IDX	(TENS_P2_IDX + TENS_P2_SIZE)
+#define TENS_P3_SIZE	2
+  [TENS_P3_IDX] = 0x0000000000000000ull, 0x0000000005f5e100ull,
+
+#define TENS_P4_IDX	(TENS_P3_IDX + TENS_P3_SIZE)
+#define TENS_P4_SIZE	2
+  [TENS_P4_IDX] = 0x0000000000000000ull, 0x002386f26fc10000ull,
+
+#define TENS_P5_IDX	(TENS_P4_IDX + TENS_P4_SIZE)
+#define TENS_P5_SIZE	3
+  [TENS_P5_IDX] = 0x0000000000000000ull, 0x85acef8100000000ull,
+  0x000004ee2d6d415bull,
+
+#define TENS_P6_IDX	(TENS_P5_IDX + TENS_P5_SIZE)
+#define TENS_P6_SIZE	5
+  [TENS_P6_IDX] = 0x0000000000000000ull, 0x0000000000000000ull,
+  0x6e38ed64bf6a1f01ull, 0xe93ff9f4daa797edull, 0x0000000000184f03ull,
+
+#define TENS_P7_IDX	(TENS_P6_IDX + TENS_P6_SIZE)
+#define TENS_P7_SIZE	8
+  [TENS_P7_IDX] = 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x03df99092e953e01ull, 0x2374e42f0f1538fdull,
+  0xc404dc08d3cff5ecull, 0xa6337f19bccdb0daull, 0x0000024ee91f2603ull,
+
+#define TENS_P8_IDX	(TENS_P7_IDX + TENS_P7_SIZE)
+#define TENS_P8_SIZE	15
+  [TENS_P8_IDX] = 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0xbed3875b982e7c01ull, 0x12152f87d8d99f72ull, 0xcf4a6e706bde50c6ull,
+  0x26b2716ed595d80full, 0x1d153624adc666b0ull, 0x63ff540e3c42d35aull,
+  0x65f9ef17cc5573c0ull, 0x80dcc7f755bc28f2ull, 0x5fdcefcef46eeddcull,
+  0x00000000000553f7ull,
+#if FLT128_MAX_EXP > 1024
+# define TENS_P9_IDX	(TENS_P8_IDX + TENS_P8_SIZE)
+# define TENS_P9_SIZE	28
+  [TENS_P9_IDX] = 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x77f27267fc6cf801ull, 0x5d96976f8f9546dcull,
+  0xc31e1ad9b83a8a97ull, 0x94e6574746c40513ull, 0x4475b579c88976c1ull,
+  0xaa1da1bf28f8733bull, 0x1e25cfea703ed321ull, 0xbc51fb2eb21a2f22ull,
+  0xbfa3edac96e14f5dull, 0xe7fc7153329c57aeull, 0x85a91924c3fc0695ull,
+  0xb2908ee0f95f635eull, 0x1366732a93abade4ull, 0x69be5b0e9449775cull,
+  0xb099bc817343afacull, 0xa269974845a71d46ull, 0x8a0b1f138cb07303ull,
+  0xc1d238d98cab8a97ull, 0x0000001c633415d4ull,
+
+# define TENS_P10_IDX	(TENS_P9_IDX + TENS_P9_SIZE)
+# define TENS_P10_SIZE	55
+  [TENS_P10_IDX] = 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0xf55b2b722919f001ull, 0x1ec29f866e7c215bull, 0x15c51a88991c4e87ull,
+  0x4c7d1e1a140ac535ull, 0x0ed1440ecc2cd819ull, 0x7de16cfb896634eeull,
+  0x9fce837d1e43f61full, 0x233e55c7231d2b9cull, 0xf451218b65dc60d7ull,
+  0xc96359861c5cd134ull, 0xa7e89431922bbb9full, 0x62be695a9f9f2a07ull,
+  0x045b7a748e1042c4ull, 0x8ad822a51abe1de3ull, 0xd814b505ba34c411ull,
+  0x8fc51a16bf3fdeb3ull, 0xf56deeecb1b896bcull, 0xb6f4654b31fb6bfdull,
+  0x6b7595fb101a3616ull, 0x80d98089dc1a47feull, 0x9a20288280bda5a5ull,
+  0xfc8f1f9031eb0f66ull, 0xe26a7b7e976a3310ull, 0x3ce3a0b8df68368aull,
+  0x75a351a28e4262ceull, 0x445975836cb0b6c9ull, 0xc356e38a31b5653full,
+  0x0190fba035faaba6ull, 0x88bc491b9fc4ed52ull, 0x005b80411640114aull,
+  0x1e8d4649f4f3235eull, 0x73c5534936a8de06ull, 0xc1a6970ca7e6bd2aull,
+  0xd2db49ef47187094ull, 0xae6209d4926c3f5bull, 0x34f4a3c62d433949ull,
+  0xd9d61a05d4305d94ull, 0x0000000000000325ull,
+
+# define TENS_P11_IDX	(TENS_P10_IDX + TENS_P10_SIZE)
+# define TENS_P11_SIZE	108
+  [TENS_P11_IDX] = 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0xe30968651333e001ull, 0x49e28dcfb27d4d3full,
+  0xee87e354ec2e4721ull, 0x368b8abbb6067584ull, 0x2ed56d55a5e5a191ull,
+  0xea50d142fd827773ull, 0x98342c9e51b78db2ull, 0x866ed6f1c850dabcull,
+  0x9279498719342c12ull, 0x66912e4ad2f869c2ull, 0x57a7842d71c7fd8full,
+  0xfb7fedcc235552ebull, 0x38209ce1f3861ce0ull, 0x34c101349713b449ull,
+  0xa7a8289c8c6c54deull, 0xe3cb64f32dbb6643ull, 0xe3892ee98074ff01ull,
+  0xa8f16f9210c17f94ull, 0x967abbb3a8281ed6ull, 0x9952fbed5a151440ull,
+  0xafe609c313b41e44ull, 0xf111821fa2bca416ull, 0x91bac974fb1264b4ull,
+  0x8e48ff35d6c7d6abull, 0xc4a656654419bd43ull, 0x33554c36685e5510ull,
+  0x0dbd21feab498697ull, 0x982da4663cfe491dull, 0x9e110c7bcbea4ca7ull,
+  0x5fc5a04779c56b8aull, 0x1aa9f44484d80e2eull, 0x6a57b1ab730f203cull,
+  0x87a7dc62d752f7a6ull, 0x40660460944545ffull, 0xc9ac375d77c1a42full,
+  0x744695f0e866d7efull, 0xa1fc6b9681428c85ull, 0x7bf03c19d7917c7bull,
+  0x5715f7915b33eb41ull, 0xdb0708fd8f6cae5full, 0x785ce6b7b125ac8eull,
+  0x6f46eadb56c6815bull, 0x195355d84eeebeeeull, 0x9d7389c0a244de3cull,
+  0xcf99d01953761abdull, 0x0d76ce39de9ec24bull, 0x2e55ecee70beb181ull,
+  0xf56d9d4bd5f86079ull, 0x13ef5a83fb8886fbull, 0x3f3389a4408f43c5ull,
+  0x58ccf45cfad37943ull, 0x415c7f3ef82df846ull, 0x8b3d5cf42915e818ull,
+  0xf8dbb57a6a445f27ull, 0x8ad803ecca8f0070ull, 0x038f9245b2e87c34ull,
+  0xc7c9dee0bedd8a6cull, 0x2ad3fa140eac7d56ull, 0xf775677ce0de0840ull,
+  0x92be221ef1bd0ad5ull, 0xce9d04a487fa1fb9ull, 0x3f6f7024d2c36fa9ull,
+  0x907855eeb028af62ull, 0x4efac5dcd83e49d6ull, 0x77cd8c6be7151aabull,
+  0x0af908b40a753b7dull, 0xe50f30278c983623ull, 0x1d08e2d694222771ull,
+  0xf2ee5ca6f7e928e6ull, 0x11eb962b1b61b93cull, 0xce2bcba19648b21cull,
+  0x7bbebe3034f77154ull, 0x8ce329ace526a319ull, 0xb5dc53d5de4a74d2ull,
+  0x000000000009e8b3ull,
+
+# define TENS_P12_IDX	(TENS_P11_IDX + TENS_P11_SIZE)
+# define TENS_P12_SIZE	214
+  [TENS_P12_IDX] = 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull,
+  0xd4724e8d2a67c001ull, 0xf89a1e908efe7ae7ull, 0x54e05154ef084117ull,
+  0x506be82913b1bb51ull, 0xe599574efb29b172ull, 0x806c0ed3f0da6146ull,
+  0x45155e93b86ae5beull, 0x7e1e7c34c0591cc2ull, 0x1d1f4cce7c4823daull,
+  0xd6bfdf759b8ba1e8ull, 0xc2dfae78e341be10ull, 0x0f237f1a016b67b2ull,
+  0xaf6a25743dbeabcdull, 0x142e0e80cab3e6d7ull, 0x2c23481161959127ull,
+  0xcb4bf98287009701ull, 0x88052f8cf8169c84ull, 0xbc13176168dde6d4ull,
+  0x54ab9c41ff0b0905ull, 0x1a1c304e7613b224ull, 0x441c2d473bfe167bull,
+  0x78f061814f6cea9cull, 0x30c7ae41eb659fb8ull, 0xa1ebcad7947e0d0eull,
+  0x2130504dd97d9556ull, 0xf2acd5071a8309cbull, 0xfd82373a3f8ec72aull,
+  0x280f4d3295a842bcull, 0x811a4f04f3618ac0ull, 0xd3967a1b6dc3a5b4ull,
+  0xdcfe388f15b8c898ull, 0x8738b909454eb2a0ull, 0x2bd9cc1110c4e996ull,
+  0x655fec303297cd0cull, 0xf4090ee8ae0725b1ull, 0x398c6fed037d19eeull,
+  0xc994a4503b9af26bull, 0x75a697b2b5341743ull, 0x3ccb5b92ac50b9c1ull,
+  0xa8329761ffe06205ull, 0xeb83cadbdfea5242ull, 0x3c20ee69e79dadf7ull,
+  0x7021b97a1e0a6817ull, 0x176ca776743074faull, 0xeca19beb77fb8af6ull,
+  0xaf63b71292baf1deull, 0xa4eb8f8cde35c88bull, 0x40b464a0e137d5e9ull,
+  0x42923bbd87d1cde8ull, 0x2e2690f3cd8f62ffull, 0x59c89f1b095edc16ull,
+  0x5138753d1fa8fd5dull, 0x80152f18390a2b29ull, 0xf984d83e2dd8d925ull,
+  0xc19e1faf7a872e74ull, 0xecf9b5d0ed4d542dull, 0xc53c0adf9462ea75ull,
+  0x37a2d4390caea134ull, 0x2181327ec8fa2e8aull, 0x2d2408206e7bb827ull,
+  0x5893d4b850be10e0ull, 0x1f2b2322ab312bb9ull, 0xbf627ede440b3f25ull,
+  0xb608b89572dac789ull, 0x86deb3f078787e2aull, 0xbb9373f46fee7aabull,
+  0xf7d8b57e27ecf57bull, 0x3d04e8d2fca26a9full, 0x3172826ac9df13cbull,
+  0xa8fcd8e0cd9e8d7cull, 0x307641d9b2c39497ull, 0x2608c4cf1cc939c1ull,
+  0x3d326a7eb6d1c7bfull, 0x8e13e25feeaf19e6ull, 0x2dfe6d97ee63302bull,
+  0xe41d3cc425971d58ull, 0xab8db59a0a80627cull, 0xe90afb779eea37c8ull,
+  0x9ee3352c90ca19cfull, 0xfe78d6823613c850ull, 0x5b060904788f6e50ull,
+  0x3fecb534b71bd1a4ull, 0x20c33857b32c450cull, 0x0239f4cea6e9cfdaull,
+  0xa19adb9548497187ull, 0x95aca6a8b492ed8aull, 0xcf1b23504dcd6cd9ull,
+  0x1a67778cfbe8b12aull, 0xc32da38338eb3accull, 0xa03f40a8fb126ab1ull,
+  0xe9ce4724ed5bf546ull, 0x73a130d84c4a74fdull, 0xa2ebd6c1d9960e2dull,
+  0x6f233b7c94ab6febull, 0x8e7b9a7349126080ull, 0xd298f9994b8c9091ull,
+  0xa96ddeff35e836b5ull, 0x6b0dd9bc96119b31ull, 0x282566fbc6cc3f8dull,
+  0xd6769f3b72b882e7ull, 0x00fc509ba674343dull, 0xd6266a3fdcbf7789ull,
+  0x4e89541bae9641fdull, 0x53400d0311953407ull, 0xe5b533458e0dd75aull,
+  0x108b89bc108f19adull, 0xe03b2b6341a4c954ull, 0x97aced8e437b3d7full,
+  0x2c5508c2cbd66670ull, 0x5c4f2ef0650ebc69ull, 0x9985a2df904ff6bfull,
+  0x5ed8d2399faddd9eull, 0xe3e51cb925585832ull, 0x56c02d9a0ff4f1d4ull,
+  0xc1a08a138c4ef804ull, 0xe6d2767113fd01c8ull, 0x9d0176cca7c234f4ull,
+  0x4d8bfa89d0d73df2ull, 0x2b17e0b2544f10cdull, 0xfd86fe49b70a5c7dull,
+  0x214495bbdf373f41ull, 0x00d313d584e857fdull, 0xa4ba47440496fcbeull,
+  0xaec29e6ee8cac982ull, 0x7000a51987ec7038ull, 0xff66e42caeee333bull,
+  0x03b4f63b8afd6b25ull, 0x5ab8d9c7bd7991dcull, 0x48741a6c2ed4684eull,
+  0x2fdc6349af06940dull, 0xe974996fb03d7ecdull, 0x52ec8721ac7867f9ull,
+  0x8edd2d00bcdd9d4aull, 0x41c759f83557de06ull, 0xa75409f23956d4b9ull,
+  0xb6100fab123cd8a1ull, 0x2e8d623b3e7b21e2ull, 0xbca35f7792959da2ull,
+  0x35fcb457200c03a5ull, 0xf74eb9281bb6c6e4ull, 0x87cc1d213d5d0b54ull,
+  0x18ae42404964046full, 0x8bd2b496d868b275ull, 0xc234d8f51c5563f4ull,
+  0xf9151ffff868e970ull, 0x271133eeae7be4a2ull, 0x25254932bb0fd922ull,
+  0x104bcd64a60a9fc0ull, 0x0000006230290145ull
+#endif
+};
+
+#else
+# error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+
+/* Each of array variable above defines one mpn integer which is a power of 10.
+   This table points to those variables, indexed by the exponent.  */
+
+const struct mp_power _fpioconst_pow10[FLT128_MAX_10_EXP_LOG + 1] =
+{
+  { TENS_P0_IDX, TENS_P0_SIZE,		4,	     },
+  { TENS_P1_IDX, TENS_P1_SIZE,		7,	   4 },
+  { TENS_P2_IDX, TENS_P2_SIZE,		14,	  10 },
+  { TENS_P3_IDX, TENS_P3_SIZE,		27,	  24 },
+  { TENS_P4_IDX, TENS_P4_SIZE,		54,	  50 },
+  { TENS_P5_IDX, TENS_P5_SIZE,		107,	 103 },
+  { TENS_P6_IDX, TENS_P6_SIZE,		213,	 210 },
+  { TENS_P7_IDX, TENS_P7_SIZE,		426,	 422 },
+  { TENS_P8_IDX, TENS_P8_SIZE,  	851,	 848 },
+#if FLT128_MAX_EXP > 1024
+  { TENS_P9_IDX, TENS_P9_SIZE,	 	1701,	1698 },
+  { TENS_P10_IDX, TENS_P10_SIZE,	3402,	3399 },
+  { TENS_P11_IDX, TENS_P11_SIZE,	6804,	6800 },
+  { TENS_P12_IDX, TENS_P12_SIZE, 	13607, 13604 }
+#endif
+};
+
+#if LAST_POW10 > _LAST_POW10
+# error "Need to expand 10^(2^i) table for i up to" LAST_POW10
+#endif
diff --git a/libquadmath/printf/fpioconst.h b/libquadmath/printf/fpioconst.h
new file mode 100644
index 000000000..418755547
--- /dev/null
+++ b/libquadmath/printf/fpioconst.h
@@ -0,0 +1,64 @@
+/* Header file for constants used in floating point <-> decimal conversions.
+   Copyright (C) 1995, 1996, 1997, 1998, 1999, 2002, 2003
+   Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#ifndef _FPIOCONST_H
+#define	_FPIOCONST_H
+
+#include <float.h>
+#include <math.h>
+
+
+/* These values are used by __printf_fp, where they are noncritical (if the
+   value is not large enough, it will just be slower); and by
+   strtof/strtod/strtold, where it is critical (it's used for overflow
+   detection).
+
+   XXX These should be defined in <float.h>.  For the time being, we have the
+   IEEE754 values here.  */
+
+#define FLT128_MAX_10_EXP_LOG	12 /* = floor(log_2(FLT128_MAX_10_EXP)) */
+
+
+/* The array with the number representation. */
+#define __tens __quadmath_tens
+extern const mp_limb_t __tens[] attribute_hidden;
+
+/* Table of powers of ten.  This is used by __printf_fp and by
+   strtof/strtod/strtold.  */
+struct mp_power
+  {
+    size_t arrayoff;		/* Offset in `__tens'.  */
+    mp_size_t arraysize;	/* Size of the array.  */
+    int p_expo;			/* Exponent of the number 10^(2^i).  */
+    int m_expo;			/* Exponent of the number 10^-(2^i-1).  */
+  };
+#define _fpioconst_pow10 __quadmath_fpioconst_pow10
+extern const struct mp_power _fpioconst_pow10[FLT128_MAX_10_EXP_LOG + 1]
+     attribute_hidden;
+
+/* The constants in the array `_fpioconst_pow10' have an offset.  */
+#if BITS_PER_MP_LIMB == 32
+# define _FPIO_CONST_OFFSET	2
+#else
+# define _FPIO_CONST_OFFSET	1
+#endif
+
+
+#endif	/* fpioconst.h */
diff --git a/libquadmath/printf/gmp-impl.h b/libquadmath/printf/gmp-impl.h
new file mode 100644
index 000000000..ca49e1966
--- /dev/null
+++ b/libquadmath/printf/gmp-impl.h
@@ -0,0 +1,181 @@
+/* Include file for internal GNU MP types and definitions.
+
+Copyright (C) 1991, 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <stdlib.h>
+#include "quadmath-imp.h"
+
+#undef alloca
+#define alloca __builtin_alloca
+
+#define ABS(x) (x >= 0 ? x : -x)
+#ifndef MIN
+#define MIN(l,o) ((l) < (o) ? (l) : (o))
+#endif
+#ifndef MAX
+#define MAX(h,i) ((h) > (i) ? (h) : (i))
+#endif
+
+#define BITS_PER_MP_LIMB (__SIZEOF_LONG__ * __CHAR_BIT__)
+#define BYTES_PER_MP_LIMB (BITS_PER_MP_LIMB / __CHAR_BIT__)
+typedef unsigned long int	mp_limb_t;
+typedef long int		mp_limb_signed_t;
+
+typedef mp_limb_t *             mp_ptr;
+typedef const mp_limb_t *	mp_srcptr;
+typedef long int                mp_size_t;
+typedef long int                mp_exp_t;
+
+/* Define stuff for longlong.h.  */
+typedef unsigned int UQItype	__attribute__ ((mode (QI)));
+typedef 	 int SItype	__attribute__ ((mode (SI)));
+typedef unsigned int USItype	__attribute__ ((mode (SI)));
+typedef		 int DItype	__attribute__ ((mode (DI)));
+typedef unsigned int UDItype	__attribute__ ((mode (DI)));
+
+typedef mp_limb_t UWtype;
+typedef unsigned int UHWtype;
+#define W_TYPE_SIZE BITS_PER_MP_LIMB
+
+#ifdef HAVE_HIDDEN_VISIBILITY
+#define attribute_hidden __attribute__((__visibility__ ("hidden")))
+#else
+#define attribute_hidden
+#endif
+
+#include "../../gcc/longlong.h"
+
+/* Copy NLIMBS *limbs* from SRC to DST.  */
+#define MPN_COPY_INCR(DST, SRC, NLIMBS) \
+  do {									\
+    mp_size_t __i;							\
+    for (__i = 0; __i < (NLIMBS); __i++)				\
+      (DST)[__i] = (SRC)[__i];						\
+  } while (0)
+#define MPN_COPY_DECR(DST, SRC, NLIMBS) \
+  do {									\
+    mp_size_t __i;							\
+    for (__i = (NLIMBS) - 1; __i >= 0; __i--)				\
+      (DST)[__i] = (SRC)[__i];						\
+  } while (0)
+#define MPN_COPY MPN_COPY_INCR
+
+/* Zero NLIMBS *limbs* AT DST.  */
+#define MPN_ZERO(DST, NLIMBS) \
+  do {									\
+    mp_size_t __i;							\
+    for (__i = 0; __i < (NLIMBS); __i++)				\
+      (DST)[__i] = 0;							\
+  } while (0)
+
+#define MPN_MUL_N_RECURSE(prodp, up, vp, size, tspace) \
+  do {									\
+    if ((size) < KARATSUBA_THRESHOLD)					\
+      impn_mul_n_basecase (prodp, up, vp, size);			\
+    else								\
+      impn_mul_n (prodp, up, vp, size, tspace);			\
+  } while (0);
+
+#define __MPN(x) __quadmath_mpn_##x
+
+/* Internal mpn calls */
+#define impn_mul_n_basecase	__MPN(impn_mul_n_basecase)
+#define impn_mul_n		__MPN(impn_mul_n)
+
+/* Prototypes for internal mpn calls.  */
+void impn_mul_n_basecase (mp_ptr prodp, mp_srcptr up, mp_srcptr vp,
+			  mp_size_t size) attribute_hidden;
+void impn_mul_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size,
+		 mp_ptr tspace) attribute_hidden;
+
+#define mpn_add_n		__MPN(add_n)
+#define mpn_addmul_1		__MPN(addmul_1)
+#define mpn_cmp			__MPN(cmp)
+#define mpn_divrem		__MPN(divrem)
+#define mpn_lshift		__MPN(lshift)
+#define mpn_mul			__MPN(mul)
+#define mpn_mul_1		__MPN(mul_1)
+#define mpn_rshift		__MPN(rshift)
+#define mpn_sub_n		__MPN(sub_n)
+#define mpn_submul_1		__MPN(submul_1)
+
+mp_limb_t mpn_add_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t)
+     attribute_hidden;
+mp_limb_t mpn_addmul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t)
+     attribute_hidden;
+int mpn_cmp (mp_srcptr, mp_srcptr, mp_size_t) attribute_hidden;
+mp_limb_t mpn_divrem (mp_ptr, mp_size_t, mp_ptr, mp_size_t, mp_srcptr,
+		      mp_size_t) attribute_hidden;
+mp_limb_t mpn_lshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int)
+     attribute_hidden;
+mp_limb_t mpn_mul (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t)
+     attribute_hidden;
+mp_limb_t mpn_mul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t)
+     attribute_hidden;
+mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int)
+     attribute_hidden;
+mp_limb_t mpn_sub_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t)
+     attribute_hidden;
+mp_limb_t mpn_submul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t)
+     attribute_hidden;
+
+#define mpn_extract_flt128 __MPN(extract_flt128)
+mp_size_t mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size, int *expt,
+			      int *is_neg, __float128 value) attribute_hidden;
+
+#define mpn_construct_float128 __MPN(construct_float128)
+__float128 mpn_construct_float128 (mp_srcptr frac_ptr, int expt, int sign)
+     attribute_hidden;
+
+#define mpn_divmod(qp,np,nsize,dp,dsize) mpn_divrem (qp,0,np,nsize,dp,dsize)
+
+static inline mp_limb_t
+mpn_add_1 (register mp_ptr res_ptr,
+	   register mp_srcptr s1_ptr,
+	   register mp_size_t s1_size,
+	   register mp_limb_t s2_limb)
+{
+  register mp_limb_t x;
+
+  x = *s1_ptr++;
+  s2_limb = x + s2_limb;
+  *res_ptr++ = s2_limb;
+  if (s2_limb < x)
+    {
+      while (--s1_size != 0)
+	{
+	  x = *s1_ptr++ + 1;
+	  *res_ptr++ = x;
+	  if (x != 0)
+	    goto fin;
+	}
+
+      return 1;
+    }
+
+ fin:
+  if (res_ptr != s1_ptr)
+    {
+      mp_size_t i;
+      for (i = 0; i < s1_size - 1; i++)
+	res_ptr[i] = s1_ptr[i];
+    }
+  return 0;
+}
diff --git a/libquadmath/printf/lshift.c b/libquadmath/printf/lshift.c
new file mode 100644
index 000000000..58aa8d464
--- /dev/null
+++ b/libquadmath/printf/lshift.c
@@ -0,0 +1,87 @@
+/* mpn_lshift -- Shift left low level.
+
+Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+/* Shift U (pointed to by UP and USIZE digits long) CNT bits to the left
+   and store the USIZE least significant digits of the result at WP.
+   Return the bits shifted out from the most significant digit.
+
+   Argument constraints:
+   1. 0 < CNT < BITS_PER_MP_LIMB
+   2. If the result is to be written over the input, WP must be >= UP.
+*/
+
+mp_limb_t
+#if __STDC__
+mpn_lshift (register mp_ptr wp,
+	    register mp_srcptr up, mp_size_t usize,
+	    register unsigned int cnt)
+#else
+mpn_lshift (wp, up, usize, cnt)
+     register mp_ptr wp;
+     register mp_srcptr up;
+     mp_size_t usize;
+     register unsigned int cnt;
+#endif
+{
+  register mp_limb_t high_limb, low_limb;
+  register unsigned sh_1, sh_2;
+  register mp_size_t i;
+  mp_limb_t retval;
+
+#ifdef DEBUG
+  if (usize == 0 || cnt == 0)
+    abort ();
+#endif
+
+  sh_1 = cnt;
+#if 0
+  if (sh_1 == 0)
+    {
+      if (wp != up)
+	{
+	  /* Copy from high end to low end, to allow specified input/output
+	     overlapping.  */
+	  for (i = usize - 1; i >= 0; i--)
+	    wp[i] = up[i];
+	}
+      return 0;
+    }
+#endif
+
+  wp += 1;
+  sh_2 = BITS_PER_MP_LIMB - sh_1;
+  i = usize - 1;
+  low_limb = up[i];
+  retval = low_limb >> sh_2;
+  high_limb = low_limb;
+  while (--i >= 0)
+    {
+      low_limb = up[i];
+      wp[i] = (high_limb << sh_1) | (low_limb >> sh_2);
+      high_limb = low_limb;
+    }
+  wp[i] = high_limb << sh_1;
+
+  return retval;
+}
diff --git a/libquadmath/printf/mul.c b/libquadmath/printf/mul.c
new file mode 100644
index 000000000..d31fa36fa
--- /dev/null
+++ b/libquadmath/printf/mul.c
@@ -0,0 +1,148 @@
+/* mpn_mul -- Multiply two natural numbers.
+
+Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+/* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
+   and v (pointed to by VP, with VSIZE limbs), and store the result at
+   PRODP.  USIZE + VSIZE limbs are always stored, but if the input
+   operands are normalized.  Return the most significant limb of the
+   result.
+
+   NOTE: The space pointed to by PRODP is overwritten before finished
+   with U and V, so overlap is an error.
+
+   Argument constraints:
+   1. USIZE >= VSIZE.
+   2. PRODP != UP and PRODP != VP, i.e. the destination
+      must be distinct from the multiplier and the multiplicand.  */
+
+/* If KARATSUBA_THRESHOLD is not already defined, define it to a
+   value which is good on most machines.  */
+#ifndef KARATSUBA_THRESHOLD
+#define KARATSUBA_THRESHOLD 32
+#endif
+
+mp_limb_t
+#if __STDC__
+mpn_mul (mp_ptr prodp,
+	 mp_srcptr up, mp_size_t usize,
+	 mp_srcptr vp, mp_size_t vsize)
+#else
+mpn_mul (prodp, up, usize, vp, vsize)
+     mp_ptr prodp;
+     mp_srcptr up;
+     mp_size_t usize;
+     mp_srcptr vp;
+     mp_size_t vsize;
+#endif
+{
+  mp_ptr prod_endp = prodp + usize + vsize - 1;
+  mp_limb_t cy;
+  mp_ptr tspace;
+
+  if (vsize < KARATSUBA_THRESHOLD)
+    {
+      /* Handle simple cases with traditional multiplication.
+
+	 This is the most critical code of the entire function.  All
+	 multiplies rely on this, both small and huge.  Small ones arrive
+	 here immediately.  Huge ones arrive here as this is the base case
+	 for Karatsuba's recursive algorithm below.  */
+      mp_size_t i;
+      mp_limb_t cy_limb;
+      mp_limb_t v_limb;
+
+      if (vsize == 0)
+	return 0;
+
+      /* Multiply by the first limb in V separately, as the result can be
+	 stored (not added) to PROD.  We also avoid a loop for zeroing.  */
+      v_limb = vp[0];
+      if (v_limb <= 1)
+	{
+	  if (v_limb == 1)
+	    MPN_COPY (prodp, up, usize);
+	  else
+	    MPN_ZERO (prodp, usize);
+	  cy_limb = 0;
+	}
+      else
+	cy_limb = mpn_mul_1 (prodp, up, usize, v_limb);
+
+      prodp[usize] = cy_limb;
+      prodp++;
+
+      /* For each iteration in the outer loop, multiply one limb from
+	 U with one limb from V, and add it to PROD.  */
+      for (i = 1; i < vsize; i++)
+	{
+	  v_limb = vp[i];
+	  if (v_limb <= 1)
+	    {
+	      cy_limb = 0;
+	      if (v_limb == 1)
+		cy_limb = mpn_add_n (prodp, prodp, up, usize);
+	    }
+	  else
+	    cy_limb = mpn_addmul_1 (prodp, up, usize, v_limb);
+
+	  prodp[usize] = cy_limb;
+	  prodp++;
+	}
+      return cy_limb;
+    }
+
+  tspace = (mp_ptr) alloca (2 * vsize * BYTES_PER_MP_LIMB);
+  MPN_MUL_N_RECURSE (prodp, up, vp, vsize, tspace);
+
+  prodp += vsize;
+  up += vsize;
+  usize -= vsize;
+  if (usize >= vsize)
+    {
+      mp_ptr tp = (mp_ptr) alloca (2 * vsize * BYTES_PER_MP_LIMB);
+      do
+	{
+	  MPN_MUL_N_RECURSE (tp, up, vp, vsize, tspace);
+	  cy = mpn_add_n (prodp, prodp, tp, vsize);
+	  mpn_add_1 (prodp + vsize, tp + vsize, vsize, cy);
+	  prodp += vsize;
+	  up += vsize;
+	  usize -= vsize;
+	}
+      while (usize >= vsize);
+    }
+
+  /* True: usize < vsize.  */
+
+  /* Make life simple: Recurse.  */
+
+  if (usize != 0)
+    {
+      mpn_mul (tspace, vp, vsize, up, usize);
+      cy = mpn_add_n (prodp, prodp, tspace, vsize);
+      mpn_add_1 (prodp + vsize, tspace + vsize, usize, cy);
+    }
+
+  return *prod_endp;
+}
diff --git a/libquadmath/printf/mul_1.c b/libquadmath/printf/mul_1.c
new file mode 100644
index 000000000..48a273f07
--- /dev/null
+++ b/libquadmath/printf/mul_1.c
@@ -0,0 +1,58 @@
+/* mpn_mul_1 -- Multiply a limb vector with a single limb and
+   store the product in a second limb vector.
+
+Copyright (C) 1991, 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+mp_limb_t
+mpn_mul_1 (res_ptr, s1_ptr, s1_size, s2_limb)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     mp_size_t s1_size;
+     register mp_limb_t s2_limb;
+{
+  register mp_limb_t cy_limb;
+  register mp_size_t j;
+  register mp_limb_t prod_high, prod_low;
+
+  /* The loop counter and index J goes from -S1_SIZE to -1.  This way
+     the loop becomes faster.  */
+  j = -s1_size;
+
+  /* Offset the base pointers to compensate for the negative indices.  */
+  s1_ptr -= j;
+  res_ptr -= j;
+
+  cy_limb = 0;
+  do
+    {
+      umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
+
+      prod_low += cy_limb;
+      cy_limb = (prod_low < cy_limb) + prod_high;
+
+      res_ptr[j] = prod_low;
+    }
+  while (++j != 0);
+
+  return cy_limb;
+}
diff --git a/libquadmath/printf/mul_n.c b/libquadmath/printf/mul_n.c
new file mode 100644
index 000000000..c4bc1bed2
--- /dev/null
+++ b/libquadmath/printf/mul_n.c
@@ -0,0 +1,220 @@
+/* mpn_mul_n -- Multiply two natural numbers of length n.
+
+Copyright (C) 1991, 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+/* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
+   both with SIZE limbs, and store the result at PRODP.  2 * SIZE limbs are
+   always stored.  Return the most significant limb.
+
+   Argument constraints:
+   1. PRODP != UP and PRODP != VP, i.e. the destination
+      must be distinct from the multiplier and the multiplicand.  */
+
+/* If KARATSUBA_THRESHOLD is not already defined, define it to a
+   value which is good on most machines.  */
+#ifndef KARATSUBA_THRESHOLD
+#define KARATSUBA_THRESHOLD 32
+#endif
+
+/* The code can't handle KARATSUBA_THRESHOLD smaller than 2.  */
+#if KARATSUBA_THRESHOLD < 2
+#undef KARATSUBA_THRESHOLD
+#define KARATSUBA_THRESHOLD 2
+#endif
+
+/* Handle simple cases with traditional multiplication.
+
+   This is the most critical code of multiplication.  All multiplies rely
+   on this, both small and huge.  Small ones arrive here immediately.  Huge
+   ones arrive here as this is the base case for Karatsuba's recursive
+   algorithm below.  */
+
+void
+#if __STDC__
+impn_mul_n_basecase (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
+#else
+impn_mul_n_basecase (prodp, up, vp, size)
+     mp_ptr prodp;
+     mp_srcptr up;
+     mp_srcptr vp;
+     mp_size_t size;
+#endif
+{
+  mp_size_t i;
+  mp_limb_t cy_limb;
+  mp_limb_t v_limb;
+
+  /* Multiply by the first limb in V separately, as the result can be
+     stored (not added) to PROD.  We also avoid a loop for zeroing.  */
+  v_limb = vp[0];
+  if (v_limb <= 1)
+    {
+      if (v_limb == 1)
+	MPN_COPY (prodp, up, size);
+      else
+	MPN_ZERO (prodp, size);
+      cy_limb = 0;
+    }
+  else
+    cy_limb = mpn_mul_1 (prodp, up, size, v_limb);
+
+  prodp[size] = cy_limb;
+  prodp++;
+
+  /* For each iteration in the outer loop, multiply one limb from
+     U with one limb from V, and add it to PROD.  */
+  for (i = 1; i < size; i++)
+    {
+      v_limb = vp[i];
+      if (v_limb <= 1)
+	{
+	  cy_limb = 0;
+	  if (v_limb == 1)
+	    cy_limb = mpn_add_n (prodp, prodp, up, size);
+	}
+      else
+	cy_limb = mpn_addmul_1 (prodp, up, size, v_limb);
+
+      prodp[size] = cy_limb;
+      prodp++;
+    }
+}
+
+void
+#if __STDC__
+impn_mul_n (mp_ptr prodp,
+	     mp_srcptr up, mp_srcptr vp, mp_size_t size, mp_ptr tspace)
+#else
+impn_mul_n (prodp, up, vp, size, tspace)
+     mp_ptr prodp;
+     mp_srcptr up;
+     mp_srcptr vp;
+     mp_size_t size;
+     mp_ptr tspace;
+#endif
+{
+  if ((size & 1) != 0)
+    {
+      /* The size is odd, the code code below doesn't handle that.
+	 Multiply the least significant (size - 1) limbs with a recursive
+	 call, and handle the most significant limb of S1 and S2
+	 separately.  */
+      /* A slightly faster way to do this would be to make the Karatsuba
+	 code below behave as if the size were even, and let it check for
+	 odd size in the end.  I.e., in essence move this code to the end.
+	 Doing so would save us a recursive call, and potentially make the
+	 stack grow a lot less.  */
+
+      mp_size_t esize = size - 1;	/* even size */
+      mp_limb_t cy_limb;
+
+      MPN_MUL_N_RECURSE (prodp, up, vp, esize, tspace);
+      cy_limb = mpn_addmul_1 (prodp + esize, up, esize, vp[esize]);
+      prodp[esize + esize] = cy_limb;
+      cy_limb = mpn_addmul_1 (prodp + esize, vp, size, up[esize]);
+
+      prodp[esize + size] = cy_limb;
+    }
+  else
+    {
+      /* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
+
+	 Split U in two pieces, U1 and U0, such that
+	 U = U0 + U1*(B**n),
+	 and V in V1 and V0, such that
+	 V = V0 + V1*(B**n).
+
+	 UV is then computed recursively using the identity
+
+		2n   n          n                     n
+	 UV = (B  + B )U V  +  B (U -U )(V -V )  +  (B + 1)U V
+			1 1        1  0   0  1              0 0
+
+	 Where B = 2**BITS_PER_MP_LIMB.  */
+
+      mp_size_t hsize = size >> 1;
+      mp_limb_t cy;
+      int negflg;
+
+      /*** Product H.	 ________________  ________________
+			|_____U1 x V1____||____U0 x V0_____|  */
+      /* Put result in upper part of PROD and pass low part of TSPACE
+	 as new TSPACE.  */
+      MPN_MUL_N_RECURSE (prodp + size, up + hsize, vp + hsize, hsize, tspace);
+
+      /*** Product M.	 ________________
+			|_(U1-U0)(V0-V1)_|  */
+      if (mpn_cmp (up + hsize, up, hsize) >= 0)
+	{
+	  mpn_sub_n (prodp, up + hsize, up, hsize);
+	  negflg = 0;
+	}
+      else
+	{
+	  mpn_sub_n (prodp, up, up + hsize, hsize);
+	  negflg = 1;
+	}
+      if (mpn_cmp (vp + hsize, vp, hsize) >= 0)
+	{
+	  mpn_sub_n (prodp + hsize, vp + hsize, vp, hsize);
+	  negflg ^= 1;
+	}
+      else
+	{
+	  mpn_sub_n (prodp + hsize, vp, vp + hsize, hsize);
+	  /* No change of NEGFLG.  */
+	}
+      /* Read temporary operands from low part of PROD.
+	 Put result in low part of TSPACE using upper part of TSPACE
+	 as new TSPACE.  */
+      MPN_MUL_N_RECURSE (tspace, prodp, prodp + hsize, hsize, tspace + size);
+
+      /*** Add/copy product H.  */
+      MPN_COPY (prodp + hsize, prodp + size, hsize);
+      cy = mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
+
+      /*** Add product M (if NEGFLG M is a negative number).  */
+      if (negflg)
+	cy -= mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
+      else
+	cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
+
+      /*** Product L.	 ________________  ________________
+			|________________||____U0 x V0_____|  */
+      /* Read temporary operands from low part of PROD.
+	 Put result in low part of TSPACE using upper part of TSPACE
+	 as new TSPACE.  */
+      MPN_MUL_N_RECURSE (tspace, up, vp, hsize, tspace + size);
+
+      /*** Add/copy Product L (twice).  */
+
+      cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
+      if (cy)
+	mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
+
+      MPN_COPY (prodp, tspace, hsize);
+      cy = mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
+      if (cy)
+	mpn_add_1 (prodp + size, prodp + size, size, 1);
+    }
+}
diff --git a/libquadmath/printf/printf_fp.c b/libquadmath/printf/printf_fp.c
new file mode 100644
index 000000000..eb663726d
--- /dev/null
+++ b/libquadmath/printf/printf_fp.c
@@ -0,0 +1,1306 @@
+/* Floating point output for `printf'.
+   Copyright (C) 1995-2003, 2006-2008, 2011 Free Software Foundation, Inc.
+
+   This file is part of the GNU C Library.
+   Written by Ulrich Drepper <drepper@gnu.ai.mit.edu>, 1995.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <config.h>
+#include <float.h>
+#include <math.h>
+#include <string.h>
+#include <unistd.h>
+#include <stdlib.h>
+#define NDEBUG
+#include <assert.h>
+#ifdef HAVE_ERRNO_H
+#include <errno.h>
+#endif
+#include <stdio.h>
+#include <stdarg.h>
+#include "quadmath-printf.h"
+#include "fpioconst.h"
+
+#ifdef USE_I18N_NUMBER_H
+#include "_i18n_number.h"
+#endif
+
+
+/* Macros for doing the actual output.  */
+
+#define outchar(ch)							      \
+  do									      \
+    {									      \
+      register const int outc = (ch);					      \
+      if (PUTC (outc, fp) == EOF)					      \
+	{								      \
+	  if (buffer_malloced)						      \
+	    free (wbuffer);						      \
+	  return -1;							      \
+	}								      \
+      ++done;								      \
+    } while (0)
+
+#define PRINT(ptr, wptr, len)						      \
+  do									      \
+    {									      \
+      register size_t outlen = (len);					      \
+      if (len > 20)							      \
+	{								      \
+	  if (PUT (fp, wide ? (const char *) wptr : ptr, outlen) != outlen)   \
+	    {								      \
+	      if (buffer_malloced)					      \
+		free (wbuffer);						      \
+	      return -1;						      \
+	    }								      \
+	  ptr += outlen;						      \
+	  done += outlen;						      \
+	}								      \
+      else								      \
+	{								      \
+	  if (wide)							      \
+	    while (outlen-- > 0)					      \
+	      outchar (*wptr++);					      \
+	  else								      \
+	    while (outlen-- > 0)					      \
+	      outchar (*ptr++);						      \
+	}								      \
+    } while (0)
+
+#define PADN(ch, len)							      \
+  do									      \
+    {									      \
+      if (PAD (fp, ch, len) != len)					      \
+	{								      \
+	  if (buffer_malloced)						      \
+	    free (wbuffer);						      \
+	  return -1;							      \
+	}								      \
+      done += len;							      \
+    }									      \
+  while (0)
+
+
+/* We use the GNU MP library to handle large numbers.
+
+   An MP variable occupies a varying number of entries in its array.  We keep
+   track of this number for efficiency reasons.  Otherwise we would always
+   have to process the whole array.  */
+#define MPN_VAR(name) mp_limb_t *name; mp_size_t name##size
+
+#define MPN_ASSIGN(dst,src)						      \
+  memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
+#define MPN_GE(u,v) \
+  (u##size > v##size || (u##size == v##size && mpn_cmp (u, v, u##size) >= 0))
+
+extern mp_size_t mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size,
+				     int *expt, int *is_neg,
+				     __float128 value) attribute_hidden;
+static unsigned int guess_grouping (unsigned int intdig_max,
+				    const char *grouping);
+
+
+static wchar_t *group_number (wchar_t *buf, wchar_t *bufend,
+			      unsigned int intdig_no, const char *grouping,
+			      wchar_t thousands_sep, int ngroups);
+
+
+int
+__quadmath_printf_fp (struct __quadmath_printf_file *fp,
+		      const struct printf_info *info,
+		      const void *const *args)
+{
+  /* The floating-point value to output.  */
+  __float128 fpnum;
+
+  /* Locale-dependent representation of decimal point.	*/
+  const char *decimal;
+  wchar_t decimalwc;
+
+  /* Locale-dependent thousands separator and grouping specification.  */
+  const char *thousands_sep = NULL;
+  wchar_t thousands_sepwc = L_('\0');
+  const char *grouping;
+
+  /* "NaN" or "Inf" for the special cases.  */
+  const char *special = NULL;
+  const wchar_t *wspecial = NULL;
+
+  /* We need just a few limbs for the input before shifting to the right
+     position.	*/
+  mp_limb_t fp_input[(FLT128_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB];
+  /* We need to shift the contents of fp_input by this amount of bits.	*/
+  int to_shift = 0;
+
+  /* The fraction of the floting-point value in question  */
+  MPN_VAR(frac);
+  /* and the exponent.	*/
+  int exponent;
+  /* Sign of the exponent.  */
+  int expsign = 0;
+  /* Sign of float number.  */
+  int is_neg = 0;
+
+  /* Scaling factor.  */
+  MPN_VAR(scale);
+
+  /* Temporary bignum value.  */
+  MPN_VAR(tmp);
+
+  /* Digit which is result of last hack_digit() call.  */
+  wchar_t digit;
+
+  /* The type of output format that will be used: 'e'/'E' or 'f'.  */
+  int type;
+
+  /* Counter for number of written characters.	*/
+  int done = 0;
+
+  /* General helper (carry limb).  */
+  mp_limb_t cy;
+
+  /* Nonzero if this is output on a wide character stream.  */
+  int wide = info->wide;
+
+  /* Buffer in which we produce the output.  */
+  wchar_t *wbuffer = NULL;
+  /* Flag whether wbuffer is malloc'ed or not.  */
+  int buffer_malloced = 0;
+
+  auto wchar_t hack_digit (void);
+
+  wchar_t hack_digit (void)
+    {
+      mp_limb_t hi;
+
+      if (expsign != 0 && type == 'f' && exponent-- > 0)
+	hi = 0;
+      else if (scalesize == 0)
+	{
+	  hi = frac[fracsize - 1];
+	  frac[fracsize - 1] = mpn_mul_1 (frac, frac, fracsize - 1, 10);
+	}
+      else
+	{
+	  if (fracsize < scalesize)
+	    hi = 0;
+	  else
+	    {
+	      hi = mpn_divmod (tmp, frac, fracsize, scale, scalesize);
+	      tmp[fracsize - scalesize] = hi;
+	      hi = tmp[0];
+
+	      fracsize = scalesize;
+	      while (fracsize != 0 && frac[fracsize - 1] == 0)
+		--fracsize;
+	      if (fracsize == 0)
+		{
+		  /* We're not prepared for an mpn variable with zero
+		     limbs.  */
+		  fracsize = 1;
+		  return L_('0') + hi;
+		}
+	    }
+
+	  mp_limb_t _cy = mpn_mul_1 (frac, frac, fracsize, 10);
+	  if (_cy != 0)
+	    frac[fracsize++] = _cy;
+	}
+
+      return L_('0') + hi;
+    }
+
+  /* Figure out the decimal point character.  */
+#ifdef USE_NL_LANGINFO
+  if (info->extra == 0)
+    decimal = nl_langinfo (DECIMAL_POINT);
+  else
+    {
+      decimal = nl_langinfo (MON_DECIMAL_POINT);
+      if (*decimal == '\0')
+	decimal = nl_langinfo (DECIMAL_POINT);
+    }
+  /* The decimal point character must never be zero.  */
+  assert (*decimal != '\0');
+#elif defined USE_LOCALECONV
+  const struct lconv *lc = localeconv ();
+  if (info->extra == 0)
+    decimal = lc->decimal_point;
+  else
+    {
+      decimal = lc->mon_decimal_point;
+      if (decimal == NULL || *decimal == '\0')
+	decimal = lc->decimal_point;
+    }
+  if (decimal == NULL || *decimal == '\0')
+    decimal = ".";
+#else
+  decimal = ".";
+#endif
+#ifdef USE_NL_LANGINFO_WC
+  if (info->extra == 0)
+    decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC);
+  else
+    {
+      decimalwc = nl_langinfo_wc (_NL_MONETARY_DECIMAL_POINT_WC);
+      if (decimalwc == L_('\0'))
+	decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC);
+    }
+  /* The decimal point character must never be zero.  */
+  assert (decimalwc != L_('\0'));
+#else
+  decimalwc = L_('.');
+#endif
+
+#if defined USE_NL_LANGINFO && defined USE_NL_LANGINFO_WC
+  if (info->group)
+    {
+      if (info->extra == 0)
+	grouping = nl_langinfo (GROUPING);
+      else
+	grouping = nl_langinfo (MON_GROUPING);
+
+      if (*grouping <= 0 || *grouping == CHAR_MAX)
+	grouping = NULL;
+      else
+	{
+	  /* Figure out the thousands separator character.  */
+	  if (wide)
+	    {
+	      if (info->extra == 0)
+		thousands_sepwc = nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC);
+	      else
+		thousands_sepwc = nl_langinfo_wc (_NL_MONETARY_THOUSANDS_SEP_WC);
+
+	      if (thousands_sepwc == L_('\0'))
+		grouping = NULL;
+	    }
+	  else
+	    {
+	      if (info->extra == 0)
+		thousands_sep = nl_langinfo (THOUSANDS_SEP);
+	      else
+		thousands_sep = nl_langinfo (MON_THOUSANDS_SEP);
+	      if (*thousands_sep == '\0')
+		grouping = NULL;
+	    }
+	}
+    }
+  else
+#elif defined USE_NL_LANGINFO
+  if (info->group && !wide)
+    {
+      if (info->extra == 0)
+	grouping = nl_langinfo (GROUPING);
+      else
+	grouping = nl_langinfo (MON_GROUPING);
+
+      if (*grouping <= 0 || *grouping == CHAR_MAX)
+	grouping = NULL;
+      else
+	{
+	  /* Figure out the thousands separator character.  */
+	  if (info->extra == 0)
+	    thousands_sep = nl_langinfo (THOUSANDS_SEP);
+	  else
+	    thousands_sep = nl_langinfo (MON_THOUSANDS_SEP);
+
+	  if (*thousands_sep == '\0')
+	    grouping = NULL;
+	}
+    }
+  else
+#elif defined USE_LOCALECONV
+  if (info->group && !wide)
+    {
+      if (info->extra == 0)
+	grouping = lc->grouping;
+      else
+	grouping = lc->mon_grouping;
+
+      if (grouping == NULL || *grouping <= 0 || *grouping == CHAR_MAX)
+	grouping = NULL;
+      else
+	{
+	  /* Figure out the thousands separator character.  */
+	  if (info->extra == 0)
+	    thousands_sep = lc->thousands_sep;
+	  else
+	    thousands_sep = lc->mon_thousands_sep;
+
+	  if (thousands_sep == NULL || *thousands_sep == '\0')
+	    grouping = NULL;
+	}
+    }
+  else
+#endif
+    grouping = NULL;
+  if (grouping != NULL && !wide)
+    /* If we are printing multibyte characters and there is a
+       multibyte representation for the thousands separator,
+       we must ensure the wide character thousands separator
+       is available, even if it is fake.  */
+    thousands_sepwc = (wchar_t) 0xfffffffe;
+
+  /* Fetch the argument value.	*/
+    {
+      fpnum = **(const __float128 **) args[0];
+
+      /* Check for special values: not a number or infinity.  */
+      if (isnanq (fpnum))
+	{
+	  ieee854_float128 u = { .value = fpnum };
+	  is_neg = u.ieee.negative != 0;
+	  if (isupper (info->spec))
+	    {
+	      special = "NAN";
+	      wspecial = L_("NAN");
+	    }
+	    else
+	      {
+		special = "nan";
+		wspecial = L_("nan");
+	      }
+	}
+      else if (isinfq (fpnum))
+	{
+	  is_neg = fpnum < 0;
+	  if (isupper (info->spec))
+	    {
+	      special = "INF";
+	      wspecial = L_("INF");
+	    }
+	  else
+	    {
+	      special = "inf";
+	      wspecial = L_("inf");
+	    }
+	}
+      else
+	{
+	  fracsize = mpn_extract_flt128 (fp_input,
+					 (sizeof (fp_input) /
+					  sizeof (fp_input[0])),
+					 &exponent, &is_neg, fpnum);
+	  to_shift = 1 + fracsize * BITS_PER_MP_LIMB - FLT128_MANT_DIG;
+	}
+    }
+
+  if (special)
+    {
+      int width = info->width;
+
+      if (is_neg || info->showsign || info->space)
+	--width;
+      width -= 3;
+
+      if (!info->left && width > 0)
+	PADN (' ', width);
+
+      if (is_neg)
+	outchar ('-');
+      else if (info->showsign)
+	outchar ('+');
+      else if (info->space)
+	outchar (' ');
+
+      PRINT (special, wspecial, 3);
+
+      if (info->left && width > 0)
+	PADN (' ', width);
+
+      return done;
+    }
+
+
+  /* We need three multiprecision variables.  Now that we have the exponent
+     of the number we can allocate the needed memory.  It would be more
+     efficient to use variables of the fixed maximum size but because this
+     would be really big it could lead to memory problems.  */
+  {
+    mp_size_t bignum_size = ((ABS (exponent) + BITS_PER_MP_LIMB - 1)
+			     / BITS_PER_MP_LIMB
+			     + (FLT128_MANT_DIG / BITS_PER_MP_LIMB > 2 ? 8 : 4))
+			    * sizeof (mp_limb_t);
+    frac = (mp_limb_t *) alloca (bignum_size);
+    tmp = (mp_limb_t *) alloca (bignum_size);
+    scale = (mp_limb_t *) alloca (bignum_size);
+  }
+
+  /* We now have to distinguish between numbers with positive and negative
+     exponents because the method used for the one is not applicable/efficient
+     for the other.  */
+  scalesize = 0;
+  if (exponent > 2)
+    {
+      /* |FP| >= 8.0.  */
+      int scaleexpo = 0;
+      int explog = FLT128_MAX_10_EXP_LOG;
+      int exp10 = 0;
+      const struct mp_power *powers = &_fpioconst_pow10[explog + 1];
+      int cnt_h, cnt_l, i;
+
+      if ((exponent + to_shift) % BITS_PER_MP_LIMB == 0)
+	{
+	  MPN_COPY_DECR (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
+			 fp_input, fracsize);
+	  fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
+	}
+      else
+	{
+	  cy = mpn_lshift (frac + (exponent + to_shift) / BITS_PER_MP_LIMB,
+			   fp_input, fracsize,
+			   (exponent + to_shift) % BITS_PER_MP_LIMB);
+	  fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB;
+	  if (cy)
+	    frac[fracsize++] = cy;
+	}
+      MPN_ZERO (frac, (exponent + to_shift) / BITS_PER_MP_LIMB);
+
+      assert (powers > &_fpioconst_pow10[0]);
+      do
+	{
+	  --powers;
+
+	  /* The number of the product of two binary numbers with n and m
+	     bits respectively has m+n or m+n-1 bits.	*/
+	  if (exponent >= scaleexpo + powers->p_expo - 1)
+	    {
+	      if (scalesize == 0)
+		{
+		  if (FLT128_MANT_DIG > _FPIO_CONST_OFFSET * BITS_PER_MP_LIMB)
+		    {
+#define _FPIO_CONST_SHIFT \
+  (((FLT128_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB) \
+   - _FPIO_CONST_OFFSET)
+		      /* 64bit const offset is not enough for
+			 IEEE quad long double.  */
+		      tmpsize = powers->arraysize + _FPIO_CONST_SHIFT;
+		      memcpy (tmp + _FPIO_CONST_SHIFT,
+			      &__tens[powers->arrayoff],
+			      tmpsize * sizeof (mp_limb_t));
+		      MPN_ZERO (tmp, _FPIO_CONST_SHIFT);
+		      /* Adjust exponent, as scaleexpo will be this much
+			 bigger too.  */
+		      exponent += _FPIO_CONST_SHIFT * BITS_PER_MP_LIMB;
+		    }
+		  else
+		    {
+		      tmpsize = powers->arraysize;
+		      memcpy (tmp, &__tens[powers->arrayoff],
+			      tmpsize * sizeof (mp_limb_t));
+		    }
+		}
+	      else
+		{
+		  cy = mpn_mul (tmp, scale, scalesize,
+				&__tens[powers->arrayoff
+					+ _FPIO_CONST_OFFSET],
+				powers->arraysize - _FPIO_CONST_OFFSET);
+		  tmpsize = scalesize + powers->arraysize - _FPIO_CONST_OFFSET;
+		  if (cy == 0)
+		    --tmpsize;
+		}
+
+	      if (MPN_GE (frac, tmp))
+		{
+		  int cnt;
+		  MPN_ASSIGN (scale, tmp);
+		  count_leading_zeros (cnt, scale[scalesize - 1]);
+		  scaleexpo = (scalesize - 2) * BITS_PER_MP_LIMB - cnt - 1;
+		  exp10 |= 1 << explog;
+		}
+	    }
+	  --explog;
+	}
+      while (powers > &_fpioconst_pow10[0]);
+      exponent = exp10;
+
+      /* Optimize number representations.  We want to represent the numbers
+	 with the lowest number of bytes possible without losing any
+	 bytes. Also the highest bit in the scaling factor has to be set
+	 (this is a requirement of the MPN division routines).  */
+      if (scalesize > 0)
+	{
+	  /* Determine minimum number of zero bits at the end of
+	     both numbers.  */
+	  for (i = 0; scale[i] == 0 && frac[i] == 0; i++)
+	    ;
+
+	  /* Determine number of bits the scaling factor is misplaced.	*/
+	  count_leading_zeros (cnt_h, scale[scalesize - 1]);
+
+	  if (cnt_h == 0)
+	    {
+	      /* The highest bit of the scaling factor is already set.	So
+		 we only have to remove the trailing empty limbs.  */
+	      if (i > 0)
+		{
+		  MPN_COPY_INCR (scale, scale + i, scalesize - i);
+		  scalesize -= i;
+		  MPN_COPY_INCR (frac, frac + i, fracsize - i);
+		  fracsize -= i;
+		}
+	    }
+	  else
+	    {
+	      if (scale[i] != 0)
+		{
+		  count_trailing_zeros (cnt_l, scale[i]);
+		  if (frac[i] != 0)
+		    {
+		      int cnt_l2;
+		      count_trailing_zeros (cnt_l2, frac[i]);
+		      if (cnt_l2 < cnt_l)
+			cnt_l = cnt_l2;
+		    }
+		}
+	      else
+		count_trailing_zeros (cnt_l, frac[i]);
+
+	      /* Now shift the numbers to their optimal position.  */
+	      if (i == 0 && BITS_PER_MP_LIMB - cnt_h > cnt_l)
+		{
+		  /* We cannot save any memory.	 So just roll both numbers
+		     so that the scaling factor has its highest bit set.  */
+
+		  (void) mpn_lshift (scale, scale, scalesize, cnt_h);
+		  cy = mpn_lshift (frac, frac, fracsize, cnt_h);
+		  if (cy != 0)
+		    frac[fracsize++] = cy;
+		}
+	      else if (BITS_PER_MP_LIMB - cnt_h <= cnt_l)
+		{
+		  /* We can save memory by removing the trailing zero limbs
+		     and by packing the non-zero limbs which gain another
+		     free one. */
+
+		  (void) mpn_rshift (scale, scale + i, scalesize - i,
+				     BITS_PER_MP_LIMB - cnt_h);
+		  scalesize -= i + 1;
+		  (void) mpn_rshift (frac, frac + i, fracsize - i,
+				     BITS_PER_MP_LIMB - cnt_h);
+		  fracsize -= frac[fracsize - i - 1] == 0 ? i + 1 : i;
+		}
+	      else
+		{
+		  /* We can only save the memory of the limbs which are zero.
+		     The non-zero parts occupy the same number of limbs.  */
+
+		  (void) mpn_rshift (scale, scale + (i - 1),
+				     scalesize - (i - 1),
+				     BITS_PER_MP_LIMB - cnt_h);
+		  scalesize -= i;
+		  (void) mpn_rshift (frac, frac + (i - 1),
+				     fracsize - (i - 1),
+				     BITS_PER_MP_LIMB - cnt_h);
+		  fracsize -= frac[fracsize - (i - 1) - 1] == 0 ? i : i - 1;
+		}
+	    }
+	}
+    }
+  else if (exponent < 0)
+    {
+      /* |FP| < 1.0.  */
+      int exp10 = 0;
+      int explog = FLT128_MAX_10_EXP_LOG;
+      const struct mp_power *powers = &_fpioconst_pow10[explog + 1];
+
+      /* Now shift the input value to its right place.	*/
+      cy = mpn_lshift (frac, fp_input, fracsize, to_shift);
+      frac[fracsize++] = cy;
+      assert (cy == 1 || (frac[fracsize - 2] == 0 && frac[0] == 0));
+
+      expsign = 1;
+      exponent = -exponent;
+
+      assert (powers != &_fpioconst_pow10[0]);
+      do
+	{
+	  --powers;
+
+	  if (exponent >= powers->m_expo)
+	    {
+	      int i, incr, cnt_h, cnt_l;
+	      mp_limb_t topval[2];
+
+	      /* The mpn_mul function expects the first argument to be
+		 bigger than the second.  */
+	      if (fracsize < powers->arraysize - _FPIO_CONST_OFFSET)
+		cy = mpn_mul (tmp, &__tens[powers->arrayoff
+					   + _FPIO_CONST_OFFSET],
+			      powers->arraysize - _FPIO_CONST_OFFSET,
+			      frac, fracsize);
+	      else
+		cy = mpn_mul (tmp, frac, fracsize,
+			      &__tens[powers->arrayoff + _FPIO_CONST_OFFSET],
+			      powers->arraysize - _FPIO_CONST_OFFSET);
+	      tmpsize = fracsize + powers->arraysize - _FPIO_CONST_OFFSET;
+	      if (cy == 0)
+		--tmpsize;
+
+	      count_leading_zeros (cnt_h, tmp[tmpsize - 1]);
+	      incr = (tmpsize - fracsize) * BITS_PER_MP_LIMB
+		     + BITS_PER_MP_LIMB - 1 - cnt_h;
+
+	      assert (incr <= powers->p_expo);
+
+	      /* If we increased the exponent by exactly 3 we have to test
+		 for overflow.	This is done by comparing with 10 shifted
+		 to the right position.	 */
+	      if (incr == exponent + 3)
+		{
+		  if (cnt_h <= BITS_PER_MP_LIMB - 4)
+		    {
+		      topval[0] = 0;
+		      topval[1]
+			= ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4 - cnt_h);
+		    }
+		  else
+		    {
+		      topval[0] = ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4);
+		      topval[1] = 0;
+		      (void) mpn_lshift (topval, topval, 2,
+					 BITS_PER_MP_LIMB - cnt_h);
+		    }
+		}
+
+	      /* We have to be careful when multiplying the last factor.
+		 If the result is greater than 1.0 be have to test it
+		 against 10.0.  If it is greater or equal to 10.0 the
+		 multiplication was not valid.  This is because we cannot
+		 determine the number of bits in the result in advance.  */
+	      if (incr < exponent + 3
+		  || (incr == exponent + 3 &&
+		      (tmp[tmpsize - 1] < topval[1]
+		       || (tmp[tmpsize - 1] == topval[1]
+			   && tmp[tmpsize - 2] < topval[0]))))
+		{
+		  /* The factor is right.  Adapt binary and decimal
+		     exponents.	 */
+		  exponent -= incr;
+		  exp10 |= 1 << explog;
+
+		  /* If this factor yields a number greater or equal to
+		     1.0, we must not shift the non-fractional digits down. */
+		  if (exponent < 0)
+		    cnt_h += -exponent;
+
+		  /* Now we optimize the number representation.	 */
+		  for (i = 0; tmp[i] == 0; ++i);
+		  if (cnt_h == BITS_PER_MP_LIMB - 1)
+		    {
+		      MPN_COPY (frac, tmp + i, tmpsize - i);
+		      fracsize = tmpsize - i;
+		    }
+		  else
+		    {
+		      count_trailing_zeros (cnt_l, tmp[i]);
+
+		      /* Now shift the numbers to their optimal position.  */
+		      if (i == 0 && BITS_PER_MP_LIMB - 1 - cnt_h > cnt_l)
+			{
+			  /* We cannot save any memory.	 Just roll the
+			     number so that the leading digit is in a
+			     separate limb.  */
+
+			  cy = mpn_lshift (frac, tmp, tmpsize, cnt_h + 1);
+			  fracsize = tmpsize + 1;
+			  frac[fracsize - 1] = cy;
+			}
+		      else if (BITS_PER_MP_LIMB - 1 - cnt_h <= cnt_l)
+			{
+			  (void) mpn_rshift (frac, tmp + i, tmpsize - i,
+					     BITS_PER_MP_LIMB - 1 - cnt_h);
+			  fracsize = tmpsize - i;
+			}
+		      else
+			{
+			  /* We can only save the memory of the limbs which
+			     are zero.	The non-zero parts occupy the same
+			     number of limbs.  */
+
+			  (void) mpn_rshift (frac, tmp + (i - 1),
+					     tmpsize - (i - 1),
+					     BITS_PER_MP_LIMB - 1 - cnt_h);
+			  fracsize = tmpsize - (i - 1);
+			}
+		    }
+		}
+	    }
+	  --explog;
+	}
+      while (powers != &_fpioconst_pow10[1] && exponent > 0);
+      /* All factors but 10^-1 are tested now.	*/
+      if (exponent > 0)
+	{
+	  int cnt_l;
+
+	  cy = mpn_mul_1 (tmp, frac, fracsize, 10);
+	  tmpsize = fracsize;
+	  assert (cy == 0 || tmp[tmpsize - 1] < 20);
+
+	  count_trailing_zeros (cnt_l, tmp[0]);
+	  if (cnt_l < MIN (4, exponent))
+	    {
+	      cy = mpn_lshift (frac, tmp, tmpsize,
+			       BITS_PER_MP_LIMB - MIN (4, exponent));
+	      if (cy != 0)
+		frac[tmpsize++] = cy;
+	    }
+	  else
+	    (void) mpn_rshift (frac, tmp, tmpsize, MIN (4, exponent));
+	  fracsize = tmpsize;
+	  exp10 |= 1;
+	  assert (frac[fracsize - 1] < 10);
+	}
+      exponent = exp10;
+    }
+  else
+    {
+      /* This is a special case.  We don't need a factor because the
+	 numbers are in the range of 1.0 <= |fp| < 8.0.  We simply
+	 shift it to the right place and divide it by 1.0 to get the
+	 leading digit.	 (Of course this division is not really made.)	*/
+      assert (0 <= exponent && exponent < 3 &&
+	      exponent + to_shift < BITS_PER_MP_LIMB);
+
+      /* Now shift the input value to its right place.	*/
+      cy = mpn_lshift (frac, fp_input, fracsize, (exponent + to_shift));
+      frac[fracsize++] = cy;
+      exponent = 0;
+    }
+
+  {
+    int width = info->width;
+    wchar_t *wstartp, *wcp;
+    size_t chars_needed;
+    int expscale;
+    int intdig_max, intdig_no = 0;
+    int fracdig_min;
+    int fracdig_max;
+    int dig_max;
+    int significant;
+    int ngroups = 0;
+    char spec = tolower (info->spec);
+
+    if (spec == 'e')
+      {
+	type = info->spec;
+	intdig_max = 1;
+	fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
+	chars_needed = 1 + 1 + (size_t) fracdig_max + 1 + 1 + 4;
+	/*	       d   .	 ddd	     e	 +-  ddd  */
+	dig_max = __INT_MAX__;		/* Unlimited.  */
+	significant = 1;		/* Does not matter here.  */
+      }
+    else if (spec == 'f')
+      {
+	type = 'f';
+	fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec;
+	dig_max = __INT_MAX__;		/* Unlimited.  */
+	significant = 1;		/* Does not matter here.  */
+	if (expsign == 0)
+	  {
+	    intdig_max = exponent + 1;
+	    /* This can be really big!	*/  /* XXX Maybe malloc if too big? */
+	    chars_needed = (size_t) exponent + 1 + 1 + (size_t) fracdig_max;
+	  }
+	else
+	  {
+	    intdig_max = 1;
+	    chars_needed = 1 + 1 + (size_t) fracdig_max;
+	  }
+      }
+    else
+      {
+	dig_max = info->prec < 0 ? 6 : (info->prec == 0 ? 1 : info->prec);
+	if ((expsign == 0 && exponent >= dig_max)
+	    || (expsign != 0 && exponent > 4))
+	  {
+	    if ('g' - 'G' == 'e' - 'E')
+	      type = 'E' + (info->spec - 'G');
+	    else
+	      type = isupper (info->spec) ? 'E' : 'e';
+	    fracdig_max = dig_max - 1;
+	    intdig_max = 1;
+	    chars_needed = 1 + 1 + (size_t) fracdig_max + 1 + 1 + 4;
+	  }
+	else
+	  {
+	    type = 'f';
+	    intdig_max = expsign == 0 ? exponent + 1 : 0;
+	    fracdig_max = dig_max - intdig_max;
+	    /* We need space for the significant digits and perhaps
+	       for leading zeros when < 1.0.  The number of leading
+	       zeros can be as many as would be required for
+	       exponential notation with a negative two-digit
+	       exponent, which is 4.  */
+	    chars_needed = (size_t) dig_max + 1 + 4;
+	  }
+	fracdig_min = info->alt ? fracdig_max : 0;
+	significant = 0;		/* We count significant digits.	 */
+      }
+
+    if (grouping)
+      {
+	/* Guess the number of groups we will make, and thus how
+	   many spaces we need for separator characters.  */
+	ngroups = guess_grouping (intdig_max, grouping);
+	/* Allocate one more character in case rounding increases the
+	   number of groups.  */
+	chars_needed += ngroups + 1;
+      }
+
+    /* Allocate buffer for output.  We need two more because while rounding
+       it is possible that we need two more characters in front of all the
+       other output.  If the amount of memory we have to allocate is too
+       large use `malloc' instead of `alloca'.  */
+    if (__builtin_expect (chars_needed >= (size_t) -1 / sizeof (wchar_t) - 2
+			  || chars_needed < fracdig_max, 0))
+      {
+	/* Some overflow occurred.  */
+#if defined HAVE_ERRNO_H && defined ERANGE
+	errno = ERANGE;
+#endif
+	return -1;
+      }
+    size_t wbuffer_to_alloc = (2 + chars_needed) * sizeof (wchar_t);
+    buffer_malloced = wbuffer_to_alloc >= 4096;
+    if (__builtin_expect (buffer_malloced, 0))
+      {
+	wbuffer = (wchar_t *) malloc (wbuffer_to_alloc);
+	if (wbuffer == NULL)
+	  /* Signal an error to the caller.  */
+	  return -1;
+      }
+    else
+      wbuffer = (wchar_t *) alloca (wbuffer_to_alloc);
+    wcp = wstartp = wbuffer + 2;	/* Let room for rounding.  */
+
+    /* Do the real work: put digits in allocated buffer.  */
+    if (expsign == 0 || type != 'f')
+      {
+	assert (expsign == 0 || intdig_max == 1);
+	while (intdig_no < intdig_max)
+	  {
+	    ++intdig_no;
+	    *wcp++ = hack_digit ();
+	  }
+	significant = 1;
+	if (info->alt
+	    || fracdig_min > 0
+	    || (fracdig_max > 0 && (fracsize > 1 || frac[0] != 0)))
+	  *wcp++ = decimalwc;
+      }
+    else
+      {
+	/* |fp| < 1.0 and the selected type is 'f', so put "0."
+	   in the buffer.  */
+	*wcp++ = L_('0');
+	--exponent;
+	*wcp++ = decimalwc;
+      }
+
+    /* Generate the needed number of fractional digits.	 */
+    int fracdig_no = 0;
+    int added_zeros = 0;
+    while (fracdig_no < fracdig_min + added_zeros
+	   || (fracdig_no < fracdig_max && (fracsize > 1 || frac[0] != 0)))
+      {
+	++fracdig_no;
+	*wcp = hack_digit ();
+	if (*wcp++ != L_('0'))
+	  significant = 1;
+	else if (significant == 0)
+	  {
+	    ++fracdig_max;
+	    if (fracdig_min > 0)
+	      ++added_zeros;
+	  }
+      }
+
+    /* Do rounding.  */
+    digit = hack_digit ();
+    if (digit > L_('4'))
+      {
+	wchar_t *wtp = wcp;
+
+	if (digit == L_('5')
+	    && ((*(wcp - 1) != decimalwc && (*(wcp - 1) & 1) == 0)
+		|| ((*(wcp - 1) == decimalwc && (*(wcp - 2) & 1) == 0))))
+	  {
+	    /* This is the critical case.	 */
+	    if (fracsize == 1 && frac[0] == 0)
+	      /* Rest of the number is zero -> round to even.
+		 (IEEE 754-1985 4.1 says this is the default rounding.)  */
+	      goto do_expo;
+	    else if (scalesize == 0)
+	      {
+		/* Here we have to see whether all limbs are zero since no
+		   normalization happened.  */
+		size_t lcnt = fracsize;
+		while (lcnt >= 1 && frac[lcnt - 1] == 0)
+		  --lcnt;
+		if (lcnt == 0)
+		  /* Rest of the number is zero -> round to even.
+		     (IEEE 754-1985 4.1 says this is the default rounding.)  */
+		  goto do_expo;
+	      }
+	  }
+
+	if (fracdig_no > 0)
+	  {
+	    /* Process fractional digits.  Terminate if not rounded or
+	       radix character is reached.  */
+	    int removed = 0;
+	    while (*--wtp != decimalwc && *wtp == L_('9'))
+	      {
+		*wtp = L_('0');
+		++removed;
+	      }
+	    if (removed == fracdig_min && added_zeros > 0)
+	      --added_zeros;
+	    if (*wtp != decimalwc)
+	      /* Round up.  */
+	      (*wtp)++;
+	    else if (__builtin_expect (spec == 'g' && type == 'f' && info->alt
+				       && wtp == wstartp + 1
+				       && wstartp[0] == L_('0'),
+				       0))
+	      /* This is a special case: the rounded number is 1.0,
+		 the format is 'g' or 'G', and the alternative format
+		 is selected.  This means the result must be "1.".  */
+	      --added_zeros;
+	  }
+
+	if (fracdig_no == 0 || *wtp == decimalwc)
+	  {
+	    /* Round the integer digits.  */
+	    if (*(wtp - 1) == decimalwc)
+	      --wtp;
+
+	    while (--wtp >= wstartp && *wtp == L_('9'))
+	      *wtp = L_('0');
+
+	    if (wtp >= wstartp)
+	      /* Round up.  */
+	      (*wtp)++;
+	    else
+	      /* It is more critical.  All digits were 9's.  */
+	      {
+		if (type != 'f')
+		  {
+		    *wstartp = '1';
+		    exponent += expsign == 0 ? 1 : -1;
+
+		    /* The above exponent adjustment could lead to 1.0e-00,
+		       e.g. for 0.999999999.  Make sure exponent 0 always
+		       uses + sign.  */
+		    if (exponent == 0)
+		      expsign = 0;
+		  }
+		else if (intdig_no == dig_max)
+		  {
+		    /* This is the case where for type %g the number fits
+		       really in the range for %f output but after rounding
+		       the number of digits is too big.	 */
+		    *--wstartp = decimalwc;
+		    *--wstartp = L_('1');
+
+		    if (info->alt || fracdig_no > 0)
+		      {
+			/* Overwrite the old radix character.  */
+			wstartp[intdig_no + 2] = L_('0');
+			++fracdig_no;
+		      }
+
+		    fracdig_no += intdig_no;
+		    intdig_no = 1;
+		    fracdig_max = intdig_max - intdig_no;
+		    ++exponent;
+		    /* Now we must print the exponent.	*/
+		    type = isupper (info->spec) ? 'E' : 'e';
+		  }
+		else
+		  {
+		    /* We can simply add another another digit before the
+		       radix.  */
+		    *--wstartp = L_('1');
+		    ++intdig_no;
+		  }
+
+		/* While rounding the number of digits can change.
+		   If the number now exceeds the limits remove some
+		   fractional digits.  */
+		if (intdig_no + fracdig_no > dig_max)
+		  {
+		    wcp -= intdig_no + fracdig_no - dig_max;
+		    fracdig_no -= intdig_no + fracdig_no - dig_max;
+		  }
+	      }
+	  }
+      }
+
+  do_expo:
+    /* Now remove unnecessary '0' at the end of the string.  */
+    while (fracdig_no > fracdig_min + added_zeros && *(wcp - 1) == L_('0'))
+      {
+	--wcp;
+	--fracdig_no;
+      }
+    /* If we eliminate all fractional digits we perhaps also can remove
+       the radix character.  */
+    if (fracdig_no == 0 && !info->alt && *(wcp - 1) == decimalwc)
+      --wcp;
+
+    if (grouping)
+      {
+	/* Rounding might have changed the number of groups.  We allocated
+	   enough memory but we need here the correct number of groups.  */
+	if (intdig_no != intdig_max)
+	  ngroups = guess_grouping (intdig_no, grouping);
+
+	/* Add in separator characters, overwriting the same buffer.  */
+	wcp = group_number (wstartp, wcp, intdig_no, grouping, thousands_sepwc,
+			    ngroups);
+      }
+
+    /* Write the exponent if it is needed.  */
+    if (type != 'f')
+      {
+	if (__builtin_expect (expsign != 0 && exponent == 4 && spec == 'g', 0))
+	  {
+	    /* This is another special case.  The exponent of the number is
+	       really smaller than -4, which requires the 'e'/'E' format.
+	       But after rounding the number has an exponent of -4.  */
+	    assert (wcp >= wstartp + 1);
+	    assert (wstartp[0] == L_('1'));
+	    memcpy (wstartp, L_("0.0001"), 6 * sizeof (wchar_t));
+	    wstartp[1] = decimalwc;
+	    if (wcp >= wstartp + 2)
+	      {
+		size_t cnt;
+		for (cnt = 0; cnt < wcp - (wstartp + 2); cnt++)
+		  wstartp[6 + cnt] = L_('0');
+		wcp += 4;
+	      }
+	    else
+	      wcp += 5;
+	  }
+	else
+	  {
+	    *wcp++ = (wchar_t) type;
+	    *wcp++ = expsign ? L_('-') : L_('+');
+
+	    /* Find the magnitude of the exponent.	*/
+	    expscale = 10;
+	    while (expscale <= exponent)
+	      expscale *= 10;
+
+	    if (exponent < 10)
+	      /* Exponent always has at least two digits.  */
+	      *wcp++ = L_('0');
+	    else
+	      do
+		{
+		  expscale /= 10;
+		  *wcp++ = L_('0') + (exponent / expscale);
+		  exponent %= expscale;
+		}
+	      while (expscale > 10);
+	    *wcp++ = L_('0') + exponent;
+	  }
+      }
+
+    /* Compute number of characters which must be filled with the padding
+       character.  */
+    if (is_neg || info->showsign || info->space)
+      --width;
+    width -= wcp - wstartp;
+
+    if (!info->left && info->pad != '0' && width > 0)
+      PADN (info->pad, width);
+
+    if (is_neg)
+      outchar ('-');
+    else if (info->showsign)
+      outchar ('+');
+    else if (info->space)
+      outchar (' ');
+
+    if (!info->left && info->pad == '0' && width > 0)
+      PADN ('0', width);
+
+    {
+      char *buffer = NULL;
+      char *buffer_end __attribute__((__unused__)) = NULL;
+      char *cp = NULL;
+      char *tmpptr;
+
+      if (! wide)
+	{
+	  /* Create the single byte string.  */
+	  size_t decimal_len;
+	  size_t thousands_sep_len;
+	  wchar_t *copywc;
+#ifdef USE_I18N_NUMBER_H
+	  size_t factor = (info->i18n
+			   ? nl_langinfo_wc (_NL_CTYPE_MB_CUR_MAX)
+			   : 1);
+#else
+	  size_t factor = 1;
+#endif
+
+	  decimal_len = strlen (decimal);
+
+	  if (thousands_sep == NULL)
+	    thousands_sep_len = 0;
+	  else
+	    thousands_sep_len = strlen (thousands_sep);
+
+	  size_t nbuffer = (2 + chars_needed * factor + decimal_len
+			    + ngroups * thousands_sep_len);
+	  if (__builtin_expect (buffer_malloced, 0))
+	    {
+	      buffer = (char *) malloc (nbuffer);
+	      if (buffer == NULL)
+		{
+		  /* Signal an error to the caller.  */
+		  free (wbuffer);
+		  return -1;
+		}
+	    }
+	  else
+	    buffer = (char *) alloca (nbuffer);
+	  buffer_end = buffer + nbuffer;
+
+	  /* Now copy the wide character string.  Since the character
+	     (except for the decimal point and thousands separator) must
+	     be coming from the ASCII range we can esily convert the
+	     string without mapping tables.  */
+	  for (cp = buffer, copywc = wstartp; copywc < wcp; ++copywc)
+	    if (*copywc == decimalwc)
+	      memcpy (cp, decimal, decimal_len), cp += decimal_len;
+	    else if (*copywc == thousands_sepwc)
+	      memcpy (cp, thousands_sep, thousands_sep_len), cp += thousands_sep_len;
+	    else
+	      *cp++ = (char) *copywc;
+	}
+
+      tmpptr = buffer;
+#if USE_I18N_NUMBER_H
+      if (__builtin_expect (info->i18n, 0))
+	{
+	  tmpptr = _i18n_number_rewrite (tmpptr, cp, buffer_end);
+	  cp = buffer_end;
+	  assert ((uintptr_t) buffer <= (uintptr_t) tmpptr);
+	  assert ((uintptr_t) tmpptr < (uintptr_t) buffer_end);
+	}
+#endif
+
+      PRINT (tmpptr, wstartp, wide ? wcp - wstartp : cp - tmpptr);
+
+      /* Free the memory if necessary.  */
+      if (__builtin_expect (buffer_malloced, 0))
+	{
+	  free (buffer);
+	  free (wbuffer);
+	}
+    }
+
+    if (info->left && width > 0)
+      PADN (info->pad, width);
+  }
+  return done;
+}
+
+/* Return the number of extra grouping characters that will be inserted
+   into a number with INTDIG_MAX integer digits.  */
+
+static unsigned int
+guess_grouping (unsigned int intdig_max, const char *grouping)
+{
+  unsigned int groups;
+
+  /* We treat all negative values like CHAR_MAX.  */
+
+  if (*grouping == CHAR_MAX || *grouping <= 0)
+    /* No grouping should be done.  */
+    return 0;
+
+  groups = 0;
+  while (intdig_max > (unsigned int) *grouping)
+    {
+      ++groups;
+      intdig_max -= *grouping++;
+
+      if (*grouping == 0)
+	{
+	  /* Same grouping repeats.  */
+	  groups += (intdig_max - 1) / grouping[-1];
+	  break;
+	}
+      else if (*grouping == CHAR_MAX || *grouping <= 0)
+	/* No more grouping should be done.  */
+	break;
+    }
+
+  return groups;
+}
+
+/* Group the INTDIG_NO integer digits of the number in [BUF,BUFEND).
+   There is guaranteed enough space past BUFEND to extend it.
+   Return the new end of buffer.  */
+
+static wchar_t *
+group_number (wchar_t *buf, wchar_t *bufend, unsigned int intdig_no,
+	      const char *grouping, wchar_t thousands_sep, int ngroups)
+{
+  wchar_t *p;
+
+  if (ngroups == 0)
+    return bufend;
+
+  /* Move the fractional part down.  */
+  memmove (buf + intdig_no + ngroups, buf + intdig_no,
+	   (bufend - (buf + intdig_no)) * sizeof (wchar_t));
+
+  p = buf + intdig_no + ngroups - 1;
+  do
+    {
+      unsigned int len = *grouping++;
+      do
+	*p-- = buf[--intdig_no];
+      while (--len > 0);
+      *p-- = thousands_sep;
+
+      if (*grouping == 0)
+	/* Same grouping repeats.  */
+	--grouping;
+      else if (*grouping == CHAR_MAX || *grouping <= 0)
+	/* No more grouping should be done.  */
+	break;
+    } while (intdig_no > (unsigned int) *grouping);
+
+  /* Copy the remaining ungrouped digits.  */
+  do
+    *p-- = buf[--intdig_no];
+  while (p > buf);
+
+  return bufend + ngroups;
+}
diff --git a/libquadmath/printf/printf_fphex.c b/libquadmath/printf/printf_fphex.c
new file mode 100644
index 000000000..941e93307
--- /dev/null
+++ b/libquadmath/printf/printf_fphex.c
@@ -0,0 +1,463 @@
+/* Print floating point number in hexadecimal notation according to ISO C99.
+   Copyright (C) 1997-2002,2004,2006 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <config.h>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#define NDEBUG
+#include <assert.h>
+#include "quadmath-printf.h"
+#include "_itoa.h"
+#include "_itowa.h"
+
+
+/* Macros for doing the actual output.  */
+
+#define outchar(ch)							      \
+  do									      \
+    {									      \
+      register const int outc = (ch);					      \
+      if (PUTC (outc, fp) == EOF)					      \
+	return -1;							      \
+      ++done;								      \
+    } while (0)
+
+#define PRINT(ptr, wptr, len)						      \
+  do									      \
+    {									      \
+      register size_t outlen = (len);					      \
+      if (wide)								      \
+	while (outlen-- > 0)						      \
+	  outchar (*wptr++);						      \
+      else								      \
+	while (outlen-- > 0)						      \
+	  outchar (*ptr++);						      \
+    } while (0)
+
+#define PADN(ch, len)							      \
+  do									      \
+    {									      \
+      if (PAD (fp, ch, len) != len)					      \
+	return -1;							      \
+      done += len;							      \
+    }									      \
+  while (0)
+
+
+
+int
+__quadmath_printf_fphex (struct __quadmath_printf_file *fp,
+			 const struct printf_info *info,
+			 const void *const *args)
+{
+  /* The floating-point value to output.  */
+  ieee854_float128 fpnum;
+
+  /* Locale-dependent representation of decimal point.	*/
+  const char *decimal;
+  wchar_t decimalwc;
+
+  /* "NaN" or "Inf" for the special cases.  */
+  const char *special = NULL;
+  const wchar_t *wspecial = NULL;
+
+  /* Buffer for the generated number string for the mantissa.  The
+     maximal size for the mantissa is 128 bits.  */
+  char numbuf[32];
+  char *numstr;
+  char *numend;
+  wchar_t wnumbuf[32];
+  wchar_t *wnumstr;
+  wchar_t *wnumend;
+  int negative;
+
+  /* The maximal exponent of two in decimal notation has 5 digits.  */
+  char expbuf[5];
+  char *expstr;
+  wchar_t wexpbuf[5];
+  wchar_t *wexpstr;
+  int expnegative;
+  int exponent;
+
+  /* Non-zero is mantissa is zero.  */
+  int zero_mantissa;
+
+  /* The leading digit before the decimal point.  */
+  char leading;
+
+  /* Precision.  */
+  int precision = info->prec;
+
+  /* Width.  */
+  int width = info->width;
+
+  /* Number of characters written.  */
+  int done = 0;
+
+  /* Nonzero if this is output on a wide character stream.  */
+  int wide = info->wide;
+
+  /* Figure out the decimal point character.  */
+#ifdef USE_NL_LANGINFO
+  if (info->extra == 0)
+    decimal = nl_langinfo (DECIMAL_POINT);
+  else
+    {
+      decimal = nl_langinfo (MON_DECIMAL_POINT);
+      if (*decimal == '\0')
+	decimal = nl_langinfo (DECIMAL_POINT);
+    }
+  /* The decimal point character must never be zero.  */
+  assert (*decimal != '\0');
+#elif defined USE_LOCALECONV
+  const struct lconv *lc = localeconv ();
+  if (info->extra == 0)
+    decimal = lc->decimal_point;
+  else
+    {
+      decimal = lc->mon_decimal_point;
+      if (decimal == NULL || *decimal == '\0')
+	decimal = lc->decimal_point;
+    }
+  if (decimal == NULL || *decimal == '\0')
+    decimal = ".";
+#else
+  decimal = ".";
+#endif
+#ifdef USE_NL_LANGINFO_WC
+  if (info->extra == 0)
+    decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC);
+  else
+    {
+      decimalwc = nl_langinfo_wc (_NL_MONETARY_DECIMAL_POINT_WC);
+      if (decimalwc == L_('\0'))
+	decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC);
+    }
+  /* The decimal point character must never be zero.  */
+  assert (decimalwc != L_('\0'));
+#else
+  decimalwc = L_('.');
+#endif
+
+  /* Fetch the argument value.	*/
+    {
+      fpnum.value = **(const __float128 **) args[0];
+
+      /* Check for special values: not a number or infinity.  */
+      if (isnanq (fpnum.value))
+	{
+	  negative = fpnum.ieee.negative != 0;
+	  if (isupper (info->spec))
+	    {
+	      special = "NAN";
+	      wspecial = L_("NAN");
+	    }
+	  else
+	    {
+	      special = "nan";
+	      wspecial = L_("nan");
+	    }
+	}
+      else
+	{
+	  if (isinfq (fpnum.value))
+	    {
+	      if (isupper (info->spec))
+		{
+		  special = "INF";
+		  wspecial = L_("INF");
+		}
+	      else
+		{
+		  special = "inf";
+		  wspecial = L_("inf");
+		}
+	    }
+
+	  negative = signbitq (fpnum.value);
+	}
+    }
+
+  if (special)
+    {
+      int width = info->width;
+
+      if (negative || info->showsign || info->space)
+	--width;
+      width -= 3;
+
+      if (!info->left && width > 0)
+	PADN (' ', width);
+
+      if (negative)
+	outchar ('-');
+      else if (info->showsign)
+	outchar ('+');
+      else if (info->space)
+	outchar (' ');
+
+      PRINT (special, wspecial, 3);
+
+      if (info->left && width > 0)
+	PADN (' ', width);
+
+      return done;
+    }
+
+    {
+      /* We have 112 bits of mantissa plus one implicit digit.  Since
+	 112 bits are representable without rest using hexadecimal
+	 digits we use only the implicit digits for the number before
+	 the decimal point.  */
+      uint64_t num0, num1;
+
+      assert (sizeof (long double) == 16);
+
+      num0 = fpnum.ieee.mant_high;
+      num1 = fpnum.ieee.mant_low;
+
+      zero_mantissa = (num0|num1) == 0;
+
+      if (sizeof (unsigned long int) > 6)
+	{
+	  numstr = _itoa_word (num1, numbuf + sizeof numbuf, 16,
+			       info->spec == 'A');
+	  wnumstr = _itowa_word (num1,
+				 wnumbuf + sizeof (wnumbuf) / sizeof (wchar_t),
+				 16, info->spec == 'A');
+	}
+      else
+	{
+	  numstr = _itoa (num1, numbuf + sizeof numbuf, 16,
+			  info->spec == 'A');
+	  wnumstr = _itowa (num1,
+			    wnumbuf + sizeof (wnumbuf) / sizeof (wchar_t),
+			    16, info->spec == 'A');
+	}
+
+      while (numstr > numbuf + (sizeof numbuf - 64 / 4))
+	{
+	  *--numstr = '0';
+	  *--wnumstr = L_('0');
+	}
+
+      if (sizeof (unsigned long int) > 6)
+	{
+	  numstr = _itoa_word (num0, numstr, 16, info->spec == 'A');
+	  wnumstr = _itowa_word (num0, wnumstr, 16, info->spec == 'A');
+	}
+      else
+	{
+	  numstr = _itoa (num0, numstr, 16, info->spec == 'A');
+	  wnumstr = _itowa (num0, wnumstr, 16, info->spec == 'A');
+	}
+
+      /* Fill with zeroes.  */
+      while (numstr > numbuf + (sizeof numbuf - 112 / 4))
+	{
+	  *--numstr = '0';
+	  *--wnumstr = L_('0');
+	}
+
+      leading = fpnum.ieee.exponent == 0 ? '0' : '1';
+
+      exponent = fpnum.ieee.exponent;
+
+      if (exponent == 0)
+	{
+	  if (zero_mantissa)
+	    expnegative = 0;
+	  else
+	    {
+	      /* This is a denormalized number.  */
+	      expnegative = 1;
+	      exponent = IEEE854_FLOAT128_BIAS - 1;
+	    }
+	}
+      else if (exponent >= IEEE854_FLOAT128_BIAS)
+	{
+	  expnegative = 0;
+	  exponent -= IEEE854_FLOAT128_BIAS;
+	}
+      else
+	{
+	  expnegative = 1;
+	  exponent = -(exponent - IEEE854_FLOAT128_BIAS);
+	}
+    }
+
+  /* Look for trailing zeroes.  */
+  if (! zero_mantissa)
+    {
+      wnumend = &wnumbuf[sizeof wnumbuf / sizeof wnumbuf[0]];
+      numend = &numbuf[sizeof numbuf / sizeof numbuf[0]];
+      while (wnumend[-1] == L_('0'))
+	{
+	  --wnumend;
+	  --numend;
+	}
+
+      if (precision == -1)
+	precision = numend - numstr;
+      else if (precision < numend - numstr
+	       && (numstr[precision] > '8'
+		   || (('A' < '0' || 'a' < '0')
+		       && numstr[precision] < '0')
+		   || (numstr[precision] == '8'
+		       && (precision + 1 < numend - numstr
+			   /* Round to even.  */
+			   || (precision > 0
+			       && ((numstr[precision - 1] & 1)
+				   ^ (isdigit (numstr[precision - 1]) == 0)))
+			   || (precision == 0
+			       && ((leading & 1)
+				   ^ (isdigit (leading) == 0)))))))
+	{
+	  /* Round up.  */
+	  int cnt = precision;
+	  while (--cnt >= 0)
+	    {
+	      char ch = numstr[cnt];
+	      /* We assume that the digits and the letters are ordered
+		 like in ASCII.  This is true for the rest of GNU, too.  */
+	      if (ch == '9')
+		{
+		  wnumstr[cnt] = (wchar_t) info->spec;
+		  numstr[cnt] = info->spec;	/* This is tricky,
+		  				   think about it!  */
+		  break;
+		}
+	      else if (tolower (ch) < 'f')
+		{
+		  ++numstr[cnt];
+		  ++wnumstr[cnt];
+		  break;
+		}
+	      else
+		{
+		  numstr[cnt] = '0';
+		  wnumstr[cnt] = L_('0');
+		}
+	    }
+	  if (cnt < 0)
+	    {
+	      /* The mantissa so far was fff...f  Now increment the
+		 leading digit.  Here it is again possible that we
+		 get an overflow.  */
+	      if (leading == '9')
+		leading = info->spec;
+	      else if (tolower (leading) < 'f')
+		++leading;
+	      else
+		{
+		  leading = '1';
+		  if (expnegative)
+		    {
+		      exponent -= 4;
+		      if (exponent <= 0)
+			{
+			  exponent = -exponent;
+			  expnegative = 0;
+			}
+		    }
+		  else
+		    exponent += 4;
+		}
+	    }
+	}
+    }
+  else
+    {
+      if (precision == -1)
+	precision = 0;
+      numend = numstr;
+      wnumend = wnumstr;
+    }
+
+  /* Now we can compute the exponent string.  */
+  expstr = _itoa_word (exponent, expbuf + sizeof expbuf, 10, 0);
+  wexpstr = _itowa_word (exponent,
+			 wexpbuf + sizeof wexpbuf / sizeof (wchar_t), 10, 0);
+
+  /* Now we have all information to compute the size.  */
+  width -= ((negative || info->showsign || info->space)
+	    /* Sign.  */
+	    + 2    + 1 + 0 + precision + 1 + 1
+	    /* 0x    h   .   hhh         P   ExpoSign.  */
+	    + ((expbuf + sizeof expbuf) - expstr));
+	    /* Exponent.  */
+
+  /* Count the decimal point.
+     A special case when the mantissa or the precision is zero and the `#'
+     is not given.  In this case we must not print the decimal point.  */
+  if (precision > 0 || info->alt)
+    width -= wide ? 1 : strlen (decimal);
+
+  if (!info->left && info->pad != '0' && width > 0)
+    PADN (' ', width);
+
+  if (negative)
+    outchar ('-');
+  else if (info->showsign)
+    outchar ('+');
+  else if (info->space)
+    outchar (' ');
+
+  outchar ('0');
+  if ('X' - 'A' == 'x' - 'a')
+    outchar (info->spec + ('x' - 'a'));
+  else
+    outchar (info->spec == 'A' ? 'X' : 'x');
+
+  if (!info->left && info->pad == '0' && width > 0)
+    PADN ('0', width);
+
+  outchar (leading);
+
+  if (precision > 0 || info->alt)
+    {
+      const wchar_t *wtmp = &decimalwc;
+      PRINT (decimal, wtmp, wide ? 1 : strlen (decimal));
+    }
+
+  if (precision > 0)
+    {
+      ssize_t tofill = precision - (numend - numstr);
+      PRINT (numstr, wnumstr, MIN (numend - numstr, precision));
+      if (tofill > 0)
+	PADN ('0', tofill);
+    }
+
+  if ('P' - 'A' == 'p' - 'a')
+    outchar (info->spec + ('p' - 'a'));
+  else
+    outchar (info->spec == 'A' ? 'P' : 'p');
+
+  outchar (expnegative ? '-' : '+');
+
+  PRINT (expstr, wexpstr, (expbuf + sizeof expbuf) - expstr);
+
+  if (info->left && info->pad != '0' && width > 0)
+    PADN (info->pad, width);
+
+  return done;
+}
diff --git a/libquadmath/printf/quadmath-printf.c b/libquadmath/printf/quadmath-printf.c
new file mode 100644
index 000000000..b70f432cc
--- /dev/null
+++ b/libquadmath/printf/quadmath-printf.c
@@ -0,0 +1,422 @@
+/* GCC Quad-Precision Math Library
+   Copyright (C) 2011 Free Software Foundation, Inc.
+   Written by Jakub Jelinek  <jakub@redhat.com>
+
+This file is part of the libquadmath library.
+Libquadmath is free software; you can redistribute it and/or
+modify it under the terms of the GNU Library General Public
+License as published by the Free Software Foundation; either
+version 2 of the License, or (at your option) any later version.
+
+Libquadmath is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with libquadmath; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include <config.h>
+#include <stdarg.h>
+#include <string.h>
+#include <stdio.h>
+#include "quadmath-printf.h"
+
+/* Read a simple integer from a string and update the string pointer.
+   It is assumed that the first character is a digit.  */
+static unsigned int
+read_int (const char **pstr)
+{
+  unsigned int retval = (unsigned char) **pstr - '0';
+
+  while (isdigit ((unsigned char) *++(*pstr)))
+    {
+      retval *= 10;
+      retval += (unsigned char) **pstr - '0';
+    }
+
+  return retval;
+}
+
+#define PADSIZE 16
+static char const blanks[PADSIZE] =
+{' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' '};
+static char const zeroes[PADSIZE] =
+{'0','0','0','0','0','0','0','0','0','0','0','0','0','0','0','0'};
+static wchar_t const wblanks[PADSIZE] =
+{
+  L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' '),
+  L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' ')
+};
+static wchar_t const wzeroes[PADSIZE] =
+{
+  L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0'),
+  L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0')
+};
+
+attribute_hidden size_t
+__quadmath_do_pad (struct __quadmath_printf_file *fp, int wide, int c,
+		   size_t n)
+{
+  ssize_t i;
+  char padbuf[PADSIZE];
+  wchar_t wpadbuf[PADSIZE];
+  const char *padstr;
+  size_t w, written = 0;
+  if (wide)
+    {
+      if (c == ' ')
+	padstr = (const char *) wblanks;
+      else if (c == '0')
+	padstr = (const char *) wzeroes;
+      else
+	{
+	  padstr = (const char *) wpadbuf;
+	  for (i = 0; i < PADSIZE; i++)
+	    wpadbuf[i] = c;
+	}
+    }
+  else
+    {
+      if (c == ' ')
+	padstr = blanks;
+      else if (c == '0')
+	padstr = zeroes;
+      else
+	{
+	  padstr = (const char *) padbuf;
+	  for (i = 0; i < PADSIZE; i++)
+	    padbuf[i] = c;
+	}
+    }
+  for (i = n; i >= PADSIZE; i -= PADSIZE)
+    {
+      w = PUT (fp, (char *) padstr, PADSIZE);
+      written += w;
+      if (w != PADSIZE)
+	return written;
+    }
+  if (i > 0)
+    {
+      w = PUT (fp, (char *) padstr, i);
+      written += w;
+    }
+  return written;
+}
+
+/* This is a stripped down version of snprintf, which just handles
+   a single %eEfFgGaA format entry with Q modifier.  % has to be
+   the first character of the format string, no $ can be used.  */
+int
+quadmath_snprintf (char *str, size_t size, const char *format, ...)
+{
+  struct printf_info info;
+  va_list ap;
+  __float128 fpnum, *fpnum_addr = &fpnum, **fpnum_addr2 = &fpnum_addr;
+  struct __quadmath_printf_file qfp;
+
+  if (*format++ != '%')
+    return -1;
+
+  /* Clear information structure.  */
+  memset (&info, '\0', sizeof info);
+  /* info.alt = 0;
+  info.space = 0;
+  info.left = 0;
+  info.showsign = 0;
+  info.group = 0;
+  info.i18n = 0;
+  info.extra = 0; */
+  info.pad = ' ';
+  /* info.wide = 0; */
+
+  /* Check for spec modifiers.  */
+  do
+    {
+      switch (*format)
+	{
+	case ' ':
+	  /* Output a space in place of a sign, when there is no sign.  */
+	  info.space = 1;
+	  continue;
+	case '+':
+	  /* Always output + or - for numbers.  */
+	  info.showsign = 1;
+	  continue;
+	case '-':
+	  /* Left-justify things.  */
+	  info.left = 1;
+	  continue;
+	case '#':
+	  /* Use the "alternate form":
+	     Hex has 0x or 0X, FP always has a decimal point.  */
+	  info.alt = 1;
+	  continue;
+	case '0':
+	  /* Pad with 0s.  */
+	  info.pad = '0';
+	  continue;
+	case '\'':
+	  /* Show grouping in numbers if the locale information
+	     indicates any.  */
+	  info.group = 1;
+	  continue;
+	case 'I':
+	  /* Use the internationalized form of the output.  Currently
+	     means to use the `outdigits' of the current locale.  */
+	  info.i18n = 1;
+	  continue;
+	default:
+	  break;
+	}
+      break;
+    }
+  while (*++format);
+
+  if (info.left)
+    info.pad = ' ';
+
+  va_start (ap, format);
+
+  /* Get the field width.  */
+  /* info.width = 0; */
+  if (*format == '*')
+    {
+      /* The field width is given in an argument.
+	 A negative field width indicates left justification.  */
+      ++format;
+      info.width = va_arg (ap, int);
+    }
+  else if (isdigit (*format))
+    /* Constant width specification.  */
+    info.width = read_int (&format);
+
+  /* Get the precision.  */
+  /* -1 means none given; 0 means explicit 0.  */
+  info.prec = -1;
+  if (*format == '.')
+    {
+      ++format;
+      if (*format == '*')
+	{
+	  /* The precision is given in an argument.  */
+	  ++format;
+
+	  info.prec = va_arg (ap, int);
+	}
+      else if (isdigit (*format))
+	info.prec = read_int (&format);
+      else
+	/* "%.?" is treated like "%.0?".  */
+	info.prec = 0;
+    }
+
+  /* Check for type modifiers.  */
+  /* info.is_long_double = 0;
+  info.is_short = 0;
+  info.is_long = 0;
+  info.is_char = 0;
+  info.user = 0; */
+
+  /* We require Q modifier.  */
+  if (*format++ != 'Q')
+    {
+      va_end (ap);
+      return -1;
+    }
+
+  /* Get the format specification.  */
+  info.spec = (wchar_t) *format++;
+  if (info.spec == L_('\0') || *format != '\0')
+    {
+      va_end (ap);
+      return -1;
+    }
+
+  switch (info.spec)
+    {
+    case L_('e'):
+    case L_('E'):
+    case L_('f'):
+    case L_('F'):
+    case L_('g'):
+    case L_('G'):
+    case L_('a'):
+    case L_('A'):
+      break;
+    default:
+      va_end (ap);
+      return -1;
+    }
+
+  fpnum = va_arg (ap, __float128);
+  va_end (ap);
+
+  qfp.fp = NULL;
+  qfp.str = str;
+  qfp.size = size ? size - 1 : 0;
+  qfp.len = 0;
+  qfp.file_p = 0;
+
+  if (info.spec == L_('a') || info.spec == L_('A'))
+    __quadmath_printf_fphex (&qfp, &info, (const void *const *)&fpnum_addr2);
+  else
+    __quadmath_printf_fp (&qfp, &info, (const void *const *)&fpnum_addr2);
+
+  if (size)
+    *qfp.str = '\0';
+
+  return qfp.len;
+}
+
+#ifdef HAVE_PRINTF_HOOKS
+static int pa_flt128;
+int mod_Q attribute_hidden;
+
+static void
+flt128_va (void *mem, va_list *ap)
+{ 
+  __float128 d = va_arg (*ap, __float128);
+  memcpy (mem, &d, sizeof (d));
+}
+
+static int
+flt128_ais (const struct printf_info *info, size_t n __attribute__ ((unused)),
+	    int *argtype, int *size)
+{
+  if (info->user & mod_Q)
+    {
+      argtype[0] = pa_flt128;
+      size[0] = sizeof (__float128);
+      return 1;
+    }
+#if __GLIBC__ < 2 || (__GLIBC__ == 2 && __GLIBC_MINOR__ <= 13)
+  /* Workaround bug in glibc printf hook handling.  */
+  size[0] = -1;
+  switch (info->spec)
+    {
+    case L_('i'):
+    case L_('d'):
+    case L_('u'):
+    case L_('o'):
+    case L_('X'):
+    case L_('x'):
+#if __LONG_MAX__ != __LONG_LONG_MAX__
+      if (info->is_long_double)
+	argtype[0] = PA_INT|PA_FLAG_LONG_LONG;
+      else
+#endif
+      if (info->is_long)
+	argtype[0] = PA_INT|PA_FLAG_LONG;
+      else if (info->is_short)
+	argtype[0] = PA_INT|PA_FLAG_SHORT;
+      else if (info->is_char)
+	argtype[0] = PA_CHAR;
+      else
+	argtype[0] = PA_INT;
+      return 1;
+    case L_('e'):
+    case L_('E'):
+    case L_('f'):
+    case L_('F'):
+    case L_('g'):
+    case L_('G'):
+    case L_('a'):
+    case L_('A'):
+      if (info->is_long_double)
+	argtype[0] = PA_DOUBLE|PA_FLAG_LONG_DOUBLE;
+      else
+	argtype[0] = PA_DOUBLE;
+      return 1;
+    case L_('c'):
+      argtype[0] = PA_CHAR;
+      return 1;
+    case L_('C'):
+      argtype[0] = PA_WCHAR;
+      return 1;
+    case L_('s'):
+      argtype[0] = PA_STRING;
+      return 1;
+    case L_('S'):
+      argtype[0] = PA_WSTRING;
+      return 1;
+    case L_('p'):
+      argtype[0] = PA_POINTER;
+      return 1;
+    case L_('n'):
+      argtype[0] = PA_INT|PA_FLAG_PTR;
+      return 1;
+
+    case L_('m'):
+    default:
+      /* An unknown spec will consume no args.  */
+      return 0;
+    }
+#endif
+  return -1;
+}
+
+static int
+flt128_printf_fp (FILE *fp, const struct printf_info *info,
+		  const void *const *args)
+{
+  struct __quadmath_printf_file qpf
+    = { .fp = fp, .str = NULL, .size = 0, .len = 0, .file_p = 1 };
+
+  if ((info->user & mod_Q) == 0)
+    return -2;
+
+  return __quadmath_printf_fp (&qpf, info, args);
+}
+
+static int
+flt128_printf_fphex (FILE *fp, const struct printf_info *info,
+		     const void *const *args)
+{
+  struct __quadmath_printf_file qpf
+    = { .fp = fp, .str = NULL, .size = 0, .len = 0, .file_p = 1 };
+
+  if ((info->user & mod_Q) == 0)
+    return -2;
+
+  return __quadmath_printf_fphex (&qpf, info, args);
+}
+
+__attribute__((constructor)) static void
+register_printf_flt128 (void)
+{
+  pa_flt128 = register_printf_type (flt128_va);
+  if (pa_flt128 == -1)
+    return;
+  mod_Q = register_printf_modifier (L_("Q"));
+  if (mod_Q == -1)
+    return;
+  register_printf_specifier ('f', flt128_printf_fp, flt128_ais);
+  register_printf_specifier ('F', flt128_printf_fp, flt128_ais);
+  register_printf_specifier ('e', flt128_printf_fp, flt128_ais);
+  register_printf_specifier ('E', flt128_printf_fp, flt128_ais);
+  register_printf_specifier ('g', flt128_printf_fp, flt128_ais);
+  register_printf_specifier ('G', flt128_printf_fp, flt128_ais);
+  register_printf_specifier ('a', flt128_printf_fphex, flt128_ais);
+  register_printf_specifier ('A', flt128_printf_fphex, flt128_ais);
+}
+
+__attribute__((destructor)) static void
+unregister_printf_flt128 (void)
+{
+  /* No way to unregister printf type and modifier currently,
+     and only one printf specifier can be registered right now.  */
+  if (pa_flt128 == -1 || mod_Q == -1)
+    return;
+  register_printf_specifier ('f', NULL, NULL);
+  register_printf_specifier ('F', NULL, NULL);
+  register_printf_specifier ('e', NULL, NULL);
+  register_printf_specifier ('E', NULL, NULL);
+  register_printf_specifier ('g', NULL, NULL);
+  register_printf_specifier ('G', NULL, NULL);
+  register_printf_specifier ('a', NULL, NULL);
+  register_printf_specifier ('A', NULL, NULL);
+}
+#endif
diff --git a/libquadmath/printf/quadmath-printf.h b/libquadmath/printf/quadmath-printf.h
new file mode 100644
index 000000000..32ebcec92
--- /dev/null
+++ b/libquadmath/printf/quadmath-printf.h
@@ -0,0 +1,186 @@
+/* GCC Quad-Precision Math Library
+   Copyright (C) 2011 Free Software Foundation, Inc.
+   Written by Jakub Jelinek  <jakub@redhat.com>
+
+This file is part of the libquadmath library.
+Libquadmath is free software; you can redistribute it and/or
+modify it under the terms of the GNU Library General Public
+License as published by the Free Software Foundation; either
+version 2 of the License, or (at your option) any later version.
+
+Libquadmath is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with libquadmath; see the file COPYING.LIB.  If
+not, write to the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+#include <stdlib.h>
+#include <stdio.h>
+#ifdef HAVE_LIMITS_H
+#include <limits.h>
+#endif
+#ifdef HAVE_LANGINFO_H
+#include <langinfo.h>
+#endif
+#ifdef HAVE_CTYPE_H
+#include <ctype.h>
+#endif
+#ifdef HAVE_WCHAR_H
+#include <wchar.h>
+#endif
+#ifdef HAVE_WCTYPE_H
+#include <wctype.h>
+#endif
+#ifdef HAVE_PRINTF_HOOKS
+#include <printf.h>
+#endif
+#ifdef HAVE_LOCALE_H
+#include <locale.h>
+#endif
+#include "quadmath-imp.h"
+#include "gmp-impl.h"
+
+#ifdef HAVE_WCHAR_H
+#define L_(x) L##x
+#else
+#define L_(x) x
+#undef wchar_t
+#undef wint_t
+#undef putwc
+#undef WEOF
+#define wchar_t char
+#define wint_t int
+#define putwc(c,f) putc(c,f)
+#define WEOF EOF
+#endif
+
+#ifndef HAVE_CTYPE_H
+/* Won't work for EBCDIC.  */
+#undef isupper
+#undef isdigit
+#undef isxdigit
+#undef tolower
+#define isupper(x) \
+  ({__typeof(x) __is_x = (x); __is_x >= 'A' && __is_x <= 'Z'; })
+#define isdigit(x) \
+  ({__typeof(x) __is_x = (x); __is_x >= '0' && __is_x <= '9'; })
+#define isxdigit(x) \
+  ({__typeof(x) __is_x = (x); \
+    (__is_x >= '0' && __is_x <= '9') \
+    || ((x) >= 'A' && (x) <= 'F') \
+    || ((x) >= 'a' && (x) <= 'f'); })
+#define tolower(x) \
+  ({__typeof(x) __is_x = (x); \
+    (__is_x >= 'A' && __is_x <= 'Z') ? __is_x - 'A' + 'a' : __is_x; })
+#endif
+
+#ifndef CHAR_MAX
+#ifdef __CHAR_UNSIGNED__
+#define CHAR_MAX (2 * __SCHAR_MAX__ + 1)
+#else
+#define CHAR_MAX __SCHAR_MAX__
+#endif
+#endif
+
+#ifndef HAVE_PRINTF_HOOKS
+#define printf_info __quadmath_printf_info
+struct printf_info
+{
+  int prec;			/* Precision.  */
+  int width;			/* Width.  */
+  wchar_t spec;			/* Format letter.  */
+  unsigned int is_long_double:1;/* L flag.  */
+  unsigned int is_short:1;	/* h flag.  */
+  unsigned int is_long:1;	/* l flag.  */
+  unsigned int alt:1;		/* # flag.  */
+  unsigned int space:1;		/* Space flag.  */
+  unsigned int left:1;		/* - flag.  */
+  unsigned int showsign:1;	/* + flag.  */
+  unsigned int group:1;		/* ' flag.  */
+  unsigned int extra:1;		/* For special use.  */
+  unsigned int is_char:1;	/* hh flag.  */
+  unsigned int wide:1;		/* Nonzero for wide character streams.  */
+  unsigned int i18n:1;		/* I flag.  */
+  unsigned short int user;	/* Bits for user-installed modifiers.  */
+  wchar_t pad;			/* Padding character.  */
+};
+#endif
+
+struct __quadmath_printf_file
+{
+  FILE *fp;
+  char *str;
+  size_t size;
+  size_t len;
+  int file_p;
+};
+
+int
+__quadmath_printf_fp (struct __quadmath_printf_file *fp,
+		      const struct printf_info *info,
+		      const void *const *args) attribute_hidden;
+int
+__quadmath_printf_fphex (struct __quadmath_printf_file *fp,
+			 const struct printf_info *info,
+			 const void *const *args) attribute_hidden;
+
+size_t __quadmath_do_pad (struct __quadmath_printf_file *fp, int wide,
+			  int c, size_t n) attribute_hidden;
+
+static inline __attribute__((__unused__)) size_t
+__quadmath_do_put (struct __quadmath_printf_file *fp, int wide,
+		   const char *s, size_t n)
+{
+  size_t len;
+  if (fp->file_p)
+    {
+      if (wide)
+	{
+	  size_t cnt;
+	  const wchar_t *ls = (const wchar_t *) s;
+	  for (cnt = 0; cnt < n; cnt++)
+	    if (putwc (ls[cnt], fp->fp) == WEOF)
+	      break;
+	  return cnt;
+	}
+      return fwrite (s, 1, n, fp->fp);
+    }
+  len = MIN (fp->size, n);
+  memcpy (fp->str, s, len);
+  fp->str += len;
+  fp->size -= len;
+  fp->len += n;
+  return n;
+}
+
+static inline __attribute__((__unused__)) int
+__quadmath_do_putc (struct __quadmath_printf_file *fp, int wide,
+		    wchar_t c)
+{
+  if (fp->file_p)
+    return wide ? (int) putwc (c, fp->fp) : putc (c, fp->fp);
+  if (fp->size)
+    {
+      *(fp->str++) = c;
+      fp->size--;
+    }
+  fp->len++;
+  return (unsigned char) c;
+}
+
+#define PUT(f, s, n) __quadmath_do_put (f, wide, s, n)
+#define PAD(f, c, n) __quadmath_do_pad (f, wide, c, n)
+#define PUTC(c, f) __quadmath_do_putc (f, wide, c)
+
+#define nl_langinfo_wc(x) \
+  ({ union { const char *mb; wchar_t wc; } u; u.mb = nl_langinfo (x); u.wc; })
+
+#undef _itoa
+#define _itoa __quadmath_itoa
+
+#undef NAN
+#define NAN __builtin_nanf ("")
diff --git a/libquadmath/printf/rshift.c b/libquadmath/printf/rshift.c
new file mode 100644
index 000000000..17fd914d9
--- /dev/null
+++ b/libquadmath/printf/rshift.c
@@ -0,0 +1,88 @@
+/* mpn_rshift -- Shift right a low-level natural-number integer.
+
+Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+/* Shift U (pointed to by UP and USIZE limbs long) CNT bits to the right
+   and store the USIZE least significant limbs of the result at WP.
+   The bits shifted out to the right are returned.
+
+   Argument constraints:
+   1. 0 < CNT < BITS_PER_MP_LIMB
+   2. If the result is to be written over the input, WP must be <= UP.
+*/
+
+mp_limb_t
+#if __STDC__
+mpn_rshift (register mp_ptr wp,
+	    register mp_srcptr up, mp_size_t usize,
+	    register unsigned int cnt)
+#else
+mpn_rshift (wp, up, usize, cnt)
+     register mp_ptr wp;
+     register mp_srcptr up;
+     mp_size_t usize;
+     register unsigned int cnt;
+#endif
+{
+  register mp_limb_t high_limb, low_limb;
+  register unsigned sh_1, sh_2;
+  register mp_size_t i;
+  mp_limb_t retval;
+
+#ifdef DEBUG
+  if (usize == 0 || cnt == 0)
+    abort ();
+#endif
+
+  sh_1 = cnt;
+
+#if 0
+  if (sh_1 == 0)
+    {
+      if (wp != up)
+	{
+	  /* Copy from low end to high end, to allow specified input/output
+	     overlapping.  */
+	  for (i = 0; i < usize; i++)
+	    wp[i] = up[i];
+	}
+      return usize;
+    }
+#endif
+
+  wp -= 1;
+  sh_2 = BITS_PER_MP_LIMB - sh_1;
+  high_limb = up[0];
+  retval = high_limb << sh_2;
+  low_limb = high_limb;
+
+  for (i = 1; i < usize; i++)
+    {
+      high_limb = up[i];
+      wp[i] = (low_limb >> sh_1) | (high_limb << sh_2);
+      low_limb = high_limb;
+    }
+  wp[i] = low_limb >> sh_1;
+
+  return retval;
+}
diff --git a/libquadmath/printf/sub_n.c b/libquadmath/printf/sub_n.c
new file mode 100644
index 000000000..92e718731
--- /dev/null
+++ b/libquadmath/printf/sub_n.c
@@ -0,0 +1,62 @@
+/* mpn_sub_n -- Subtract two limb vectors of equal, non-zero length.
+
+Copyright (C) 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+mp_limb_t
+#if __STDC__
+mpn_sub_n (mp_ptr res_ptr, mp_srcptr s1_ptr, mp_srcptr s2_ptr, mp_size_t size)
+#else
+mpn_sub_n (res_ptr, s1_ptr, s2_ptr, size)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     register mp_srcptr s2_ptr;
+     mp_size_t size;
+#endif
+{
+  register mp_limb_t x, y, cy;
+  register mp_size_t j;
+
+  /* The loop counter and index J goes from -SIZE to -1.  This way
+     the loop becomes faster.  */
+  j = -size;
+
+  /* Offset the base pointers to compensate for the negative indices.  */
+  s1_ptr -= j;
+  s2_ptr -= j;
+  res_ptr -= j;
+
+  cy = 0;
+  do
+    {
+      y = s2_ptr[j];
+      x = s1_ptr[j];
+      y += cy;			/* add previous carry to subtrahend */
+      cy = (y < cy);		/* get out carry from that addition */
+      y = x - y;		/* main subtract */
+      cy = (y > x) + cy;	/* get out carry from the subtract, combine */
+      res_ptr[j] = y;
+    }
+  while (++j != 0);
+
+  return cy;
+}
diff --git a/libquadmath/printf/submul_1.c b/libquadmath/printf/submul_1.c
new file mode 100644
index 000000000..31903c628
--- /dev/null
+++ b/libquadmath/printf/submul_1.c
@@ -0,0 +1,64 @@
+/* mpn_submul_1 -- multiply the S1_SIZE long limb vector pointed to by S1_PTR
+   by S2_LIMB, subtract the S1_SIZE least significant limbs of the product
+   from the limb vector pointed to by RES_PTR.  Return the most significant
+   limb of the product, adjusted for carry-out from the subtraction.
+
+Copyright (C) 1992, 1993, 1994, 1996, 2005 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <config.h>
+#include "gmp-impl.h"
+
+mp_limb_t
+mpn_submul_1 (res_ptr, s1_ptr, s1_size, s2_limb)
+     register mp_ptr res_ptr;
+     register mp_srcptr s1_ptr;
+     mp_size_t s1_size;
+     register mp_limb_t s2_limb;
+{
+  register mp_limb_t cy_limb;
+  register mp_size_t j;
+  register mp_limb_t prod_high, prod_low;
+  register mp_limb_t x;
+
+  /* The loop counter and index J goes from -SIZE to -1.  This way
+     the loop becomes faster.  */
+  j = -s1_size;
+
+  /* Offset the base pointers to compensate for the negative indices.  */
+  res_ptr -= j;
+  s1_ptr -= j;
+
+  cy_limb = 0;
+  do
+    {
+      umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
+
+      prod_low += cy_limb;
+      cy_limb = (prod_low < cy_limb) + prod_high;
+
+      x = res_ptr[j];
+      prod_low = x - prod_low;
+      cy_limb += (prod_low > x);
+      res_ptr[j] = prod_low;
+    }
+  while (++j != 0);
+
+  return cy_limb;
+}