From: nsz Date: Sun, 11 Mar 2012 14:24:11 +0000 (+0100) Subject: extend cmath (some of the functions are dummy) X-Git-Url: http://nsz.repo.hu/git/?p=libm;a=commitdiff_plain;h=5718e964d8a5c273a91b4d86d16926f54151c58f;ds=sidebyside extend cmath (some of the functions are dummy) --- diff --git a/src/cmath/cabsl.c b/src/cmath/cabsl.c index 31ac460..40a067c 100644 --- a/src/cmath/cabsl.c +++ b/src/cmath/cabsl.c @@ -1,6 +1,13 @@ #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)); } +#endif diff --git a/src/cmath/cacos.c b/src/cmath/cacos.c index b8a2d71..3aca051 100644 --- a/src/cmath/cacos.c +++ b/src/cmath/cacos.c @@ -1,52 +1,11 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/s_cacos.c */ -/* - * Copyright (c) 2008 Stephen L. Moshier - * - * 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" +// 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 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 index 0000000..563766e --- /dev/null +++ b/src/cmath/cacosf.c @@ -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)); +} diff --git a/src/cmath/cacosh.c b/src/cmath/cacosh.c index 89cbb1f..c2dfc1b 100644 --- a/src/cmath/cacosh.c +++ b/src/cmath/cacosh.c @@ -1,48 +1,9 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/s_cacosh.c */ -/* - * Copyright (c) 2008 Stephen L. Moshier - * - * 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" +/* acosh(z) = i acos(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 index 0000000..37ff880 --- /dev/null +++ b/src/cmath/cacoshf.c @@ -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 index 0000000..2a04e27 --- /dev/null +++ b/src/cmath/cacoshl.c @@ -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 index 0000000..5992e05 --- /dev/null +++ b/src/cmath/cacosl.c @@ -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 diff --git a/src/cmath/cargl.c b/src/cmath/cargl.c index 8aacdfa..e0d5047 100644 --- a/src/cmath/cargl.c +++ b/src/cmath/cargl.c @@ -1,6 +1,13 @@ #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)); } +#endif diff --git a/src/cmath/casin.c b/src/cmath/casin.c index 3428457..79aff27 100644 --- a/src/cmath/casin.c +++ b/src/cmath/casin.c @@ -1,120 +1,16 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/s_casin.c */ -/* - * Copyright (c) 2008 Stephen L. Moshier - * - * 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" +// FIXME + +/* asin(z) = -i log(i z + sqrt(1 - z*z)) */ + 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; + 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 index 0000000..cb9863f --- /dev/null +++ b/src/cmath/casinf.c @@ -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)); +} diff --git a/src/cmath/casinh.c b/src/cmath/casinh.c index 8c8d9fa..f2b3fef 100644 --- a/src/cmath/casinh.c +++ b/src/cmath/casinh.c @@ -1,47 +1,9 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/s_casinh.c */ -/* - * Copyright (c) 2008 Stephen L. Moshier - * - * 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" +/* asinh(z) = -i asin(i 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 index 0000000..ed4af64 --- /dev/null +++ b/src/cmath/casinhf.c @@ -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 index 0000000..e5d80ce --- /dev/null +++ b/src/cmath/casinhl.c @@ -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 index 0000000..f9aa8de --- /dev/null +++ b/src/cmath/casinl.c @@ -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 index 0000000..8533bde --- /dev/null +++ b/src/cmath/catanf.c @@ -0,0 +1,115 @@ +/* origin: OpenBSD /usr/src/lib/libm/src/s_catanf.c */ +/* + * Copyright (c) 2008 Stephen L. Moshier + * + * 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; +} diff --git a/src/cmath/catanh.c b/src/cmath/catanh.c index 5545747..b162802 100644 --- a/src/cmath/catanh.c +++ b/src/cmath/catanh.c @@ -1,47 +1,9 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/s_catanh.c */ -/* - * Copyright (c) 2008 Stephen L. Moshier - * - * 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" +/* atanh = -i atan(i 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 index 0000000..e1d1e64 --- /dev/null +++ b/src/cmath/catanhf.c @@ -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 index 0000000..0a9374a --- /dev/null +++ b/src/cmath/catanhl.c @@ -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 index 0000000..5ace770 --- /dev/null +++ b/src/cmath/catanl.c @@ -0,0 +1,126 @@ +/* origin: OpenBSD /usr/src/lib/libm/src/s_catanl.c */ +/* + * Copyright (c) 2008 Stephen L. Moshier + * + * 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 +#include +#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 diff --git a/src/cmath/ccos.c b/src/cmath/ccos.c index 9757e6a..5754c23 100644 --- a/src/cmath/ccos.c +++ b/src/cmath/ccos.c @@ -1,5 +1,7 @@ #include "libm.h" +/* cos(z) = cosh(i 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 index 0000000..9b72c4f --- /dev/null +++ b/src/cmath/ccosf.c @@ -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 index 0000000..683e77f --- /dev/null +++ b/src/cmath/ccoshf.c @@ -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 index 0000000..e37825a --- /dev/null +++ b/src/cmath/ccosl.c @@ -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 index 0000000..0cf13a3 --- /dev/null +++ b/src/cmath/cexpf.c @@ -0,0 +1,83 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/s_cexpf.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, 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 index 0000000..99fffc5 --- /dev/null +++ b/src/cmath/cimagf.c @@ -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 index 0000000..d9a0780 --- /dev/null +++ b/src/cmath/cimagl.c @@ -0,0 +1,7 @@ +#include "libm.h" + +long double (cimagl)(long double complex z) +{ + union lcomplex u = {z}; + return u.a[1]; +} diff --git a/src/cmath/clog.c b/src/cmath/clog.c index 4f8a428..6f10a11 100644 --- a/src/cmath/clog.c +++ b/src/cmath/clog.c @@ -1,5 +1,9 @@ #include "libm.h" +// FIXME + +/* log(z) = log(|z|) + i arg(z) */ + 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 index 0000000..f3aec54 --- /dev/null +++ b/src/cmath/clogf.c @@ -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 index 0000000..5b84ba5 --- /dev/null +++ b/src/cmath/clogl.c @@ -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 diff --git a/src/cmath/conj.c b/src/cmath/conj.c index 8720feb..4aceea7 100644 --- a/src/cmath/conj.c +++ b/src/cmath/conj.c @@ -1,23 +1,6 @@ #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; -} -*/ diff --git a/src/cmath/conjf.c b/src/cmath/conjf.c new file mode 100644 index 0000000..3155680 --- /dev/null +++ b/src/cmath/conjf.c @@ -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 index 0000000..0133226 --- /dev/null +++ b/src/cmath/conjl.c @@ -0,0 +1,6 @@ +#include "libm.h" + +long double complex conjl(long double complex z) +{ + return cpackl(creall(z), -cimagl(z)); +} diff --git a/src/cmath/cpow.c b/src/cmath/cpow.c index 0bc36d3..f863588 100644 --- a/src/cmath/cpow.c +++ b/src/cmath/cpow.c @@ -1,61 +1,8 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/s_cpow.c */ -/* - * Copyright (c) 2008 Stephen L. Moshier - * - * 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; +/* 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 index 0000000..53c65dc --- /dev/null +++ b/src/cmath/cpowf.c @@ -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 index 0000000..c1a80a7 --- /dev/null +++ b/src/cmath/cpowl.c @@ -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 diff --git a/src/cmath/cproj.c b/src/cmath/cproj.c index cc12e46..1cf9bb9 100644 --- a/src/cmath/cproj.c +++ b/src/cmath/cproj.c @@ -2,9 +2,7 @@ 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; } diff --git a/src/cmath/cprojf.c b/src/cmath/cprojf.c new file mode 100644 index 0000000..7112974 --- /dev/null +++ b/src/cmath/cprojf.c @@ -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 index 0000000..72e50cf --- /dev/null +++ b/src/cmath/cprojl.c @@ -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 diff --git a/src/cmath/csin.c b/src/cmath/csin.c index ab72055..246a459 100644 --- a/src/cmath/csin.c +++ b/src/cmath/csin.c @@ -1,9 +1,9 @@ #include "libm.h" +/* sin(z) = -i sinh(i 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 index 0000000..3aabe8f --- /dev/null +++ b/src/cmath/csinf.c @@ -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 index 0000000..1754186 --- /dev/null +++ b/src/cmath/csinhf.c @@ -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 +__FBSDID("$FreeBSD$"); + +#include +#include + +#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 index 0000000..4ad8674 --- /dev/null +++ b/src/cmath/csinl.c @@ -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 index 0000000..16487c2 --- /dev/null +++ b/src/cmath/csqrtf.c @@ -0,0 +1,82 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/s_csqrtf.c */ +/*- + * Copyright (c) 2007 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, 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)); + } +} diff --git a/src/cmath/ctan.c b/src/cmath/ctan.c index 4462ef0..4741a4d 100644 --- a/src/cmath/ctan.c +++ b/src/cmath/ctan.c @@ -1,9 +1,9 @@ #include "libm.h" +/* tan(z) = -i tanh(i 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 index 0000000..9bbeb05 --- /dev/null +++ b/src/cmath/ctanf.c @@ -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 index 0000000..7d74613 --- /dev/null +++ b/src/cmath/ctanhf.c @@ -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 index 0000000..4b4c99b --- /dev/null +++ b/src/cmath/ctanl.c @@ -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