Update.

42 files changed, 3202 insertions, 0 deletions

diff --git a/sysdeps/ieee754/ldbl-128/e_acoshl.c b/sysdeps/ieee754/ldbl-128/e_acoshl.c
new file mode 100644
index 0000000..7f79340
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/e_acoshl.c
@@ -0,0 +1,68 @@
+/* e_acoshl.c -- long double version of e_acosh.c.
+ * Conversion to 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.
+ * ====================================================
+ */
+
+/* __ieee754_acoshl(x)
+ * Method :
+ *	Based on
+ *		acoshl(x) = logl [ x + sqrtl(x*x-1) ]
+ *	we have
+ *		acoshl(x) := logl(x)+ln2,	if x is large; else
+ *		acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
+ *		acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
+ *
+ * Special cases:
+ *	acoshl(x) is NaN with signal if x<1.
+ *	acoshl(NaN) is NaN without signal.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double
+#else
+static long double
+#endif
+one	= 1.0,
+ln2	= 0.6931471805599453094172321214581766L;
+
+#ifdef __STDC__
+	long double __ieee754_acoshl(long double x)
+#else
+	long double __ieee754_acoshl(x)
+	long double x;
+#endif
+{
+	long double t;
+	u_int64_t lx;
+	int64_t hx;
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	if(hx<0x3fff000000000000LL) {		/* x < 1 */
+	    return (x-x)/(x-x);
+	} else if(hx >=0x401b000000000000LL) {	/* x > 2**28 */
+	    if(hx >=0x7fff000000000000LL) {	/* x is inf of NaN */
+	        return x+x;
+	    } else
+		return __ieee754_logl(x)+ln2;	/* acoshl(huge)=logl(2x) */
+	} else if(((hx-0x3fff000000000000LL)|lx)==0) {
+	    return 0.0L;			/* acosh(1) = 0 */
+	} else if (hx > 0x4000000000000000LL) {	/* 2**28 > x > 2 */
+	    t=x*x;
+	    return __ieee754_logl(2.0L*x-one/(x+__ieee754_sqrtl(t-one)));
+	} else {			/* 1<x<2 */
+	    t = x-one;
+	    return __log1pl(t+__sqrtl(2.0L*t+t*t));
+	}
+}
diff --git a/sysdeps/ieee754/ldbl-128/e_atan2l.c b/sysdeps/ieee754/ldbl-128/e_atan2l.c
new file mode 100644
index 0000000..2e081d3
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/e_atan2l.c
@@ -0,0 +1,129 @@
+/* e_atan2l.c -- long double version of e_atan2.c.
+ * Conversion to 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.
+ * ====================================================
+ */
+
+/* __ieee754_atan2l(y,x)
+ * Method :
+ *	1. Reduce y to positive by atan2l(y,x)=-atan2l(-y,x).
+ *	2. Reduce x to positive by (if x and y are unexceptional):
+ *		ARG (x+iy) = arctan(y/x)   	   ... if x > 0,
+ *		ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
+ *
+ * Special cases:
+ *
+ *	ATAN2((anything), NaN ) is NaN;
+ *	ATAN2(NAN , (anything) ) is NaN;
+ *	ATAN2(+-0, +(anything but NaN)) is +-0  ;
+ *	ATAN2(+-0, -(anything but NaN)) is +-pi ;
+ *	ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
+ *	ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
+ *	ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
+ *	ATAN2(+-INF,+INF ) is +-pi/4 ;
+ *	ATAN2(+-INF,-INF ) is +-3pi/4;
+ *	ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double
+#else
+static long double
+#endif
+tiny  = 1.0e-4900L,
+zero  = 0.0,
+pi_o_4  = 7.85398163397448309615660845819875699e-01L, /* 3ffe921fb54442d18469898cc51701b8 */
+pi_o_2  = 1.57079632679489661923132169163975140e+00L, /* 3fff921fb54442d18469898cc51701b8 */
+pi      = 3.14159265358979323846264338327950280e+00L, /* 4000921fb54442d18469898cc51701b8 */
+pi_lo   = 8.67181013012378102479704402604335225e-35L; /* 3f8dcd129024e088a67cc74020bbea64 */
+
+#ifdef __STDC__
+	long double __ieee754_atan2l(long double y, long double x)
+#else
+	long double __ieee754_atan2l(y,x)
+	long double  y,x;
+#endif
+{
+	long double z;
+	int64_t k,m,hx,hy,ix,iy;
+	u_int64_t lx,ly;
+
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	ix = hx&0x7fffffffffffffffLL;
+	GET_LDOUBLE_WORDS64(hy,ly,y);
+	iy = hy&0x7fffffffffffffffLL;
+	if(((ix|((lx|-lx)>>63))>0x7fff000000000000LL)||
+	   ((iy|((ly|-ly)>>63))>0x7fff000000000000LL))	/* x or y is NaN */
+	   return x+y;
+	if((hx-0x3fff000000000000LL|lx)==0) return __atanl(y);   /* x=1.0L */
+	m = ((hy>>63)&1)|((hx>>62)&2);	/* 2*sign(x)+sign(y) */
+
+    /* when y = 0 */
+	if((iy|ly)==0) {
+	    switch(m) {
+		case 0:
+		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
+		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
+		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
+	    }
+	}
+    /* when x = 0 */
+	if((ix|lx)==0) return (hy<0)?  -pi_o_2-tiny: pi_o_2+tiny;
+
+    /* when x is INF */
+	if(ix==0x7fff000000000000LL) {
+	    if(iy==0x7fff000000000000LL) {
+		switch(m) {
+		    case 0: return  pi_o_4+tiny;/* atan(+INF,+INF) */
+		    case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
+		    case 2: return  3.0L*pi_o_4+tiny;/*atan(+INF,-INF)*/
+		    case 3: return -3.0L*pi_o_4-tiny;/*atan(-INF,-INF)*/
+		}
+	    } else {
+		switch(m) {
+		    case 0: return  zero  ;	/* atan(+...,+INF) */
+		    case 1: return -zero  ;	/* atan(-...,+INF) */
+		    case 2: return  pi+tiny  ;	/* atan(+...,-INF) */
+		    case 3: return -pi-tiny  ;	/* atan(-...,-INF) */
+		}
+	    }
+	}
+    /* when y is INF */
+	if(iy==0x7fff000000000000LL) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
+
+    /* compute y/x */
+	k = (iy-ix)>>48;
+	if(k > 120) z=pi_o_2+0.5L*pi_lo; 	/* |y/x| >  2**120 */
+	else if(hx<0&&k<-120) z=0.0L; 		/* |y|/x < -2**120 */
+	else z=__atanl(fabsl(y/x));		/* safe to do y/x */
+	switch (m) {
+	    case 0: return       z  ;	/* atan(+,+) */
+	    case 1: {
+	    	      u_int64_t zh;
+		      GET_LDOUBLE_MSW64(zh,z);
+		      SET_LDOUBLE_MSW64(z,zh ^ 0x8000000000000000ULL);
+		    }
+		    return       z  ;	/* atan(-,+) */
+	    case 2: return  pi-(z-pi_lo);/* atan(+,-) */
+	    default: /* case 3 */
+	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
+	}
+}
diff --git a/sysdeps/ieee754/ldbl-128/e_fmodl.c b/sysdeps/ieee754/ldbl-128/e_fmodl.c
new file mode 100644
index 0000000..1043f69
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/e_fmodl.c
@@ -0,0 +1,138 @@
+/* e_fmodl.c -- long double version of e_fmod.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.
+ * ====================================================
+ */
+
+/* 
+ * __ieee754_fmodl(x,y)
+ * Return x mod y in exact arithmetic
+ * Method: shift and subtract
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double one = 1.0, Zero[] = {0.0, -0.0,};
+#else
+static long double one = 1.0, Zero[] = {0.0, -0.0,};
+#endif
+
+#ifdef __STDC__
+	long double __ieee754_fmodl(long double x, long double y)
+#else
+	long double __ieee754_fmodl(x,y)
+	long double x,y;
+#endif
+{
+	int64_t n,hx,hy,hz,ix,iy,sx,i;
+	u_int64_t lx,ly,lz;
+
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	GET_LDOUBLE_WORDS64(hy,ly,y);
+	sx = hx&0x8000000000000000ULL;		/* sign of x */
+	hx ^=sx;				/* |x| */
+	hy &= 0x7fffffffffffffffLL;		/* |y| */
+
+    /* purge off exception values */
+	if((hy|ly)==0||(hx>=0x7fff000000000000LL)|| /* y=0,or x not finite */
+	  ((hy|((ly|-ly)>>63))>0x7fff000000000000LL))	/* or y is NaN */
+	    return (x*y)/(x*y);
+	if(hx<=hy) {
+	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
+	    if(lx==ly) 
+		return Zero[(u_int64_t)sx>>63];	/* |x|=|y| return x*0*/
+	}
+
+    /* determine ix = ilogb(x) */
+	if(hx<0x0001000000000000LL) {	/* subnormal x */
+	    if(hx==0) {
+		for (ix = -16431, i=lx; i>0; i<<=1) ix -=1;
+	    } else {
+		for (ix = -16382, i=hx<<15; i>0; i<<=1) ix -=1;
+	    }
+	} else ix = (hx>>48)-0x3fff;
+
+    /* determine iy = ilogb(y) */
+	if(hy<0x0001000000000000LL) {	/* subnormal y */
+	    if(hy==0) {
+		for (iy = -16431, i=ly; i>0; i<<=1) iy -=1;
+	    } else {
+		for (iy = -16382, i=hy<<15; i>0; i<<=1) iy -=1;
+	    }
+	} else iy = (hy>>48)-0x3fff;
+
+    /* set up {hx,lx}, {hy,ly} and align y to x */
+	if(ix >= -16382) 
+	    hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
+	else {		/* subnormal x, shift x to normal */
+	    n = -16382-ix;
+	    if(n<=63) {
+	        hx = (hx<<n)|(lx>>(64-n));
+	        lx <<= n;
+	    } else {
+		hx = lx<<(n-64);
+		lx = 0;
+	    }
+	}
+	if(iy >= -16382) 
+	    hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
+	else {		/* subnormal y, shift y to normal */
+	    n = -16382-iy;
+	    if(n<=63) {
+	        hy = (hy<<n)|(ly>>(64-n));
+	        ly <<= n;
+	    } else {
+		hy = ly<<(n-64);
+		ly = 0;
+	    }
+	}
+
+    /* fix point fmod */
+	n = ix - iy;
+	while(n--) {
+	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	    if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
+	    else {
+	    	if((hz|lz)==0) 		/* return sign(x)*0 */
+		    return Zero[(u_int64_t)sx>>63];
+	    	hx = hz+hz+(lz>>63); lx = lz+lz;
+	    }
+	}
+	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	if(hz>=0) {hx=hz;lx=lz;}
+
+    /* convert back to floating value and restore the sign */
+	if((hx|lx)==0) 			/* return sign(x)*0 */
+	    return Zero[(u_int64_t)sx>>63];	
+	while(hx<0x0001000000000000LL) {	/* normalize x */
+	    hx = hx+hx+(lx>>63); lx = lx+lx;
+	    iy -= 1;
+	}
+	if(iy>= -16382) {	/* normalize output */
+	    hx = ((hx-0x0001000000000000LL)|((iy+16383)<<48));
+	    SET_LDOUBLE_WORDS64(x,hx|sx,lx);
+	} else {		/* subnormal output */
+	    n = -16382 - iy;
+	    if(n<=48) {
+		lx = (lx>>n)|((u_int64_t)hx<<(64-n));
+		hx >>= n;
+	    } else if (n<=63) {
+		lx = (hx<<(64-n))|(lx>>n); hx = sx;
+	    } else {
+		lx = hx>>(n-64); hx = sx;
+	    }
+	    SET_LDOUBLE_WORDS64(x,hx|sx,lx);
+	    x *= one;		/* create necessary signal */
+	}
+	return x;		/* exact output */
+}
diff --git a/sysdeps/ieee754/ldbl-128/e_gammal_r.c b/sysdeps/ieee754/ldbl-128/e_gammal_r.c
new file mode 100644
index 0000000..f77350f
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/e_gammal_r.c
@@ -0,0 +1,52 @@
+/* Implementation of gamma function according to ISO C.
+   Copyright (C) 1997, 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+#include <math_private.h>
+
+
+long double
+__ieee754_gammal_r (long double x, int *signgamp)
+{
+  /* We don't have a real gamma implementation now.  We'll use lgamma
+     and the exp function.  But due to the required boundary
+     conditions we must check some values separately.  */
+  int64_t hx;
+  u_int64_t lx;
+
+  GET_LDOUBLE_WORDS64 (hx, lx, x);
+
+  if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
+    {
+      /* Return value for x == 0 is NaN with invalid exception.  */
+      *signgamp = 0;
+      return x / x;
+    }
+  if (hx < 0 && (u_int64_t) hx < 0xffff000000000000ULL && __rintl (x) == x)
+    {
+      /* Return value for integer x < 0 is NaN with invalid exception.  */
+      *signgamp = 0;
+      return (x - x) / (x - x);
+    }
+
+  /* XXX FIXME.  */
+  return __ieee754_expl (__ieee754_lgammal_r (x, signgamp));
+}
diff --git a/sysdeps/ieee754/ldbl-128/e_remainderl.c b/sysdeps/ieee754/ldbl-128/e_remainderl.c
new file mode 100644
index 0000000..81af247
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/e_remainderl.c
@@ -0,0 +1,78 @@
+/* e_fmodl.c -- long double version of e_fmod.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.
+ * ====================================================
+ */
+
+/* __ieee754_remainderl(x,p)
+ * Return :                  
+ * 	returns  x REM p  =  x - [x/p]*p as if in infinite 
+ * 	precise arithmetic, where [x/p] is the (infinite bit) 
+ *	integer nearest x/p (in half way case choose the even one).
+ * Method : 
+ *	Based on fmodl() return x-[x/p]chopped*p exactlp.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double zero = 0.0L;
+#else
+static long double zero = 0.0L;
+#endif
+
+
+#ifdef __STDC__
+	long double __ieee754_remainderl(long double x, long double p)
+#else
+	long double __ieee754_remainderl(x,p)
+	long double x,p;
+#endif
+{
+	int64_t hx,hp;
+	u_int64_t sx,lx,lp;
+	long double p_half;
+
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	GET_LDOUBLE_WORDS64(hp,lp,p);
+	sx = hx&0x8000000000000000ULL;
+	hp &= 0x7fffffffffffffffLL;
+	hx &= 0x7fffffffffffffffLL;
+
+    /* purge off exception values */
+	if((hp|lp)==0) return (x*p)/(x*p); 	/* p = 0 */
+	if((hx>=0x7fff000000000000LL)||			/* x not finite */
+	  ((hp>=0x7fff000000000000LL)&&			/* p is NaN */
+	  (((hp-0x7fff000000000000LL)|lp)!=0)))
+	    return (x*p)/(x*p);
+
+
+	if (hp<=0x7ffdffffffffffffLL) x = __ieee754_fmodl(x,p+p);	/* now x < 2p */
+	if (((hx-hp)|(lx-lp))==0) return zero*x;
+	x  = fabsl(x);
+	p  = fabsl(p);
+	if (hp<0x0002000000000000LL) {
+	    if(x+x>p) {
+		x-=p;
+		if(x+x>=p) x -= p;
+	    }
+	} else {
+	    p_half = 0.5L*p;
+	    if(x>p_half) {
+		x-=p;
+		if(x>=p_half) x -= p;
+	    }
+	}
+	GET_LDOUBLE_MSW64(hx,x);
+	SET_LDOUBLE_MSW64(x,hx^sx);
+	return x;
+}
diff --git a/sysdeps/ieee754/ldbl-128/ieee754.h b/sysdeps/ieee754/ldbl-128/ieee754.h
new file mode 100644
index 0000000..c3b2f38
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/ieee754.h
@@ -0,0 +1,171 @@
+/* Copyright (C) 1992, 1995, 1996, 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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.  */
+
+#ifndef _IEEE754_H
+
+#define _IEEE754_H 1
+#include <features.h>
+
+#include <endian.h>
+
+__BEGIN_DECLS
+
+union ieee754_float
+  {
+    float f;
+
+    /* This is the IEEE 754 single-precision format.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:8;
+	unsigned int mantissa:23;
+#endif				/* Big endian.  */
+#if	__BYTE_ORDER == __LITTLE_ENDIAN
+	unsigned int mantissa:23;
+	unsigned int exponent:8;
+	unsigned int negative:1;
+#endif				/* Little endian.  */
+      } ieee;
+
+    /* This format makes it easier to see if a NaN is a signalling NaN.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:8;
+	unsigned int quiet_nan:1;
+	unsigned int mantissa:22;
+#endif				/* Big endian.  */
+#if	__BYTE_ORDER == __LITTLE_ENDIAN
+	unsigned int mantissa:22;
+	unsigned int quiet_nan:1;
+	unsigned int exponent:8;
+	unsigned int negative:1;
+#endif				/* Little endian.  */
+      } ieee_nan;
+  };
+
+#define IEEE754_FLOAT_BIAS	0x7f /* Added to exponent.  */
+
+
+union ieee754_double
+  {
+    double d;
+
+    /* This is the IEEE 754 double-precision format.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:11;
+	/* Together these comprise the mantissa.  */
+	unsigned int mantissa0:20;
+	unsigned int mantissa1:32;
+#endif				/* Big endian.  */
+#if	__BYTE_ORDER == __LITTLE_ENDIAN
+	/* Together these comprise the mantissa.  */
+	unsigned int mantissa1:32;
+	unsigned int mantissa0:20;
+	unsigned int exponent:11;
+	unsigned int negative:1;
+#endif				/* Little endian.  */
+      } ieee;
+
+    /* This format makes it easier to see if a NaN is a signalling NaN.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:11;
+	unsigned int quiet_nan:1;
+	/* Together these comprise the mantissa.  */
+	unsigned int mantissa0:19;
+	unsigned int mantissa1:32;
+#else
+	/* Together these comprise the mantissa.  */
+	unsigned int mantissa1:32;
+	unsigned int mantissa0:19;
+	unsigned int quiet_nan:1;
+	unsigned int exponent:11;
+	unsigned int negative:1;
+#endif
+      } ieee_nan;
+  };
+
+#define IEEE754_DOUBLE_BIAS	0x3ff /* Added to exponent.  */
+
+
+union ieee854_long_double
+  {
+    long double d;
+
+    /* This is the IEEE 854 quad-precision format.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:15;
+	/* Together these comprise the mantissa.  */
+	unsigned int mantissa0:16;
+	unsigned int mantissa1:32;
+	unsigned int mantissa2:32;
+	unsigned int mantissa3:32;
+#endif				/* Big endian.  */
+#if	__BYTE_ORDER == __LITTLE_ENDIAN
+	/* Together these comprise the mantissa.  */
+	unsigned int mantissa3:32;
+	unsigned int mantissa2:32;
+	unsigned int mantissa1:32;
+	unsigned int mantissa0:16;
+	unsigned int exponent:15;
+	unsigned int negative:1;
+#endif				/* Little endian.  */
+      } ieee;
+
+    /* This format makes it easier to see if a NaN is a signalling NaN.  */
+    struct
+      {
+#if	__BYTE_ORDER == __BIG_ENDIAN
+	unsigned int negative:1;
+	unsigned int exponent:15;
+	unsigned int quiet_nan:1;
+	/* Together these comprise the mantissa.  */
+	unsigned int mantissa0:15;
+	unsigned int mantissa1:32;
+	unsigned int mantissa2:32;
+	unsigned int mantissa3:32;
+#else
+	/* Together these comprise the mantissa.  */
+	unsigned int mantissa3:32;
+	unsigned int mantissa2:32;
+	unsigned int mantissa1:32;
+	unsigned int mantissa0:15;
+	unsigned int quiet_nan:1;
+	unsigned int exponent:15;
+	unsigned int negative:1;
+#endif
+      } ieee_nan;
+  };
+
+#define IEEE854_LONG_DOUBLE_BIAS 0x3fff /* Added to exponent.  */
+
+__END_DECLS
+
+#endif /* ieee754.h */
diff --git a/sysdeps/ieee754/ldbl-128/ldbl2mpn.c b/sysdeps/ieee754/ldbl-128/ldbl2mpn.c
new file mode 100644
index 0000000..68ecea8
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/ldbl2mpn.c
@@ -0,0 +1,137 @@
+/* Copyright (C) 1995, 1996, 1997, 1998, 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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"
+#include <ieee754.h>
+#include <float.h>
+#include <math.h>
+#include <stdlib.h>
+
+/* Convert a `long double' 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_long_double (mp_ptr res_ptr, mp_size_t size,
+			   int *expt, int *is_neg,
+			   long double value)
+{
+  union ieee854_long_double u;
+  u.d = value;
+
+  *is_neg = u.ieee.negative;
+  *expt = (int) u.ieee.exponent - IEEE854_LONG_DOUBLE_BIAS;
+
+#if BITS_PER_MP_LIMB == 32
+  res_ptr[0] = u.ieee.mantissa3; /* Low-order 32 bits of fraction.  */
+  res_ptr[1] = u.ieee.mantissa2;
+  res_ptr[2] = u.ieee.mantissa1;
+  res_ptr[3] = u.ieee.mantissa0; /* High-order 32 bits.  */
+  #define N 4
+#elif BITS_PER_MP_LIMB == 64
+  /* Hopefully the compiler will combine the two bitfield extracts
+     and this composition into just the original quadword extract.  */
+  res_ptr[0] = ((unsigned long int) u.ieee.mantissa2 << 32) | u.ieee.mantissa3;
+  res_ptr[1] = ((unsigned long int) u.ieee.mantissa0 << 32) | u.ieee.mantissa1;
+  #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 \
+			   - (LDBL_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 = LDBL_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 = LDBL_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 = LDBL_MIN_EXP - 1 - 3 * BITS_PER_MP_LIMB - cnt;
+#endif
+	}
+    }
+  else
+    /* Add the implicit leading one bit for a normalized number.  */
+    res_ptr[N - 1] |= 1L << (LDBL_MANT_DIG - 1 - ((N - 1) * BITS_PER_MP_LIMB));
+
+  return N;
+}
diff --git a/sysdeps/ieee754/ldbl-128/math_ldbl.h b/sysdeps/ieee754/ldbl-128/math_ldbl.h
new file mode 100644
index 0000000..aecb20a
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/math_ldbl.h
@@ -0,0 +1,82 @@
+#ifndef _MATH_PRIVATE_H_
+#error "Never use <math_ldbl.h> directly; include <math_private.h> instead."
+#endif
+
+/* A union which permits us to convert between a long double and
+   four 32 bit ints or two 64 bit ints.  */
+
+#if __FLOAT_WORD_ORDER == BIG_ENDIAN
+
+typedef union
+{
+  long double value;
+  struct
+  {
+    u_int64_t msw;
+    u_int64_t lsw;
+  } parts64;
+  struct
+  {
+    u_int32_t w0, w1, w2, w3;
+  } parts32;
+} ieee854_long_double_shape_type;
+
+#endif
+
+#if __FLOAT_WORD_ORDER == LITTLE_ENDIAN
+
+typedef union
+{
+  long double value;
+  struct
+  {
+    u_int64_t lsw;
+    u_int64_t msw;
+  } parts64;
+  struct
+  {
+    u_int32_t w3, w2, w1, w0;
+  } parts32;
+} ieee854_long_double_shape_type;
+
+#endif
+
+/* Get two 64 bit ints from a long double.  */
+
+#define GET_LDOUBLE_WORDS64(ix0,ix1,d)				\
+do {								\
+  ieee854_long_double_shape_type qw_u;				\
+  qw_u.value = (d);						\
+  (ix0) = qw_u.parts64.msw;					\
+  (ix1) = qw_u.parts64.lsw;					\
+} while (0)
+
+/* Set a long double from two 64 bit ints.  */
+
+#define SET_LDOUBLE_WORDS64(d,ix0,ix1)				\
+do {								\
+  ieee854_long_double_shape_type qw_u;				\
+  qw_u.parts64.msw = (ix0);					\
+  qw_u.parts64.lsw = (ix1);					\
+  (d) = qw_u.value;						\
+} while (0)
+
+/* Get the more significant 64 bits of a long double mantissa.  */
+
+#define GET_LDOUBLE_MSW64(v,d)					\
+do {								\
+  ieee854_long_double_shape_type sh_u;				\
+  sh_u.value = (d);						\
+  (v) = sh_u.parts64.msw;					\
+} while (0)
+
+/* Set the more significant 64 bits of a long double mantissa from an int.  */
+
+#define SET_LDOUBLE_MSW64(d,v)					\
+do {								\
+  ieee854_long_double_shape_type sh_u;				\
+  sh_u.value = (d);						\
+  sh_u.parts64.msw = (v);					\
+  (d) = sh_u.value;						\
+} while (0)
+
diff --git a/sysdeps/ieee754/ldbl-128/mpn2ldbl.c b/sysdeps/ieee754/ldbl-128/mpn2ldbl.c
new file mode 100644
index 0000000..d9cc452
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/mpn2ldbl.c
@@ -0,0 +1,51 @@
+/* Copyright (C) 1995, 1996, 1997, 1998, 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <ieee754.h>
+#include <float.h>
+#include <math.h>
+
+/* Convert a multi-precision integer of the needed number of bits (113 for
+   long double) and an integral power of two to a `long double' in IEEE854
+   quad-precision format.  */
+
+long double
+__mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
+{
+  union ieee854_long_double u;
+
+  u.ieee.negative = sign;
+  u.ieee.exponent = expt + IEEE854_LONG_DOUBLE_BIAS;
+#if BITS_PER_MP_LIMB == 32
+  u.ieee.mantissa3 = frac_ptr[0];
+  u.ieee.mantissa2 = frac_ptr[1];
+  u.ieee.mantissa1 = frac_ptr[2];
+  u.ieee.mantissa0 = frac_ptr[3] & ((1 << (LDBL_MANT_DIG - 96)) - 1);
+#elif BITS_PER_MP_LIMB == 64
+  u.ieee.mantissa3 = frac_ptr[0] & ((1L << 32) - 1);
+  u.ieee.mantissa2 = frac_ptr[0] >> 32;
+  u.ieee.mantissa1 = frac_ptr[1] & ((1L << 32) - 1);
+  u.ieee.mantissa0 = (frac_ptr[1] >> 32) & ((1 << (LDBL_MANT_DIG - 96)) - 1);
+#else
+  #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+
+  return u.d;
+}
diff --git a/sysdeps/ieee754/ldbl-128/printf_fphex.c b/sysdeps/ieee754/ldbl-128/printf_fphex.c
new file mode 100644
index 0000000..e25d668
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/printf_fphex.c
@@ -0,0 +1,84 @@
+/* Print floating point number in hexadecimal notation according to
+   ISO C 9X.
+   Copyright (C) 1997, 1998, 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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.  */
+
+#define PRINT_FPHEX_LONG_DOUBLE \
+do {									      \
+      /* 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.  */						      \
+      unsigned long long int num0, num1;				      \
+									      \
+      assert (sizeof (long double) == 16);				      \
+									      \
+      num0 = (((unsigned long long int) fpnum.ldbl.ieee.mantissa0) << 32      \
+	     | fpnum.ldbl.ieee.mantissa1);				      \
+      num1 = (((unsigned long long int) fpnum.ldbl.ieee.mantissa2) << 32      \
+	     | fpnum.ldbl.ieee.mantissa3);				      \
+									      \
+      zero_mantissa = (num0|num1) == 0;					      \
+									      \
+      if (sizeof (unsigned long int) > 6)				      \
+	numstr = _itoa_word (num1, numbuf + sizeof numbuf, 16,		      \
+			     info->spec == 'A');			      \
+      else								      \
+	numstr = _itoa (num1, numbuf + sizeof numbuf, 16,		      \
+			info->spec == 'A');				      \
+									      \
+      while (numstr > numbuf + (sizeof numbuf - 64 / 4))		      \
+	*--numstr = '0';						      \
+									      \
+      if (sizeof (unsigned long int) > 6)				      \
+	numstr = _itoa_word (num0, numstr, 16, info->spec == 'A');	      \
+      else								      \
+	numstr = _itoa (num0, numstr, 16, info->spec == 'A');		      \
+									      \
+      /* Fill with zeroes.  */						      \
+      while (numstr > numbuf + (sizeof numbuf - 112 / 4))		      \
+	*--numstr = '0';						      \
+									      \
+      leading = fpnum.ldbl.ieee.exponent == 0 ? '0' : '1';		      \
+									      \
+      exponent = fpnum.ldbl.ieee.exponent;				      \
+									      \
+      if (exponent == 0)						      \
+	{								      \
+	  if (zero_mantissa)						      \
+	    expnegative = 0;						      \
+	  else								      \
+	    {								      \
+	      /* This is a denormalized number.  */			      \
+	      expnegative = 1;						      \
+	      exponent = IEEE854_LONG_DOUBLE_BIAS - 1;			      \
+	    }								      \
+	}								      \
+      else if (exponent >= IEEE854_LONG_DOUBLE_BIAS)			      \
+	{								      \
+	  expnegative = 0;						      \
+	  exponent -= IEEE854_LONG_DOUBLE_BIAS;				      \
+	}								      \
+      else								      \
+	{								      \
+	  expnegative = 1;						      \
+	  exponent = -(exponent - IEEE854_LONG_DOUBLE_BIAS);		      \
+	}								      \
+} while (0)
+
+#include <sysdeps/generic/printf_fphex.c>
diff --git a/sysdeps/ieee754/ldbl-128/s_ceill.c b/sysdeps/ieee754/ldbl-128/s_ceill.c
new file mode 100644
index 0000000..f241554
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_ceill.c
@@ -0,0 +1,84 @@
+/* s_ceill.c -- long double version of s_ceil.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
+
+/*
+ * ceill(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to ceil(x).
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double huge = 1.0e4930;
+#else
+static long double huge = 1.0e4930;
+#endif
+
+#ifdef __STDC__
+	long double __ceill(long double x)
+#else
+	long double __ceill(x)
+	long double x;
+#endif
+{
+	int64_t i0,i1,j0;
+	u_int64_t i,j;
+	GET_LDOUBLE_WORDS64(i0,i1,x);
+	j0 = ((i0>>48)&0x7fff)-0x3fff;
+	if(j0<48) {
+	    if(j0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
+		    if(i0<0) {i0=0x8000000000000000ULL;i1=0;}
+		    else if((i0|i1)!=0) { i0=0x3fff000000000000ULL;i1=0;}
+		}
+	    } else {
+		i = (0x7fffffffffffffffULL)>>j0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if(huge+x>0.0) {	/* raise inexact flag */
+		    if(i0>0) i0 += (0x0001000000000000LL)>>j0;
+		    i0 &= (~i); i1=0;
+		}
+	    }
+	} 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 */
+	    if(huge+x>0.0) { 		/* raise inexact flag */
+		if(i0>0) {
+		    if(j0==48) i0+=1;
+		    else {
+			j = i1+(1LL<<(112-j0));
+			if(j<i1) i0 +=1 ; 	/* got a carry */
+			i1=j;
+		    }
+		}
+		i1 &= (~i);
+	    }
+	}
+	SET_LDOUBLE_WORDS64(x,i0,i1);
+	return x;
+}
+weak_alias (__ceill, ceill)
diff --git a/sysdeps/ieee754/ldbl-128/s_copysignl.c b/sysdeps/ieee754/ldbl-128/s_copysignl.c
new file mode 100644
index 0000000..cece4f2
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_copysignl.c
@@ -0,0 +1,43 @@
+/* s_copysignl.c -- long double version of s_copysign.c.
+ * Conversion to 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
+
+/*
+ * copysignl(long double x, long double y)
+ * copysignl(x,y) returns a value with the magnitude of x and
+ * with the sign bit of y.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	long double __copysignl(long double x, long double y)
+#else
+	long double __copysignl(x,y)
+	long double x,y;
+#endif
+{
+	u_int64_t hx,hy;
+	GET_LDOUBLE_MSW64(hx,x);
+	GET_LDOUBLE_MSW64(hy,y);
+	SET_LDOUBLE_MSW64(x,(hx&0x7fffffffffffffffULL)
+			    |(hy&0x8000000000000000ULL));
+        return x;
+}
+weak_alias (__copysignl, copysignl)
diff --git a/sysdeps/ieee754/ldbl-128/s_cosl.c b/sysdeps/ieee754/ldbl-128/s_cosl.c
new file mode 100644
index 0000000..d1258b2
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_cosl.c
@@ -0,0 +1,83 @@
+/* s_cosl.c -- long double version of s_cos.c.
+ * Conversion to 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.
+ * ====================================================
+ */
+
+/* cosl(x)
+ * Return cosine function of x.
+ *
+ * kernel function:
+ *	__kernel_sinl		... sine function on [-pi/4,pi/4]
+ *	__kernel_cosl		... cosine function on [-pi/4,pi/4]
+ *	__ieee754_rem_pio2l	... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *	in [-pi/4 , +pi/4], and let n = k mod 4.
+ *	We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *	    0	       S	   C		 T
+ *	    1	       C	  -S		-1/T
+ *	    2	      -S	  -C		 T
+ *	    3	      -C	   S		-1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *	TRIG(x) returns trig(x) nearly rounded
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	long double __cosl(long double x)
+#else
+	long double __cosl(x)
+	long double x;
+#endif
+{
+	long double y[2],z=0.0L;
+	int64_t n, ix;
+
+    /* High word of x. */
+	GET_LDOUBLE_MSW64(ix,x);
+
+    /* |x| ~< pi/4 */
+	ix &= 0x7fffffffffffffffLL;
+	if(ix <= 0x3ffe921fb54442d1LL)
+	  return __kernel_cosl(x,z);
+
+    /* cos(Inf or NaN) is NaN */
+	else if (ix>=0x7fff000000000000LL) return x-x;
+
+    /* argument reduction needed */
+	else {
+	    n = __ieee754_rem_pio2l(x,y);
+	    switch(n&3) {
+		case 0: return  __kernel_cosl(y[0],y[1]);
+		case 1: return -__kernel_sinl(y[0],y[1],1);
+		case 2: return -__kernel_cosl(y[0],y[1]);
+		default:
+		        return  __kernel_sinl(y[0],y[1],1);
+	    }
+	}
+}
+weak_alias (__cosl, cosl)
diff --git a/sysdeps/ieee754/ldbl-128/s_fabsl.c b/sysdeps/ieee754/ldbl-128/s_fabsl.c
new file mode 100644
index 0000000..c0fd05a
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_fabsl.c
@@ -0,0 +1,39 @@
+/* s_fabsl.c -- long double version of s_fabs.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
+
+/*
+ * fabsl(x) returns the absolute value of x.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	long double __fabsl(long double x)
+#else
+	long double __fabsl(x)
+	long double x;
+#endif
+{
+	u_int64_t hx;
+	GET_LDOUBLE_MSW64(hx,x);
+	SET_LDOUBLE_MSW64(x,hx&0x7fffffffffffffffLL);
+        return x;
+}
+weak_alias (__fabsl, fabsl)
diff --git a/sysdeps/ieee754/ldbl-128/s_finitel.c b/sysdeps/ieee754/ldbl-128/s_finitel.c
new file mode 100644
index 0000000..dd176c1
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_finitel.c
@@ -0,0 +1,40 @@
+/* s_finitel.c -- long double version of s_finite.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
+
+/*
+ * finitel(x) returns 1 is x is finite, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	int __finitel(long double x)
+#else
+	int __finitel(x)
+	long double x;
+#endif
+{
+	int64_t hx;
+	GET_LDOUBLE_MSW64(hx,x);
+	return (int)((u_int64_t)((hx&0x7fffffffffffffffLL)
+				 -0x7fff000000000000LL)>>63);
+}
+weak_alias (__finitel, finitel)
diff --git a/sysdeps/ieee754/ldbl-128/s_floorl.c b/sysdeps/ieee754/ldbl-128/s_floorl.c
new file mode 100644
index 0000000..c9b8b70
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_floorl.c
@@ -0,0 +1,85 @@
+/* s_floorl.c -- long double version of s_floor.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
+
+/*
+ * floorl(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to floor(x).
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double huge = 1.0e4930;
+#else
+static long double huge = 1.0e4930;
+#endif
+
+#ifdef __STDC__
+	long double __floorl(long double x)
+#else
+	long double __floorl(x)
+	long double x;
+#endif
+{
+	int64_t i0,i1,j0;
+	u_int64_t i,j;
+	GET_LDOUBLE_WORDS64(i0,i1,x);
+	j0 = ((i0>>48)&0x7fff)-0x3fff;
+	if(j0<48) {
+	    if(j0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
+		    if(i0>=0) {i0=i1=0;}
+		    else if(((i0&0x7fffffffffffffffLL)|i1)!=0)
+			{ i0=0xbfff000000000000ULL;i1=0;}
+		}
+	    } else {
+		i = (0x7fffffffffffffffULL)>>j0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if(huge+x>0.0) {	/* raise inexact flag */
+		    if(i0<0) i0 += (0x0001000000000000LL)>>j0;
+		    i0 &= (~i); i1=0;
+		}
+	    }
+	} 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 */
+	    if(huge+x>0.0) { 		/* raise inexact flag */
+		if(i0<0) {
+		    if(j0==48) i0+=1;
+		    else {
+			j = i1+(1LL<<(112-j0));
+			if(j<i1) i0 +=1 ; 	/* got a carry */
+			i1=j;
+		    }
+		}
+		i1 &= (~i);
+	    }
+	}
+	SET_LDOUBLE_WORDS64(x,i0,i1);
+	return x;
+}
+weak_alias (__floorl, floorl)
diff --git a/sysdeps/ieee754/ldbl-128/s_fpclassifyl.c b/sysdeps/ieee754/ldbl-128/s_fpclassifyl.c
new file mode 100644
index 0000000..868c3c4
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_fpclassifyl.c
@@ -0,0 +1,44 @@
+/* Return classification value corresponding to argument.
+   Copyright (C) 1997, 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+#include "math_private.h"
+
+
+int
+__fpclassifyl (long double x)
+{
+  u_int64_t hx, lx;
+  int retval = FP_NORMAL;
+
+  GET_LDOUBLE_WORDS64 (hx, lx, x);
+  lx |= (hx & 0x0000ffffffffffffLL);
+  hx &= 0x7fff000000000000LL;
+  if ((hx | lx) == 0)
+    retval = FP_ZERO;
+  else if (hx == 0)
+    retval = FP_SUBNORMAL;
+  else if (hx == 0x7fff000000000000LL)
+    retval = lx != 0 ? FP_NAN : FP_INFINITE;
+
+  return retval;
+}
diff --git a/sysdeps/ieee754/ldbl-128/s_frexpl.c b/sysdeps/ieee754/ldbl-128/s_frexpl.c
new file mode 100644
index 0000000..6dbb60e
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_frexpl.c
@@ -0,0 +1,63 @@
+/* s_frexpl.c -- long double version of s_frexp.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
+
+/*
+ * for non-zero x
+ *	x = frexpl(arg,&exp);
+ * return a long double fp quantity x such that 0.5 <= |x| <1.0
+ * and the corresponding binary exponent "exp". That is
+ *	arg = x*2^exp.
+ * If arg is inf, 0.0, or NaN, then frexpl(arg,&exp) returns arg
+ * with *exp=0.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double
+#else
+static long double
+#endif
+two114 = 2.0769187434139310514121985316880384E+34L; /* 0x4071000000000000, 0 */
+
+#ifdef __STDC__
+	long double __frexpl(long double x, int *eptr)
+#else
+	long double __frexpl(x, eptr)
+	long double x; int *eptr;
+#endif
+{
+	u_int64_t hx, lx, ix;
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	ix = 0x7fffffffffffffffULL&hx;
+	*eptr = 0;
+	if(ix>=0x7fff000000000000ULL||((ix|lx)==0)) return x;	/* 0,inf,nan */
+	if (ix<0x0001000000000000ULL) {		/* subnormal */
+	    x *= two114;
+	    GET_LDOUBLE_MSW64(hx,x);
+	    ix = hx&0x7fffffffffffffffULL;
+	    *eptr = -114;
+	}
+	*eptr += (ix>>48)-16382;
+	hx = (hx&0x8000ffffffffffffULL) | 0x3ffe000000000000ULL;
+	SET_LDOUBLE_MSW64(x,hx);
+	return x;
+}
+weak_alias (__frexpl, frexpl)
diff --git a/sysdeps/ieee754/ldbl-128/s_ilogbl.c b/sysdeps/ieee754/ldbl-128/s_ilogbl.c
new file mode 100644
index 0000000..d2acfd3
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_ilogbl.c
@@ -0,0 +1,55 @@
+/* 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) = 0x80000001
+ * ilogbl(inf/NaN) = 0x7fffffff (no signal is raised)
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	int __ilogbl(long double x)
+#else
+	int __ilogbl(x)
+	long double x;
+#endif
+{
+	int64_t hx,lx;
+	int ix;
+
+	GET_LDOUBLE_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 return FP_ILOGBNAN;
+}
+weak_alias (__ilogbl, ilogbl)
diff --git a/sysdeps/ieee754/ldbl-128/s_isinfl.c b/sysdeps/ieee754/ldbl-128/s_isinfl.c
new file mode 100644
index 0000000..038c294
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_isinfl.c
@@ -0,0 +1,28 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Change for long double by Jakub Jelinek <jj@ultra.linux.cz>
+ * Public domain.
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: $";
+#endif
+
+/*
+ * isinfl(x) returns 1 if x is inf, -1 if x is -inf, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+int
+__isinfl (long double x)
+{
+	int64_t hx,lx;
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	lx |= (hx & 0x7fffffffffffffffLL) ^ 0x7fff000000000000LL;
+	lx |= -lx;
+	return ~(lx >> 63) & (hx >> 62);
+}
+weak_alias (__isinfl, isinfl)
diff --git a/sysdeps/ieee754/ldbl-128/s_isnanl.c b/sysdeps/ieee754/ldbl-128/s_isnanl.c
new file mode 100644
index 0000000..d2fb403
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_isnanl.c
@@ -0,0 +1,42 @@
+/* s_isnanl.c -- long double version of s_isnan.c.
+ * Conversion to 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
+
+/*
+ * isnanl(x) returns 1 is x is nan, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	int __isnanl(long double x)
+#else
+	int __isnanl(x)
+	long double x;
+#endif
+{
+	int64_t hx,lx;
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	hx &= 0x7fffffffffffffffLL;
+	hx |= (u_int64_t)(lx|(-lx))>>63;
+	hx = 0x7fff000000000000LL - hx;
+	return (int)((u_int64_t)hx>>63);
+}
+weak_alias (__isnanl, isnanl)
diff --git a/sysdeps/ieee754/ldbl-128/s_llrintl.c b/sysdeps/ieee754/ldbl-128/s_llrintl.c
new file mode 100644
index 0000000..389a65d
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_llrintl.c
@@ -0,0 +1,75 @@
+/* Round argument to nearest integral value according to current rounding
+   direction.
+   Copyright (C) 1997, 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+#include "math_private.h"
+
+static const long double two112[2] =
+{
+  5.19229685853482762853049632922009600E+33L, /* 0x406F000000000000, 0 */
+ -5.19229685853482762853049632922009600E+33L  /* 0xC06F000000000000, 0 */
+};
+
+long long int
+__llrintl (long double x)
+{
+  int32_t j0;
+  u_int64_t i0,i1;
+  volatile long double w;
+  long double t;
+  long long int result;
+  int sx;
+
+  GET_LDOUBLE_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)
+    {
+      if (j0 < -1)
+	return 0;
+      w = two112[sx] + x;
+      t = w - two112[sx];
+      GET_LDOUBLE_WORDS64 (i0, i1, x);
+      j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+      i0 &= 0x0000ffffffffffffLL;
+      i0 |= 0x0001000000000000LL;
+
+      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;
+}
+
+weak_alias (__llrintl, llrintl)
diff --git a/sysdeps/ieee754/ldbl-128/s_llroundl.c b/sysdeps/ieee754/ldbl-128/s_llroundl.c
new file mode 100644
index 0000000..82395a7
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_llroundl.c
@@ -0,0 +1,74 @@
+/* Round long double value to long long int.
+   Copyright (C) 1997, 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+#include "math_private.h"
+
+
+long long int
+__llroundl (long double x)
+{
+  int64_t j0;
+  u_int64_t i1, i0;
+  long long int result;
+  int sign;
+
+  GET_LDOUBLE_WORDS64 (i0, i1, x);
+  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+  sign = (i0 & 0x8000000000000000ULL) != 0 ? -1 : 1;
+  i0 &= 0x0000ffffffffffffLL;
+  i0 |= 0x0001000000000000LL;
+
+  if (j0 < 48)
+    {
+      if (j0 < 0)
+	return j0 < -1 ? 0 : sign;
+      else
+	{
+	  i0 += 0x0000800000000000LL >> j0;
+	  result = i0 >> (48 - j0);
+	}
+    }
+  else if (j0 < (int32_t) (8 * sizeof (long long int)) - 1)
+    {
+      if (j0 >= 112)
+	result = ((long long int) i0 << (j0 - 48)) | (i1 << (j0 - 112));
+      else
+	{
+	  u_int64_t j = i1 + (0x8000000000000000ULL >> (j0 - 48));
+	  if (j < i1)
+	    ++i0;
+
+	  result = ((long long int) i0 << (j0 - 48)) | (j >> (112 - j0));
+	}
+    }
+  else
+    {
+      /* The number is too large.  It is left implementation defined
+	 what happens.  */
+      return (long long int) x;
+    }
+
+  return sign * result;
+}
+
+weak_alias (__llroundl, llroundl)
diff --git a/sysdeps/ieee754/ldbl-128/s_logbl.c b/sysdeps/ieee754/ldbl-128/s_logbl.c
new file mode 100644
index 0000000..1fda289
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_logbl.c
@@ -0,0 +1,46 @@
+/* s_logbl.c -- long double version of s_logb.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
+
+/*
+ * long double logbl(x)
+ * IEEE 754 logb. Included to pass IEEE test suite. Not recommend.
+ * Use ilogb instead.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	long double __logbl(long double x)
+#else
+	long double __logbl(x)
+	long double x;
+#endif
+{
+	int64_t lx,hx;
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	hx &= 0x7fffffffffffffffLL;		/* high |x| */
+	if((hx|lx)==0) return -1.0/fabs(x);
+	if(hx>=0x7fff000000000000LL) return x*x;
+	if((hx>>=48)==0) 			/* IEEE 754 logb */
+		return -16382.0;
+	else
+		return (long double) (hx-0x3fff);
+}
+weak_alias (__logbl, logbl)
diff --git a/sysdeps/ieee754/ldbl-128/s_lrintl.c b/sysdeps/ieee754/ldbl-128/s_lrintl.c
new file mode 100644
index 0000000..434aa31
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_lrintl.c
@@ -0,0 +1,91 @@
+/* Round argument to nearest integral value according to current rounding
+   direction.
+   Copyright (C) 1997, 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+#include "math_private.h"
+
+static const long double two112[2] =
+{
+  5.19229685853482762853049632922009600E+33L, /* 0x406F000000000000, 0 */
+ -5.19229685853482762853049632922009600E+33L  /* 0xC06F000000000000, 0 */
+};
+
+long int
+__lrintl (long double x)
+{
+  int32_t j0;
+  u_int64_t i0,i1;
+  volatile long double w;
+  long double t;
+  long int result;
+  int sx;
+
+  GET_LDOUBLE_WORDS64 (i0, i1, x);
+  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+  sx = i0 >> 63;
+  i0 &= 0x0000ffffffffffffLL;
+  i0 |= 0x0001000000000000LL;
+
+  if (j0 < 48)
+    {
+      if (j0 < -1)
+	return 0;
+      else
+	{
+	  w = two112[sx] + x;
+	  t = w - two112[sx];
+	  GET_LDOUBLE_WORDS64 (i0, i1, x);
+	  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+	  i0 &= 0x0000ffffffffffffLL;
+	  i0 |= 0x0001000000000000LL;
+
+	  result = 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_LDOUBLE_WORDS64 (i0, i1, x);
+	  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+	  i0 &= 0x0000ffffffffffffLL;
+	  i0 |= 0x0001000000000000LL;
+
+	  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;
+}
+
+weak_alias (__lrintl, lrintl)
diff --git a/sysdeps/ieee754/ldbl-128/s_lroundl.c b/sysdeps/ieee754/ldbl-128/s_lroundl.c
new file mode 100644
index 0000000..ea82f4c
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_lroundl.c
@@ -0,0 +1,74 @@
+/* Round long double value to long int.
+   Copyright (C) 1997, 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+#include "math_private.h"
+
+
+long int
+__lroundl (long double x)
+{
+  int64_t j0;
+  u_int64_t i1, i0;
+  long int result;
+  int sign;
+
+  GET_LDOUBLE_WORDS64 (i0, i1, x);
+  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+  sign = (i0 & 0x8000000000000000ULL) != 0 ? -1 : 1;
+  i0 &= 0x0000ffffffffffffLL;
+  i0 |= 0x0001000000000000LL;
+
+  if (j0 < 48)
+    {
+      if (j0 < 0)
+	return j0 < -1 ? 0 : sign;
+      else
+	{
+	  i0 += 0x0000800000000000LL >> j0;
+	  result = 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
+	{
+	  u_int64_t j = i1 + (0x8000000000000000ULL >> (j0 - 48));
+	  if (j < i1)
+	    ++i0;
+
+	  result = ((long int) i0 << (j0 - 48)) | (j >> (112 - j0));
+	}
+    }
+  else
+    {
+      /* The number is too large.  It is left implementation defined
+	 what happens.  */
+      return (long int) x;
+    }
+
+  return sign * result;
+}
+
+weak_alias (__lroundl, lroundl)
diff --git a/sysdeps/ieee754/ldbl-128/s_modfl.c b/sysdeps/ieee754/ldbl-128/s_modfl.c
new file mode 100644
index 0000000..63d66e7
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_modfl.c
@@ -0,0 +1,88 @@
+/* s_modfl.c -- long double version of s_modf.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
+
+/*
+ * modfl(long double x, long double *iptr)
+ * return fraction part of x, and return x's integral part in *iptr.
+ * Method:
+ *	Bit twiddling.
+ *
+ * Exception:
+ *	No exception.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double one = 1.0;
+#else
+static long double one = 1.0;
+#endif
+
+#ifdef __STDC__
+	long double __modfl(long double x, long double *iptr)
+#else
+	long double __modfl(x, iptr)
+	long double x,*iptr;
+#endif
+{
+	int64_t i0,i1,j0;
+	u_int64_t i;
+	GET_LDOUBLE_WORDS64(i0,i1,x);
+	j0 = ((i0>>48)&0x7fff)-0x3fff;	/* exponent of x */
+	if(j0<48) {			/* integer part in high x */
+	    if(j0<0) {			/* |x|<1 */
+		/* *iptr = +-0 */
+	        SET_LDOUBLE_WORDS64(*iptr,i0&0x8000000000000000ULL,0);
+		return x;
+	    } else {
+		i = (0x0000ffffffffffffLL)>>j0;
+		if(((i0&i)|i1)==0) {		/* x is integral */
+		    *iptr = x;
+		    /* return +-0 */
+		    SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
+		    return x;
+		} else {
+		    SET_LDOUBLE_WORDS64(*iptr,i0&(~i),0);
+		    return x - *iptr;
+		}
+	    }
+	} else if (j0>111) {		/* no fraction part */
+	    *iptr = x*one;
+	    /* We must handle NaNs separately.  */
+	    if (j0 == 0x4000 && ((i0 & 0x0000ffffffffffffLL) | i1))
+	      return x*one;
+	    /* return +-0 */
+	    SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
+	    return x;
+	} else {			/* fraction part in low x */
+	    i = -1ULL>>(j0-48);
+	    if((i1&i)==0) { 		/* x is integral */
+		*iptr = x;
+		/* return +-0 */
+		SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
+		return x;
+	    } else {
+		SET_LDOUBLE_WORDS64(*iptr,i0,i1&(~i));
+		return x - *iptr;
+	    }
+	}
+}
+weak_alias (__modfl, modfl)
diff --git a/sysdeps/ieee754/ldbl-128/s_nearbyintl.c b/sysdeps/ieee754/ldbl-128/s_nearbyintl.c
new file mode 100644
index 0000000..bea3183
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_nearbyintl.c
@@ -0,0 +1,93 @@
+/* s_nearbyintl.c -- long double 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.
+ * ====================================================
+ */
+
+/*
+ * nearbyintl(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 rintl(x).
+ */
+
+#include <fenv.h>
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double
+#else
+static long double
+#endif
+TWO112[2]={
+  5.19229685853482762853049632922009600E+33L, /* 0x406F000000000000, 0 */
+ -5.19229685853482762853049632922009600E+33L  /* 0xC06F000000000000, 0 */
+};
+
+#ifdef __STDC__
+	long double __nearbyintl(long double x)
+#else
+	long double __nearbyintl(x)
+	long double x;
+#endif
+{
+	fenv_t env;
+	int64_t i0,j0,sx;
+	u_int64_t i,i1;
+	long double w,t;
+	GET_LDOUBLE_WORDS64(i0,i1,x);
+	sx = (((u_int64_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_LDOUBLE_MSW64(x,i0);
+		feholdexcept (&env);
+	        w = TWO112[sx]+x;
+	        t = w-TWO112[sx];
+	        fesetenv (&env);
+		GET_LDOUBLE_MSW64(i0,t);
+		SET_LDOUBLE_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_LDOUBLE_WORDS64(x,i0,i1);
+	feholdexcept (&env);
+	w = TWO112[sx]+x;
+	t = w-TWO112[sx];
+	fesetenv (&env);
+	return t;
+}
+weak_alias (__nearbyintl, nearbyintl)
diff --git a/sysdeps/ieee754/ldbl-128/s_nextafterl.c b/sysdeps/ieee754/ldbl-128/s_nextafterl.c
new file mode 100644
index 0000000..d3df668
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_nextafterl.c
@@ -0,0 +1,85 @@
+/* s_nextafterl.c -- long double version of s_nextafter.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
+
+/* IEEE functions
+ *	nextafterl(x,y)
+ *	return the next machine floating-point number of x in the
+ *	direction toward y.
+ *   Special cases:
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	long double __nextafterl(long double x, long double y)
+#else
+	long double __nextafterl(x,y)
+	long double x,y;
+#endif
+{
+	int64_t hx,hy,ix,iy;
+	u_int64_t lx,ly;
+
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	GET_LDOUBLE_WORDS64(hy,ly,y);
+	ix = hx&0x7fffffffffffffffLL;		/* |x| */
+	iy = hy&0x7fffffffffffffffLL;		/* |y| */
+
+	if(((ix>=0x7fff000000000000LL)&&((ix-0x7fff000000000000LL)|lx)!=0) ||   /* x is nan */
+	   ((iy>=0x7fff000000000000LL)&&((iy-0x7fff000000000000LL)|ly)!=0))     /* y is nan */
+	   return x+y;
+	if(x==y) return y;		/* x=y, return y */
+	if((ix|lx)==0) {			/* x == 0 */
+	    SET_LDOUBLE_WORDS64(x,hy&0x8000000000000000ULL,1);/* return +-minsubnormal */
+	    y = x*x;
+	    if(y==x) return y; else return x;	/* raise underflow flag */
+	}
+	if(hx>=0) {			/* x > 0 */
+	    if(hx>hy||((hx==hy)&&(lx>ly))) {	/* x > y, x -= ulp */
+		if(lx==0) hx--;
+		lx--;
+	    } else {				/* x < y, x += ulp */
+		lx++;
+		if(lx==0) hx++;
+	    }
+	} else {				/* x < 0 */
+	    if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){/* x < y, x -= ulp */
+		if(lx==0) hx--;
+		lx--;
+	    } else {				/* x > y, x += ulp */
+		lx++;
+		if(lx==0) hx++;
+	    }
+	}
+	hy = hx&0x7fff000000000000LL;
+	if(hy==0x7fff000000000000LL) return x+x;/* overflow  */
+	if(hy==0) {				/* underflow */
+	    y = x*x;
+	    if(y!=x) {		/* raise underflow flag */
+	        SET_LDOUBLE_WORDS64(y,hx,lx);
+		return y;
+	    }
+	}
+	SET_LDOUBLE_WORDS64(x,hx,lx);
+	return x;
+}
+weak_alias (__nextafterl, nextafterl)
+strong_alias (__nextafterl, __nexttowardl)
+weak_alias (__nextafterl, nexttowardl)
diff --git a/sysdeps/ieee754/ldbl-128/s_nexttoward.c b/sysdeps/ieee754/ldbl-128/s_nexttoward.c
new file mode 100644
index 0000000..f121be2
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_nexttoward.c
@@ -0,0 +1,97 @@
+/* s_nexttoward.c
+ * Conversion from s_nextafter.c by Ulrich Drepper, Cygnus Support,
+ * drepper@cygnus.com and 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
+
+/* IEEE functions
+ *	nexttoward(x,y)
+ *	return the next machine floating-point number of x in the
+ *	direction toward y.
+ *   Special cases:
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	double __nexttoward(double x, long double y)
+#else
+	double __nexttoward(x,y)
+	double x;
+	long double y;
+#endif
+{
+	int32_t hx,ix;
+	int64_t hy,iy;
+	u_int32_t lx;
+	u_int64_t ly;
+
+	EXTRACT_WORDS(hx,lx,x);
+	GET_LDOUBLE_WORDS64(hy,ly,y);
+	ix = hx&0x7fffffff;		/* |x| */
+	iy = hy&0x7fffffffffffffffLL;	/* |y| */
+
+	if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||   /* x is nan */
+	   ((iy>=0x7fff000000000000LL)&&((iy-0x7fff000000000000LL)|ly)!=0))
+	   						    /* y is nan */
+	   return x+y;
+	if((long double) x==y) return x;	/* x=y, return x */
+	if((ix|lx)==0) {			/* x == 0 */
+	    double x2;
+	    INSERT_WORDS(x,(u_int32_t)((hy>>32)&0x80000000),1);/* return +-minsub */
+	    x2 = x*x;
+	    if(x2==x) return x2; else return x;	/* raise underflow flag */
+	}
+	if(hx>=0) {				/* x > 0 */
+	    if (hy<0||(ix>>20)>(iy>>48)-0x3c00
+		|| ((ix>>20)==(iy>>48)-0x3c00
+		    && (((((int64_t)hx)<<28)|(lx>>4))>(hy&0x0000ffffffffffffLL)
+			|| (((((int64_t)hx)<<28)|(lx>>4))==(hy&0x0000ffffffffffffLL)
+			    && (lx&0xf)>(ly>>60))))) {	/* x > y, x -= ulp */
+		if(lx==0) hx -= 1;
+		lx -= 1;
+	    } else {				/* x < y, x += ulp */
+		lx += 1;
+		if(lx==0) hx += 1;
+	    }
+	} else {				/* x < 0 */
+	    if (hy>=0||(ix>>20)>(iy>>48)-0x3c00
+		|| ((ix>>20)==(iy>>48)-0x3c00
+		    && (((((int64_t)hx)<<28)|(lx>>4))>(hy&0x0000ffffffffffffLL)
+			|| (((((int64_t)hx)<<28)|(lx>>4))==(hy&0x0000ffffffffffffLL)
+			    && (lx&0xf)>(ly>>60))))) {	/* x < y, x -= ulp */
+		if(lx==0) hx -= 1;
+		lx -= 1;
+	    } else {				/* x > y, x += ulp */
+		lx += 1;
+		if(lx==0) hx += 1;
+	    }
+	}
+	hy = hx&0x7ff00000;
+	if(hy>=0x7ff00000) return x+x;	/* overflow  */
+	if(hy<0x00100000) {		/* underflow */
+	    double x2 = x*x;
+	    if(x2!=x) {		/* raise underflow flag */
+	        INSERT_WORDS(x2,hx,lx);
+		return x2;
+	    }
+	}
+	INSERT_WORDS(x,hx,lx);
+	return x;
+}
+weak_alias (__nexttoward, nexttoward)
diff --git a/sysdeps/ieee754/ldbl-128/s_nexttowardf.c b/sysdeps/ieee754/ldbl-128/s_nexttowardf.c
new file mode 100644
index 0000000..1a22e01
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_nexttowardf.c
@@ -0,0 +1,81 @@
+/* s_nexttowardf.c -- float version of s_nextafter.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com
+ * and 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
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	float __nexttowardf(float x, long double y)
+#else
+	float __nexttowardf(x,y)
+	float x;
+	long double y;
+#endif
+{
+	int32_t hx,ix;
+	int64_t hy,iy;
+	u_int64_t ly;
+                                
+	GET_FLOAT_WORD(hx,x);
+	GET_LDOUBLE_WORDS64(hy,ly,y);
+	ix = hx&0x7fffffff;		/* |x| */
+	iy = hy&0x7fffffffffffffffLL;	/* |y| */
+
+	if((ix>0x7f800000) ||   /* x is nan */
+	   ((iy>=0x7fff000000000000LL)&&((iy-0x7fff000000000000LL)|ly)!=0))
+				/* y is nan */
+	   return x+y;
+	if((long double) x==y) return y;	/* x=y, return y */
+	if(ix==0) {				/* x == 0 */
+	    float x2;
+	    SET_FLOAT_WORD(x,(u_int32_t)((hy>>32)&0x80000000)|1);/* return +-minsub*/
+	    x2 = x*x;
+	    if(x2==x) return x2; else return x;	/* raise underflow flag */
+	}
+	if(hx>=0) {				/* x > 0 */
+	    if(hy<0||(ix>>23)>(iy>>48)-0x3f80
+	       || ((ix>>23)==(iy>>48)-0x3f80
+		   && (ix&0x7fffff)>((hy>>25)&0x7fffff))) {/* x > y, x -= ulp */
+		hx -= 1;
+	    } else {				/* x < y, x += ulp */
+		hx += 1;
+	    }
+	} else {				/* x < 0 */
+	    if(hy>=0||(ix>>23)>(iy>>48)-0x3f80
+	       || ((ix>>23)==(iy>>48)-0x3f80
+		   && (ix&0x7fffff)>((hy>>25)&0x7fffff))) {/* x < y, x -= ulp */
+		hx -= 1;
+	    } else {				/* x > y, x += ulp */
+		hx += 1;
+	    }
+	}
+	hy = hx&0x7f800000;
+	if(hy>=0x7f800000) return x+x;	/* overflow  */
+	if(hy<0x00800000) {		/* underflow */
+	    float x2 = x*x;
+	    if(x2!=x) {		/* raise underflow flag */
+	        SET_FLOAT_WORD(x2,hx);
+		return x2;
+	    }
+	}
+	SET_FLOAT_WORD(x,hx);
+	return x;
+}
+weak_alias (__nexttowardf, nexttowardf)
diff --git a/sysdeps/ieee754/ldbl-128/s_remquol.c b/sysdeps/ieee754/ldbl-128/s_remquol.c
new file mode 100644
index 0000000..0d26958
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_remquol.c
@@ -0,0 +1,110 @@
+/* Compute remainder and a congruent to the quotient.
+   Copyright (C) 1997, 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+#include "math_private.h"
+
+
+static const long double zero = 0.0;
+
+
+long double
+__remquol (long double x, long double y, int *quo)
+{
+  int64_t hx,hy;
+  u_int64_t sx,lx,ly;
+  int cquo,qs;
+
+  GET_LDOUBLE_WORDS64 (hx, lx, x);
+  GET_LDOUBLE_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 = __ieee754_fmodl (x, 8 * y);              /* now x < 8y */
+      
+  if (((hx - hy) | (lx - ly)) == 0)
+    {
+      *quo = qs ? -1 : 1;
+      return zero * x;
+    }
+
+  x  = fabsl (x);
+  y  = fabsl (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
+    {
+      long double y_half = 0.5L * 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;
+}
+weak_alias (__remquol, remquol)
diff --git a/sysdeps/ieee754/ldbl-128/s_rintl.c b/sysdeps/ieee754/ldbl-128/s_rintl.c
new file mode 100644
index 0000000..c3fc3ba
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_rintl.c
@@ -0,0 +1,90 @@
+/* s_rintl.c -- long double version of s_rint.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
+
+/*
+ * rintl(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 rintl(x).
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double
+#else
+static long double
+#endif
+TWO112[2]={
+  5.19229685853482762853049632922009600E+33L, /* 0x406F000000000000, 0 */
+ -5.19229685853482762853049632922009600E+33L  /* 0xC06F000000000000, 0 */
+};
+
+#ifdef __STDC__
+	long double __rintl(long double x)
+#else
+	long double __rintl(x)
+	long double x;
+#endif
+{
+	int64_t i0,j0,sx;
+	u_int64_t i,i1;
+	long double w,t;
+	GET_LDOUBLE_WORDS64(i0,i1,x);
+	sx = (((u_int64_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_LDOUBLE_MSW64(x,i0);
+	        w = TWO112[sx]+x;
+	        t = w-TWO112[sx];
+		GET_LDOUBLE_MSW64(i0,t);
+		SET_LDOUBLE_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_LDOUBLE_WORDS64(x,i0,i1);
+	w = TWO112[sx]+x;
+	return w-TWO112[sx];
+}
+weak_alias (__rintl, rintl)
diff --git a/sysdeps/ieee754/ldbl-128/s_roundl.c b/sysdeps/ieee754/ldbl-128/s_roundl.c
new file mode 100644
index 0000000..f9fff48
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_roundl.c
@@ -0,0 +1,94 @@
+/* Round long double to integer away from zero.
+   Copyright (C) 1997, 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+#include "math_private.h"
+
+
+static const long double huge = 1.0E4930L;
+
+
+long double
+__roundl (long double x)
+{
+  int32_t j0;
+  u_int64_t se, i1, i0;
+
+  GET_LDOUBLE_WORDS64 (i0, i1, x);
+  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+  if (j0 < 31)
+    {
+      if (j0 < 0)
+	{
+	  if (huge + x > 0.0)
+	    {
+	      i0 &= 0x8000000000000000ULL;
+	      if (j0 == -1)
+		i0 |= 0x3fff000000000000LL;
+	      i1 = 0;
+	    }
+	}
+      else
+	{
+	  u_int64_t i = 0x0000ffffffffffffLL >> j0;
+	  if (((i0 & i) | i1) == 0)
+	    /* X is integral.  */
+	    return x;
+	  if (huge + x > 0.0)
+	    {
+	      /* Raise inexact if x != 0.  */
+	      i0 += 0x0000800000000000LL >> j0;
+	      i0 &= ~i;
+	      i1 = 0;
+	    }
+	}
+    }
+  else if (j0 > 111)
+    {
+      if (j0 == 0x4000)
+	/* Inf or NaN.  */
+	return x + x;
+      else
+	return x;
+    }
+  else
+    {
+      u_int64_t i = -1ULL >> (j0 - 48);
+      if ((i1 & i) == 0)
+	/* X is integral.  */
+	return x;
+
+      if (huge + x > 0.0)
+	{
+	  /* Raise inexact if x != 0.  */
+	  u_int64_t j = i1 + (1LL << (111 - j0));
+	  if (j < i1)
+	    i0 += 1;
+	  i1 = j;
+	}
+      i1 &= ~i;
+    }
+
+  SET_LDOUBLE_WORDS64 (x, i0, i1);
+  return x;
+}
+weak_alias (__roundl, roundl)
diff --git a/sysdeps/ieee754/ldbl-128/s_scalblnl.c b/sysdeps/ieee754/ldbl-128/s_scalblnl.c
new file mode 100644
index 0000000..5e8b58b
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_scalblnl.c
@@ -0,0 +1,70 @@
+/* s_scalblnl.c -- long double version of s_scalbn.c.
+ * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+   
+/* @(#)s_scalbn.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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
+
+/*
+ * scalblnl (long double x, long int n)
+ * scalblnl(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
+ * exponentiation or a multiplication.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double
+#else
+static long double
+#endif
+two114 = 2.0769187434139310514121985316880384E+34L, /* 0x4071000000000000, 0 */
+twom114 = 4.8148248609680896326399448564623183E-35L, /* 0x3F8D000000000000, 0 */
+huge   = 1.0E+4900L,
+tiny   = 1.0E-4900L;
+
+#ifdef __STDC__
+	long double __scalblnl (long double x, long int n)
+#else
+	long double __scalblnl (x,n)
+	long double x; long int n;
+#endif
+{
+	int64_t k,hx,lx;
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+        k = (hx>>48)&0x7fff;		/* extract exponent */
+        if (k==0) {				/* 0 or subnormal x */
+            if ((lx|(hx&0x7fffffffffffffffULL))==0) return x; /* +-0 */
+	    x *= two114;
+	    GET_LDOUBLE_MSW64(hx,x);
+	    k = ((hx>>48)&0x7fff) - 114;
+	}
+        if (k==0x7fff) return x+x;		/* NaN or Inf */
+        k = k+n;
+        if (n> 50000 || k > 0x7ffe)
+	  return huge*__copysignl(huge,x); /* overflow  */
+	if (n< -50000) return tiny*__copysignl(tiny,x); /*underflow*/
+        if (k > 0) 				/* normal result */
+	    {SET_LDOUBLE_MSW64(x,(hx&0x8000ffffffffffffULL)|(k<<48)); return x;}
+        if (k <= -114)
+	  return tiny*__copysignl(tiny,x); 	/*underflow*/
+        k += 114;				/* subnormal result */
+	SET_LDOUBLE_MSW64(x,(hx&0x8000ffffffffffffULL)|(k<<48));
+        return x*twom114;
+}
+weak_alias (__scalblnl, scalblnl)
diff --git a/sysdeps/ieee754/ldbl-128/s_scalbnl.c b/sysdeps/ieee754/ldbl-128/s_scalbnl.c
new file mode 100644
index 0000000..c54f064
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_scalbnl.c
@@ -0,0 +1,70 @@
+/* s_scalbnl.c -- long double version of s_scalbn.c.
+ * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+   
+/* @(#)s_scalbn.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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
+
+/*
+ * scalbnl (long double x, int n)
+ * scalbnl(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
+ * exponentiation or a multiplication.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+static const long double
+#else
+static long double
+#endif
+two114 = 2.0769187434139310514121985316880384E+34L, /* 0x4071000000000000, 0 */
+twom114 = 4.8148248609680896326399448564623183E-35L, /* 0x3F8D000000000000, 0 */
+huge   = 1.0E+4900L,
+tiny   = 1.0E-4900L;
+
+#ifdef __STDC__
+	long double __scalbnl (long double x, int n)
+#else
+	long double __scalbnl (x,n)
+	long double x; int n;
+#endif
+{
+	int64_t k,hx,lx;
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+        k = (hx>>48)&0x7fff;		/* extract exponent */
+        if (k==0) {				/* 0 or subnormal x */
+            if ((lx|(hx&0x7fffffffffffffffULL))==0) return x; /* +-0 */
+	    x *= two114;
+	    GET_LDOUBLE_MSW64(hx,x);
+	    k = ((hx>>48)&0x7fff) - 114;
+	}
+        if (k==0x7fff) return x+x;		/* NaN or Inf */
+        k = k+n;
+        if (n> 50000 || k > 0x7ffe)
+	  return huge*__copysignl(huge,x); /* overflow  */
+	if (n< -50000) return tiny*__copysignl(tiny,x); /*underflow*/
+        if (k > 0) 				/* normal result */
+	    {SET_LDOUBLE_MSW64(x,(hx&0x8000ffffffffffffULL)|(k<<48)); return x;}
+        if (k <= -114)
+	  return tiny*__copysignl(tiny,x); 	/*underflow*/
+        k += 114;				/* subnormal result */
+	SET_LDOUBLE_MSW64(x,(hx&0x8000ffffffffffffULL)|(k<<48));
+        return x*twom114;
+}
+weak_alias (__scalbnl, scalbnl)
diff --git a/sysdeps/ieee754/ldbl-128/s_signbitl.c b/sysdeps/ieee754/ldbl-128/s_signbitl.c
new file mode 100644
index 0000000..52a0afb
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_signbitl.c
@@ -0,0 +1,32 @@
+/* Return nonzero value if number is negative.
+   Copyright (C) 1997, 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+#include "math_private.h"
+
+int
+__signbitl (long double x)
+{
+  int64_t e;
+
+  GET_LDOUBLE_MSW64 (e, x);
+  return e < 0;
+}
diff --git a/sysdeps/ieee754/ldbl-128/s_sincosl.c b/sysdeps/ieee754/ldbl-128/s_sincosl.c
new file mode 100644
index 0000000..2d74e72
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_sincosl.c
@@ -0,0 +1,77 @@
+/* Compute sine and cosine of argument.
+   Copyright (C) 1997, 1999 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>.
+
+   The GNU C Library 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.
+
+   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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+#include "math_private.h"
+
+/* Note: We should probably introduce __kernel_sincosl to speed things up,
+   because __kernel_{cos,sin}l sometimes compute both sine and cosine.  */
+
+void
+__sincosl (long double x, long double *sinx, long double *cosx)
+{
+  int64_t ix;
+
+  /* High word of x. */
+  GET_LDOUBLE_MSW64 (ix, x);
+
+  /* |x| ~< pi/4 */
+  ix &= 0x7fffffffffffffffLL;
+  if (ix <= 0x3ffe921fb54442d1LL)
+    {
+      *sinx = __kernel_sinl (x, 0.0, 0);
+      *cosx = __kernel_cosl (x, 0.0);
+    }
+  else if (ix >= 0x7fff000000000000LL)
+    {
+      /* sin(Inf or NaN) is NaN */
+      *sinx = *cosx = x - x;
+    }
+  else
+    {
+      /* Argument reduction needed.  */
+      long double y[2];
+      int n;
+
+      n = __ieee754_rem_pio2l (x, y);
+      switch (n & 3)
+	{
+	case 0:
+	  *sinx = __kernel_sinl (y[0], y[1], 1);
+	  *cosx = __kernel_cosl (y[0], y[1]);
+	  break;
+	case 1:
+	  *sinx = __kernel_cosl (y[0], y[1]);
+	  *cosx = -__kernel_sinl (y[0], y[1], 1);
+	  break;
+	case 2:
+	  *sinx = -__kernel_sinl (y[0], y[1], 1);
+	  *cosx = -__kernel_cosl (y[0], y[1]);
+	  break;
+	default:
+	  *sinx = -__kernel_cosl (y[0], y[1]);
+	  *cosx = __kernel_sinl (y[0], y[1], 1);
+	  break;
+	}
+    }
+}
+weak_alias (__sincosl, sincosl)
diff --git a/sysdeps/ieee754/ldbl-128/s_sinl.c b/sysdeps/ieee754/ldbl-128/s_sinl.c
new file mode 100644
index 0000000..446a75f
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_sinl.c
@@ -0,0 +1,83 @@
+/* s_sinl.c -- long double version of s_sin.c.
+ * Conversion to 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.
+ * ====================================================
+ */
+
+/* sinl(x)
+ * Return sine function of x.
+ *
+ * kernel function:
+ *	__kernel_sinl		... sine function on [-pi/4,pi/4]
+ *	__kernel_cosl		... cose function on [-pi/4,pi/4]
+ *	__ieee754_rem_pio2l	... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *	in [-pi/4 , +pi/4], and let n = k mod 4.
+ *	We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *	    0	       S	   C		 T
+ *	    1	       C	  -S		-1/T
+ *	    2	      -S	  -C		 T
+ *	    3	      -C	   S		-1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *	TRIG(x) returns trig(x) nearly rounded
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	long double __sinl(long double x)
+#else
+	long double __sinl(x)
+	long double x;
+#endif
+{
+	long double y[2],z=0.0L;
+	int64_t n, ix;
+
+    /* High word of x. */
+	GET_LDOUBLE_MSW64(ix,x);
+
+    /* |x| ~< pi/4 */
+	ix &= 0x7fffffffffffffffLL;
+	if(ix <= 0x3ffe921fb54442d1LL)
+	  return __kernel_sinl(x,z,0);
+
+    /* sin(Inf or NaN) is NaN */
+	else if (ix>=0x7fff000000000000LL) return x-x;
+
+    /* argument reduction needed */
+	else {
+	    n = __ieee754_rem_pio2l(x,y);
+	    switch(n&3) {
+		case 0: return  __kernel_sinl(y[0],y[1],1);
+		case 1: return  __kernel_cosl(y[0],y[1]);
+		case 2: return -__kernel_sinl(y[0],y[1],1);
+		default:
+			return -__kernel_cosl(y[0],y[1]);
+	    }
+	}
+}
+weak_alias (__sinl, sinl)
diff --git a/sysdeps/ieee754/ldbl-128/s_tanl.c b/sysdeps/ieee754/ldbl-128/s_tanl.c
new file mode 100644
index 0000000..ea9d053
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_tanl.c
@@ -0,0 +1,77 @@
+/* s_tanl.c -- long double version of s_tan.c.
+ * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
+ */
+  
+/* @(#)s_tan.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/* tanl(x)
+ * Return tangent function of x.
+ *
+ * kernel function:
+ *	__kernel_tanl		... tangent function on [-pi/4,pi/4]
+ *	__ieee754_rem_pio2l	... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *	in [-pi/4 , +pi/4], and let n = k mod 4.
+ *	We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *	    0	       S	   C		 T
+ *	    1	       C	  -S		-1/T
+ *	    2	      -S	  -C		 T
+ *	    3	      -C	   S		-1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *	TRIG(x) returns trig(x) nearly rounded
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __STDC__
+	long double __tanl(long double x)
+#else
+	long double __tanl(x)
+	long double x;
+#endif
+{
+	long double y[2],z=0.0L;
+	int64_t n, ix;
+
+    /* High word of x. */
+	GET_LDOUBLE_MSW64(ix,x);
+
+    /* |x| ~< pi/4 */
+	ix &= 0x7fffffffffffffffLL;
+	if(ix <= 0x3ffe921fb54442d1LL) return __kernel_tanl(x,z,1);
+
+    /* tanl(Inf or NaN) is NaN */
+	else if (ix>=0x7fff000000000000LL) return x-x;		/* NaN */
+
+    /* argument reduction needed */
+	else {
+	    n = __ieee754_rem_pio2l(x,y);
+	    return __kernel_tanl(y[0],y[1],1-((n&1)<<1)); /*   1 -- n even
+							-1 -- n odd */
+	}
+}
+weak_alias (__tanl, tanl)
diff --git a/sysdeps/ieee754/ldbl-128/s_truncl.c b/sysdeps/ieee754/ldbl-128/s_truncl.c
new file mode 100644
index 0000000..bbff5a4
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/s_truncl.c
@@ -0,0 +1,57 @@
+/* Truncate argument to nearest integral value not larger than the argument.
+   Copyright (C) 1997, 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+#include "math_private.h"
+
+
+long double
+__truncl (long double x)
+{
+  int32_t j0;
+  u_int64_t i0, i1, sx;
+
+  GET_LDOUBLE_WORDS64 (i0, i1, x);
+  sx = i0 & 0x8000000000000000ULL;
+  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
+  if (j0 < 48)
+    {
+      if (j0 < 0)
+	/* The magnitude of the number is < 1 so the result is +-0.  */
+	SET_LDOUBLE_WORDS64 (x, sx, 0);
+      else
+	SET_LDOUBLE_WORDS64 (x, i0 & ~(0x0000ffffffffffffLL >> j0), 0);
+    }
+  else if (j0 > 111)
+    {
+      if (j0 == 0x4000)
+	/* x is inf or NaN.  */
+	return x + x;
+    }
+  else
+    {
+      SET_LDOUBLE_WORDS64 (x, i0, i1 & ~(0xffffffffffffffffULL >> (j0 - 48)));
+    }
+
+  return x;
+}
+weak_alias (__truncl, truncl)
diff --git a/sysdeps/ieee754/ldbl-128/strtold.c b/sysdeps/ieee754/ldbl-128/strtold.c
new file mode 100644
index 0000000..32049fc
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/strtold.c
@@ -0,0 +1,42 @@
+/* Copyright (C) 1999 Free Software Foundation, Inc.
+
+   The GNU C Library 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.
+
+   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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C 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 <math.h>
+
+/* The actual implementation for all floating point sizes is in strtod.c.
+   These macros tell it to produce the `long double' version, `strtold'.  */
+
+# define FLOAT		long double
+# define FLT		LDBL
+# ifdef USE_IN_EXTENDED_LOCALE_MODEL
+#  define STRTOF	__strtold_l
+# else
+#  define STRTOF	strtold
+# endif
+# define MPN2FLOAT	__mpn_construct_long_double
+# define FLOAT_HUGE_VAL	HUGE_VALL
+# define SET_MANTISSA(flt, mant) \
+  do { union ieee854_long_double u;					      \
+       u.d = (flt);							      \
+       u.ieee.mantissa0 = 0x8000;					      \
+       u.ieee.mantissa1 = 0;						      \
+       u.ieee.mantissa2 = ((mant) >> 32);	      			      \
+       u.ieee.mantissa3 = (mant) & 0xffffffff;				      \
+       (flt) = u.d;							      \
+  } while (0)
+
+# include "strtod.c"