extend cmath (some of the functions are dummy)
authornsz <nsz@port70.net>
Sun, 11 Mar 2012 14:24:11 +0000 (15:24 +0100)
committernsz <nsz@port70.net>
Sun, 11 Mar 2012 14:24:11 +0000 (15:24 +0100)
47 files changed:
src/cmath/cabsl.c
src/cmath/cacos.c
src/cmath/cacosf.c [new file with mode: 0644]
src/cmath/cacosh.c
src/cmath/cacoshf.c [new file with mode: 0644]
src/cmath/cacoshl.c [new file with mode: 0644]
src/cmath/cacosl.c [new file with mode: 0644]
src/cmath/cargl.c
src/cmath/casin.c
src/cmath/casinf.c [new file with mode: 0644]
src/cmath/casinh.c
src/cmath/casinhf.c [new file with mode: 0644]
src/cmath/casinhl.c [new file with mode: 0644]
src/cmath/casinl.c [new file with mode: 0644]
src/cmath/catanf.c [new file with mode: 0644]
src/cmath/catanh.c
src/cmath/catanhf.c [new file with mode: 0644]
src/cmath/catanhl.c [new file with mode: 0644]
src/cmath/catanl.c [new file with mode: 0644]
src/cmath/ccos.c
src/cmath/ccosf.c [new file with mode: 0644]
src/cmath/ccoshf.c [new file with mode: 0644]
src/cmath/ccosl.c [new file with mode: 0644]
src/cmath/cexpf.c [new file with mode: 0644]
src/cmath/cimagf.c [new file with mode: 0644]
src/cmath/cimagl.c [new file with mode: 0644]
src/cmath/clog.c
src/cmath/clogf.c [new file with mode: 0644]
src/cmath/clogl.c [new file with mode: 0644]
src/cmath/conj.c
src/cmath/conjf.c [new file with mode: 0644]
src/cmath/conjl.c [new file with mode: 0644]
src/cmath/cpow.c
src/cmath/cpowf.c [new file with mode: 0644]
src/cmath/cpowl.c [new file with mode: 0644]
src/cmath/cproj.c
src/cmath/cprojf.c [new file with mode: 0644]
src/cmath/cprojl.c [new file with mode: 0644]
src/cmath/csin.c
src/cmath/csinf.c [new file with mode: 0644]
src/cmath/csinhf.c [new file with mode: 0644]
src/cmath/csinl.c [new file with mode: 0644]
src/cmath/csqrtf.c [new file with mode: 0644]
src/cmath/ctan.c
src/cmath/ctanf.c [new file with mode: 0644]
src/cmath/ctanhf.c [new file with mode: 0644]
src/cmath/ctanl.c [new file with mode: 0644]

index 31ac460..40a067c 100644 (file)
@@ -1,6 +1,13 @@
 #include "libm.h"
 
 #include "libm.h"
 
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double cabsl(long double complex z)
+{
+       return cabs(z);
+}
+#else
 long double cabsl(long double complex z)
 {
        return hypotl(creall(z), cimagl(z));
 }
 long double cabsl(long double complex z)
 {
        return hypotl(creall(z), cimagl(z));
 }
+#endif
index b8a2d71..3aca051 100644 (file)
@@ -1,52 +1,11 @@
-/* origin: OpenBSD /usr/src/lib/libm/src/s_cacos.c */
-/*
- * 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.
- */
-/*
- *      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 "libm.h"
 
 #include "libm.h"
 
+// FIXME: Hull et al. "Implementing the complex arcsine and arccosine functions using exception handling" 1997
+
+/* acos(z) = pi/2 - asin(z) */
+
 double complex cacos(double complex z)
 {
 double complex cacos(double complex z)
 {
-       double complex w;
-
-       w = casin(z);
-       w = (M_PI_2 - creal(w)) - cimag(w) * I;
-       return w;
+       z = casin(z);
+       return cpack(M_PI_2 - creal(z), -cimag(z));
 }
 }
diff --git a/src/cmath/cacosf.c b/src/cmath/cacosf.c
new file mode 100644 (file)
index 0000000..563766e
--- /dev/null
@@ -0,0 +1,9 @@
+#include "libm.h"
+
+// FIXME
+
+float complex cacosf(float complex z)
+{
+       z = casinf(z);
+       return cpackf((float)M_PI_2 - crealf(z), -cimagf(z));
+}
index 89cbb1f..c2dfc1b 100644 (file)
@@ -1,48 +1,9 @@
-/* origin: OpenBSD /usr/src/lib/libm/src/s_cacosh.c */
-/*
- * 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.
- */
-/*
- *      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
- *
- */
-
-// FIXME: reeval
-
 #include "libm.h"
 
 #include "libm.h"
 
+/* acosh(z) = i acos(z) */
+
 double complex cacosh(double complex z)
 {
 double complex cacosh(double complex z)
 {
-       return I * cacos(z);
+       z = cacos(z);
+       return cpack(-cimag(z), creal(z));
 }
 }
diff --git a/src/cmath/cacoshf.c b/src/cmath/cacoshf.c
new file mode 100644 (file)
index 0000000..37ff880
--- /dev/null
@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float complex cacoshf(float complex z)
+{
+       z = cacosf(z);
+       return cpackf(-cimagf(z), crealf(z));
+}
diff --git a/src/cmath/cacoshl.c b/src/cmath/cacoshl.c
new file mode 100644 (file)
index 0000000..2a04e27
--- /dev/null
@@ -0,0 +1,14 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex cacoshl(long double complex z)
+{
+       return cacosh(z);
+}
+#else
+long double complex cacoshl(long double complex z)
+{
+       z = cacosl(z);
+       return cpackl(-cimagl(z), creall(z));
+}
+#endif
diff --git a/src/cmath/cacosl.c b/src/cmath/cacosl.c
new file mode 100644 (file)
index 0000000..5992e05
--- /dev/null
@@ -0,0 +1,16 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex cacosl(long double complex z)
+{
+       return cacos(z);
+}
+#else
+// FIXME
+#define PI_2 1.57079632679489661923132169163975144L
+long double complex cacosl(long double complex z)
+{
+       z = casinl(z);
+       return cpackl(PI_2 - creall(z), -cimagl(z));
+}
+#endif
index 8aacdfa..e0d5047 100644 (file)
@@ -1,6 +1,13 @@
 #include "libm.h"
 
 #include "libm.h"
 
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double cargl(long double complex z)
+{
+       return carg(z);
+}
+#else
 long double cargl(long double complex z)
 {
        return atan2l(cimagl(z), creall(z));
 }
 long double cargl(long double complex z)
 {
        return atan2l(cimagl(z), creall(z));
 }
+#endif
index 3428457..79aff27 100644 (file)
-/* origin: OpenBSD /usr/src/lib/libm/src/s_casin.c */
-/*
- * 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.
- */
-/*
- *      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 "libm.h"
 
 #include "libm.h"
 
+// FIXME
+
+/* asin(z) = -i log(i z + sqrt(1 - z*z)) */
+
 double complex casin(double complex z)
 {
        double complex w;
 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);
        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;
-               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 */
-       // FIXME: 1 - ()() - () I
-       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 */
-       // FIXME: -I
-       w = zz * (-1.0 * I);
-       return w;
+       w = cpack(1.0 - (x - y)*(x + y), -2.0*x*y);
+       return clog(cpack(-y, x) + csqrt(w));
 }
 }
diff --git a/src/cmath/casinf.c b/src/cmath/casinf.c
new file mode 100644 (file)
index 0000000..cb9863f
--- /dev/null
@@ -0,0 +1,14 @@
+#include "libm.h"
+
+// FIXME
+
+float complex casinf(float complex z)
+{
+       float complex w;
+       float x, y;
+
+       x = crealf(z);
+       y = cimagf(z);
+       w = cpackf(1.0 - (x - y)*(x + y), -2.0*x*y);
+       return clogf(cpackf(-y, x) + csqrtf(w));
+}
index 8c8d9fa..f2b3fef 100644 (file)
@@ -1,47 +1,9 @@
-/* origin: OpenBSD /usr/src/lib/libm/src/s_casinh.c */
-/*
- * 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.
- */
-/*
- *      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 "libm.h"
 
 #include "libm.h"
 
+/* asinh(z) = -i asin(i z) */
+
 double complex casinh(double complex z)
 {
 double complex casinh(double complex z)
 {
-       // FIXME: -I, reeval
-       return -1.0 * I * casin(z * I);
+       z = casin(cpack(-cimag(z), creal(z)));
+       return cpack(cimag(z), -creal(z));
 }
 }
diff --git a/src/cmath/casinhf.c b/src/cmath/casinhf.c
new file mode 100644 (file)
index 0000000..ed4af64
--- /dev/null
@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float complex casinhf(float complex z)
+{
+       z = casinf(cpackf(-cimagf(z), crealf(z)));
+       return cpackf(cimagf(z), -crealf(z));
+}
diff --git a/src/cmath/casinhl.c b/src/cmath/casinhl.c
new file mode 100644 (file)
index 0000000..e5d80ce
--- /dev/null
@@ -0,0 +1,14 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex casinhl(long double complex z)
+{
+       return casinh(z);
+}
+#else
+long double complex casinhl(long double complex z)
+{
+       z = casinl(cpackl(-cimagl(z), creall(z)));
+       return cpackl(cimagl(z), -creall(z));
+}
+#endif
diff --git a/src/cmath/casinl.c b/src/cmath/casinl.c
new file mode 100644 (file)
index 0000000..f9aa8de
--- /dev/null
@@ -0,0 +1,20 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex casinl(long double complex z)
+{
+       return casin(z);
+}
+#else
+// FIXME
+long double complex casinl(long double complex z)
+{
+       long double complex w;
+       long double x, y;
+
+       x = creall(z);
+       y = cimagl(z);
+       w = cpackl(1.0 - (x - y)*(x + y), -2.0*x*y);
+       return clogl(cpackl(-y, x) + csqrtl(w));
+}
+#endif
diff --git a/src/cmath/catanf.c b/src/cmath/catanf.c
new file mode 100644 (file)
index 0000000..8533bde
--- /dev/null
@@ -0,0 +1,115 @@
+/* origin: OpenBSD /usr/src/lib/libm/src/s_catanf.c */
+/*
+ * 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.
+ */
+/*
+ *      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 "libm.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.0f)
+               t += 0.5f;
+       else
+               t -= 0.5f;
+
+       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:
+       // FIXME
+       w = MAXNUMF + MAXNUMF * I;
+       return w;
+}
index 5545747..b162802 100644 (file)
@@ -1,47 +1,9 @@
-/* origin: OpenBSD /usr/src/lib/libm/src/s_catanh.c */
-/*
- * 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.
- */
-/*
- *      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 "libm.h"
 
 #include "libm.h"
 
+/* atanh = -i atan(i z) */
+
 double complex catanh(double complex z)
 {
 double complex catanh(double complex z)
 {
-       // FIXME: -i, reeval
-       return -1.0 * I * catan(z * I);
+       z = catan(cpack(-cimag(z), creal(z)));
+       return cpack(cimag(z), -creal(z));
 }
 }
diff --git a/src/cmath/catanhf.c b/src/cmath/catanhf.c
new file mode 100644 (file)
index 0000000..e1d1e64
--- /dev/null
@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float complex catanhf(float complex z)
+{
+       z = catanf(cpackf(-cimagf(z), crealf(z)));
+       return cpackf(cimagf(z), -crealf(z));
+}
diff --git a/src/cmath/catanhl.c b/src/cmath/catanhl.c
new file mode 100644 (file)
index 0000000..0a9374a
--- /dev/null
@@ -0,0 +1,14 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex catanhl(long double complex z)
+{
+       return catanh(z);
+}
+#else
+long double complex catanhl(long double complex z)
+{
+       z = catanl(cpackl(-cimagl(z), creall(z)));
+       return cpackl(cimagl(z), -creall(z));
+}
+#endif
diff --git a/src/cmath/catanl.c b/src/cmath/catanl.c
new file mode 100644 (file)
index 0000000..5ace770
--- /dev/null
@@ -0,0 +1,126 @@
+/* origin: OpenBSD /usr/src/lib/libm/src/s_catanl.c */
+/*
+ * 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.
+ */
+/*
+ *      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 "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex catanl(long double complex z)
+{
+       return catan(z);
+}
+#else
+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:
+       // FIXME
+       w = LDBL_MAX + LDBL_MAX * I;
+       return w;
+}
+#endif
index 9757e6a..5754c23 100644 (file)
@@ -1,5 +1,7 @@
 #include "libm.h"
 
 #include "libm.h"
 
+/* cos(z) = cosh(i z) */
+
 double complex ccos(double complex z)
 {
        return ccosh(cpack(-cimag(z), creal(z)));
 double complex ccos(double complex z)
 {
        return ccosh(cpack(-cimag(z), creal(z)));
diff --git a/src/cmath/ccosf.c b/src/cmath/ccosf.c
new file mode 100644 (file)
index 0000000..9b72c4f
--- /dev/null
@@ -0,0 +1,6 @@
+#include "libm.h"
+
+float complex ccosf(float complex z)
+{
+       return ccoshf(cpackf(-cimagf(z), crealf(z)));
+}
diff --git a/src/cmath/ccoshf.c b/src/cmath/ccoshf.c
new file mode 100644 (file)
index 0000000..683e77f
--- /dev/null
@@ -0,0 +1,90 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_ccoshf.c */
+/*-
+ * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
+ * 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 unmodified, 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.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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.
+ */
+/*
+ * Hyperbolic cosine of a complex argument.  See s_ccosh.c for details.
+ */
+
+#include "libm.h"
+
+static const float huge = 0x1p127;
+
+float complex ccoshf(float complex z)
+{
+       float x, y, h;
+       int32_t hx, hy, ix, iy;
+
+       x = crealf(z);
+       y = cimagf(z);
+
+       GET_FLOAT_WORD(hx, x);
+       GET_FLOAT_WORD(hy, y);
+
+       ix = 0x7fffffff & hx;
+       iy = 0x7fffffff & hy;
+
+       if (ix < 0x7f800000 && iy < 0x7f800000) {
+               if (iy == 0)
+                       return cpackf(coshf(x), x * y);
+               if (ix < 0x41100000)    /* small x: normal case */
+                       return cpackf(coshf(x) * cosf(y), sinhf(x) * sinf(y));
+
+               /* |x| >= 9, so cosh(x) ~= exp(|x|) */
+               if (ix < 0x42b17218) {
+                       /* x < 88.7: expf(|x|) won't overflow */
+                       h = expf(fabsf(x)) * 0.5f;
+                       return cpackf(h * cosf(y), copysignf(h, x) * sinf(y));
+               } else if (ix < 0x4340b1e7) {
+                       /* x < 192.7: scale to avoid overflow */
+                       z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
+                       return cpackf(crealf(z), cimagf(z) * copysignf(1, x));
+               } else {
+                       /* x >= 192.7: the result always overflows */
+                       h = huge * x;
+                       return cpackf(h * h * cosf(y), h * sinf(y));
+               }
+       }
+
+       if (ix == 0 && iy >= 0x7f800000)
+               return cpackf(y - y, copysignf(0, x * (y - y)));
+
+       if (iy == 0 && ix >= 0x7f800000) {
+               if ((hx & 0x7fffff) == 0)
+                       return cpackf(x * x, copysignf(0, x) * y);
+               return cpackf(x * x, copysignf(0, (x + x) * y));
+       }
+
+       if (ix < 0x7f800000 && iy >= 0x7f800000)
+               return cpackf(y - y, x * (y - y));
+
+       if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
+               if (iy >= 0x7f800000)
+                       return cpackf(x * x, x * (y - y));
+               return cpackf((x * x) * cosf(y), x * sinf(y));
+       }
+
+       return cpackf((x * x) * (y - y), (x + x) * (y - y));
+}
diff --git a/src/cmath/ccosl.c b/src/cmath/ccosl.c
new file mode 100644 (file)
index 0000000..e37825a
--- /dev/null
@@ -0,0 +1,13 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex ccosl(long double complex z)
+{
+       return ccos(z);
+}
+#else
+long double complex ccosl(long double complex z)
+{
+       return ccoshl(cpackl(-cimagl(z), creall(z)));
+}
+#endif
diff --git a/src/cmath/cexpf.c b/src/cmath/cexpf.c
new file mode 100644 (file)
index 0000000..0cf13a3
--- /dev/null
@@ -0,0 +1,83 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_cexpf.c */
+/*-
+ * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
+ * 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.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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 AUTHOR 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.
+ */
+
+#include "libm.h"
+
+static const uint32_t
+exp_ovfl  = 0x42b17218,  /* MAX_EXP * ln2 ~= 88.722839355 */
+cexp_ovfl = 0x43400074;  /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
+
+float complex cexpf(float complex z)
+{
+       float x, y, exp_x;
+       uint32_t hx, hy;
+
+       x = crealf(z);
+       y = cimagf(z);
+
+       GET_FLOAT_WORD(hy, y);
+       hy &= 0x7fffffff;
+
+       /* cexp(x + I 0) = exp(x) + I 0 */
+       if (hy == 0)
+               return cpackf(expf(x), y);
+       GET_FLOAT_WORD(hx, x);
+       /* cexp(0 + I y) = cos(y) + I sin(y) */
+       if ((hx & 0x7fffffff) == 0)
+               return cpackf(cosf(y), sinf(y));
+
+       if (hy >= 0x7f800000) {
+               if ((hx & 0x7fffffff) != 0x7f800000) {
+                       /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
+                       return cpackf(y - y, y - y);
+               } else if (hx & 0x80000000) {
+                       /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
+                       return cpackf(0.0, 0.0);
+               } else {
+                       /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
+                       return cpackf(x, y - y);
+               }
+       }
+
+       if (hx >= exp_ovfl && hx <= cexp_ovfl) {
+               /*
+                * x is between 88.7 and 192, so we must scale to avoid
+                * overflow in expf(x).
+                */
+               return __ldexp_cexpf(z, 0);
+       } else {
+               /*
+                * Cases covered here:
+                *  -  x < exp_ovfl and exp(x) won't overflow (common case)
+                *  -  x > cexp_ovfl, so exp(x) * s overflows for all s > 0
+                *  -  x = +-Inf (generated by exp())
+                *  -  x = NaN (spurious inexact exception from y)
+                */
+               exp_x = expf(x);
+               return cpackf(exp_x * cosf(y), exp_x * sinf(y));
+       }
+}
diff --git a/src/cmath/cimagf.c b/src/cmath/cimagf.c
new file mode 100644 (file)
index 0000000..99fffc5
--- /dev/null
@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float (cimagf)(float complex z)
+{
+       union fcomplex u = {z};
+       return u.a[1];
+}
diff --git a/src/cmath/cimagl.c b/src/cmath/cimagl.c
new file mode 100644 (file)
index 0000000..d9a0780
--- /dev/null
@@ -0,0 +1,7 @@
+#include "libm.h"
+
+long double (cimagl)(long double complex z)
+{
+       union lcomplex u = {z};
+       return u.a[1];
+}
index 4f8a428..6f10a11 100644 (file)
@@ -1,5 +1,9 @@
 #include "libm.h"
 
 #include "libm.h"
 
+// FIXME
+
+/* log(z) = log(|z|) + i arg(z) */
+
 double complex clog(double complex z)
 {
        double r, phi;
 double complex clog(double complex z)
 {
        double r, phi;
diff --git a/src/cmath/clogf.c b/src/cmath/clogf.c
new file mode 100644 (file)
index 0000000..f3aec54
--- /dev/null
@@ -0,0 +1,12 @@
+#include "libm.h"
+
+// FIXME
+
+float complex clogf(float complex z)
+{
+       float r, phi;
+
+       r = cabsf(z);
+       phi = cargf(z);
+       return cpackf(logf(r), phi);
+}
diff --git a/src/cmath/clogl.c b/src/cmath/clogl.c
new file mode 100644 (file)
index 0000000..5b84ba5
--- /dev/null
@@ -0,0 +1,18 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex clogl(long double complex z)
+{
+       return clog(z);
+}
+#else
+// FIXME
+long double complex clogl(long double complex z)
+{
+       long double r, phi;
+
+       r = cabsl(z);
+       phi = cargl(z);
+       return cpackl(logl(r), phi);
+}
+#endif
index 8720feb..4aceea7 100644 (file)
@@ -1,23 +1,6 @@
 #include "libm.h"
 
 #include "libm.h"
 
-// nice, but internally ugly:
-// macros generate gratuitous code the compiler has to optimize away
 double complex conj(double complex z)
 {
        return cpack(creal(z), -cimag(z));
 }
 double complex conj(double complex z)
 {
        return cpack(creal(z), -cimag(z));
 }
-
-/*
-// always works, but ugly union
-double complex conj(double complex z) {
-       union dcomplex u = {z};
-
-       u.a[1] = -u.a[1];
-       return u.z;
-}
-
-// reasonable, needs clever compiler that understands *I
-double complex conj(double complex z) {
-       return creal(z) - cimag(z)*I;
-}
-*/
diff --git a/src/cmath/conjf.c b/src/cmath/conjf.c
new file mode 100644 (file)
index 0000000..3155680
--- /dev/null
@@ -0,0 +1,6 @@
+#include "libm.h"
+
+float complex conjf(float complex z)
+{
+       return cpackf(crealf(z), -cimagf(z));
+}
diff --git a/src/cmath/conjl.c b/src/cmath/conjl.c
new file mode 100644 (file)
index 0000000..0133226
--- /dev/null
@@ -0,0 +1,6 @@
+#include "libm.h"
+
+long double complex conjl(long double complex z)
+{
+       return cpackl(creall(z), -cimagl(z));
+}
index 0bc36d3..f863588 100644 (file)
@@ -1,61 +1,8 @@
-/* origin: OpenBSD /usr/src/lib/libm/src/s_cpow.c */
-/*
- * 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.
- */
-/*
- *      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 "libm.h"
 
 #include "libm.h"
 
-double complex cpow(double complex a, double complex z)
-{
-       double x, y, r, theta, absa, arga;
+/* pow(z, c) = exp(c log(z)), See C99 G.6.4.1 */
 
 
-       x = creal(z);
-       y = cimag(z);
-       absa = cabs(a);
-       if (absa == 0.0)
-               return 0.0;
-       arga = carg(a);
-       r = pow(absa, x);
-       theta = x * arga;
-       if (y != 0.0) {
-               r = r * exp(-y * arga);
-               theta += y * log(absa);
-       }
-       return cpack(r * cos(theta), (r * sin(theta)));
+double complex cpow(double complex z, double complex c)
+{
+       return cexp(c * clog(z));
 }
 }
diff --git a/src/cmath/cpowf.c b/src/cmath/cpowf.c
new file mode 100644 (file)
index 0000000..53c65dc
--- /dev/null
@@ -0,0 +1,6 @@
+#include "libm.h"
+
+float complex cpowf(float complex z, float complex c)
+{
+       return cexpf(c * clogf(z));
+}
diff --git a/src/cmath/cpowl.c b/src/cmath/cpowl.c
new file mode 100644 (file)
index 0000000..c1a80a7
--- /dev/null
@@ -0,0 +1,13 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex cpowl(long double complex z, long double complex c)
+{
+       return cpow(z, c);
+}
+#else
+long double complex cpowl(long double complex z, long double complex c)
+{
+       return cexpl(c * clogl(z));
+}
+#endif
index cc12e46..1cf9bb9 100644 (file)
@@ -2,9 +2,7 @@
 
 double complex cproj(double complex z)
 {
 
 double complex cproj(double complex z)
 {
-       union dcomplex u = {z};
-
-       if (isinf(u.a[0]) || isinf(u.a[1]))
-               return cpack(INFINITY, copysign(0.0, u.a[1]));
+       if (isinf(creal(z)) || isinf(cimag(z)))
+               return cpack(INFINITY, copysign(0.0, creal(z)));
        return z;
 }
        return z;
 }
diff --git a/src/cmath/cprojf.c b/src/cmath/cprojf.c
new file mode 100644 (file)
index 0000000..7112974
--- /dev/null
@@ -0,0 +1,8 @@
+#include "libm.h"
+
+float complex cprojf(float complex z)
+{
+       if (isinf(crealf(z)) || isinf(cimagf(z)))
+               return cpackf(INFINITY, copysignf(0.0, crealf(z)));
+       return z;
+}
diff --git a/src/cmath/cprojl.c b/src/cmath/cprojl.c
new file mode 100644 (file)
index 0000000..72e50cf
--- /dev/null
@@ -0,0 +1,15 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex cprojl(long double complex z)
+{
+       return cproj(z);
+}
+#else
+long double complex cprojl(long double complex z)
+{
+       if (isinf(creall(z)) || isinf(cimagl(z)))
+               return cpackl(INFINITY, copysignl(0.0, creall(z)));
+       return z;
+}
+#endif
index ab72055..246a459 100644 (file)
@@ -1,9 +1,9 @@
 #include "libm.h"
 
 #include "libm.h"
 
+/* sin(z) = -i sinh(i z) */
+
 double complex csin(double complex z)
 {
 double complex csin(double complex z)
 {
-       double complex r;
-
-       r = csinh(cpack(-cimag(z), creal(z)));
-       return cpack(cimag(r), -creal(r));
+       z = csinh(cpack(-cimag(z), creal(z)));
+       return cpack(cimag(z), -creal(z));
 }
 }
diff --git a/src/cmath/csinf.c b/src/cmath/csinf.c
new file mode 100644 (file)
index 0000000..3aabe8f
--- /dev/null
@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float complex csinf(float complex z)
+{
+       z = csinhf(cpackf(-cimagf(z), crealf(z)));
+       return cpackf(cimagf(z), -crealf(z));
+}
diff --git a/src/cmath/csinhf.c b/src/cmath/csinhf.c
new file mode 100644 (file)
index 0000000..1754186
--- /dev/null
@@ -0,0 +1,105 @@
+/*-
+ * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
+ * 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 unmodified, 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.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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.
+ */
+
+/*
+ * Hyperbolic sine of a complex argument z.  See s_csinh.c for details.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <complex.h>
+#include <math.h>
+
+#include "math_private.h"
+
+static const float huge = 0x1p127;
+
+float complex
+csinhf(float complex z)
+{
+       float x, y, h;
+       int32_t hx, hy, ix, iy;
+
+       x = crealf(z);
+       y = cimagf(z);
+
+       GET_FLOAT_WORD(hx, x);
+       GET_FLOAT_WORD(hy, y);
+
+       ix = 0x7fffffff & hx;
+       iy = 0x7fffffff & hy;
+
+       if (ix < 0x7f800000 && iy < 0x7f800000) {
+               if (iy == 0)
+                       return (cpackf(sinhf(x), y));
+               if (ix < 0x41100000)    /* small x: normal case */
+                       return (cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y)));
+
+               /* |x| >= 9, so cosh(x) ~= exp(|x|) */
+               if (ix < 0x42b17218) {
+                       /* x < 88.7: expf(|x|) won't overflow */
+                       h = expf(fabsf(x)) * 0.5f;
+                       return (cpackf(copysignf(h, x) * cosf(y), h * sinf(y)));
+               } else if (ix < 0x4340b1e7) {
+                       /* x < 192.7: scale to avoid overflow */
+                       z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
+                       return (cpackf(crealf(z) * copysignf(1, x), cimagf(z)));
+               } else {
+                       /* x >= 192.7: the result always overflows */
+                       h = huge * x;
+                       return (cpackf(h * cosf(y), h * h * sinf(y)));
+               }
+       }
+
+       if (ix == 0 && iy >= 0x7f800000)
+               return (cpackf(copysignf(0, x * (y - y)), y - y));
+
+       if (iy == 0 && ix >= 0x7f800000) {
+               if ((hx & 0x7fffff) == 0)
+                       return (cpackf(x, y));
+               return (cpackf(x, copysignf(0, y)));
+       }
+
+       if (ix < 0x7f800000 && iy >= 0x7f800000)
+               return (cpackf(y - y, x * (y - y)));
+
+       if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
+               if (iy >= 0x7f800000)
+                       return (cpackf(x * x, x * (y - y)));
+               return (cpackf(x * cosf(y), INFINITY * sinf(y)));
+       }
+
+       return (cpackf((x * x) * (y - y), (x + x) * (y - y)));
+}
+
+float complex
+csinf(float complex z)
+{
+
+       z = csinhf(cpackf(-cimagf(z), crealf(z)));
+       return (cpackf(cimagf(z), -crealf(z)));
+}
diff --git a/src/cmath/csinl.c b/src/cmath/csinl.c
new file mode 100644 (file)
index 0000000..4ad8674
--- /dev/null
@@ -0,0 +1,14 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex csinl(long double complex z)
+{
+       return csin(z);
+}
+#else
+long double complex csinl(long double complex z)
+{
+       z = csinhl(cpackl(-cimagl(z), creall(z)));
+       return cpackl(cimagl(z), -creall(z));
+}
+#endif
diff --git a/src/cmath/csqrtf.c b/src/cmath/csqrtf.c
new file mode 100644 (file)
index 0000000..16487c2
--- /dev/null
@@ -0,0 +1,82 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_csqrtf.c */
+/*-
+ * Copyright (c) 2007 David Schultz <das@FreeBSD.ORG>
+ * 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.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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 AUTHOR 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.
+ */
+
+#include "libm.h"
+
+/*
+ * gcc doesn't implement complex multiplication or division correctly,
+ * so we need to handle infinities specially. We turn on this pragma to
+ * notify conforming c99 compilers that the fast-but-incorrect code that
+ * gcc generates is acceptable, since the special cases have already been
+ * handled.
+ */
+#pragma STDC CX_LIMITED_RANGE ON
+
+float complex csqrtf(float complex z)
+{
+       float a = crealf(z), b = cimagf(z);
+       double t;
+
+       /* Handle special cases. */
+       if (z == 0)
+               return cpackf(0, b);
+       if (isinf(b))
+               return cpackf(INFINITY, b);
+       if (isnan(a)) {
+               t = (b - b) / (b - b);  /* raise invalid if b is not a NaN */
+               return cpackf(a, t);  /* return NaN + NaN i */
+       }
+       if (isinf(a)) {
+               /*
+                * csqrtf(inf + NaN i)  = inf +  NaN i
+                * csqrtf(inf + y i)    = inf +  0 i
+                * csqrtf(-inf + NaN i) = NaN +- inf i
+                * csqrtf(-inf + y i)   = 0   +  inf i
+                */
+               if (signbit(a))
+                       return cpackf(fabsf(b - b), copysignf(a, b));
+               else
+                       return cpackf(a, copysignf(b - b, b));
+       }
+       /*
+        * The remaining special case (b is NaN) is handled just fine by
+        * the normal code path below.
+        */
+
+       /*
+        * We compute t in double precision to avoid overflow and to
+        * provide correct rounding in nearly all cases.
+        * This is Algorithm 312, CACM vol 10, Oct 1967.
+        */
+       if (a >= 0) {
+               t = sqrt((a + hypot(a, b)) * 0.5);
+               return cpackf(t, b / (2.0 * t));
+       } else {
+               t = sqrt((-a + hypot(a, b)) * 0.5);
+               return cpackf(fabsf(b) / (2.0 * t), copysignf(t, b));
+       }
+}
index 4462ef0..4741a4d 100644 (file)
@@ -1,9 +1,9 @@
 #include "libm.h"
 
 #include "libm.h"
 
+/* tan(z) = -i tanh(i z) */
+
 double complex ctan(double complex z)
 {
 double complex ctan(double complex z)
 {
-       double complex r;
-
-       r = ctanh(cpack(-cimag(z), creal(z)));
-       return cpack(cimag(r), -creal(r));
+       z = ctanh(cpack(-cimag(z), creal(z)));
+       return cpack(cimag(z), -creal(z));
 }
 }
diff --git a/src/cmath/ctanf.c b/src/cmath/ctanf.c
new file mode 100644 (file)
index 0000000..9bbeb05
--- /dev/null
@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float complex ctanf(float complex z)
+{
+       z = ctanhf(cpackf(-cimagf(z), crealf(z)));
+       return cpackf(cimagf(z), -crealf(z));
+}
diff --git a/src/cmath/ctanhf.c b/src/cmath/ctanhf.c
new file mode 100644 (file)
index 0000000..7d74613
--- /dev/null
@@ -0,0 +1,66 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_ctanhf.c */
+/*-
+ * Copyright (c) 2011 David Schultz
+ * 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 unmodified, 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.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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.
+ */
+/*
+ * Hyperbolic tangent of a complex argument z.  See s_ctanh.c for details.
+ */
+
+#include "libm.h"
+
+float complex ctanhf(float complex z)
+{
+       float x, y;
+       float t, beta, s, rho, denom;
+       uint32_t hx, ix;
+
+       x = crealf(z);
+       y = cimagf(z);
+
+       GET_FLOAT_WORD(hx, x);
+       ix = hx & 0x7fffffff;
+
+       if (ix >= 0x7f800000) {
+               if (ix & 0x7fffff)
+                       return cpackf(x, (y == 0 ? y : x * y));
+               SET_FLOAT_WORD(x, hx - 0x40000000);
+               return cpackf(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)));
+       }
+
+       if (!isfinite(y))
+               return cpackf(y - y, y - y);
+
+       if (ix >= 0x41300000) { /* x >= 11 */
+               float exp_mx = expf(-fabsf(x));
+               return cpackf(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx);
+       }
+
+       t = tanf(y);
+       beta = 1.0 + t * t;
+       s = sinhf(x);
+       rho = sqrtf(1 + s * s);
+       denom = 1 + beta * s * s;
+       return cpackf((beta * rho * s) / denom, t / denom);
+}
diff --git a/src/cmath/ctanl.c b/src/cmath/ctanl.c
new file mode 100644 (file)
index 0000000..4b4c99b
--- /dev/null
@@ -0,0 +1,14 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex ctanl(long double complex z)
+{
+       return ctan(z);
+}
+#else
+long double complex ctanl(long double complex z)
+{
+       z = ctanhl(cpackl(-cimagl(z), creall(z)));
+       return cpackl(cimagl(z), -creall(z));
+}
+#endif