PR fortran/46416

	* quadmath.h (cbrtq, finiteq, isnanq, signbitq, sqrtq): Remove
	const from prototype argument.
	(cimagq, conjq, cprojq, crealq, fdimq, fmaxq, fminq, ilogbq,
	llrintq, log2q, lrintq, nearbyintq, remquoq): New prototypes.
	(__quadmath_extern_inline): Define.
	(cimagq, conjq, crealq): New inlines.
	* Makefile.am (libquadmath_la_SOURCES): Add math/cimagq.c,
	math/conjq.c, math/cprojq.c, math/crealq.c, math/fdimq.c,
	math/fmaxq.c, math/fminq.c, math/ilogbq.c, math/llrintq.c,
	math/log2q.c, math/lrintq.c, math/nearbyintq.c and math/remquoq.c.
	* Makefile.in: Regenerated.
	* quadmath_weak.h (cimagq, conjq, cprojq, crealq, fdimq, fmaxq,
	fminq, ilogbq, llrintq, log2q, lrintq, nearbyintq, remquoq): Add.
	* quadmath-imp.h (__LITTLE_ENDIAN__): Don't define.
	(ieee854_float128): Use __BYTE_ORDER == __ORDER_BIG_ENDIAN__
	tests instead of __BIG_ENDIAN__.
	* quadmath.map (QUADMATH_1.0): Add cimagq, conjq, cprojq, crealq,
	fdimq, fmaxq, fminq, ilogbq, llrintq, log2q, lrintq, nearbyintq
	and remquoq.
	* libquadmath.texi (cimagq, conjq, cprojq, crealq, fdimq, fmaxq,
	fminq, ilogbq, llrintq, log2q, lrintq, nearbyintq, remquoq): Add.
	* math/cprojq.c: New file.
	* math/ilogbq.c: New file.
	* math/fminq.c: New file.
	* math/llrintq.c: New file.
	* math/log2q.c: New file.
	* math/lrintq.c: New file.
	* math/crealq.c: New file.
	* math/nearbyintq.c: New file.
	* math/fmaxq.c: New file.
	* math/conjq.c: New file.
	* math/remquoq.c: New file.
	* math/cimagq.c: New file.
	* math/fdimq.c: New file.

	* math/ldexpq.c: Include errno.h.  Set errno to ERANGE if needed.


git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@168854 138bc75d-0d04-0410-961f-82ee72b054a4

14 files changed, 895 insertions, 1 deletions

diff --git a/libquadmath/math/cimagq.c b/libquadmath/math/cimagq.c
new file mode 100644
index 00000000000..9b3bf70a10d
--- /dev/null
+++ b/libquadmath/math/cimagq.c
@@ -0,0 +1,27 @@
+/* Return imaginary part of complex __float128 value.
+   Copyright (C) 1997, 1998 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 "quadmath-imp.h"
+
+__float128
+cimagq (__complex128 z)
+{
+  return __imag__ z;
+}
diff --git a/libquadmath/math/conjq.c b/libquadmath/math/conjq.c
new file mode 100644
index 00000000000..8587d12193b
--- /dev/null
+++ b/libquadmath/math/conjq.c
@@ -0,0 +1,27 @@
+/* Return complex conjugate of complex __float128 value.
+   Copyright (C) 1997, 1998 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 "quadmath-imp.h"
+
+__complex128
+conjq (__complex128 z)
+{
+  return ~z;
+}
diff --git a/libquadmath/math/cprojq.c b/libquadmath/math/cprojq.c
new file mode 100644
index 00000000000..6092c732503
--- /dev/null
+++ b/libquadmath/math/cprojq.c
@@ -0,0 +1,40 @@
+/* Compute projection of complex __float128 value to Riemann sphere.
+   Copyright (C) 1997, 1999, 2010 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 "quadmath-imp.h"
+
+
+__complex128
+cprojq (__complex128 x)
+{
+  if (isnanq (__real__ x) && isnanq (__imag__ x))
+    return x;
+  else if (!finiteq (__real__ x) || !finiteq (__imag__ x))
+    {
+      __complex128 res;
+
+      __real__ res = __builtin_inf ();
+      __imag__ res = copysignq (0.0, __imag__ x);
+
+      return res;
+    }
+
+  return x;
+}
diff --git a/libquadmath/math/crealq.c b/libquadmath/math/crealq.c
new file mode 100644
index 00000000000..71f4a4405aa
--- /dev/null
+++ b/libquadmath/math/crealq.c
@@ -0,0 +1,27 @@
+/* Return real part of complex __float128 value.
+   Copyright (C) 1997, 1998 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 "quadmath-imp.h"
+
+__float128
+crealq (__complex128 z)
+{
+  return __real__ z;
+}
diff --git a/libquadmath/math/fdimq.c b/libquadmath/math/fdimq.c
new file mode 100644
index 00000000000..539fb08c641
--- /dev/null
+++ b/libquadmath/math/fdimq.c
@@ -0,0 +1,43 @@
+/* Return positive difference between arguments.
+   Copyright (C) 1997, 2004, 2009 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 <errno.h>
+#include "quadmath-imp.h"
+
+__float128
+fdimq (__float128 x, __float128 y)
+{
+  int clsx = fpclassifyq (x);
+  int clsy = fpclassifyq (y);
+
+  if (clsx == QUADFP_NAN || clsy == QUADFP_NAN
+      || (y < 0 && clsx == QUADFP_INFINITE && clsy == QUADFP_INFINITE))
+    /* Raise invalid flag.  */
+    return x - y;
+
+  if (x <= y)
+    return 0.0Q;
+
+  __float128 r = x - y;
+  if (isinfq (r))
+    errno = ERANGE;
+
+  return r;
+}
diff --git a/libquadmath/math/fmaxq.c b/libquadmath/math/fmaxq.c
new file mode 100644
index 00000000000..e8ed6f440aa
--- /dev/null
+++ b/libquadmath/math/fmaxq.c
@@ -0,0 +1,28 @@
+/* Return maximum numeric value of X and Y.
+   Copyright (C) 1997 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 "quadmath-imp.h"
+
+
+__float128
+fmaxq (__float128 x, __float128 y)
+{
+  return (__builtin_isgreaterequal (x, y) || isnanq (y)) ? x : y;
+}
diff --git a/libquadmath/math/fminq.c b/libquadmath/math/fminq.c
new file mode 100644
index 00000000000..2fbe4116a6d
--- /dev/null
+++ b/libquadmath/math/fminq.c
@@ -0,0 +1,28 @@
+/* Return minimum numeric value of X and Y.
+   Copyright (C) 1997 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 "quadmath-imp.h"
+
+
+__float128
+fminq (__float128 x, __float128 y)
+{
+  return (__builtin_islessequal (x, y) || isnanq (y)) ? x : y;
+}
diff --git a/libquadmath/math/ilogbq.c b/libquadmath/math/ilogbq.c
new file mode 100644
index 00000000000..47986f5b311
--- /dev/null
+++ b/libquadmath/math/ilogbq.c
@@ -0,0 +1,64 @@
+/* s_ilogbl.c -- long double version of s_ilogb.c.
+ * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: $";
+#endif
+
+/* ilogbl(long double x)
+ * return the binary exponent of non-zero x
+ * ilogbl(0) = FP_ILOGB0
+ * ilogbl(NaN) = FP_ILOGBNAN (no signal is raised)
+ * ilogbl(+-Inf) = INT_MAX (no signal is raised)
+ */
+
+#include <limits.h>
+#include <math.h>
+#include "quadmath-imp.h"
+
+#ifndef FP_ILOGB0
+# define FP_ILOGB0 INT_MIN
+#endif
+#ifndef FP_ILOGBNAN
+# define FP_ILOGBNAN INT_MAX
+#endif
+
+int
+ilogbq (__float128 x)
+{
+	int64_t hx,lx;
+	int ix;
+
+	GET_FLT128_WORDS64(hx,lx,x);
+	hx &= 0x7fffffffffffffffLL;
+	if(hx <= 0x0001000000000000LL) {
+	    if((hx|lx)==0)
+		return FP_ILOGB0;	/* ilogbl(0) = FP_ILOGB0 */
+	    else			/* subnormal x */
+		if(hx==0) {
+		    for (ix = -16431; lx>0; lx<<=1) ix -=1;
+		} else {
+		    for (ix = -16382, hx<<=15; hx>0; hx<<=1) ix -=1;
+		}
+	    return ix;
+	}
+	else if (hx<0x7fff000000000000LL) return (hx>>48)-0x3fff;
+	else if (FP_ILOGBNAN != INT_MAX) {
+	    /* ISO C99 requires ilogbl(+-Inf) == INT_MAX.  */
+	    if (((hx^0x7fff000000000000LL)|lx) == 0)
+		return INT_MAX;
+	}
+	return FP_ILOGBNAN;
+}
diff --git a/libquadmath/math/ldexpq.c b/libquadmath/math/ldexpq.c
index 73f52e509bf..394d4590cca 100644
--- a/libquadmath/math/ldexpq.c
+++ b/libquadmath/math/ldexpq.c
@@ -14,6 +14,7 @@
  * ====================================================
  */
 
+#include <errno.h>
 #include "quadmath-imp.h"
 
 __float128
@@ -21,6 +22,6 @@ ldexpq (__float128 value, int exp)
 {
   if(!finiteq(value)||value==0.0Q) return value;
   value = scalbnq(value,exp);
-  /* if(!__finitel(value)||value==0.0Q) __set_errno (ERANGE); */
+  if(!finiteq(value)||value==0.0Q) errno = ERANGE;
   return value;
 }
diff --git a/libquadmath/math/llrintq.c b/libquadmath/math/llrintq.c
new file mode 100644
index 00000000000..eef31d823b6
--- /dev/null
+++ b/libquadmath/math/llrintq.c
@@ -0,0 +1,71 @@
+/* Round argument to nearest integral value according to current rounding
+   direction.
+   Copyright (C) 1997, 1999, 2006 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
+   		  Jakub Jelinek <jj@ultra.linux.cz>, 1999.
+
+   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 "quadmath-imp.h"
+
+static const __float128 two112[2] =
+{
+  5.19229685853482762853049632922009600E+33Q, /* 0x406F000000000000, 0 */
+ -5.19229685853482762853049632922009600E+33Q  /* 0xC06F000000000000, 0 */
+};
+
+long long int
+llrintq (__float128 x)
+{
+  int32_t j0;
+  uint64_t i0,i1;
+  volatile __float128 w;
+  __float128 t;
+  long long int result;
+  int sx;
+
+  GET_FLT128_WORDS64 (i0, i1, x);
+  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+  sx = i0 >> 63;
+  i0 &= 0x0000ffffffffffffLL;
+  i0 |= 0x0001000000000000LL;
+
+  if (j0 < (int32_t) (8 * sizeof (long long int)) - 1)
+    {
+      w = two112[sx] + x;
+      t = w - two112[sx];
+      GET_FLT128_WORDS64 (i0, i1, t);
+      j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+      i0 &= 0x0000ffffffffffffLL;
+      i0 |= 0x0001000000000000LL;
+
+      if (j0 < 0)
+	result = 0;
+      else if (j0 <= 48)
+	result = i0 >> (48 - j0);
+      else
+	result = ((long long int) i0 << (j0 - 48)) | (i1 >> (112 - j0));
+    }
+  else
+    {
+      /* The number is too large.  It is left implementation defined
+	 what happens.  */
+      return (long long int) x;
+    }
+
+  return sx ? -result : result;
+}
diff --git a/libquadmath/math/log2q.c b/libquadmath/math/log2q.c
new file mode 100644
index 00000000000..963b38c8483
--- /dev/null
+++ b/libquadmath/math/log2q.c
@@ -0,0 +1,248 @@
+/*                                                      log2l.c
+ *      Base 2 logarithm, 128-bit long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, log2l();
+ *
+ * y = log2l( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base 2 logarithm of x.
+ *
+ * The argument is separated into its exponent and fractional
+ * parts.  If the exponent is between -1 and +1, the (natural)
+ * logarithm of the fraction is approximated by
+ *
+ *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
+ *
+ * Otherwise, setting  z = 2(x-1)/x+1),
+ *
+ *     log(x) = z + z^3 P(z)/Q(z).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0.5, 2.0     100,000    2.6e-34     4.9e-35
+ *    IEEE     exp(+-10000)  100,000    9.6e-35     4.0e-35
+ *
+ * In the tests over the interval exp(+-10000), the logarithms
+ * of the random arguments were uniformly distributed over
+ * [-10000, +10000].
+ *
+ */
+
+/*
+   Cephes Math Library Release 2.2:  January, 1991
+   Copyright 1984, 1991 by Stephen L. Moshier
+   Adapted for glibc November, 2001
+
+    This 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.
+
+    This 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 this library; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA 
+ */
+
+#include "quadmath-imp.h"
+
+/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 5.3e-37,
+ * relative peak error spread = 2.3e-14
+ */
+static const __float128 P[13] =
+{
+  1.313572404063446165910279910527789794488E4Q,
+  7.771154681358524243729929227226708890930E4Q,
+  2.014652742082537582487669938141683759923E5Q,
+  3.007007295140399532324943111654767187848E5Q,
+  2.854829159639697837788887080758954924001E5Q,
+  1.797628303815655343403735250238293741397E5Q,
+  7.594356839258970405033155585486712125861E4Q,
+  2.128857716871515081352991964243375186031E4Q,
+  3.824952356185897735160588078446136783779E3Q,
+  4.114517881637811823002128927449878962058E2Q,
+  2.321125933898420063925789532045674660756E1Q,
+  4.998469661968096229986658302195402690910E-1Q,
+  1.538612243596254322971797716843006400388E-6Q
+};
+static const __float128 Q[12] =
+{
+  3.940717212190338497730839731583397586124E4Q,
+  2.626900195321832660448791748036714883242E5Q,
+  7.777690340007566932935753241556479363645E5Q,
+  1.347518538384329112529391120390701166528E6Q,
+  1.514882452993549494932585972882995548426E6Q,
+  1.158019977462989115839826904108208787040E6Q,
+  6.132189329546557743179177159925690841200E5Q,
+  2.248234257620569139969141618556349415120E5Q,
+  5.605842085972455027590989944010492125825E4Q,
+  9.147150349299596453976674231612674085381E3Q,
+  9.104928120962988414618126155557301584078E2Q,
+  4.839208193348159620282142911143429644326E1Q
+/* 1.000000000000000000000000000000000000000E0Q, */
+};
+
+/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
+ * where z = 2(x-1)/(x+1)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 1.1e-35,
+ * relative peak error spread 1.1e-9
+ */
+static const __float128 R[6] =
+{
+  1.418134209872192732479751274970992665513E5Q,
+ -8.977257995689735303686582344659576526998E4Q,
+  2.048819892795278657810231591630928516206E4Q,
+ -2.024301798136027039250415126250455056397E3Q,
+  8.057002716646055371965756206836056074715E1Q,
+ -8.828896441624934385266096344596648080902E-1Q
+};
+static const __float128 S[6] =
+{
+  1.701761051846631278975701529965589676574E6Q,
+ -1.332535117259762928288745111081235577029E6Q,
+  4.001557694070773974936904547424676279307E5Q,
+ -5.748542087379434595104154610899551484314E4Q,
+  3.998526750980007367835804959888064681098E3Q,
+ -1.186359407982897997337150403816839480438E2Q
+/* 1.000000000000000000000000000000000000000E0Q, */
+};
+
+static const __float128
+/* log2(e) - 1 */
+LOG2EA = 4.4269504088896340735992468100189213742664595E-1Q,
+/* sqrt(2)/2 */
+SQRTH = 7.071067811865475244008443621048490392848359E-1Q;
+
+
+/* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
+
+static __float128
+neval (__float128 x, const __float128 *p, int n)
+{
+  __float128 y;
+
+  p += n;
+  y = *p--;
+  do
+    {
+      y = y * x + *p--;
+    }
+  while (--n > 0);
+  return y;
+}
+
+
+/* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
+
+static __float128
+deval (__float128 x, const __float128 *p, int n)
+{
+  __float128 y;
+
+  p += n;
+  y = x + *p--;
+  do
+    {
+      y = y * x + *p--;
+    }
+  while (--n > 0);
+  return y;
+}
+
+
+
+__float128
+log2q (__float128 x)
+{
+  __float128 z;
+  __float128 y;
+  int e;
+  int64_t hx, lx;
+
+/* Test for domain */
+  GET_FLT128_WORDS64 (hx, lx, x);
+  if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
+    return (-1.0Q / (x - x));
+  if (hx < 0)
+    return (x - x) / (x - x);
+  if (hx >= 0x7fff000000000000LL)
+    return (x + x);
+
+/* separate mantissa from exponent */
+
+/* Note, frexp is used so that denormal numbers
+ * will be handled properly.
+ */
+  x = frexpq (x, &e);
+
+
+/* logarithm using log(x) = z + z**3 P(z)/Q(z),
+ * where z = 2(x-1)/x+1)
+ */
+  if ((e > 2) || (e < -2))
+    {
+      if (x < SQRTH)
+	{			/* 2( 2x-1 )/( 2x+1 ) */
+	  e -= 1;
+	  z = x - 0.5Q;
+	  y = 0.5Q * z + 0.5Q;
+	}
+      else
+	{			/*  2 (x-1)/(x+1)   */
+	  z = x - 0.5Q;
+	  z -= 0.5Q;
+	  y = 0.5Q * x + 0.5Q;
+	}
+      x = z / y;
+      z = x * x;
+      y = x * (z * neval (z, R, 5) / deval (z, S, 5));
+      goto done;
+    }
+
+
+/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
+
+  if (x < SQRTH)
+    {
+      e -= 1;
+      x = 2.0 * x - 1.0Q;	/*  2x - 1  */
+    }
+  else
+    {
+      x = x - 1.0Q;
+    }
+  z = x * x;
+  y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
+  y = y - 0.5 * z;
+
+done:
+
+/* Multiply log of fraction by log2(e)
+ * and base 2 exponent by 1
+ */
+  z = y * LOG2EA;
+  z += x * LOG2EA;
+  z += y;
+  z += x;
+  z += e;
+  return (z);
+}
diff --git a/libquadmath/math/lrintq.c b/libquadmath/math/lrintq.c
new file mode 100644
index 00000000000..d1497ae3882
--- /dev/null
+++ b/libquadmath/math/lrintq.c
@@ -0,0 +1,85 @@
+/* Round argument to nearest integral value according to current rounding
+   direction.
+   Copyright (C) 1997, 1999, 2004, 2006 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
+   		  Jakub Jelinek <jj@ultra.linux.cz>, 1999.
+
+   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 "quadmath-imp.h"
+
+static const __float128 two112[2] =
+{
+  5.19229685853482762853049632922009600E+33Q, /* 0x406F000000000000, 0 */
+ -5.19229685853482762853049632922009600E+33Q  /* 0xC06F000000000000, 0 */
+};
+
+long int
+lrintq (__float128 x)
+{
+  int32_t j0;
+  uint64_t i0,i1;
+  volatile __float128 w;
+  __float128 t;
+  long int result;
+  int sx;
+
+  GET_FLT128_WORDS64 (i0, i1, x);
+  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+  sx = i0 >> 63;
+  i0 &= 0x0000ffffffffffffLL;
+  i0 |= 0x0001000000000000LL;
+
+  if (j0 < 48)
+    {
+      w = two112[sx] + x;
+      t = w - two112[sx];
+      GET_FLT128_WORDS64 (i0, i1, t);
+      j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+      i0 &= 0x0000ffffffffffffLL;
+      i0 |= 0x0001000000000000LL;
+
+      result = (j0 < 0 ? 0 : i0 >> (48 - j0));
+    }
+  else if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
+    {
+      if (j0 >= 112)
+	result = ((long int) i0 << (j0 - 48)) | (i1 << (j0 - 112));
+      else
+	{
+	  w = two112[sx] + x;
+	  t = w - two112[sx];
+	  GET_FLT128_WORDS64 (i0, i1, t);
+	  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+	  i0 &= 0x0000ffffffffffffLL;
+	  i0 |= 0x0001000000000000LL;
+
+	  if (j0 == 48)
+	    result = (long int) i0;
+	  else
+	    result = ((long int) i0 << (j0 - 48)) | (i1 >> (112 - j0));
+	}
+    }
+  else
+    {
+      /* The number is too large.  It is left implementation defined
+	 what happens.  */
+      return (long int) x;
+    }
+
+  return sx ? -result : result;
+}
diff --git a/libquadmath/math/nearbyintq.c b/libquadmath/math/nearbyintq.c
new file mode 100644
index 00000000000..8e92c5afd4b
--- /dev/null
+++ b/libquadmath/math/nearbyintq.c
@@ -0,0 +1,98 @@
+/* nearbyintq.c -- __float128 version of s_nearbyint.c.
+ * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * nearbyintq(x)
+ * Return x rounded to integral value according to the prevailing
+ * rounding mode.
+ * Method:
+ *	Using floating addition.
+ * Exception:
+ *	Inexact flag raised if x not equal to rintq(x).
+ */
+
+#include "quadmath-imp.h"
+#ifdef HAVE_FENV_H
+# include <fenv.h>
+# if defined HAVE_FEHOLDEXCEPT && defined HAVE_FESETENV
+#   define USE_FENV_H
+# endif
+#endif
+
+static const __float128
+TWO112[2]={
+  5.19229685853482762853049632922009600E+33Q, /* 0x406F000000000000, 0 */
+ -5.19229685853482762853049632922009600E+33Q  /* 0xC06F000000000000, 0 */
+};
+
+__float128
+nearbyintq(__float128 x)
+{
+#ifdef USE_FENV_H
+	fenv_t env;
+#endif
+	int64_t i0,j0,sx;
+	uint64_t i,i1;
+	__float128 w,t;
+	GET_FLT128_WORDS64(i0,i1,x);
+	sx = (((uint64_t)i0)>>63);
+	j0 = ((i0>>48)&0x7fff)-0x3fff;
+	if(j0<48) {
+	    if(j0<0) {
+		if(((i0&0x7fffffffffffffffLL)|i1)==0) return x;
+		i1 |= (i0&0x0000ffffffffffffLL);
+		i0 &= 0xffffe00000000000ULL;
+		i0 |= ((i1|-i1)>>16)&0x0000800000000000LL;
+		SET_FLT128_MSW64(x,i0);
+#ifdef USE_FENV_H
+		feholdexcept (&env);
+#endif
+	        w = TWO112[sx]+x;
+	        t = w-TWO112[sx];
+#ifdef USE_FENV_H
+	        fesetenv (&env);
+#endif
+		GET_FLT128_MSW64(i0,t);
+		SET_FLT128_MSW64(t,(i0&0x7fffffffffffffffLL)|(sx<<63));
+	        return t;
+	    } else {
+		i = (0x0000ffffffffffffLL)>>j0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		i>>=1;
+		if(((i0&i)|i1)!=0) {
+		    if(j0==47) i1 = 0x4000000000000000ULL; else
+		    i0 = (i0&(~i))|((0x0000200000000000LL)>>j0);
+		}
+	    }
+	} else if (j0>111) {
+	    if(j0==0x4000) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	} else {
+	    i = -1ULL>>(j0-48);
+	    if((i1&i)==0) return x;	/* x is integral */
+	    i>>=1;
+	    if((i1&i)!=0) i1 = (i1&(~i))|((0x4000000000000000LL)>>(j0-48));
+	}
+	SET_FLT128_WORDS64(x,i0,i1);
+#ifdef USE_FENV_H
+	feholdexcept (&env);
+#endif
+	w = TWO112[sx]+x;
+	t = w-TWO112[sx];
+#ifdef USE_FENV_H	
+	fesetenv (&env);
+#endif
+	return t;
+}
diff --git a/libquadmath/math/remquoq.c b/libquadmath/math/remquoq.c
new file mode 100644
index 00000000000..3e3b4f68ce4
--- /dev/null
+++ b/libquadmath/math/remquoq.c
@@ -0,0 +1,107 @@
+/* Compute remainder and a congruent to the quotient.
+   Copyright (C) 1997, 1999, 2002 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
+		  Jakub Jelinek <jj@ultra.linux.cz>, 1999.
+
+   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 "quadmath-imp.h"
+
+
+static const __float128 zero = 0.0;
+
+
+__float128
+remquoq (__float128 x, __float128 y, int *quo)
+{
+  int64_t hx,hy;
+  u_int64_t sx,lx,ly,qs;
+  int cquo;
+
+  GET_FLT128_WORDS64 (hx, lx, x);
+  GET_FLT128_WORDS64 (hy, ly, y);
+  sx = hx & 0x8000000000000000ULL;
+  qs = sx ^ (hy & 0x8000000000000000ULL);
+  hy &= 0x7fffffffffffffffLL;
+  hx &= 0x7fffffffffffffffLL;
+
+  /* Purge off exception values.  */
+  if ((hy | ly) == 0)
+    return (x * y) / (x * y); 			/* y = 0 */
+  if ((hx >= 0x7fff000000000000LL)		/* x not finite */
+      || ((hy >= 0x7fff000000000000LL)		/* y is NaN */
+	  && (((hy - 0x7fff000000000000LL) | ly) != 0)))
+    return (x * y) / (x * y);
+
+  if (hy <= 0x7ffbffffffffffffLL)
+    x = fmodq (x, 8 * y);              /* now x < 8y */
+      
+  if (((hx - hy) | (lx - ly)) == 0)
+    {
+      *quo = qs ? -1 : 1;
+      return zero * x;
+    }
+
+  x  = fabsq (x);
+  y  = fabsq (y);
+  cquo = 0;
+
+  if (x >= 4 * y)
+    {
+      x -= 4 * y;
+      cquo += 4;
+    }
+  if (x >= 2 * y)
+    {
+      x -= 2 * y;
+      cquo += 2;
+    }
+
+  if (hy < 0x0002000000000000LL)
+    {
+      if (x + x > y)
+	{
+	  x -= y;
+	  ++cquo;
+	  if (x + x >= y)
+	    {
+	      x -= y;
+	      ++cquo;
+	    }
+	}
+    }
+  else
+    {
+      __float128 y_half = 0.5Q * y;
+      if (x > y_half)
+	{
+	  x -= y;
+	  ++cquo;
+	  if (x >= y_half)
+	    {
+	      x -= y;
+	      ++cquo;
+	    }
+	}
+    }
+
+  *quo = qs ? -cquo : cquo;
+
+  if (sx)
+    x = -x;
+  return x;
+}