--- /dev/null
+#include "libm.h"
+
+double cabs(double complex z)
+{
+ return hypot(creal(z), cimag(z));
+}
--- /dev/null
+#include "libm.h"
+
+float cabsf(float complex z)
+{
+ return hypotf(crealf(z), cimagf(z));
+}
--- /dev/null
+#include "libm.h"
+
+long double cabsl(long double complex z)
+{
+ return hypotl(creall(z), cimagl(z));
+}
--- /dev/null
+/* 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"
+
+double complex cacos(double complex z)
+{
+ double complex w;
+
+ w = casin(z);
+ w = (M_PI_2 - creal(w)) - cimag(w) * I;
+ return w;
+}
--- /dev/null
+/* 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"
+
+double complex cacosh(double complex z)
+{
+ return I * cacos(z);
+}
--- /dev/null
+#include "libm.h"
+
+double carg(double complex z)
+{
+ return atan2(cimag(z), creal(z));
+}
--- /dev/null
+#include "libm.h"
+
+float cargf(float complex z)
+{
+ return atan2f(cimagf(z), crealf(z));
+}
--- /dev/null
+#include "libm.h"
+
+long double cargl(long double complex z)
+{
+ return atan2l(cimagl(z), creall(z));
+}
--- /dev/null
+/* 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"
+
+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;
+ 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;
+}
--- /dev/null
+/* 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"
+
+double complex casinh(double complex z)
+{
+ // FIXME: -I, reeval
+ return -1.0 * I * casin(z * I);
+}
--- /dev/null
+/* origin: OpenBSD /usr/src/lib/libm/src/s_catan.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:
+ *
+ * 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 "libm.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:
+ // FIXME
+ w = MAXNUM + MAXNUM * I;
+ return w;
+}
--- /dev/null
+/* 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"
+
+double complex catanh(double complex z)
+{
+ // FIXME: -i, reeval
+ return -1.0 * I * catan(z * I);
+}
--- /dev/null
+#include "libm.h"
+
+double complex ccos(double complex z)
+{
+ return ccosh(cpack(-cimag(z), creal(z)));
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_ccosh.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 z = x + i y.
+ *
+ * cosh(z) = cosh(x+iy)
+ * = cosh(x) cos(y) + i sinh(x) sin(y).
+ *
+ * Exceptional values are noted in the comments within the source code.
+ * These values and the return value were taken from n1124.pdf.
+ */
+
+#include "libm.h"
+
+static const double huge = 0x1p1023;
+
+double complex ccosh(double complex z)
+{
+ double x, y, h;
+ int32_t hx, hy, ix, iy, lx, ly;
+
+ x = creal(z);
+ y = cimag(z);
+
+ EXTRACT_WORDS(hx, lx, x);
+ EXTRACT_WORDS(hy, ly, y);
+
+ ix = 0x7fffffff & hx;
+ iy = 0x7fffffff & hy;
+
+ /* Handle the nearly-non-exceptional cases where x and y are finite. */
+ if (ix < 0x7ff00000 && iy < 0x7ff00000) {
+ if ((iy | ly) == 0)
+ return cpack(cosh(x), x * y);
+ if (ix < 0x40360000) /* small x: normal case */
+ return cpack(cosh(x) * cos(y), sinh(x) * sin(y));
+
+ /* |x| >= 22, so cosh(x) ~= exp(|x|) */
+ if (ix < 0x40862e42) {
+ /* x < 710: exp(|x|) won't overflow */
+ h = exp(fabs(x)) * 0.5;
+ return cpack(h * cos(y), copysign(h, x) * sin(y));
+ } else if (ix < 0x4096bbaa) {
+ /* x < 1455: scale to avoid overflow */
+ z = __ldexp_cexp(cpack(fabs(x), y), -1);
+ return cpack(creal(z), cimag(z) * copysign(1, x));
+ } else {
+ /* x >= 1455: the result always overflows */
+ h = huge * x;
+ return cpack(h * h * cos(y), h * sin(y));
+ }
+ }
+
+ /*
+ * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0.
+ * The sign of 0 in the result is unspecified. Choice = normally
+ * the same as dNaN. Raise the invalid floating-point exception.
+ *
+ * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0.
+ * The sign of 0 in the result is unspecified. Choice = normally
+ * the same as d(NaN).
+ */
+ if ((ix | lx) == 0 && iy >= 0x7ff00000)
+ return cpack(y - y, copysign(0, x * (y - y)));
+
+ /*
+ * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
+ *
+ * cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0.
+ * The sign of 0 in the result is unspecified.
+ */
+ if ((iy | ly) == 0 && ix >= 0x7ff00000) {
+ if (((hx & 0xfffff) | lx) == 0)
+ return cpack(x * x, copysign(0, x) * y);
+ return cpack(x * x, copysign(0, (x + x) * y));
+ }
+
+ /*
+ * cosh(x +- I Inf) = dNaN + I dNaN.
+ * Raise the invalid floating-point exception for finite nonzero x.
+ *
+ * cosh(x + I NaN) = d(NaN) + I d(NaN).
+ * Optionally raises the invalid floating-point exception for finite
+ * nonzero x. Choice = don't raise (except for signaling NaNs).
+ */
+ if (ix < 0x7ff00000 && iy >= 0x7ff00000)
+ return cpack(y - y, x * (y - y));
+
+ /*
+ * cosh(+-Inf + I NaN) = +Inf + I d(NaN).
+ *
+ * cosh(+-Inf +- I Inf) = +Inf + I dNaN.
+ * The sign of Inf in the result is unspecified. Choice = always +.
+ * Raise the invalid floating-point exception.
+ *
+ * cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y)
+ */
+ if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
+ if (iy >= 0x7ff00000)
+ return cpack(x * x, x * (y - y));
+ return cpack((x * x) * cos(y), x * sin(y));
+ }
+
+ /*
+ * cosh(NaN + I NaN) = d(NaN) + I d(NaN).
+ *
+ * cosh(NaN +- I Inf) = d(NaN) + I d(NaN).
+ * Optionally raises the invalid floating-point exception.
+ * Choice = raise.
+ *
+ * cosh(NaN + I y) = d(NaN) + I d(NaN).
+ * Optionally raises the invalid floating-point exception for finite
+ * nonzero y. Choice = don't raise (except for signaling NaNs).
+ */
+ return cpack((x * x) * (y - y), (x + x) * (y - y));
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_cexp.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 = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */
+cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
+
+double complex cexp(double complex z)
+{
+ double x, y, exp_x;
+ uint32_t hx, hy, lx, ly;
+
+ x = creal(z);
+ y = cimag(z);
+
+ EXTRACT_WORDS(hy, ly, y);
+ hy &= 0x7fffffff;
+
+ /* cexp(x + I 0) = exp(x) + I 0 */
+ if ((hy | ly) == 0)
+ return cpack(exp(x), y);
+ EXTRACT_WORDS(hx, lx, x);
+ /* cexp(0 + I y) = cos(y) + I sin(y) */
+ if (((hx & 0x7fffffff) | lx) == 0)
+ return cpack(cos(y), sin(y));
+
+ if (hy >= 0x7ff00000) {
+ if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
+ /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
+ return cpack(y - y, y - y);
+ } else if (hx & 0x80000000) {
+ /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
+ return cpack(0.0, 0.0);
+ } else {
+ /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
+ return cpack(x, y - y);
+ }
+ }
+
+ if (hx >= exp_ovfl && hx <= cexp_ovfl) {
+ /*
+ * x is between 709.7 and 1454.3, so we must scale to avoid
+ * overflow in exp(x).
+ */
+ return __ldexp_cexp(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 = exp(x);
+ return cpack(exp_x * cos(y), exp_x * sin(y));
+ }
+}
--- /dev/null
+#include "libm.h"
+
+double (cimag)(double complex z)
+{
+ union dcomplex u = {z};
+ return u.a[1];
+}
--- /dev/null
+#include "libm.h"
+
+double complex clog(double complex z)
+{
+ double r, phi;
+
+ r = cabs(z);
+ phi = carg(z);
+ return cpack(log(r), phi);
+}
--- /dev/null
+#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));
+}
+
+/*
+// 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;
+}
+*/
--- /dev/null
+/* 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"
+
+double complex cpow(double complex a, double complex z)
+{
+ double x, y, r, theta, absa, arga;
+
+ 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)));
+}
--- /dev/null
+#include "libm.h"
+
+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]));
+ return z;
+}
--- /dev/null
+#include <complex.h>
+
+double creal(double complex z)
+{
+ return z;
+}
--- /dev/null
+#include <complex.h>
+
+float crealf(float complex z)
+{
+ return z;
+}
--- /dev/null
+#include <complex.h>
+
+long double creall(long double complex z)
+{
+ return z;
+}
--- /dev/null
+#include "libm.h"
+
+double complex csin(double complex z)
+{
+ double complex r;
+
+ r = csinh(cpack(-cimag(z), creal(z)));
+ return cpack(cimag(r), -creal(r));
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_csinh.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 sine of a complex argument z = x + i y.
+ *
+ * sinh(z) = sinh(x+iy)
+ * = sinh(x) cos(y) + i cosh(x) sin(y).
+ *
+ * Exceptional values are noted in the comments within the source code.
+ * These values and the return value were taken from n1124.pdf.
+ */
+
+#include "libm.h"
+
+static const double huge = 0x1p1023;
+
+double complex csinh(double complex z)
+{
+ double x, y, h;
+ int32_t hx, hy, ix, iy, lx, ly;
+
+ x = creal(z);
+ y = cimag(z);
+
+ EXTRACT_WORDS(hx, lx, x);
+ EXTRACT_WORDS(hy, ly, y);
+
+ ix = 0x7fffffff & hx;
+ iy = 0x7fffffff & hy;
+
+ /* Handle the nearly-non-exceptional cases where x and y are finite. */
+ if (ix < 0x7ff00000 && iy < 0x7ff00000) {
+ if ((iy | ly) == 0)
+ return cpack(sinh(x), y);
+ if (ix < 0x40360000) /* small x: normal case */
+ return cpack(sinh(x) * cos(y), cosh(x) * sin(y));
+
+ /* |x| >= 22, so cosh(x) ~= exp(|x|) */
+ if (ix < 0x40862e42) {
+ /* x < 710: exp(|x|) won't overflow */
+ h = exp(fabs(x)) * 0.5;
+ return (cpack(copysign(h, x) * cos(y), h * sin(y)));
+ } else if (ix < 0x4096bbaa) {
+ /* x < 1455: scale to avoid overflow */
+ z = __ldexp_cexp(cpack(fabs(x), y), -1);
+ return cpack(creal(z) * copysign(1, x), cimag(z));
+ } else {
+ /* x >= 1455: the result always overflows */
+ h = huge * x;
+ return cpack(h * cos(y), h * h * sin(y));
+ }
+ }
+
+ /*
+ * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN.
+ * The sign of 0 in the result is unspecified. Choice = normally
+ * the same as dNaN. Raise the invalid floating-point exception.
+ *
+ * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN).
+ * The sign of 0 in the result is unspecified. Choice = normally
+ * the same as d(NaN).
+ */
+ if ((ix | lx) == 0 && iy >= 0x7ff00000)
+ return cpack(copysign(0, x * (y - y)), y - y);
+
+ /*
+ * sinh(+-Inf +- I 0) = +-Inf + I +-0.
+ *
+ * sinh(NaN +- I 0) = d(NaN) + I +-0.
+ */
+ if ((iy | ly) == 0 && ix >= 0x7ff00000) {
+ if (((hx & 0xfffff) | lx) == 0)
+ return (cpack(x, y));
+ return cpack(x, copysign(0, y));
+ }
+
+ /*
+ * sinh(x +- I Inf) = dNaN + I dNaN.
+ * Raise the invalid floating-point exception for finite nonzero x.
+ *
+ * sinh(x + I NaN) = d(NaN) + I d(NaN).
+ * Optionally raises the invalid floating-point exception for finite
+ * nonzero x. Choice = don't raise (except for signaling NaNs).
+ */
+ if (ix < 0x7ff00000 && iy >= 0x7ff00000)
+ return cpack(y - y, x * (y - y));
+
+ /*
+ * sinh(+-Inf + I NaN) = +-Inf + I d(NaN).
+ * The sign of Inf in the result is unspecified. Choice = normally
+ * the same as d(NaN).
+ *
+ * sinh(+-Inf +- I Inf) = +Inf + I dNaN.
+ * The sign of Inf in the result is unspecified. Choice = always +.
+ * Raise the invalid floating-point exception.
+ *
+ * sinh(+-Inf + I y) = +-Inf cos(y) + I Inf sin(y)
+ */
+ if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
+ if (iy >= 0x7ff00000)
+ return cpack(x * x, x * (y - y));
+ return cpack(x * cos(y), INFINITY * sin(y));
+ }
+
+ /*
+ * sinh(NaN + I NaN) = d(NaN) + I d(NaN).
+ *
+ * sinh(NaN +- I Inf) = d(NaN) + I d(NaN).
+ * Optionally raises the invalid floating-point exception.
+ * Choice = raise.
+ *
+ * sinh(NaN + I y) = d(NaN) + I d(NaN).
+ * Optionally raises the invalid floating-point exception for finite
+ * nonzero y. Choice = don't raise (except for signaling NaNs).
+ */
+ return cpack((x * x) * (y - y), (x + x) * (y - y));
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_csqrt.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
+
+/* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */
+#define THRESH 0x1.a827999fcef32p+1022
+
+double complex csqrt(double complex z)
+{
+ double complex result;
+ double a, b;
+ double t;
+ int scale;
+
+ a = creal(z);
+ b = cimag(z);
+
+ /* Handle special cases. */
+ if (z == 0)
+ return cpack(0, b);
+ if (isinf(b))
+ return cpack(INFINITY, b);
+ if (isnan(a)) {
+ t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
+ return cpack(a, t); /* return NaN + NaN i */
+ }
+ if (isinf(a)) {
+ /*
+ * csqrt(inf + NaN i) = inf + NaN i
+ * csqrt(inf + y i) = inf + 0 i
+ * csqrt(-inf + NaN i) = NaN +- inf i
+ * csqrt(-inf + y i) = 0 + inf i
+ */
+ if (signbit(a))
+ return cpack(fabs(b - b), copysign(a, b));
+ else
+ return cpack(a, copysign(b - b, b));
+ }
+ /*
+ * The remaining special case (b is NaN) is handled just fine by
+ * the normal code path below.
+ */
+
+ /* Scale to avoid overflow. */
+ if (fabs(a) >= THRESH || fabs(b) >= THRESH) {
+ a *= 0.25;
+ b *= 0.25;
+ scale = 1;
+ } else {
+ scale = 0;
+ }
+
+ /* Algorithm 312, CACM vol 10, Oct 1967. */
+ if (a >= 0) {
+ t = sqrt((a + hypot(a, b)) * 0.5);
+ result = cpack(t, b / (2 * t));
+ } else {
+ t = sqrt((-a + hypot(a, b)) * 0.5);
+ result = cpack(fabs(b) / (2 * t), copysign(t, b));
+ }
+
+ /* Rescale. */
+ if (scale)
+ result *= 2;
+ return result;
+}
--- /dev/null
+#include "libm.h"
+
+double complex ctan(double complex z)
+{
+ double complex r;
+
+ r = ctanh(cpack(-cimag(z), creal(z)));
+ return cpack(cimag(r), -creal(r));
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_ctanh.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 = x + i y.
+ *
+ * The algorithm is from:
+ *
+ * W. Kahan. Branch Cuts for Complex Elementary Functions or Much
+ * Ado About Nothing's Sign Bit. In The State of the Art in
+ * Numerical Analysis, pp. 165 ff. Iserles and Powell, eds., 1987.
+ *
+ * Method:
+ *
+ * Let t = tan(x)
+ * beta = 1/cos^2(y)
+ * s = sinh(x)
+ * rho = cosh(x)
+ *
+ * We have:
+ *
+ * tanh(z) = sinh(z) / cosh(z)
+ *
+ * sinh(x) cos(y) + i cosh(x) sin(y)
+ * = ---------------------------------
+ * cosh(x) cos(y) + i sinh(x) sin(y)
+ *
+ * cosh(x) sinh(x) / cos^2(y) + i tan(y)
+ * = -------------------------------------
+ * 1 + sinh^2(x) / cos^2(y)
+ *
+ * beta rho s + i t
+ * = ----------------
+ * 1 + beta s^2
+ *
+ * Modifications:
+ *
+ * I omitted the original algorithm's handling of overflow in tan(x) after
+ * verifying with nearpi.c that this can't happen in IEEE single or double
+ * precision. I also handle large x differently.
+ */
+
+#include "libm.h"
+
+double complex ctanh(double complex z)
+{
+ double x, y;
+ double t, beta, s, rho, denom;
+ uint32_t hx, ix, lx;
+
+ x = creal(z);
+ y = cimag(z);
+
+ EXTRACT_WORDS(hx, lx, x);
+ ix = hx & 0x7fffffff;
+
+ /*
+ * ctanh(NaN + i 0) = NaN + i 0
+ *
+ * ctanh(NaN + i y) = NaN + i NaN for y != 0
+ *
+ * The imaginary part has the sign of x*sin(2*y), but there's no
+ * special effort to get this right.
+ *
+ * ctanh(+-Inf +- i Inf) = +-1 +- 0
+ *
+ * ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite
+ *
+ * The imaginary part of the sign is unspecified. This special
+ * case is only needed to avoid a spurious invalid exception when
+ * y is infinite.
+ */
+ if (ix >= 0x7ff00000) {
+ if ((ix & 0xfffff) | lx) /* x is NaN */
+ return cpack(x, (y == 0 ? y : x * y));
+ SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */
+ return cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)));
+ }
+
+ /*
+ * ctanh(x + i NAN) = NaN + i NaN
+ * ctanh(x +- i Inf) = NaN + i NaN
+ */
+ if (!isfinite(y))
+ return cpack(y - y, y - y);
+
+ /*
+ * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the
+ * approximation sinh^2(huge) ~= exp(2*huge) / 4.
+ * We use a modified formula to avoid spurious overflow.
+ */
+ if (ix >= 0x40360000) { /* x >= 22 */
+ double exp_mx = exp(-fabs(x));
+ return cpack(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx);
+ }
+
+ /* Kahan's algorithm */
+ t = tan(y);
+ beta = 1.0 + t * t; /* = 1 / cos^2(y) */
+ s = sinh(x);
+ rho = sqrt(1 + s * s); /* = cosh(x) */
+ denom = 1 + beta * s * s;
+ return cpack((beta * rho * s) / denom, t / denom);
+}
} while (0)
/* fdlibm kernel functions */
+
int __rem_pio2_large(double*,double*,int,int,int);
int __rem_pio2(double,double*);
double __cos(double,double);
double __tan(double,double,int);
double __ldexp_exp(double,int);
-#if 0
double complex __ldexp_cexp(double complex,int);
-#endif
int __rem_pio2f(float,double*);
float __sindf(double);
float __cosdf(double);
float __tandf(double,int);
float __ldexp_expf(float,int);
-#if 0
float complex __ldexp_cexpf(float complex,int);
-#endif
-/* long double precision kernel functions */
long double __sinl(long double, long double, int);
long double __cosl(long double, long double);
long double __tanl(long double, long double, int);
long double __polevll(long double, long double *, int);
long double __p1evll(long double, long double *, int);
-// FIXME: nan
-/*
- * Common routine to process the arguments to nan(), nanf(), and nanl().
- */
-void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
-
// FIXME: not needed when -fexcess-precision=standard is supported (>=gcc4.5)
/*
* Attempt to get strict C99 semantics for assignment with non-C99 compilers.
return exp_x * scale;
}
-// FIXME
-#if 0
double complex __ldexp_cexp(double complex z, int expt)
{
double x, y, exp_x, scale1, scale2;
half_expt = expt - half_expt;
INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
- return (cpack(cos(y) * exp_x * scale1 * scale2,
- sin(y) * exp_x * scale1 * scale2));
+ return cpack(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2);
}
-#endif