diff --git a/library/changes b/library/changes
index 16daba6..91fc3be 100644
--- a/library/changes
+++ b/library/changes
@@ -6,6 +6,12 @@
   in Ghostscript. Wish I managed to track this down yesterday for 1.196
   release... <aantonijevic>
 
+- Implemented atanh() and atanhf() which were not listed in the TODO
+  file but were still unimplemented up until now.
+
+- Replaced ldexp() and modf().
+
+
 c.lib 1.196 (11.10.2005)
 
 - Removed the various workarounds associated with <unistd.h>, required
diff --git a/library/include/math.h b/library/include/math.h
index 6cb9a50..4946753 100644
--- a/library/include/math.h
+++ b/library/include/math.h
@@ -1,5 +1,5 @@
 /*
- * $Id: math.h,v 1.17 2005-10-09 10:38:56 obarthel Exp $
+ * $Id: math.h,v 1.18 2005-10-16 09:05:03 obarthel Exp $
  *
  * :ts=4
  *
@@ -244,6 +244,7 @@ extern float tanhf(float x);
 
 extern float acoshf(float x);
 extern float asinhf(float x);
+extern float atanhf(float x);
 extern float cbrtf(float x);
 extern float copysignf(float x, float y);
 extern float erfcf(float x);
@@ -275,6 +276,7 @@ extern int ilogbf(float x);
 
 extern double acosh(double x);
 extern double asinh(double x);
+extern double atanh(double x);
 extern double cbrt(double x);
 extern double copysign(double x, double y);
 extern double erf(double x);
diff --git a/library/math_atanh.c b/library/math_atanh.c
index ed16784..6c93847 100644
--- a/library/math_atanh.c
+++ b/library/math_atanh.c
@@ -1,5 +1,5 @@
 /*
- * $Id: math_atanh.c,v 1.1 2005-05-29 11:19:00 obarthel Exp $
+ * $Id: math_atanh.c,v 1.2 2005-10-16 09:05:02 obarthel Exp $
  *
  * :ts=4
  *
@@ -29,6 +29,15 @@
  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  * POSSIBILITY OF SUCH DAMAGE.
+ *
+ *
+ * PowerPC math library based in part on work by Sun Microsystems
+ * 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.
  */
 
 #ifndef _MATH_HEADERS_H
@@ -41,11 +50,31 @@
 
 /****************************************************************************/
 
+static const double zero = 0.0;
+static const double one = 1.0, huge = 1e300;
+
+/****************************************************************************/
+
 double
 atanh(double x)
 {
-	/* ZZZ unimplemented */
-	return(0);
+	double t;
+	LONG hx,ix;
+	ULONG lx;
+	EXTRACT_WORDS(hx,lx,x);
+	ix = hx&0x7fffffff;
+	if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
+	    return (x-x)/(x-x);
+	if(ix==0x3ff00000) 
+	    return x/zero;
+	if(ix<0x3e300000&&(huge+x)>zero) return x;	/* x<2**-28 */
+	SET_HIGH_WORD(x,ix);
+	if(ix<0x3fe00000) {		/* x < 0.5 */
+	    t = x+x;
+	    t = 0.5*log1p(t+t*x/(one-x));
+	} else 
+	    t = 0.5*log1p((x+x)/(one-x));
+	if(hx>=0) return t; else return -t;
 }
 
 /****************************************************************************/
diff --git a/library/math_atanhf.c b/library/math_atanhf.c
index 0eb97aa..baf759c 100644
--- a/library/math_atanhf.c
+++ b/library/math_atanhf.c
@@ -1,5 +1,5 @@
 /*
- * $Id: math_atanhf.c,v 1.1 2005-05-29 11:19:00 obarthel Exp $
+ * $Id: math_atanhf.c,v 1.2 2005-10-16 09:05:02 obarthel Exp $
  *
  * :ts=4
  *
@@ -29,6 +29,17 @@
  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  * POSSIBILITY OF SUCH DAMAGE.
+ *
+ *
+ * PowerPC math library based in part on work by Sun Microsystems
+ * 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.
+ *
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
  */
 
 #ifndef _MATH_HEADERS_H
@@ -41,11 +52,29 @@
 
 /****************************************************************************/
 
-float
-atanhf(float x)
+static const float zero = 0.0;
+static const float one = 1.0, huge = 1e30;
+
+/****************************************************************************/
+
+float atanhf(float x)
 {
-	/* ZZZ unimplemented */
-	return(0);
+	float t;
+	LONG hx,ix;
+	GET_FLOAT_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if (ix>0x3f800000) 		/* |x|>1 */
+	    return (x-x)/(x-x);
+	if(ix==0x3f800000) 
+	    return x/zero;
+	if(ix<0x31800000&&(huge+x)>zero) return x;	/* x<2**-28 */
+	SET_FLOAT_WORD(x,ix);
+	if(ix<0x3f000000) {		/* x < 0.5 */
+	    t = x+x;
+	    t = (float)0.5*log1pf(t+t*x/(one-x));
+	} else 
+	    t = (float)0.5*log1pf((x+x)/(one-x));
+	if(hx>=0) return t; else return -t;
 }
 
 /****************************************************************************/
diff --git a/library/math_ldexp.c b/library/math_ldexp.c
index 5d0745b..42f8b79 100644
--- a/library/math_ldexp.c
+++ b/library/math_ldexp.c
@@ -1,5 +1,5 @@
 /*
- * $Id: math_ldexp.c,v 1.3 2005-02-25 10:14:21 obarthel Exp $
+ * $Id: math_ldexp.c,v 1.4 2005-10-16 09:05:02 obarthel Exp $
  *
  * :ts=4
  *
@@ -50,122 +50,22 @@
 
 /****************************************************************************/
 
-#if defined(IEEE_FLOATING_POINT_SUPPORT)
-
-#define MANT_MASK	0x800FFFFF	/* Mantissa extraction mask */
-#define ZPOS_MASK	0x3FF00000	/* Positive # mask for exp = 0 */
-#define ZNEG_MASK	0x3FF00000	/* Negative # mask for exp = 0 */
-
-#define EXP_MASK	0x7FF00000	/* Mask for exponent */
-#define EXP_SHIFTS	20			/* Shifts to get into LSB's */
-#define EXP_BIAS	1023		/* Exponent bias */
-
-union dtol
-{
-	double	dval;
-	long	ival[2];
-};
-
-INLINE STATIC const double
-__ldexp(double x,int n)
-{
-	union dtol number;
-	long *iptr, cn;
-
-	number.dval = x;
-
-	iptr = &number.ival[0];
-
-	cn = (((*iptr) & EXP_MASK) >> EXP_SHIFTS) - EXP_BIAS;
-
-	(*iptr) &= ~EXP_MASK;
-
-	n += EXP_BIAS;
-
-	/* ZZZ we can't just muck with the exponent, we
-	 * have to check for underflow and overflow, too!
-	 */
-	(*iptr) |= ((n + cn) << EXP_SHIFTS) & EXP_MASK;
-
-	return(number.dval);
-}
-
-#endif /* IEEE_FLOATING_POINT_SUPPORT */
-
-/****************************************************************************/
-
-#if defined(M68881_FLOATING_POINT_SUPPORT)
-
-INLINE STATIC const double
-__ldexp(double x,int n)
-{
-	double result;
-
-	__asm ("fscale%.l %2,%0"
-	       : "=f" (result)
-	       : "0" (x),
-	       "dmi" (n));
-
-	return(result);
-}
-
-#endif /* M68881_FLOATING_POINT_SUPPORT */
-
-/****************************************************************************/
-
-#if defined(PPC_FLOATING_POINT_SUPPORT)
-
-#define MANT_MASK	0x800FFFFF	/* Mantissa extraction mask */
-#define ZPOS_MASK	0x3FF00000	/* Positive # mask for exp = 0 */
-#define ZNEG_MASK	0x3FF00000	/* Negative # mask for exp = 0 */
-
-#define EXP_MASK	0x7FF00000	/* Mask for exponent */
-#define EXP_SHIFTS	20			/* Shifts to get into LSB's */
-#define EXP_BIAS	1023		/* Exponent bias */
-
-union dtol
-{
-	double	dval;
-	long	ival[2];
-};
-
-INLINE STATIC const double
-__ldexp(double x,int n)
-{
-	union dtol number;
-	long *iptr, cn;
-
-	number.dval = x;
-
-	iptr = &number.ival[0];
-
-	cn = (((*iptr) & EXP_MASK) >> EXP_SHIFTS) - EXP_BIAS;
-
-	(*iptr) &= ~EXP_MASK;
-
-	n += EXP_BIAS;
-
-	/* ZZZ we can't just muck with the exponent, we
-	 * have to check for underflow and overflow, too!
-	 */
-	(*iptr) |= ((n + cn) << EXP_SHIFTS) & EXP_MASK;
-
-	return(number.dval);
-}
-
-#endif /* PPC_FLOATING_POINT_SUPPORT */
-
-/****************************************************************************/
-
 double
-ldexp(double x,int n)
+ldexp(double value, int exp)
 {
 	double result;
 
-	if(x != 0.0)
-		result = __ldexp(x,n);
+	if(isinf(x) || fpclassify(x) == FP_ZERO)
+	{
+		result = x;
+	}
 	else
-		result = 0.0;
+	{
+		result = scalbn(x,exp);
+
+		if(isinf(result) || (result < DBL_MIN || result > -DBL_MIN))
+			__set_errno(ERANGE);
+	}
 
 	return(result);
 }
diff --git a/library/math_modf.c b/library/math_modf.c
index 983d1cd..26e02d9 100644
--- a/library/math_modf.c
+++ b/library/math_modf.c
@@ -1,5 +1,5 @@
 /*
- * $Id: math_modf.c,v 1.5 2005-02-25 10:14:21 obarthel Exp $
+ * $Id: math_modf.c,v 1.6 2005-10-16 09:05:02 obarthel Exp $
  *
  * :ts=4
  *
@@ -56,38 +56,10 @@
 
 /****************************************************************************/
 
-#if defined(IEEE_FLOATING_POINT_SUPPORT)
-
-INLINE STATIC const double
-__modf(double x,double *nptr)
-{
-	double int_n;
-	double result;
-
-	if(x < 0.0)
-	{
-		int_n = ceil(x);
-
-		result = int_n - x;
-	}
-	else
-	{
-		int_n = floor(x);
-
-		result = x - int_n;
-	}
-
-	(*nptr) = int_n;
-
-	return(result);
-}
-
-#endif /* IEEE_FLOATING_POINT_SUPPORT */
+#if defined(M68881_FLOATING_POINT_SUPPORT)
 
 /****************************************************************************/
 
-#if defined(M68881_FLOATING_POINT_SUPPORT)
-
 INLINE STATIC const double
 __modf(double x,double *nptr)
 {
@@ -105,37 +77,62 @@ __modf(double x,double *nptr)
 	return(result);
 }
 
-#endif /* M68881_FLOATING_POINT_SUPPORT */
+/****************************************************************************/
+
+#else
 
 /****************************************************************************/
 
-#if defined(__PPC__)
+static const double one = 1.0;
 
-INLINE STATIC const double
-__modf(double x,double *nptr)
+INLINE STATIC double
+__modf(double x, double *iptr)
 {
-	double int_n;
-	double result;
-
-	if(x < 0.0)
-	{
-		int_n = ceil(x);
-
-		result = int_n - x;
+	LONG i0,i1,j0;
+	ULONG i;
+	EXTRACT_WORDS(i0,i1,x);
+	j0 = ((i0>>20)&0x7ff)-0x3ff;	/* exponent of x */
+	if(j0<20) {			/* integer part in high x */
+	    if(j0<0) {			/* |x|<1 */
+	        INSERT_WORDS(*iptr,i0&0x80000000,0);	/* *iptr = +-0 */
+		return x;
+	    } else {
+		i = (0x000fffff)>>j0;
+		if(((i0&i)|i1)==0) {		/* x is integral */
+		    ULONG high;
+		    *iptr = x;
+		    GET_HIGH_WORD(high,x);
+		    INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
+		    return x;
+		} else {
+		    INSERT_WORDS(*iptr,i0&(~i),0);
+		    return x - *iptr;
+		}
+	    }
+	} else if (j0>51) {		/* no fraction part */
+	    ULONG high;
+	    *iptr = x*one;
+	    GET_HIGH_WORD(high,x);
+	    INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
+	    return x;
+	} else {			/* fraction part in low x */
+	    i = ((ULONG)(0xffffffff))>>(j0-20);
+	    if((i1&i)==0) { 		/* x is integral */
+	        ULONG high;
+		*iptr = x;
+		GET_HIGH_WORD(high,x);
+		INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
+		return x;
+	    } else {
+	        INSERT_WORDS(*iptr,i0,i1&(~i));
+		return x - *iptr;
+	    }
 	}
-	else
-	{
-		int_n = floor(x);
-
-		result = x - int_n;
-	}
-
-	(*nptr) = int_n;
-
-	return(result);
 }
 
-#endif /* PPC_FLOATING_POINT_SUPPORT */
+/****************************************************************************/
+
+#endif /* M68881_FLOATING_POINT_SUPPORT */
 
 /****************************************************************************/