Add various complex math routines from OpenBSD.

src/s_cabs.c

@@ -0,0 +1,30 @@
+/*	$OpenBSD: s_cabs.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.h>
+#include <math.h>
+
+double
+cabs(double complex z)
+{
+	return hypot(__real__ z, __imag__ z);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(cabsl, cabs);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_cabsf.c

@@ -0,0 +1,25 @@
+/*	$OpenBSD: s_cabsf.c,v 1.1 2008/09/07 20:36:09 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
+cabsf(float complex z)
+{
+	return hypotf(__real__ z, __imag__ z);
+}

src/s_cabsl.c

@@ -0,0 +1,26 @@
+/*	$OpenBSD: s_cabsl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2011 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>
+
+long double
+cabsl(long double complex z)
+{
+	return hypotl(__real__ z, __imag__ z);
+}

src/s_cacos.c

@@ -0,0 +1,65 @@
+/*	$OpenBSD: s_cacos.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.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);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(cacosl, cacos);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_cacosf.c

@@ -0,0 +1,60 @@
+/*	$OpenBSD: s_cacosf.c,v 1.2 2011/07/20 19:28:33 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.
+ */
+
+/*							cacosf()
+ *
+ *	Complex circular arc cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * void cacosf();
+ * cmplxf z, w;
+ *
+ * cacosf( &z, &w );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ *
+ * w = arccos z  =  PI/2 - arcsin z.
+ *
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       9.2e-6       1.2e-6
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+float complex
+cacosf(float complex z)
+{
+	float complex w;
+
+	w = casinf( z );
+	w = ((float)M_PI_2 - crealf (w)) - cimagf (w) * I;
+	return (w);
+}

src/s_cacosh.c

@@ -0,0 +1,60 @@
+/*	$OpenBSD: s_cacosh.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.h>
+#include <math.h>
+
+double complex
+cacosh(double complex z)
+{
+	double complex w;
+
+	w = I * cacos (z);
+	return (w);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(cacoshl, cacosh);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_cacoshf.c

@@ -0,0 +1,55 @@
+/*	$OpenBSD: s_cacoshf.c,v 1.1 2008/09/07 20:36:09 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.
+ */
+
+/*							cacoshf
+ *
+ *	Complex inverse hyperbolic cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float complex cacoshf();
+ * float complex z, w;
+ *
+ * w = cacoshf (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>
+
+float complex
+cacoshf(float complex z)
+{
+	float complex w;
+
+	w = I * cacosf (z);
+	return (w);
+}

src/s_cacoshl.c

@@ -0,0 +1,56 @@
+/*	$OpenBSD: s_cacoshl.c,v 1.1 2011/07/08 19:25:31 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.
+ */
+
+/*							cacoshl
+ *
+ *	Complex inverse hyperbolic cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex cacoshl();
+ * long double complex z, w;
+ *
+ * w = cacoshl (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>
+
+long double complex
+cacoshl(long double complex z)
+{
+	long double complex w;
+
+	w = I * cacosl(z);
+	return (w);
+}

src/s_cacosl.c

@@ -0,0 +1,63 @@
+/*	$OpenBSD: s_cacosl.c,v 1.3 2011/07/20 21:02:51 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.
+ */
+
+/*							cacosl()
+ *
+ *	Complex circular arc cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex cacosl();
+ * long double complex z, w;
+ *
+ * w = cacosl( 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>
+
+static const long double PIO2L = 1.570796326794896619231321691639751442098585L;
+
+long double complex
+cacosl(long double complex z)
+{
+	long double complex w;
+
+	w = casinl(z);
+	w = (PIO2L - creall(w)) - cimagl(w) * I;
+	return (w);
+}

src/s_casin.c

@@ -0,0 +1,134 @@
+/*	$OpenBSD: s_casin.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.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);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(casinl, casin);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_casinf.c

@@ -0,0 +1,132 @@
+/*	$OpenBSD: s_casinf.c,v 1.3 2011/07/20 19:28:33 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.
+ */
+
+/*							casinf()
+ *
+ *	Complex circular arc sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * void casinf();
+ * cmplxf z, w;
+ *
+ * casinf( &z, &w );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Inverse complex sine:
+ *
+ *                               2
+ * w = -i clog( iz + csqrt( 1 - z ) ).
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       1.1e-5      1.5e-6
+ * Larger relative error can be observed for z near zero.
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+float complex
+casinf(float complex z)
+{
+	float complex w;
+	float x, y;
+	static float complex ca, ct, zz, z2;
+	/*
+	float cn, n;
+	static float a, b, s, t, u, v, y2;
+	static cmplxf sum;
+	*/
+
+	x = crealf(z);
+	y = cimagf(z);
+
+	if(y == 0.0f) {
+		if(fabsf(x) > 1.0f) {
+			w = (float)M_PI_2 + 0.0f * I;
+			/*mtherr( "casinf", DOMAIN );*/
+		}
+		else {
+			w = asinf (x) + 0.0f * I;
+		}
+		return (w);
+	}
+
+	/* Power series expansion */
+	/*
+	b = cabsf(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 = fabsf(ct.r) + fabsf(ct.i);
+		}
+		while(b > MACHEPF);
+		w->r = sum.r;
+		w->i = sum.i;
+		return;
+	}
+	*/
+
+
+	ca = x + y * I;
+	ct = ca * I;	/* iz */
+	/* sqrt( 1 - z*z) */
+	/* cmul( &ca, &ca, &zz ) */
+	/*x * x  -  y * y */
+	zz = (x - y) * (x + y) + (2.0f * x * y) * I;
+	zz = 1.0f - crealf(zz) - cimagf(zz) * I;
+	z2 = csqrtf (zz);
+
+	zz = ct + z2;
+	zz = clogf (zz);
+	/* multiply by 1/i = -i */
+	w = zz * (-1.0f * I);
+	return (w);
+}

src/s_casinh.c

@@ -0,0 +1,60 @@
+/*	$OpenBSD: s_casinh.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.h>
+#include <math.h>
+
+double complex
+casinh(double complex z)
+{
+	double complex w;
+
+	w = -1.0 * I * casin (z * I);
+	return (w);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(casinhl, casinh);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_casinhf.c

@@ -0,0 +1,55 @@
+/*	$OpenBSD: s_casinhf.c,v 1.1 2008/09/07 20:36:09 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.
+ */
+
+/*							casinhf
+ *
+ *	Complex inverse hyperbolic sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float complex casinhf();
+ * float complex z, w;
+ *
+ * w = casinhf (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>
+
+float complex
+casinhf(float complex z)
+{
+	float complex w;
+
+	w = -1.0f * I * casinf (z * I);
+	return (w);
+}

src/s_casinhl.c

@@ -0,0 +1,56 @@
+/*	$OpenBSD: s_casinhl.c,v 1.1 2011/07/08 19:25:31 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.
+ */
+
+/*							casinhl
+ *
+ *	Complex inverse hyperbolic sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex casinhf();
+ * long double complex z, w;
+ *
+ * w = casinhl (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>
+
+long double complex
+casinhl(long double complex z)
+{
+	long double complex w;
+
+	w = -1.0L * I * casinl(z * I);
+	return (w);
+}

src/s_casinl.c

@@ -0,0 +1,130 @@
+/*	$OpenBSD: s_casinl.c,v 1.3 2011/07/20 21:02:51 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.
+ */
+
+/*							casinl()
+ *
+ *	Complex circular arc sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex casinl();
+ * long double complex z, w;
+ *
+ * w = casinl( z );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Inverse complex sine:
+ *
+ *                               2
+ * w = -i clog( iz + csqrt( 1 - z ) ).
+ *
+ *
+ * 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 <float.h>
+#include <math.h>
+
+#if	LDBL_MANT_DIG == 64
+static const long double MACHEPL= 5.42101086242752217003726400434970855712890625E-20L;
+#elif	LDBL_MANT_DIG == 113
+static const long double MACHEPL = 9.629649721936179265279889712924636592690508e-35L;
+#endif
+
+static const long double PIO2L = 1.570796326794896619231321691639751442098585L;
+
+long double complex
+casinl(long double complex z)
+{
+	long double complex w;
+	long double x, y, b;
+	static long double complex ca, ct, zz, z2;
+
+	x = creall(z);
+	y = cimagl(z);
+
+	if (y == 0.0L) {
+		if (fabsl(x) > 1.0L) {
+			w = PIO2L + 0.0L * I;
+			/*mtherr( "casinl", DOMAIN );*/
+		}
+		else {
+			w = asinl(x) + 0.0L * I;
+		}
+		return (w);
+	}
+
+	/* Power series expansion */
+	b = cabsl(z);
+	if (b < 0.125L) {
+		long double complex sum;
+		long double n, cn;
+
+		z2 = (x - y) * (x + y) + (2.0L * x * y) * I;
+		cn = 1.0L;
+		n = 1.0L;
+		ca = x + y * I;
+		sum = x + y * I;
+		do {
+			ct = z2 * ca;
+			ca = ct;
+
+			cn *= n;
+			n += 1.0L;
+			cn /= n;
+			n += 1.0L;
+			b = cn/n;
+
+			ct *= b;
+			sum += ct;
+			b = cabsl(ct);
+		}
+
+		while (b > MACHEPL);
+		w = sum;
+		return w;
+	}
+
+	ca = x + y * I;
+	ct = ca * I;	/* iz */
+	/* sqrt(1 - z*z) */
+	/* cmul(&ca, &ca, &zz) */
+	/* x * x  -  y * y */
+	zz = (x - y) * (x + y) + (2.0L * x * y) * I;
+	zz = 1.0L - creall(zz) - cimagl(zz) * I;
+	z2 = csqrtl(zz);
+
+	zz = ct + z2;
+	zz = clogl(zz);
+	/* multiply by 1/i = -i */
+	w = zz * (-1.0L * I);
+	return (w);
+}

src/s_catan.c

@@ -0,0 +1,131 @@
+/*	$OpenBSD: s_catan.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.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);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(catanl, catan);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_catanf.c

@@ -0,0 +1,124 @@
+/*	$OpenBSD: s_catanf.c,v 1.2 2010/07/18 18:42:26 guenther 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.
+ */
+
+/*							catanf()
+ *
+ *	Complex circular arc tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float complex catanf();
+ * float complex z, w;
+ *
+ * w = catanf( 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.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000        2.3e-6      5.2e-8
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+#define MAXNUMF 1.0e38F
+
+static const double DP1 = 3.140625;
+static const double DP2 = 9.67502593994140625E-4;
+static const double DP3 = 1.509957990978376432E-7;
+
+static float
+_redupif(float xx)
+{
+	float x, t;
+	long i;
+
+	x = xx;
+	t = x/(float)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);
+}
+
+float complex
+catanf(float complex z)
+{
+	float complex w;
+	float a, t, x, x2, y;
+
+	x = crealf(z);
+	y = cimagf(z);
+
+	if((x == 0.0f) && (y > 1.0f))
+		goto ovrf;
+
+	x2 = x * x;
+	a = 1.0f - x2 - (y * y);
+	if (a == 0.0f)
+		goto ovrf;
+
+	t = 0.5f * atan2f(2.0f * x, a);
+	w = _redupif(t);
+
+	t = y - 1.0f;
+	a = x2 + (t * t);
+	if(a == 0.0f)
+		goto ovrf;
+
+	t = y + 1.0f;
+	a = (x2 + (t * t))/a;
+	w = w + (0.25f * logf (a)) * I;
+	return (w);
+
+ovrf:
+	/*mtherr( "catanf", OVERFLOW );*/
+	w = MAXNUMF + MAXNUMF * I;
+	return (w);
+}

src/s_catanh.c

@@ -0,0 +1,60 @@
+/*	$OpenBSD: s_catanh.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.h>
+#include <math.h>
+
+double complex
+catanh(double complex z)
+{
+	double complex w;
+
+	w = -1.0 * I * catan (z * I);
+	return (w);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(catanhl, catanh);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_catanhf.c

@@ -0,0 +1,55 @@
+/*	$OpenBSD: s_catanhf.c,v 1.1 2008/09/07 20:36:09 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.
+ */
+
+/*							catanhf
+ *
+ *	Complex inverse hyperbolic tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float complex catanhf();
+ * float complex z, w;
+ *
+ * w = catanhf (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>
+
+float complex
+catanhf(float complex z)
+{
+	float complex w;
+
+	w = -1.0f * I * catanf (z * I);
+	return (w);
+}

src/s_catanhl.c

@@ -0,0 +1,56 @@
+/*	$OpenBSD: s_catanhl.c,v 1.1 2011/07/08 19:25:31 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.
+ */
+
+/*							catanhl
+ *
+ *	Complex inverse hyperbolic tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex catanhl();
+ * long double complex z, w;
+ *
+ * w = catanhl (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>
+
+long double complex
+catanhl(long double complex z)
+{
+	long double complex w;
+
+	w = -1.0L * I * catanl(z * I);
+	return (w);
+}

src/s_catanl.c

@@ -0,0 +1,127 @@
+/*	$OpenBSD: s_catanl.c,v 1.3 2011/07/20 21:02:51 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.
+ */
+
+/*							catanl()
+ *
+ *	Complex circular arc tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex catanl();
+ * long double complex z, w;
+ *
+ * w = catanl( 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.
+ *
+ *
+ *
+ * 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 <float.h>
+#include <math.h>
+
+static const long double PIL = 3.141592653589793238462643383279502884197169L;
+static const long double DP1 = 3.14159265358979323829596852490908531763125L;
+static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L;
+static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L;
+
+static long double
+redupil(long double x)
+{
+	long double t;
+	long i;
+
+	t = x / PIL;
+	if (t >= 0.0L)
+		t += 0.5L;
+	else
+		t -= 0.5L;
+
+	i = t;	/* the multiple */
+	t = i;
+	t = ((x - t * DP1) - t * DP2) - t * DP3;
+	return (t);
+}
+
+long double complex
+catanl(long double complex z)
+{
+	long double complex w;
+	long double a, t, x, x2, y;
+
+	x = creall(z);
+	y = cimagl(z);
+
+	if ((x == 0.0L) && (y > 1.0L))
+		goto ovrf;
+
+	x2 = x * x;
+	a = 1.0L - x2 - (y * y);
+	if (a == 0.0L)
+		goto ovrf;
+
+	t = atan2l(2.0L * x, a) * 0.5L;
+	w = redupil(t);
+
+	t = y - 1.0L;
+	a = x2 + (t * t);
+	if (a == 0.0L)
+		goto ovrf;
+
+	t = y + 1.0L;
+	a = (x2 + (t * t)) / a;
+	w = w + (0.25L * logl(a)) * I;
+	return (w);
+
+ovrf:
+	/*mtherr( "catanl", OVERFLOW );*/
+	w = LDBL_MAX + LDBL_MAX * I;
+	return (w);
+}

src/s_ccos.c

@@ -0,0 +1,89 @@
+/*	$OpenBSD: s_ccos.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.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);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(ccosl, ccos);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_ccosf.c

@@ -0,0 +1,84 @@
+/*	$OpenBSD: s_ccosf.c,v 1.2 2010/07/18 18:42:26 guenther 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.
+ */
+
+/*							ccosf()
+ *
+ *	Complex circular cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * void ccosf();
+ * cmplxf z, w;
+ *
+ * ccosf( &z, &w );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * If
+ *     z = x + iy,
+ *
+ * then
+ *
+ *     w = cos x  cosh y  -  i sin x sinh y.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       1.8e-7       5.5e-8
+ */
+
+#include <complex.h>
+#include <math.h>
+
+/* calculate cosh and sinh */
+
+static void
+_cchshf(float xx, float *c, float *s)
+{
+	float x, e, ei;
+
+	x = xx;
+	if(fabsf(x) <= 0.5f) {
+		*c = coshf(x);
+		*s = sinhf(x);
+	}
+	else {
+		e = expf(x);
+		ei = 0.5f/e;
+		e = 0.5f * e;
+		*s = e - ei;
+		*c = e + ei;
+	}
+}
+
+float complex
+ccosf(float complex z)
+{
+	float complex w;
+	float ch, sh;
+
+	_cchshf( cimagf(z), &ch, &sh );
+	w = cosf( crealf(z) ) * ch + ( -sinf( crealf(z) ) * sh) * I;
+	return (w);
+}

src/s_ccoshl.c

@@ -0,0 +1,59 @@
+/*	$OpenBSD: s_ccoshl.c,v 1.2 2011/07/20 19:28:33 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.
+ */
+
+/*							ccoshl
+ *
+ *	Complex hyperbolic cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex ccoshl();
+ * long double complex z, w;
+ *
+ * w = ccoshl (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>
+
+long double complex
+ccoshl(long double complex z)
+{
+	long double complex w;
+	long double x, y;
+
+	x = creall(z);
+	y = cimagl(z);
+	w = coshl(x) * cosl(y) + (sinhl(x) * sinl(y)) * I;
+	return (w);
+}

src/s_ccosl.c

@@ -0,0 +1,82 @@
+/*	$OpenBSD: s_ccosl.c,v 1.2 2011/07/20 19:28:33 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.
+ */
+
+/*							ccosl()
+ *
+ *	Complex circular cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex ccosl();
+ * long double complex z, w;
+ *
+ * w = ccosl( 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>
+
+static void
+cchshl(long double x, long double *c, long double *s)
+{
+	long double e, ei;
+
+	if(fabsl(x) <= 0.5L) {
+		*c = coshl(x);
+		*s = sinhl(x);
+	} else {
+		e = expl(x);
+		ei = 0.5L/e;
+		e = 0.5L * e;
+		*s = e - ei;
+		*c = e + ei;
+	}
+}
+
+long double complex
+ccosl(long double complex z)
+{
+	long double complex w;
+	long double ch, sh;
+
+	cchshl(cimagl(z), &ch, &sh);
+	w = cosl(creall(z)) * ch + (-sinl(creall(z)) * sh) * I;
+	return (w);
+}

src/s_cexpl.c

@@ -0,0 +1,69 @@
+/*	$OpenBSD: s_cexpl.c,v 1.2 2011/07/20 19:28:33 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.
+ */
+
+/*							cexpl()
+ *
+ *	Complex exponential function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex cexpl();
+ * long double complex z, w;
+ *
+ * w = cexpl( 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>
+
+long double complex
+cexpl(long double complex z)
+{
+	long double complex w;
+	long double r;
+
+	r = expl(creall(z));
+	w = r * cosl(cimagl(z)) + (r * sinl(cimagl(z))) * I;
+	return (w);
+}

src/s_clog.c

@@ -0,0 +1,77 @@
+/*	$OpenBSD: s_clog.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.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);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(clogl, clog);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_clogf.c

@@ -0,0 +1,72 @@
+/*	$OpenBSD: s_clogf.c,v 1.2 2010/07/18 18:42:26 guenther 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.
+ */
+
+/*							clogf.c
+ *
+ *	Complex natural logarithm
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * void clogf();
+ * cmplxf z, w;
+ *
+ * clogf( &z, &w );
+ *
+ *
+ *
+ * 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
+ *    IEEE      -10,+10     30000       1.9e-6       6.2e-8
+ *
+ * Larger relative error can be observed for z near 1 +i0.
+ * In IEEE arithmetic the peak absolute error is 3.1e-7.
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+float complex
+clogf(float complex z)
+{
+	float complex w;
+	float p, rr, x, y;
+
+	x = crealf(z);
+	y = cimagf(z);
+	rr = atan2f(y, x);
+	p = cabsf(z);
+	p = logf(p);
+	w = p + rr * I;
+	return (w);
+}

src/s_clogl.c

@@ -0,0 +1,73 @@
+/*	$OpenBSD: s_clogl.c,v 1.2 2011/07/20 19:28:33 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.
+ */
+
+/*							clogl.c
+ *
+ *	Complex natural logarithm
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex clogl();
+ * long double complex z, w;
+ *
+ * w = clogl( 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>
+
+long double complex
+clogl(long double complex z)
+{
+	long double complex w;
+	long double p, rr;
+
+	/*rr = sqrt(z->r * z->r  +  z->i * z->i);*/
+	p = cabsl(z);
+	p = logl(p);
+	rr = atan2l(cimagl(z), creall(z));
+	w = p + rr * I;
+	return (w);
+}

src/s_csin.c

@@ -0,0 +1,91 @@
+/*	$OpenBSD: s_csin.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.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);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(csinl, csin);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_csinf.c

@@ -0,0 +1,85 @@
+/*	$OpenBSD: s_csinf.c,v 1.2 2010/07/18 18:42:26 guenther 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.
+ */
+
+/*							csinf()
+ *
+ *	Complex circular sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * void csinf();
+ * cmplxf z, w;
+ *
+ * csinf( &z, &w );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * If
+ *     z = x + iy,
+ *
+ * then
+ *
+ *     w = sin x  cosh y  +  i cos x sinh y.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       1.9e-7      5.5e-8
+ *
+ */
+
+#include <complex.h>
+#include <math.h>
+
+/* calculate cosh and sinh */
+
+static void
+cchshf(float xx, float *c, float *s)
+{
+	float x, e, ei;
+
+	x = xx;
+	if(fabsf(x) <= 0.5f) {
+		*c = coshf(x);
+		*s = sinhf(x);
+	}
+	else {
+		e = expf(x);
+		ei = 0.5f/e;
+		e = 0.5f * e;
+		*s = e - ei;
+		*c = e + ei;
+	}
+}
+
+float complex
+csinf(float complex z)
+{
+	float complex w;
+	float ch, sh;
+
+	cchshf(cimagf(z), &ch, &sh);
+	w = sinf(crealf(z)) * ch  + (cosf(crealf(z)) * sh) * I;
+	return (w);
+}

src/s_csinhl.c

@@ -0,0 +1,58 @@
+/*	$OpenBSD: s_csinhl.c,v 1.2 2011/07/20 19:28:33 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.
+ */
+
+/*							csinhl
+ *
+ *	Complex hyperbolic sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex csinhl();
+ * long double complex z, w;
+ *
+ * w = csinhl (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>
+
+long double complex
+csinhl(long double complex z)
+{
+	long double complex w;
+	long double x, y;
+
+	x = creall(z);
+	y = cimagl(z);
+	w = sinhl(x) * cosl(y) + (coshl(x) * sinl(y)) * I;
+	return (w);
+}

src/s_csinl.c

@@ -0,0 +1,84 @@
+/*	$OpenBSD: s_csinl.c,v 1.2 2011/07/20 19:28:33 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.
+ */
+
+/*							csinl()
+ *
+ *	Complex circular sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex csinl();
+ * long double complex z, w;
+ *
+ * w = csinl( z );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * If
+ *     z = x + iy,
+ *
+ * then
+ *
+ *     w = sin x  cosh y  +  i cos x sinh y.
+ *
+ *
+ *
+ * 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>
+
+static void
+cchshl(long double x, long double *c, long double *s)
+{
+	long double e, ei;
+
+	if(fabsl(x) <= 0.5L) {
+		*c = coshl(x);
+		*s = sinhl(x);
+	} else {
+		e = expl(x);
+		ei = 0.5L/e;
+		e = 0.5L * e;
+		*s = e - ei;
+		*c = e + ei;
+	}
+}
+
+long double complex
+csinl(long double complex z)
+{
+	long double complex w;
+	long double ch, sh;
+
+	cchshl(cimagl(z), &ch, &sh);
+	w = sinl(creall(z)) * ch + (cosl(creall(z)) * sh) * I;
+	return (w);
+}

src/s_ctan.c

@@ -0,0 +1,157 @@
+/*	$OpenBSD: s_ctan.c,v 1.6 2013/07/03 04:46:36 espie 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 <float.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);
+}
+
+#if	LDBL_MANT_DIG == DBL_MANT_DIG
+__strong_alias(ctanl, ctan);
+#endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */

src/s_ctanf.c

@@ -0,0 +1,148 @@
+/*	$OpenBSD: s_ctanf.c,v 1.2 2011/07/20 19:28:33 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.
+ */
+
+/*							ctanf()
+ *
+ *	Complex circular tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * void ctanf();
+ * cmplxf z, w;
+ *
+ * ctanf( &z, &w );
+ *
+ *
+ *
+ * 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.
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000       3.3e-7       5.1e-8
+ */
+
+#include <complex.h>
+#include <math.h>
+
+#define MACHEPF 3.0e-8
+#define MAXNUMF 1.0e38f
+
+static const double DP1 = 3.140625;
+static const double DP2 = 9.67502593994140625E-4;
+static const double DP3 = 1.509957990978376432E-7;
+
+static float
+_redupif(float xx)
+{
+	float x, t;
+	long i;
+
+	x = xx;
+	t = x/(float)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 float
+_ctansf(float complex z)
+{
+	float f, x, x2, y, y2, rn, t, d;
+
+	x = fabsf(2.0f * crealf(z));
+	y = fabsf(2.0f * cimagf(z));
+
+	x = _redupif(x);
+
+	x = x * x;
+	y = y * y;
+	x2 = 1.0f;
+	y2 = 1.0f;
+	f = 1.0f;
+	rn = 0.0f;
+	d = 0.0f;
+	do {
+		rn += 1.0f;
+		f *= rn;
+		rn += 1.0f;
+		f *= rn;
+		x2 *= x;
+		y2 *= y;
+		t = y2 + x2;
+		t /= f;
+		d += t;
+
+		rn += 1.0f;
+		f *= rn;
+		rn += 1.0f;
+		f *= rn;
+		x2 *= x;
+		y2 *= y;
+		t = y2 - x2;
+		t /= f;
+		d += t;
+	}
+	while (fabsf(t/d) > MACHEPF)
+		;
+	return(d);
+}
+
+float complex
+ctanf(float complex z)
+{
+	float complex w;
+	float d;
+
+	d = cosf( 2.0f * crealf(z) ) + coshf( 2.0f * cimagf(z) );
+
+	if(fabsf(d) < 0.25f)
+		d = _ctansf(z);
+
+	if (d == 0.0f) {
+		/*mtherr( "ctanf", OVERFLOW );*/
+		w = MAXNUMF + MAXNUMF * I;
+		return (w);
+	}
+	w = sinf (2.0f * crealf(z)) / d + (sinhf (2.0f * cimagf(z)) / d) * I;
+	return (w);
+}

src/s_ctanhl.c

@@ -0,0 +1,60 @@
+/*	$OpenBSD: s_ctanhl.c,v 1.2 2011/07/20 19:28:33 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.
+ */
+
+/*							ctanhl
+ *
+ *	Complex hyperbolic tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex ctanhl();
+ * long double complex z, w;
+ *
+ * w = ctanhl (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>
+
+long double complex
+ctanhl(long double complex z)
+{
+	long double complex w;
+	long double x, y, d;
+
+	x = creall(z);
+	y = cimagl(z);
+	d = coshl(2.0L * x) + cosl(2.0L * y);
+	w = sinhl(2.0L * x) / d + (sinl(2.0L * y) / d) * I;
+	return (w);
+}

src/s_ctanl.c

@@ -0,0 +1,157 @@
+/*	$OpenBSD: s_ctanl.c,v 1.3 2011/07/20 21:02:51 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.
+ */
+
+/*							ctanl()
+ *
+ *	Complex circular tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex ctanl();
+ * long double complex z, w;
+ *
+ * w = ctanl( 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.
+ *
+ *
+ * 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 <float.h>
+#include <math.h>
+
+#if	LDBL_MANT_DIG == 64
+static const long double MACHEPL= 5.42101086242752217003726400434970855712890625E-20L;
+#elif	LDBL_MANT_DIG == 113
+static const long double MACHEPL = 9.629649721936179265279889712924636592690508e-35L;
+#endif
+
+static const long double PIL = 3.141592653589793238462643383279502884197169L;
+static const long double DP1 = 3.14159265358979323829596852490908531763125L;
+static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L;
+static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L;
+
+static long double
+redupil(long double x)
+{
+	long double t;
+	long i;
+
+	t = x / PIL;
+	if (t >= 0.0L)
+		t += 0.5L;
+	else
+		t -= 0.5L;
+
+	i = t;	/* the multiple */
+	t = i;
+	t = ((x - t * DP1) - t * DP2) - t * DP3;
+	return (t);
+}
+
+static long double
+ctansl(long double complex z)
+{
+	long double f, x, x2, y, y2, rn, t;
+	long double d;
+
+	x = fabsl(2.0L * creall(z));
+	y = fabsl(2.0L * cimagl(z));
+
+	x = redupil(x);
+
+	x = x * x;
+	y = y * y;
+	x2 = 1.0L;
+	y2 = 1.0L;
+	f = 1.0L;
+	rn = 0.0L;
+	d = 0.0L;
+	do {
+		rn += 1.0L;
+		f *= rn;
+		rn += 1.0L;
+		f *= rn;
+		x2 *= x;
+		y2 *= y;
+		t = y2 + x2;
+		t /= f;
+		d += t;
+
+		rn += 1.0L;
+		f *= rn;
+		rn += 1.0L;
+		f *= rn;
+		x2 *= x;
+		y2 *= y;
+		t = y2 - x2;
+		t /= f;
+		d += t;
+	}
+	while (fabsl(t/d) > MACHEPL);
+	return(d);
+}
+
+long double complex
+ctanl(long double complex z)
+{
+	long double complex w;
+	long double d, x, y;
+
+	x = creall(z);
+	y = cimagl(z);
+	d = cosl(2.0L * x) + coshl(2.0L * y);
+
+	if (fabsl(d) < 0.25L) {
+		d = fabsl(d);
+		d = ctansl(z);
+	}
+	if (d == 0.0L) {
+		/*mtherr( "ctan", OVERFLOW );*/
+		w = LDBL_MAX + LDBL_MAX * I;
+		return (w);
+	}
+
+	w = sinl(2.0L * x) / d + (sinhl(2.0L * y) / d) * I;
+	return (w);
+}