- noieee_src: adapt complex versions of the functions it already

supports (which is all, except the float ones)
ok millert@

29 files changed, 1794 insertions, 229 deletions

diff --git a/lib/libm/Makefile b/lib/libm/Makefile
index a1d109a8392..90e669abc91 100644
--- a/lib/libm/Makefile
+++ b/lib/libm/Makefile
@@ -1,5 +1,5 @@
 #  $NetBSD: Makefile,v 1.28 1995/11/20 22:06:19 jtc Exp $
-#  $OpenBSD: Makefile,v 1.57 2008/09/16 22:06:24 martynas Exp $
+#  $OpenBSD: Makefile,v 1.58 2008/10/07 22:25:53 martynas Exp $
 #
 #  @(#)Makefile 5.1beta 93/09/24 
 # 
@@ -54,8 +54,8 @@ ARCH_SRCS = e_sqrt.c e_sqrtf.c e_remainder.c e_remainderf.c \
 	    s_ceil.c s_ceilf.c s_floor.c s_floorf.c s_rint.c s_rintf.c
 .elif (${MACHINE_ARCH} == "vax")
 .PATH:	${.CURDIR}/arch/vax
-NOIEEE_ARCH=n_infnan.S n_argred.S n_sqrt.S
-ARCH_SRCS = n_atan2.S n_cabs.S n_cbrt.S n_sincos.S n_tan.S n_support.S 
+NOIEEE_ARCH = n_argred.S n_infnan.S n_sqrt.S
+ARCH_SRCS = n_atan2.S n_cbrt.S n_hypot.S n_sincos.S n_support.S n_tan.S
 .endif
 
 .PATH:	${.CURDIR}/man
@@ -100,14 +100,17 @@ COMMON_SRCS = b_exp__D.c b_log__D.c b_tgamma.c \
 	w_gammaf_r.c w_lgamma.c w_lgammaf.c
 
 # math routines for non-IEEE architectures.
-NOIEEE_SRCS = n_asincos.c n_acosh.c n_asinh.c n_atan.c n_atanh.c n_cosh.c \
-	n_erf.c n_exp.c n_exp__E.c n_expm1.c n_fdim.c n_floor.c n_fmax.c \
-	n_fmaxf.c n_fmin.c n_fminf.c n_fmod.c n_tgamma.c \
-	n_lgamma.c n_j0.c n_j1.c n_jn.c n_log.c n_log10.c n_log1p.c \
-	n_log__L.c n_pow.c n_sinh.c n_tanh.c n_atan2.c n_cabs.c n_cbrt.c \
-	n_sqrt.c n_sincos.c n_tan.c n_argred.c n_support.c n_infnan.c \
-	n_round.c
-
+NOIEEE_SRCS = n_acosh.c n_argred.c n_asincos.c n_asinh.c n_atan.c \
+	n_atan2.c n_atanh.c n_cabs.c n_cacos.c n_cacosh.c n_carg.c \
+	n_casin.c n_casinh.c n_catan.c n_catanh.c n_cbrt.c n_ccos.c \
+	n_ccosh.c n_cexp.c n_cimag.c n_cimagf.c n_clog.c n_conj.c \
+	n_conjf.c n_cosh.c n_cpow.c n_cproj.c n_cprojf.c n_creal.c \
+	n_crealf.c n_csin.c n_csinh.c n_csqrt.c n_ctan.c n_ctanh.c \
+	n_erf.c n_exp.c n_exp__E.c n_expm1.c n_fdim.c n_floor.c \
+	n_fmax.c n_fmaxf.c n_fmin.c n_fminf.c n_fmod.c n_hypot.c \
+	n_infnan.c n_j0.c n_j1.c n_jn.c n_lgamma.c n_log.c n_log10.c \
+	n_log1p.c n_log__L.c n_pow.c n_round.c n_sincos.c n_sinh.c \
+	n_sqrt.c n_support.c n_tan.c n_tanh.c n_tgamma.c
 
 # OpenBSD's C library supplies these functions:
 #COMMON_SRCS+=	s_fabs.c s_frexp.c s_ldexp.c s_modf.c
diff --git a/lib/libm/arch/vax/n_cabs.S b/lib/libm/arch/vax/n_hypot.S
index 02f01ae8277..bd5cbf49df1 100644
--- a/lib/libm/arch/vax/n_cabs.S
+++ b/lib/libm/arch/vax/n_hypot.S
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_cabs.S,v 1.4 2008/09/13 21:25:40 martynas Exp $	*/
+/*	$OpenBSD: n_hypot.S,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
 /*	$NetBSD: n_cabs.S,v 1.1 1995/10/10 23:40:26 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -45,8 +45,7 @@
 /*	entry for c functions cabs and hypot */
 	.text
 	_ALIGN_TEXT
-ALTENTRY(hypot)
-ENTRY(cabs, 0x8000|R2|R3|R4|R5|R6)	# enable floating overflow
+ENTRY(hypot, 0x8000|R2|R3|R4|R5|R6)	# enable floating overflow
 	movq	4(ap),r0	# r0:1 = x
 	movq	12(ap),r2	# r2:3 = y
 	bicw3	$0x7f,r0,r4	# r4 has signed biased exp of x
diff --git a/lib/libm/noieee_src/n_cabs.c b/lib/libm/noieee_src/n_cabs.c
index db0326c5271..c07c1891945 100644
--- a/lib/libm/noieee_src/n_cabs.c
+++ b/lib/libm/noieee_src/n_cabs.c
@@ -1,222 +1,25 @@
-/*	$OpenBSD: n_cabs.c,v 1.11 2008/09/13 21:25:40 martynas Exp $	*/
-/*	$NetBSD: n_cabs.c,v 1.1 1995/10/10 23:36:39 ragge Exp $	*/
+/*	$OpenBSD: n_cabs.c,v 1.12 2008/10/07 22:25:53 martynas Exp $	*/
 /*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#ifndef lint
-static char sccsid[] = "@(#)cabs.c	8.1 (Berkeley) 6/4/93";
-#endif /* not lint */
-
-/* HYPOT(X,Y)
- * RETURN THE SQUARE ROOT OF X^2 + Y^2  WHERE Z=X+iY
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 11/28/84;
- * REVISED BY K.C. NG, 7/12/85.
- *
- * Required system supported functions :
- *	copysign(x,y)
- *	finite(x)
- *	scalbn(x,N)
- *	sqrt(x)
- *
- * Method :
- *	1. replace x by |x| and y by |y|, and swap x and
- *	   y if y > x (hence x is never smaller than y).
- *	2. Hypot(x,y) is computed by:
- *	   Case I, x/y > 2
- *
- *				       y
- *		hypot = x + -----------------------------
- *			 		    2
- *			    sqrt ( 1 + [x/y]  )  +  x/y
- *
- *	   Case II, x/y <= 2
- *				                   y
- *		hypot = x + --------------------------------------------------
- *				          		     2
- *				     			[x/y]   -  2
- *			   (sqrt(2)+1) + (x-y)/y + -----------------------------
- *			 		    			  2
- *			    			  sqrt ( 1 + [x/y]  )  + sqrt(2)
- *
- *
- *
- * Special cases:
- *	hypot(x,y) is INF if x or y is +INF or -INF; else
- *	hypot(x,y) is NAN if x or y is NAN.
- *
- * Accuracy:
- * 	hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units
- *	in the last place). See Kahan's "Interval Arithmetic Options in the
- *	Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
- *      1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
- *	code follows in	comments.) In a test run with 500,000 random arguments
- *	on a VAX, the maximum observed error was .959 ulps.
- *
- * 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.
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  */
 
-#include "math.h"
-#include "mathimpl.h"
-
-vc(r2p1hi, 2.4142135623730950345E0   ,8279,411a,ef32,99fc,   2, .9A827999FCEF32)
-vc(r2p1lo, 1.4349369327986523769E-17 ,597d,2484,754b,89b3, -55, .84597D89B3754B)
-vc(sqrt2,  1.4142135623730950622E0   ,04f3,40b5,de65,33f9,   1, .B504F333F9DE65)
-
-ic(r2p1hi, 2.4142135623730949234E0   ,   1, 1.3504F333F9DE6)
-ic(r2p1lo, 1.2537167179050217666E-16 , -53, 1.21165F626CDD5)
-ic(sqrt2,  1.4142135623730951455E0   ,   0, 1.6A09E667F3BCD)
-
-#ifdef vccast
-#define	r2p1hi	vccast(r2p1hi)
-#define	r2p1lo	vccast(r2p1lo)
-#define	sqrt2	vccast(sqrt2)
-#endif
+#include <complex.h>
+#include <math.h>
 
 double
-hypot(double x, double y)
+cabs(double complex z)
 {
-	static const double zero=0, one=1,
-		      small=1.0E-18;	/* fl(1+small)==1 */
-	static const ibig=30;	/* fl(1+2**(2*ibig))==1 */
-	double t,r;
-	int exp;
-
-	if(finite(x))
-	    if(finite(y))
-	    {
-		x=copysign(x,one);
-		y=copysign(y,one);
-		if(y > x)
-		    { t=x; x=y; y=t; }
-		if(x == zero) return(zero);
-		if(y == zero) return(x);
-		exp= logb(x);
-		if (exp - (int)logb(y) > ibig) {
-			if (one + small >= 1.0)	/* raise inexact flag */
-				return(x);	/* return |x| */
-		}
-
-	    /* start computing sqrt(x^2 + y^2) */
-		r=x-y;
-		if(r>y) { 	/* x/y > 2 */
-		    r=x/y;
-		    r=r+sqrt(one+r*r); }
-		else {		/* 1 <= x/y <= 2 */
-		    r/=y; t=r*(r+2.0);
-		    r+=t/(sqrt2+sqrt(2.0+t));
-		    r+=r2p1lo; r+=r2p1hi; }
-
-		r=y/r;
-		return(x+r);
-
-	    }
-
-	    else if(isinf(y))		/* y is +-INF */
-		     return(copysign(y,one));
-	    else
-		     return(y);		/* y is NaN and x is finite */
-
-	else if(isinf(x))		/* x is +-INF */
-	         return (copysign(x,one));
-	else if(finite(y))
-	         return(x);		/* x is NaN, y is finite */
-	else if (isnan(y))
-		return (y);
-	else return(copysign(y,one));	/* y is INF */
-}
-
-/* CABS(Z)
- * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER  Z = X + iY
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 11/28/84.
- * REVISED BY K.C. NG, 7/12/85.
- *
- * Required kernel function :
- *	hypot(x,y)
- *
- * Method :
- *	cabs(z) = hypot(x,y) .
- */
-
-struct complex { double x, y; };
-
-double
-cabs(struct complex z)
-{
-	return hypot(z.x,z.y);
-}
-
-/* A faster but less accurate version of cabs(x,y) */
-#if 0
-double
-hypot(double x, double y)
-{
-	static const double zero=0, one=1;
-		      small=1.0E-18;	/* fl(1+small)==1 */
-	static const ibig=30;	/* fl(1+2**(2*ibig))==1 */
-	double temp;
-	int exp;
-
-	if(finite(x))
-	    if(finite(y))
-	    {
-		x=copysign(x,one);
-		y=copysign(y,one);
-		if(y > x)
-		    { temp=x; x=y; y=temp; }
-		if(x == zero) return(zero);
-		if(y == zero) return(x);
-		exp= logb(x);
-		x=scalbn(x,-exp);
-		if (exp - (int)logb(y) > ibig) {
-			if (one + small >= 1.0)		/* raise inexact flag */
-				return(scalbn(x,exp));	/* return |x| */
-		}
-		else y=scalbn(y,-exp);
-		return(scalbn(sqrt(x*x+y*y),exp));
-	    }
-
-	    else if(isinf(y))		/* y is +-INF */
-		     return(copysign(y,one));
-	    else
-		     return(y);		/* y is NaN and x is finite */
-
-	else if(isinf(x))		/* x is +-INF */
-	         return (copysign(x,one));
-	else if(finite(y))
-	         return(x);		/* x is NaN, y is finite */
-	else if(isnan(y)) return(y);	/* x and y is NaN */
-	else return(copysign(y,one));	/* y is INF */
+	return hypot(__real__ z, __imag__ z);
 }
-#endif
diff --git a/lib/libm/noieee_src/n_cacos.c b/lib/libm/noieee_src/n_cacos.c
new file mode 100644
index 00000000000..929aa0fa8b3
--- /dev/null
+++ b/lib/libm/noieee_src/n_cacos.c
@@ -0,0 +1,60 @@
+/*	$OpenBSD: n_cacos.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							cacos()
+ *
+ *	Complex circular arc cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex cacos();
+ * double complex z, w;
+ *
+ * w = cacos (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ *
+ * w = arccos z  =  PI/2 - arcsin z.
+ *
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10      5200      1.6e-15      2.8e-16
+ *    IEEE      -10,+10     30000      1.8e-14      2.2e-15
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+cacos(double complex z)
+{
+	double complex w;
+
+	w = casin (z);
+	w = (M_PI_2 - creal (w)) - cimag (w) * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_cacosh.c b/lib/libm/noieee_src/n_cacosh.c
new file mode 100644
index 00000000000..246aadab4ed
--- /dev/null
+++ b/lib/libm/noieee_src/n_cacosh.c
@@ -0,0 +1,55 @@
+/*	$OpenBSD: n_cacosh.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							cacosh
+ *
+ *	Complex inverse hyperbolic cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex cacosh();
+ * double complex z, w;
+ *
+ * w = cacosh (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * acosh z = i acos z .
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       1.6e-14     2.1e-15
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+cacosh(double complex z)
+{
+	double complex w;
+
+	w = I * cacos (z);
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_carg.c b/lib/libm/noieee_src/n_carg.c
new file mode 100644
index 00000000000..2d8033b0f07
--- /dev/null
+++ b/lib/libm/noieee_src/n_carg.c
@@ -0,0 +1,25 @@
+/*	$OpenBSD: n_carg.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double
+carg(double complex z)
+{
+	return atan2 (__imag__ z, __real__ z);
+}
diff --git a/lib/libm/noieee_src/n_casin.c b/lib/libm/noieee_src/n_casin.c
new file mode 100644
index 00000000000..cebc2ffe9bf
--- /dev/null
+++ b/lib/libm/noieee_src/n_casin.c
@@ -0,0 +1,129 @@
+/*	$OpenBSD: n_casin.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							casin()
+ *
+ *	Complex circular arc sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex casin();
+ * double complex z, w;
+ *
+ * w = casin (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Inverse complex sine:
+ *
+ *                               2
+ * w = -i clog( iz + csqrt( 1 - z ) ).
+ *
+ * casin(z) = -i casinh(iz)
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10     10100       2.1e-15     3.4e-16
+ *    IEEE      -10,+10     30000       2.2e-14     2.7e-15
+ * Larger relative error can be observed for z near zero.
+ * Also tested by csin(casin(z)) = z.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+casin(double complex z)
+{
+	double complex w;
+	static double complex ca, ct, zz, z2;
+	double x, y;
+
+	x = creal (z);
+	y = cimag (z);
+
+	if (y == 0.0) {
+		if (fabs(x) > 1.0) {
+			w = M_PI_2 + 0.0 * I;
+			/*mtherr ("casin", DOMAIN);*/
+		}
+		else {
+			w = asin (x) + 0.0 * I;
+		}
+		return (w);
+	}
+
+	/* Power series expansion */
+	/*
+	b = cabs(z);
+	if( b < 0.125 ) {
+		z2.r = (x - y) * (x + y);
+		z2.i = 2.0 * x * y;
+
+		cn = 1.0;
+		n = 1.0;
+		ca.r = x;
+		ca.i = y;
+		sum.r = x;
+		sum.i = y;
+		do {
+			ct.r = z2.r * ca.r  -  z2.i * ca.i;
+			ct.i = z2.r * ca.i  +  z2.i * ca.r;
+			ca.r = ct.r;
+			ca.i = ct.i;
+
+			cn *= n;
+			n += 1.0;
+			cn /= n;
+			n += 1.0;
+			b = cn/n;
+
+			ct.r *= b;
+			ct.i *= b;
+			sum.r += ct.r;
+			sum.i += ct.i;
+			b = fabs(ct.r) + fabs(ct.i);
+		}
+		while( b > MACHEP );
+		w->r = sum.r;
+		w->i = sum.i;
+		return;
+	}
+	*/
+
+	ca = x + y * I;
+	ct = ca * I;
+	/* sqrt( 1 - z*z) */
+	/* cmul( &ca, &ca, &zz ) */
+	/*x * x  -  y * y */
+	zz = (x - y) * (x + y) + (2.0 * x * y) * I;
+
+	zz = 1.0 - creal(zz) - cimag(zz) * I;
+	z2 = csqrt (zz);
+
+	zz = ct + z2;
+	zz = clog (zz);
+	/* multiply by 1/i = -i */
+	w = zz * (-1.0 * I);
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_casinh.c b/lib/libm/noieee_src/n_casinh.c
new file mode 100644
index 00000000000..ee3df8af99a
--- /dev/null
+++ b/lib/libm/noieee_src/n_casinh.c
@@ -0,0 +1,55 @@
+/*	$OpenBSD: n_casinh.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							casinh
+ *
+ *	Complex inverse hyperbolic sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex casinh();
+ * double complex z, w;
+ *
+ * w = casinh (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * casinh z = -i casin iz .
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       1.8e-14     2.6e-15
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+casinh(double complex z)
+{
+	double complex w;
+
+	w = -1.0 * I * casin (z * I);
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_catan.c b/lib/libm/noieee_src/n_catan.c
new file mode 100644
index 00000000000..431fe0a209d
--- /dev/null
+++ b/lib/libm/noieee_src/n_catan.c
@@ -0,0 +1,126 @@
+/*	$OpenBSD: n_catan.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							catan()
+ *
+ *	Complex circular arc tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex catan();
+ * double complex z, w;
+ *
+ * w = catan (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * If
+ *     z = x + iy,
+ *
+ * then
+ *          1       (    2x     )
+ * Re w  =  - arctan(-----------)  +  k PI
+ *          2       (     2    2)
+ *                  (1 - x  - y )
+ *
+ *               ( 2         2)
+ *          1    (x  +  (y+1) )
+ * Im w  =  - log(------------)
+ *          4    ( 2         2)
+ *               (x  +  (y-1) )
+ *
+ * Where k is an arbitrary integer.
+ *
+ * catan(z) = -i catanh(iz).
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10      5900       1.3e-16     7.8e-18
+ *    IEEE      -10,+10     30000       2.3e-15     8.5e-17
+ * The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,
+ * had peak relative error 1.5e-16, rms relative error
+ * 2.9e-17.  See also clog().
+ */
+
+#include <complex.h>
+#include <math.h>
+
+#define MAXNUM 1.0e308
+
+static const double DP1 = 3.14159265160560607910E0;
+static const double DP2 = 1.98418714791870343106E-9;
+static const double DP3 = 1.14423774522196636802E-17;
+
+static double
+_redupi(double x)
+{
+	double t;
+	long i;
+
+	t = x/M_PI;
+	if(t >= 0.0)
+		t += 0.5;
+	else
+		t -= 0.5;
+
+	i = t;	/* the multiple */
+	t = i;
+	t = ((x - t * DP1) - t * DP2) - t * DP3;
+	return (t);
+}
+
+double complex
+catan(double complex z)
+{
+	double complex w;
+	double a, t, x, x2, y;
+
+	x = creal (z);
+	y = cimag (z);
+
+	if ((x == 0.0) && (y > 1.0))
+		goto ovrf;
+
+	x2 = x * x;
+	a = 1.0 - x2 - (y * y);
+	if (a == 0.0)
+		goto ovrf;
+
+	t = 0.5 * atan2 (2.0 * x, a);
+	w = _redupi (t);
+
+	t = y - 1.0;
+	a = x2 + (t * t);
+	if (a == 0.0)
+	goto ovrf;
+
+	t = y + 1.0;
+	a = (x2 + (t * t))/a;
+	w = w + (0.25 * log (a)) * I;
+	return (w);
+
+ovrf:
+	/*mtherr ("catan", OVERFLOW);*/
+	w = MAXNUM + MAXNUM * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_catanh.c b/lib/libm/noieee_src/n_catanh.c
new file mode 100644
index 00000000000..c1c1f8cff03
--- /dev/null
+++ b/lib/libm/noieee_src/n_catanh.c
@@ -0,0 +1,55 @@
+/*	$OpenBSD: n_catanh.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							catanh
+ *
+ *	Complex inverse hyperbolic tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex catanh();
+ * double complex z, w;
+ *
+ * w = catanh (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Inverse tanh, equal to  -i catan (iz);
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       2.3e-16     6.2e-17
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+catanh(double complex z)
+{
+	double complex w;
+
+	w = -1.0 * I * catan (z * I);
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_ccos.c b/lib/libm/noieee_src/n_ccos.c
new file mode 100644
index 00000000000..21ff994418b
--- /dev/null
+++ b/lib/libm/noieee_src/n_ccos.c
@@ -0,0 +1,84 @@
+/*	$OpenBSD: n_ccos.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							ccos()
+ *
+ *	Complex circular cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex ccos();
+ * double complex z, w;
+ *
+ * w = ccos (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * If
+ *     z = x + iy,
+ *
+ * then
+ *
+ *     w = cos x  cosh y  -  i sin x sinh y.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10      8400       4.5e-17     1.3e-17
+ *    IEEE      -10,+10     30000       3.8e-16     1.0e-16
+ */
+
+#include <complex.h>
+#include <math.h>
+
+/* calculate cosh and sinh */
+
+static void
+_cchsh(double x, double *c, double *s)
+{
+	double e, ei;
+
+	if (fabs(x) <= 0.5) {
+		*c = cosh(x);
+		*s = sinh(x);
+	}
+	else {
+		e = exp(x);
+		ei = 0.5/e;
+		e = 0.5 * e;
+		*s = e - ei;
+		*c = e + ei;
+	}
+}
+
+double complex
+ccos(double complex z)
+{
+	double complex w;
+	double ch, sh;
+
+	_cchsh( cimag(z), &ch, &sh );
+	w = cos(creal (z)) * ch - (sin (creal (z)) * sh) * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_ccosh.c b/lib/libm/noieee_src/n_ccosh.c
new file mode 100644
index 00000000000..7151977977e
--- /dev/null
+++ b/lib/libm/noieee_src/n_ccosh.c
@@ -0,0 +1,58 @@
+/*	$OpenBSD: n_ccosh.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							ccosh
+ *
+ *	Complex hyperbolic cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex ccosh();
+ * double complex z, w;
+ *
+ * w = ccosh (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * ccosh(z) = cosh x  cos y + i sinh x sin y .
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       2.9e-16     8.1e-17
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+ccosh(double complex z)
+{
+	double complex w;
+	double x, y;
+
+	x = creal(z);
+	y = cimag(z);
+	w = cosh (x) * cos (y)  +  (sinh (x) * sin (y)) * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_cexp.c b/lib/libm/noieee_src/n_cexp.c
new file mode 100644
index 00000000000..8dbaf613e5e
--- /dev/null
+++ b/lib/libm/noieee_src/n_cexp.c
@@ -0,0 +1,70 @@
+/*	$OpenBSD: n_cexp.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							cexp()
+ *
+ *	Complex exponential function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex cexp ();
+ * double complex z, w;
+ *
+ * w = cexp (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the exponential of the complex argument z
+ * into the complex result w.
+ *
+ * If
+ *     z = x + iy,
+ *     r = exp(x),
+ *
+ * then
+ *
+ *     w = r cos y + i r sin y.
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10      8700       3.7e-17     1.1e-17
+ *    IEEE      -10,+10     30000       3.0e-16     8.7e-17
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+cexp(double complex z)
+{
+	double complex w;
+	double r, x, y;
+
+	x = creal (z);
+	y = cimag (z);
+	r = exp (x);
+	w = r * cos (y) + r * sin (y) * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_cimag.c b/lib/libm/noieee_src/n_cimag.c
new file mode 100644
index 00000000000..5bd9d882bbf
--- /dev/null
+++ b/lib/libm/noieee_src/n_cimag.c
@@ -0,0 +1,25 @@
+/*	$OpenBSD: n_cimag.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double
+cimag(double complex z)
+{
+	return __imag__ z;
+}
diff --git a/lib/libm/noieee_src/n_cimagf.c b/lib/libm/noieee_src/n_cimagf.c
new file mode 100644
index 00000000000..245de7e827f
--- /dev/null
+++ b/lib/libm/noieee_src/n_cimagf.c
@@ -0,0 +1,25 @@
+/*	$OpenBSD: n_cimagf.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+float
+cimagf(float complex z)
+{
+	return __imag__ z;
+}
diff --git a/lib/libm/noieee_src/n_clog.c b/lib/libm/noieee_src/n_clog.c
new file mode 100644
index 00000000000..a3cd4b59ead
--- /dev/null
+++ b/lib/libm/noieee_src/n_clog.c
@@ -0,0 +1,72 @@
+/*	$OpenBSD: n_clog.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							clog.c
+ *
+ *	Complex natural logarithm
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex clog();
+ * double complex z, w;
+ *
+ * w = clog (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns complex logarithm to the base e (2.718...) of
+ * the complex argument x.
+ *
+ * If z = x + iy, r = sqrt( x**2 + y**2 ),
+ * then
+ *       w = log(r) + i arctan(y/x).
+ * 
+ * The arctangent ranges from -PI to +PI.
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10      7000       8.5e-17     1.9e-17
+ *    IEEE      -10,+10     30000       5.0e-15     1.1e-16
+ *
+ * Larger relative error can be observed for z near 1 +i0.
+ * In IEEE arithmetic the peak absolute error is 5.2e-16, rms
+ * absolute error 1.0e-16.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+clog(double complex z)
+{
+	double complex w;
+	double p, rr;
+
+	/*rr = sqrt( z->r * z->r  +  z->i * z->i );*/
+	rr = cabs(z);
+	p = log(rr);
+	rr = atan2 (cimag (z), creal (z));
+	w = p + rr * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_conj.c b/lib/libm/noieee_src/n_conj.c
new file mode 100644
index 00000000000..b776d2de957
--- /dev/null
+++ b/lib/libm/noieee_src/n_conj.c
@@ -0,0 +1,25 @@
+/*	$OpenBSD: n_conj.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+conj(double complex z)
+{
+	return ~z;
+}
diff --git a/lib/libm/noieee_src/n_conjf.c b/lib/libm/noieee_src/n_conjf.c
new file mode 100644
index 00000000000..f0589a986c1
--- /dev/null
+++ b/lib/libm/noieee_src/n_conjf.c
@@ -0,0 +1,25 @@
+/*	$OpenBSD: n_conjf.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+float complex
+conjf(float complex z)
+{
+	return ~z;
+}
diff --git a/lib/libm/noieee_src/n_cpow.c b/lib/libm/noieee_src/n_cpow.c
new file mode 100644
index 00000000000..85cb93c228a
--- /dev/null
+++ b/lib/libm/noieee_src/n_cpow.c
@@ -0,0 +1,71 @@
+/*	$OpenBSD: n_cpow.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							cpow
+ *
+ *	Complex power function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex cpow();
+ * double complex a, z, w;
+ *
+ * w = cpow (a, z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Raises complex A to the complex Zth power.
+ * Definition is per AMS55 # 4.2.8,
+ * analytically equivalent to cpow(a,z) = cexp(z clog(a)).
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       9.4e-15     1.5e-15
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+cpow(double complex a, double complex z)
+{
+	double complex w;
+	double x, y, r, theta, absa, arga;
+
+	x = creal (z);
+	y = cimag (z);
+	absa = cabs (a);
+	if (absa == 0.0) {
+	return (0.0 + 0.0 * I);
+	}
+	arga = carg (a);
+	r = pow (absa, x);
+	theta = x * arga;
+	if (y != 0.0) {
+		r = r * exp (-y * arga);
+		theta = theta + y * log (absa);
+	}
+	w = r * cos (theta) + (r * sin (theta)) * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_cproj.c b/lib/libm/noieee_src/n_cproj.c
new file mode 100644
index 00000000000..373cef865fb
--- /dev/null
+++ b/lib/libm/noieee_src/n_cproj.c
@@ -0,0 +1,32 @@
+/*	$OpenBSD: n_cproj.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+cproj(double complex z)
+{
+	double complex res;
+
+	if (isinf(__real__ z) || isinf(__imag__ z)) {
+		__real__ res = INFINITY;
+		__imag__ res = copysign(0.0, __imag__ z);
+	}
+
+	return res;
+}
diff --git a/lib/libm/noieee_src/n_cprojf.c b/lib/libm/noieee_src/n_cprojf.c
new file mode 100644
index 00000000000..978d36f9395
--- /dev/null
+++ b/lib/libm/noieee_src/n_cprojf.c
@@ -0,0 +1,32 @@
+/*	$OpenBSD: n_cprojf.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+float complex
+cprojf(float complex z)
+{
+	float complex res;
+
+	if (isinf(__real__ z) || isinf(__imag__ z)) {
+		__real__ res = INFINITY;
+		__imag__ res = copysign(0.0, __imag__ z);
+	}
+
+	return res;
+}
diff --git a/lib/libm/noieee_src/n_creal.c b/lib/libm/noieee_src/n_creal.c
new file mode 100644
index 00000000000..3bceac5719e
--- /dev/null
+++ b/lib/libm/noieee_src/n_creal.c
@@ -0,0 +1,25 @@
+/*	$OpenBSD: n_creal.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double
+creal(double complex z)
+{
+	return __real__ z;
+}
diff --git a/lib/libm/noieee_src/n_crealf.c b/lib/libm/noieee_src/n_crealf.c
new file mode 100644
index 00000000000..a31c5393696
--- /dev/null
+++ b/lib/libm/noieee_src/n_crealf.c
@@ -0,0 +1,25 @@
+/*	$OpenBSD: n_crealf.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+float
+crealf(float complex z)
+{
+	return __real__ z;
+}
diff --git a/lib/libm/noieee_src/n_csin.c b/lib/libm/noieee_src/n_csin.c
new file mode 100644
index 00000000000..017f6d6068b
--- /dev/null
+++ b/lib/libm/noieee_src/n_csin.c
@@ -0,0 +1,86 @@
+/*	$OpenBSD: n_csin.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							csin()
+ *
+ *	Complex circular sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex csin();
+ * double complex z, w;
+ *
+ * w = csin (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * If
+ *     z = x + iy,
+ *
+ * then
+ *
+ *     w = sin x  cosh y  +  i cos x sinh y.
+ *
+ * csin(z) = -i csinh(iz).
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10      8400       5.3e-17     1.3e-17
+ *    IEEE      -10,+10     30000       3.8e-16     1.0e-16
+ * Also tested by csin(casin(z)) = z.
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+/* calculate cosh and sinh */
+
+static void
+cchsh(double x, double *c, double *s)
+{
+	double e, ei;
+
+	if (fabs(x) <= 0.5) {
+		*c = cosh(x);
+		*s = sinh(x);
+	}
+	else {
+		e = exp(x);
+		ei = 0.5/e;
+		e = 0.5 * e;
+		*s = e - ei;
+		*c = e + ei;
+	}
+}
+
+double complex
+csin(double complex z)
+{
+	double complex w;
+	double ch, sh;
+
+	cchsh( cimag (z), &ch, &sh );
+	w = sin (creal(z)) * ch + (cos (creal(z)) * sh) * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_csinh.c b/lib/libm/noieee_src/n_csinh.c
new file mode 100644
index 00000000000..122985612bb
--- /dev/null
+++ b/lib/libm/noieee_src/n_csinh.c
@@ -0,0 +1,57 @@
+/*	$OpenBSD: n_csinh.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							csinh
+ *
+ *	Complex hyperbolic sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex csinh();
+ * double complex z, w;
+ *
+ * w = csinh (z);
+ *
+ * DESCRIPTION:
+ *
+ * csinh z = (cexp(z) - cexp(-z))/2
+ *         = sinh x * cos y  +  i cosh x * sin y .
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       3.1e-16     8.2e-17
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+csinh(double complex z)
+{
+	double complex w;
+	double x, y;
+
+	x = creal(z);
+	y = cimag(z);
+	w = sinh (x) * cos (y)  +  (cosh (x) * sin (y)) * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_csqrt.c b/lib/libm/noieee_src/n_csqrt.c
new file mode 100644
index 00000000000..13f0903fcb2
--- /dev/null
+++ b/lib/libm/noieee_src/n_csqrt.c
@@ -0,0 +1,131 @@
+/*	$OpenBSD: n_csqrt.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							csqrt()
+ *
+ *	Complex square root
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex csqrt();
+ * double complex z, w;
+ *
+ * w = csqrt (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ *
+ * If z = x + iy,  r = |z|, then
+ *
+ *                       1/2
+ * Re w  =  [ (r + x)/2 ]   ,
+ *
+ *                       1/2
+ * Im w  =  [ (r - x)/2 ]   .
+ *
+ * Cancellation error in r-x or r+x is avoided by using the
+ * identity  2 Re w Im w  =  y.
+ *
+ * Note that -w is also a square root of z.  The root chosen
+ * is always in the right half plane and Im w has the same sign as y.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10     25000       3.2e-17     9.6e-18
+ *    IEEE      -10,+10   1,000,000     2.9e-16     6.1e-17
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+csqrt(double complex z)
+{
+	double complex w;
+	double x, y, r, t, scale;
+
+	x = creal (z);
+	y = cimag (z);
+
+	if (y == 0.0) {
+		if (x == 0.0) {
+			w = 0.0 + y * I;
+		}
+		else {
+			r = fabs (x);
+			r = sqrt (r);
+			if (x < 0.0) {
+				w = 0.0 + r * I;
+			}
+			else {
+				w = r + y * I;
+			}
+		}
+		return (w);
+	}
+	if (x == 0.0) {
+		r = fabs (y);
+		r = sqrt (0.5*r);
+		if (y > 0)
+			w = r + r * I;
+		else
+			w = r - r * I;
+		return (w);
+	}
+	/* Rescale to avoid internal overflow or underflow.  */
+	if ((fabs(x) > 4.0) || (fabs(y) > 4.0)) {
+		x *= 0.25;
+		y *= 0.25;
+		scale = 2.0;
+	}
+	else {
+		x *= 1.8014398509481984e16;  /* 2^54 */
+		y *= 1.8014398509481984e16;
+		scale = 7.450580596923828125e-9; /* 2^-27 */
+#if 0
+		x *= 4.0;
+		y *= 4.0;
+		scale = 0.5;
+#endif
+	}
+	w = x + y * I;
+	r = cabs(w);
+	if (x > 0) {
+		t = sqrt(0.5 * r + 0.5 * x);
+		r = scale * fabs((0.5 * y) / t);
+		t *= scale;
+	}
+	else {
+		r = sqrt( 0.5 * r - 0.5 * x );
+		t = scale * fabs( (0.5 * y) / r );
+		r *= scale;
+	}
+	if (y < 0)
+		w = t - r * I;
+	else
+		w = t + r * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_ctan.c b/lib/libm/noieee_src/n_ctan.c
new file mode 100644
index 00000000000..6ae9cca8d30
--- /dev/null
+++ b/lib/libm/noieee_src/n_ctan.c
@@ -0,0 +1,152 @@
+/*	$OpenBSD: n_ctan.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							ctan()
+ *
+ *	Complex circular tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex ctan();
+ * double complex z, w;
+ *
+ * w = ctan (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * If
+ *     z = x + iy,
+ *
+ * then
+ *
+ *           sin 2x  +  i sinh 2y
+ *     w  =  --------------------.
+ *            cos 2x  +  cosh 2y
+ *
+ * On the real axis the denominator is zero at odd multiples
+ * of PI/2.  The denominator is evaluated by its Taylor
+ * series near these points.
+ *
+ * ctan(z) = -i ctanh(iz).
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10      5200       7.1e-17     1.6e-17
+ *    IEEE      -10,+10     30000       7.2e-16     1.2e-16
+ * Also tested by ctan * ccot = 1 and catan(ctan(z))  =  z.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+#define MACHEP 1.1e-16
+#define MAXNUM 1.0e308
+
+static const double DP1 = 3.14159265160560607910E0;
+static const double DP2 = 1.98418714791870343106E-9;
+static const double DP3 = 1.14423774522196636802E-17;
+
+static double
+_redupi(double x)
+{
+	double t;
+	long i;
+
+	t = x/M_PI;
+	if (t >= 0.0)
+		t += 0.5;
+	else
+		t -= 0.5;
+
+	i = t;	/* the multiple */
+	t = i;
+	t = ((x - t * DP1) - t * DP2) - t * DP3;
+	return (t);
+}
+
+/*  Taylor series expansion for cosh(2y) - cos(2x)	*/
+
+static double
+_ctans(double complex z)
+{
+	double f, x, x2, y, y2, rn, t;
+	double d;
+
+	x = fabs (2.0 * creal (z));
+	y = fabs (2.0 * cimag(z));
+
+	x = _redupi(x);
+
+	x = x * x;
+	y = y * y;
+	x2 = 1.0;
+	y2 = 1.0;
+	f = 1.0;
+	rn = 0.0;
+	d = 0.0;
+	do {
+		rn += 1.0;
+		f *= rn;
+		rn += 1.0;
+		f *= rn;
+		x2 *= x;
+		y2 *= y;
+		t = y2 + x2;
+		t /= f;
+		d += t;
+
+		rn += 1.0;
+		f *= rn;
+		rn += 1.0;
+		f *= rn;
+		x2 *= x;
+		y2 *= y;
+		t = y2 - x2;
+		t /= f;
+		d += t;
+	}
+	while (fabs(t/d) > MACHEP)
+		;
+	return (d);
+}
+
+double complex
+ctan(double complex z)
+{
+	double complex w;
+	double d;
+
+	d = cos (2.0 * creal (z)) + cosh (2.0 * cimag (z));
+
+	if (fabs(d) < 0.25)
+		d = _ctans (z);
+
+	if (d == 0.0) {
+		/*mtherr ("ctan", OVERFLOW);*/
+		w = MAXNUM + MAXNUM * I;
+		return (w);
+	}
+
+	w = sin (2.0 * creal(z)) / d + (sinh (2.0 * cimag(z)) / d) * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_ctanh.c b/lib/libm/noieee_src/n_ctanh.c
new file mode 100644
index 00000000000..4d3536f9986
--- /dev/null
+++ b/lib/libm/noieee_src/n_ctanh.c
@@ -0,0 +1,59 @@
+/*	$OpenBSD: n_ctanh.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							ctanh
+ *
+ *	Complex hyperbolic tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex ctanh();
+ * double complex z, w;
+ *
+ * w = ctanh (z);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * tanh z = (sinh 2x  +  i sin 2y) / (cosh 2x + cos 2y) .
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       1.7e-14     2.4e-16
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+double complex
+ctanh(double complex z)
+{
+	double complex w;
+	double x, y, d;
+
+	x = creal(z);
+	y = cimag(z);
+	d = cosh (2.0 * x) + cos (2.0 * y);
+	w = sinh (2.0 * x) / d  +  (sin (2.0 * y) / d) * I;
+	return (w);
+}
diff --git a/lib/libm/noieee_src/n_hypot.c b/lib/libm/noieee_src/n_hypot.c
new file mode 100644
index 00000000000..4f613607400
--- /dev/null
+++ b/lib/libm/noieee_src/n_hypot.c
@@ -0,0 +1,201 @@
+/*	$OpenBSD: n_hypot.c,v 1.1 2008/10/07 22:25:53 martynas Exp $	*/
+/*	$NetBSD: n_cabs.c,v 1.1 1995/10/10 23:36:39 ragge Exp $	*/
+/*
+ * Copyright (c) 1985, 1993
+ *	The Regents of the University of California.  All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 3. Neither the name of the University nor the names of its contributors
+ *    may be used to endorse or promote products derived from this software
+ *    without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#ifndef lint
+static char sccsid[] = "@(#)cabs.c	8.1 (Berkeley) 6/4/93";
+#endif /* not lint */
+
+/* HYPOT(X,Y)
+ * RETURN THE SQUARE ROOT OF X^2 + Y^2  WHERE Z=X+iY
+ * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
+ * CODED IN C BY K.C. NG, 11/28/84;
+ * REVISED BY K.C. NG, 7/12/85.
+ *
+ * Required system supported functions :
+ *	copysign(x,y)
+ *	finite(x)
+ *	scalbn(x,N)
+ *	sqrt(x)
+ *
+ * Method :
+ *	1. replace x by |x| and y by |y|, and swap x and
+ *	   y if y > x (hence x is never smaller than y).
+ *	2. Hypot(x,y) is computed by:
+ *	   Case I, x/y > 2
+ *
+ *				       y
+ *		hypot = x + -----------------------------
+ *			 		    2
+ *			    sqrt ( 1 + [x/y]  )  +  x/y
+ *
+ *	   Case II, x/y <= 2
+ *				                   y
+ *		hypot = x + --------------------------------------------------
+ *				          		     2
+ *				     			[x/y]   -  2
+ *			   (sqrt(2)+1) + (x-y)/y + -----------------------------
+ *			 		    			  2
+ *			    			  sqrt ( 1 + [x/y]  )  + sqrt(2)
+ *
+ *
+ *
+ * Special cases:
+ *	hypot(x,y) is INF if x or y is +INF or -INF; else
+ *	hypot(x,y) is NAN if x or y is NAN.
+ *
+ * Accuracy:
+ * 	hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units
+ *	in the last place). See Kahan's "Interval Arithmetic Options in the
+ *	Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
+ *      1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
+ *	code follows in	comments.) In a test run with 500,000 random arguments
+ *	on a VAX, the maximum observed error was .959 ulps.
+ *
+ * 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 "mathimpl.h"
+
+vc(r2p1hi, 2.4142135623730950345E0   ,8279,411a,ef32,99fc,   2, .9A827999FCEF32)
+vc(r2p1lo, 1.4349369327986523769E-17 ,597d,2484,754b,89b3, -55, .84597D89B3754B)
+vc(sqrt2,  1.4142135623730950622E0   ,04f3,40b5,de65,33f9,   1, .B504F333F9DE65)
+
+ic(r2p1hi, 2.4142135623730949234E0   ,   1, 1.3504F333F9DE6)
+ic(r2p1lo, 1.2537167179050217666E-16 , -53, 1.21165F626CDD5)
+ic(sqrt2,  1.4142135623730951455E0   ,   0, 1.6A09E667F3BCD)
+
+#ifdef vccast
+#define	r2p1hi	vccast(r2p1hi)
+#define	r2p1lo	vccast(r2p1lo)
+#define	sqrt2	vccast(sqrt2)
+#endif
+
+double
+hypot(double x, double y)
+{
+	static const double zero=0, one=1,
+		      small=1.0E-18;	/* fl(1+small)==1 */
+	static const ibig=30;	/* fl(1+2**(2*ibig))==1 */
+	double t,r;
+	int exp;
+
+	if(finite(x))
+	    if(finite(y))
+	    {
+		x=copysign(x,one);
+		y=copysign(y,one);
+		if(y > x)
+		    { t=x; x=y; y=t; }
+		if(x == zero) return(zero);
+		if(y == zero) return(x);
+		exp= logb(x);
+		if (exp - (int)logb(y) > ibig) {
+			if (one + small >= 1.0)	/* raise inexact flag */
+				return(x);	/* return |x| */
+		}
+
+	    /* start computing sqrt(x^2 + y^2) */
+		r=x-y;
+		if(r>y) { 	/* x/y > 2 */
+		    r=x/y;
+		    r=r+sqrt(one+r*r); }
+		else {		/* 1 <= x/y <= 2 */
+		    r/=y; t=r*(r+2.0);
+		    r+=t/(sqrt2+sqrt(2.0+t));
+		    r+=r2p1lo; r+=r2p1hi; }
+
+		r=y/r;
+		return(x+r);
+
+	    }
+
+	    else if(isinf(y))		/* y is +-INF */
+		     return(copysign(y,one));
+	    else
+		     return(y);		/* y is NaN and x is finite */
+
+	else if(isinf(x))		/* x is +-INF */
+	         return (copysign(x,one));
+	else if(finite(y))
+	         return(x);		/* x is NaN, y is finite */
+	else if (isnan(y))
+		return (y);
+	else return(copysign(y,one));	/* y is INF */
+}
+
+/* A faster but less accurate version of cabs(x,y) */
+#if 0
+double
+hypot(double x, double y)
+{
+	static const double zero=0, one=1;
+		      small=1.0E-18;	/* fl(1+small)==1 */
+	static const ibig=30;	/* fl(1+2**(2*ibig))==1 */
+	double temp;
+	int exp;
+
+	if(finite(x))
+	    if(finite(y))
+	    {
+		x=copysign(x,one);
+		y=copysign(y,one);
+		if(y > x)
+		    { temp=x; x=y; y=temp; }
+		if(x == zero) return(zero);
+		if(y == zero) return(x);
+		exp= logb(x);
+		x=scalbn(x,-exp);
+		if (exp - (int)logb(y) > ibig) {
+			if (one + small >= 1.0)		/* raise inexact flag */
+				return(scalbn(x,exp));	/* return |x| */
+		}
+		else y=scalbn(y,-exp);
+		return(scalbn(sqrt(x*x+y*y),exp));
+	    }
+
+	    else if(isinf(y))		/* y is +-INF */
+		     return(copysign(y,one));
+	    else
+		     return(y);		/* y is NaN and x is finite */
+
+	else if(isinf(x))		/* x is +-INF */
+	         return (copysign(x,one));
+	else if(finite(y))
+	         return(x);		/* x is NaN, y is finite */
+	else if(isnan(y)) return(y);	/* x and y is NaN */
+	else return(copysign(y,one));	/* y is INF */
+}
+#endif