2005-12-13 Ulrich Drepper <drepper@redhat.com>

9 files changed, 896 insertions, 2 deletions

diff --git a/stdlib/abort.c b/stdlib/abort.c
new file mode 100644
index 0000000000..00788f22c7
--- /dev/null
+++ b/stdlib/abort.c
@@ -0,0 +1,140 @@
+/* Copyright (C) 1991,93,95,96,97,98,2001,02 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 <bits/libc-lock.h>
+#include <signal.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+
+/* Try to get a machine dependent instruction which will make the
+   program crash.  This is used in case everything else fails.  */
+#include <abort-instr.h>
+#ifndef ABORT_INSTRUCTION
+/* No such instruction is available.  */
+# define ABORT_INSTRUCTION
+#endif
+
+#ifdef USE_IN_LIBIO
+# include <libio/libioP.h>
+# define fflush(s) _IO_flush_all_lockp (0)
+#endif
+
+/* We must avoid to run in circles.  Therefore we remember how far we
+   already got.  */
+static int stage;
+
+/* We should be prepared for multiple threads trying to run abort.  */
+__libc_lock_define_initialized_recursive (static, lock);
+
+
+/* Cause an abnormal program termination with core-dump.  */
+void
+abort (void)
+{
+  struct sigaction act;
+  sigset_t sigs;
+
+  /* First acquire the lock.  */
+  __libc_lock_lock_recursive (lock);
+
+  /* Now it's for sure we are alone.  But recursive calls are possible.  */
+
+  /* Unlock SIGABRT.  */
+  if (stage == 0)
+    {
+      ++stage;
+      if (__sigemptyset (&sigs) == 0 &&
+	  __sigaddset (&sigs, SIGABRT) == 0)
+	__sigprocmask (SIG_UNBLOCK, &sigs, (sigset_t *) NULL);
+    }
+
+  /* Flush all streams.  We cannot close them now because the user
+     might have registered a handler for SIGABRT.  */
+  if (stage == 1)
+    {
+      ++stage;
+      fflush (NULL);
+    }
+
+  /* Send signal which possibly calls a user handler.  */
+  if (stage == 2)
+    {
+      /* This stage is special: we must allow repeated calls of
+	 `abort' when a user defined handler for SIGABRT is installed.
+	 This is risky since the `raise' implementation might also
+	 fail but I don't see another possibility.  */
+      int save_stage = stage;
+
+      stage = 0;
+      __libc_lock_unlock_recursive (lock);
+
+      raise (SIGABRT);
+
+      __libc_lock_lock_recursive (lock);
+      stage = save_stage + 1;
+    }
+
+  /* There was a handler installed.  Now remove it.  */
+  if (stage == 3)
+    {
+      ++stage;
+      memset (&act, '\0', sizeof (struct sigaction));
+      act.sa_handler = SIG_DFL;
+      __sigfillset (&act.sa_mask);
+      act.sa_flags = 0;
+      __sigaction (SIGABRT, &act, NULL);
+    }
+
+  /* Now close the streams which also flushes the output the user
+     defined handler might has produced.  */
+  if (stage == 4)
+    {
+      ++stage;
+      __fcloseall ();
+    }
+
+  /* Try again.  */
+  if (stage == 5)
+    {
+      ++stage;
+      raise (SIGABRT);
+    }
+
+  /* Now try to abort using the system specific command.  */
+  if (stage == 6)
+    {
+      ++stage;
+      ABORT_INSTRUCTION;
+    }
+
+  /* If we can't signal ourselves and the abort instruction failed, exit.  */
+  if (stage == 7)
+    {
+      ++stage;
+      _exit (127);
+    }
+
+  /* If even this fails try to use the provided instruction to crash
+     or otherwise make sure we never return.  */
+  while (1)
+    /* Try for ever and ever.  */
+    ABORT_INSTRUCTION;
+}
+libc_hidden_def (abort)
diff --git a/stdlib/add_n.c b/stdlib/add_n.c
new file mode 100644
index 0000000000..280e30545a
--- /dev/null
+++ b/stdlib/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 <gmp.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/stdlib/addmul_1.c b/stdlib/addmul_1.c
new file mode 100644
index 0000000000..6ae1e57ad9
--- /dev/null
+++ b/stdlib/addmul_1.c
@@ -0,0 +1,65 @@
+/* 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 <gmp.h>
+#include "gmp-impl.h"
+#include "longlong.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/stdlib/cmp.c b/stdlib/cmp.c
new file mode 100644
index 0000000000..e7661702b6
--- /dev/null
+++ b/stdlib/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 <gmp.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/stdlib/dbl2mpn.c b/stdlib/dbl2mpn.c
new file mode 100644
index 0000000000..4444467946
--- /dev/null
+++ b/stdlib/dbl2mpn.c
@@ -0,0 +1,32 @@
+/* Copyright (C) 1993, 1997 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 <gmp.h>
+#include "gmp-impl.h"
+
+/* Convert a `double' to a multi-precision integer representing the
+   significand scaled up by the highest possible number of significant bits
+   of fraction (DBL_MANT_DIG), and an integral power of two (MPN frexp). */
+
+mp_size_t
+__mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
+		      int *expt, int *is_neg,
+		      double value)
+{
+#error "__mpn_extract_double is not implemented for this floating point format"
+}
diff --git a/stdlib/div.c b/stdlib/div.c
new file mode 100644
index 0000000000..5268f4c494
--- /dev/null
+++ b/stdlib/div.c
@@ -0,0 +1,86 @@
+/* Copyright (C) 1992, 1997, 1999 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.  */
+
+/*
+ * Copyright (c) 1990 Regents of the University of California.
+ * All rights reserved.
+ *
+ * This code is derived from software contributed to Berkeley by
+ * Chris Torek.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 4. Neither the name of the University nor the names of its contributors
+ *    may be used to endorse or promote products derived from this software
+ *    without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include <stdlib.h>
+
+/* Return the `div_t' representation of NUMER over DENOM.  */
+div_t
+div (numer, denom)
+     int numer, denom;
+{
+  div_t result;
+
+  result.quot = numer / denom;
+  result.rem = numer % denom;
+
+  /* The ANSI standard says that |QUOT| <= |NUMER / DENOM|, where
+     NUMER / DENOM is to be computed in infinite precision.  In
+     other words, we should always truncate the quotient towards
+     zero, never -infinity.  Machine division and remainer may
+     work either way when one or both of NUMER or DENOM is
+     negative.  If only one is negative and QUOT has been
+     truncated towards -infinity, REM will have the same sign as
+     DENOM and the opposite sign of NUMER; if both are negative
+     and QUOT has been truncated towards -infinity, REM will be
+     positive (will have the opposite sign of NUMER).  These are
+     considered `wrong'.  If both are NUM and DENOM are positive,
+     RESULT will always be positive.  This all boils down to: if
+     NUMER >= 0, but REM < 0, we got the wrong answer.  In that
+     case, to get the right answer, add 1 to QUOT and subtract
+     DENOM from REM.  */
+
+  if (numer >= 0 && result.rem < 0)
+    {
+      ++result.quot;
+      result.rem -= denom;
+    }
+
+  return result;
+}
diff --git a/stdlib/divmod_1.c b/stdlib/divmod_1.c
new file mode 100644
index 0000000000..51a47d85d3
--- /dev/null
+++ b/stdlib/divmod_1.c
@@ -0,0 +1,208 @@
+/* mpn_divmod_1(quot_ptr, dividend_ptr, dividend_size, divisor_limb) --
+   Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
+   Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
+   Return the single-limb remainder.
+   There are no constraints on the value of the divisor.
+
+   QUOT_PTR and DIVIDEND_PTR might point to the same limb.
+
+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 <gmp.h>
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifndef UMUL_TIME
+#define UMUL_TIME 1
+#endif
+
+#ifndef UDIV_TIME
+#define UDIV_TIME UMUL_TIME
+#endif
+
+/* FIXME: We should be using invert_limb (or invert_normalized_limb)
+   here (not udiv_qrnnd).  */
+
+mp_limb_t
+#if __STDC__
+mpn_divmod_1 (mp_ptr quot_ptr,
+	      mp_srcptr dividend_ptr, mp_size_t dividend_size,
+	      mp_limb_t divisor_limb)
+#else
+mpn_divmod_1 (quot_ptr, dividend_ptr, dividend_size, divisor_limb)
+     mp_ptr quot_ptr;
+     mp_srcptr dividend_ptr;
+     mp_size_t dividend_size;
+     mp_limb_t divisor_limb;
+#endif
+{
+  mp_size_t i;
+  mp_limb_t n1, n0, r;
+  int dummy;
+
+  /* ??? Should this be handled at all?  Rely on callers?  */
+  if (dividend_size == 0)
+    return 0;
+
+  /* If multiplication is much faster than division, and the
+     dividend is large, pre-invert the divisor, and use
+     only multiplications in the inner loop.  */
+
+  /* This test should be read:
+       Does it ever help to use udiv_qrnnd_preinv?
+	 && Does what we save compensate for the inversion overhead?  */
+  if (UDIV_TIME > (2 * UMUL_TIME + 6)
+      && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME)
+    {
+      int normalization_steps;
+
+      count_leading_zeros (normalization_steps, divisor_limb);
+      if (normalization_steps != 0)
+	{
+	  mp_limb_t divisor_limb_inverted;
+
+	  divisor_limb <<= normalization_steps;
+
+	  /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
+	     result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
+	     most significant bit (with weight 2**N) implicit.  */
+
+	  /* Special case for DIVISOR_LIMB == 100...000.  */
+	  if (divisor_limb << 1 == 0)
+	    divisor_limb_inverted = ~(mp_limb_t) 0;
+	  else
+	    udiv_qrnnd (divisor_limb_inverted, dummy,
+			-divisor_limb, 0, divisor_limb);
+
+	  n1 = dividend_ptr[dividend_size - 1];
+	  r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
+
+	  /* Possible optimization:
+	     if (r == 0
+	     && divisor_limb > ((n1 << normalization_steps)
+			     | (dividend_ptr[dividend_size - 2] >> ...)))
+	     ...one division less... */
+
+	  for (i = dividend_size - 2; i >= 0; i--)
+	    {
+	      n0 = dividend_ptr[i];
+	      udiv_qrnnd_preinv (quot_ptr[i + 1], r, r,
+				 ((n1 << normalization_steps)
+				  | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
+				 divisor_limb, divisor_limb_inverted);
+	      n1 = n0;
+	    }
+	  udiv_qrnnd_preinv (quot_ptr[0], r, r,
+			     n1 << normalization_steps,
+			     divisor_limb, divisor_limb_inverted);
+	  return r >> normalization_steps;
+	}
+      else
+	{
+	  mp_limb_t divisor_limb_inverted;
+
+	  /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
+	     result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
+	     most significant bit (with weight 2**N) implicit.  */
+
+	  /* Special case for DIVISOR_LIMB == 100...000.  */
+	  if (divisor_limb << 1 == 0)
+	    divisor_limb_inverted = ~(mp_limb_t) 0;
+	  else
+	    udiv_qrnnd (divisor_limb_inverted, dummy,
+			-divisor_limb, 0, divisor_limb);
+
+	  i = dividend_size - 1;
+	  r = dividend_ptr[i];
+
+	  if (r >= divisor_limb)
+	    r = 0;
+	  else
+	    {
+	      quot_ptr[i] = 0;
+	      i--;
+	    }
+
+	  for (; i >= 0; i--)
+	    {
+	      n0 = dividend_ptr[i];
+	      udiv_qrnnd_preinv (quot_ptr[i], r, r,
+				 n0, divisor_limb, divisor_limb_inverted);
+	    }
+	  return r;
+	}
+    }
+  else
+    {
+      if (UDIV_NEEDS_NORMALIZATION)
+	{
+	  int normalization_steps;
+
+	  count_leading_zeros (normalization_steps, divisor_limb);
+	  if (normalization_steps != 0)
+	    {
+	      divisor_limb <<= normalization_steps;
+
+	      n1 = dividend_ptr[dividend_size - 1];
+	      r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
+
+	      /* Possible optimization:
+		 if (r == 0
+		 && divisor_limb > ((n1 << normalization_steps)
+				 | (dividend_ptr[dividend_size - 2] >> ...)))
+		 ...one division less... */
+
+	      for (i = dividend_size - 2; i >= 0; i--)
+		{
+		  n0 = dividend_ptr[i];
+		  udiv_qrnnd (quot_ptr[i + 1], r, r,
+			      ((n1 << normalization_steps)
+			       | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
+			      divisor_limb);
+		  n1 = n0;
+		}
+	      udiv_qrnnd (quot_ptr[0], r, r,
+			  n1 << normalization_steps,
+			  divisor_limb);
+	      return r >> normalization_steps;
+	    }
+	}
+      /* No normalization needed, either because udiv_qrnnd doesn't require
+	 it, or because DIVISOR_LIMB is already normalized.  */
+
+      i = dividend_size - 1;
+      r = dividend_ptr[i];
+
+      if (r >= divisor_limb)
+	r = 0;
+      else
+	{
+	  quot_ptr[i] = 0;
+	  i--;
+	}
+
+      for (; i >= 0; i--)
+	{
+	  n0 = dividend_ptr[i];
+	  udiv_qrnnd (quot_ptr[i], r, r, n0, divisor_limb);
+	}
+      return r;
+    }
+}
diff --git a/stdlib/divrem.c b/stdlib/divrem.c
new file mode 100644
index 0000000000..c97d01ec17
--- /dev/null
+++ b/stdlib/divrem.c
@@ -0,0 +1,245 @@
+/* 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 <gmp.h>
+#include "gmp-impl.h"
+#include "longlong.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/stdlib/strtod_l.c b/stdlib/strtod_l.c
index 3a1c1ebd4b..c7901c29fe 100644
--- a/stdlib/strtod_l.c
+++ b/stdlib/strtod_l.c
@@ -68,8 +68,8 @@ extern unsigned long long int ____strtoull_l_internal (const char *, char **,
    and _LONG_LONG_LIMB in it can take effect into gmp.h.  */
 #include <gmp-mparam.h>
 #include <gmp.h>
-#include <gmp-impl.h>
-#include <longlong.h>
+#include "gmp-impl.h"
+#include "longlong.h"
 #include "fpioconst.h"
 
 #define NDEBUG 1