diff options
| author |   martynas <martynas@openbsd.org> | 2011-07-08 19:25:31 +0000 |
|---|---|---|
| committer | martynas <martynas@openbsd.org> | 2011-07-08 19:25:31 +0000 |
| commit | de3697aa9479c3be5ff81d594e19ca819a9fb9ec (patch) | |
| tree | 836d295c7bed312c6f332a7545bcc8399e6cc16e /lib/libm/src | |
| parent | COMPAT_43 does not do all it used to. COMPAT_AOUT is wholly useless. (diff) | |
| download | wireguard-openbsd-de3697aa9479c3be5ff81d594e19ca819a9fb9ec.tar.xz wireguard-openbsd-de3697aa9479c3be5ff81d594e19ca819a9fb9ec.zip | |
Finalize work on complex math routines, now that we have the
extended-precision support. Mostly from Cephes.
Diffstat (limited to 'lib/libm/src')
44 files changed, 1805 insertions, 23 deletions
diff --git a/lib/libm/src/s_cabs.c b/lib/libm/src/s_cabs.c index d2f3f942c3f..3d28fd4de64 100644 --- a/lib/libm/src/s_cabs.c +++ b/lib/libm/src/s_cabs.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_cabs.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_cabs.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org> * @@ -15,7 +15,11 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double @@ -23,3 +27,12 @@ cabs(double complex z) { return hypot(__real__ z, __imag__ z); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double cabsl(long double complex); +#else /* lint */ +__weak_alias(cabsl, cabs); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_cabsl.c b/lib/libm/src/s_cabsl.c new file mode 100644 index 00000000000..2e68ec0e693 --- /dev/null +++ b/lib/libm/src/s_cabsl.c @@ -0,0 +1,26 @@ +/* $OpenBSD: s_cabsl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2011 Martynas Venckus <martynas@openbsd.org> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +#include <complex.h> +#include <math.h> + +long double +cabsl(long double complex z) +{ + return hypotl(__real__ z, __imag__ z); +} diff --git a/lib/libm/src/s_cacos.c b/lib/libm/src/s_cacos.c index 96d26626b98..1ff40d4f858 100644 --- a/lib/libm/src/s_cacos.c +++ b/lib/libm/src/s_cacos.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_cacos.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_cacos.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* cacos() * * Complex circular arc cosine @@ -46,7 +48,9 @@ * IEEE -10,+10 30000 1.8e-14 2.2e-15 */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -58,3 +62,12 @@ cacos(double complex z) w = (M_PI_2 - creal (w)) - cimag (w) * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex cacosl(long double complex); +#else /* lint */ +__weak_alias(cacosl, cacos); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_cacosh.c b/lib/libm/src/s_cacosh.c index 43ba99ad6f9..8db79ddc0b2 100644 --- a/lib/libm/src/s_cacosh.c +++ b/lib/libm/src/s_cacosh.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_cacosh.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_cacosh.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* cacosh * * Complex inverse hyperbolic cosine @@ -42,7 +44,9 @@ * */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -53,3 +57,12 @@ cacosh(double complex z) w = I * cacos (z); return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex cacoshl(long double complex); +#else /* lint */ +__weak_alias(cacoshl, cacosh); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_cacoshl.c b/lib/libm/src/s_cacoshl.c new file mode 100644 index 00000000000..67a2362c273 --- /dev/null +++ b/lib/libm/src/s_cacoshl.c @@ -0,0 +1,56 @@ +/* $OpenBSD: s_cacoshl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* cacoshl + * + * Complex inverse hyperbolic cosine + * + * + * + * SYNOPSIS: + * + * long double complex cacoshl(); + * long double complex z, w; + * + * w = cacoshl (z); + * + * + * + * DESCRIPTION: + * + * acosh z = i acos z . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.6e-14 2.1e-15 + * + */ + +#include <complex.h> +#include <math.h> + +long double complex +cacoshl(long double complex z) +{ + long double complex w; + + w = I * cacosl(z); + return (w); +} diff --git a/lib/libm/src/s_cacosl.c b/lib/libm/src/s_cacosl.c new file mode 100644 index 00000000000..54fd29a5664 --- /dev/null +++ b/lib/libm/src/s_cacosl.c @@ -0,0 +1,63 @@ +/* $OpenBSD: s_cacosl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* cacosl() + * + * Complex circular arc cosine + * + * + * + * SYNOPSIS: + * + * long double complex cacosl(); + * long double complex z, w; + * + * w = cacosl( z ); + * + * + * + * DESCRIPTION: + * + * + * w = arccos z = PI/2 - arcsin z. + * + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 5200 1.6e-15 2.8e-16 + * IEEE -10,+10 30000 1.8e-14 2.2e-15 + */ + +#include <complex.h> +#include <math.h> + +static long double PIO2L = 1.570796326794896619231321691639751442098585L; + +long double complex +cacosl(long double complex z) +{ + long double complex w; + + w = casinl(z); + w = (PIO2L - creal(w)) - cimag(w) * I; + return (w); +} diff --git a/lib/libm/src/s_carg.c b/lib/libm/src/s_carg.c index db51662a1d9..08f321c36ab 100644 --- a/lib/libm/src/s_carg.c +++ b/lib/libm/src/s_carg.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_carg.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_carg.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org> * @@ -15,7 +15,11 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double @@ -23,3 +27,12 @@ carg(double complex z) { return atan2 (__imag__ z, __real__ z); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double cargl(long double complex); +#else /* lint */ +__weak_alias(cargl, carg); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_cargl.c b/lib/libm/src/s_cargl.c new file mode 100644 index 00000000000..ba2fb5cdabe --- /dev/null +++ b/lib/libm/src/s_cargl.c @@ -0,0 +1,26 @@ +/* $OpenBSD: s_cargl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2011 Martynas Venckus <martynas@openbsd.org> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +#include <complex.h> +#include <math.h> + +long double +cargl(long double complex z) +{ + return atan2l(__imag__ z, __real__ z); +} diff --git a/lib/libm/src/s_casin.c b/lib/libm/src/s_casin.c index 1c1b53541bb..64ebe6bb7a8 100644 --- a/lib/libm/src/s_casin.c +++ b/lib/libm/src/s_casin.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_casin.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_casin.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* casin() * * Complex circular arc sine @@ -49,7 +51,9 @@ * Also tested by csin(casin(z)) = z. */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -127,3 +131,12 @@ casin(double complex z) w = zz * (-1.0 * I); return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex casinl(long double complex); +#else /* lint */ +__weak_alias(casinl, casin); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_casinh.c b/lib/libm/src/s_casinh.c index a249cf15ede..8f1f02fa54a 100644 --- a/lib/libm/src/s_casinh.c +++ b/lib/libm/src/s_casinh.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_casinh.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_casinh.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* casinh * * Complex inverse hyperbolic sine @@ -42,7 +44,9 @@ * */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -53,3 +57,12 @@ casinh(double complex z) w = -1.0 * I * casin (z * I); return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex casinhl(long double complex); +#else /* lint */ +__weak_alias(casinhl, casinh); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_casinhl.c b/lib/libm/src/s_casinhl.c new file mode 100644 index 00000000000..7a8a31ce27c --- /dev/null +++ b/lib/libm/src/s_casinhl.c @@ -0,0 +1,56 @@ +/* $OpenBSD: s_casinhl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* casinhl + * + * Complex inverse hyperbolic sine + * + * + * + * SYNOPSIS: + * + * long double complex casinhf(); + * long double complex z, w; + * + * w = casinhl (z); + * + * + * + * DESCRIPTION: + * + * casinh z = -i casin iz . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.8e-14 2.6e-15 + * + */ + +#include <complex.h> +#include <math.h> + +long double complex +casinhl(long double complex z) +{ + long double complex w; + + w = -1.0L * I * casinl(z * I); + return (w); +} diff --git a/lib/libm/src/s_casinl.c b/lib/libm/src/s_casinl.c new file mode 100644 index 00000000000..8b50af7ef2e --- /dev/null +++ b/lib/libm/src/s_casinl.c @@ -0,0 +1,130 @@ +/* $OpenBSD: s_casinl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* casinl() + * + * Complex circular arc sine + * + * + * + * SYNOPSIS: + * + * long double complex casinl(); + * long double complex z, w; + * + * w = casinl( z ); + * + * + * + * DESCRIPTION: + * + * Inverse complex sine: + * + * 2 + * w = -i clog( iz + csqrt( 1 - z ) ). + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 10100 2.1e-15 3.4e-16 + * IEEE -10,+10 30000 2.2e-14 2.7e-15 + * Larger relative error can be observed for z near zero. + * Also tested by csin(casin(z)) = z. + */ + +#include <complex.h> +#include <float.h> +#include <math.h> + +#if LDBL_MANT_DIG == 64 +static long double MACHEPL= 5.42101086242752217003726400434970855712890625E-20L; +#elif LDBL_MANT_DIG == 113 +static long double MACHEPL = 9.629649721936179265279889712924636592690508e-35L; +#endif + +static long double PIO2L = 1.570796326794896619231321691639751442098585L; + +long double complex +casinl(long double complex z) +{ + long double complex w; + long double x, y, b; + static long double complex ca, ct, zz, z2; + + x = creal(z); + y = cimag(z); + + if (y == 0.0L) { + if (fabsl(x) > 1.0L) { + w = PIO2L + 0.0L * I; + /*mtherr( "casinl", DOMAIN );*/ + } + else { + w = asinl(x) + 0.0L * I; + } + return (w); + } + + /* Power series expansion */ + b = cabsl(z); + if (b < 0.125L) { + long double complex sum; + long double n, cn; + + z2 = (x - y) * (x + y) + (2.0L * x * y) * I; + cn = 1.0L; + n = 1.0L; + ca = x + y * I; + sum = x + y * I; + do { + ct = z2 * ca; + ca = ct; + + cn *= n; + n += 1.0L; + cn /= n; + n += 1.0L; + b = cn/n; + + ct *= b; + sum += ct; + b = cabsl(ct); + } + + while (b > MACHEPL); + w = sum; + return w; + } + + ca = x + y * I; + ct = ca * I; /* iz */ + /* sqrt(1 - z*z) */ + /* cmul(&ca, &ca, &zz) */ + /* x * x - y * y */ + zz = (x - y) * (x + y) + (2.0L * x * y) * I; + zz = 1.0L - creal(zz) - cimag(zz) * I; + z2 = csqrtl(zz); + + zz = ct + z2; + zz = clogl(zz); + /* multiply by 1/i = -i */ + w = zz * (-1.0L * I); + return (w); +} diff --git a/lib/libm/src/s_catan.c b/lib/libm/src/s_catan.c index ddd7940456e..78861b2f929 100644 --- a/lib/libm/src/s_catan.c +++ b/lib/libm/src/s_catan.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_catan.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_catan.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* catan() * * Complex circular arc tangent @@ -62,7 +64,9 @@ * 2.9e-17. See also clog(). */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> #define MAXNUM 1.0e308 @@ -124,3 +128,12 @@ ovrf: w = MAXNUM + MAXNUM * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex catanl(long double complex); +#else /* lint */ +__weak_alias(catanl, catan); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_catanh.c b/lib/libm/src/s_catanh.c index 7a6d6f99646..7b23934b0fe 100644 --- a/lib/libm/src/s_catanh.c +++ b/lib/libm/src/s_catanh.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_catanh.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_catanh.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* catanh * * Complex inverse hyperbolic tangent @@ -42,7 +44,9 @@ * */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -53,3 +57,12 @@ catanh(double complex z) w = -1.0 * I * catan (z * I); return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex catanhl(long double complex); +#else /* lint */ +__weak_alias(catanhl, catanh); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_catanhl.c b/lib/libm/src/s_catanhl.c new file mode 100644 index 00000000000..e2134f93261 --- /dev/null +++ b/lib/libm/src/s_catanhl.c @@ -0,0 +1,56 @@ +/* $OpenBSD: s_catanhl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* catanhl + * + * Complex inverse hyperbolic tangent + * + * + * + * SYNOPSIS: + * + * long double complex catanhl(); + * long double complex z, w; + * + * w = catanhl (z); + * + * + * + * DESCRIPTION: + * + * Inverse tanh, equal to -i catan (iz); + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 2.3e-16 6.2e-17 + * + */ + +#include <complex.h> +#include <math.h> + +long double complex +catanhl(long double complex z) +{ + long double complex w; + + w = -1.0L * I * catanl(z * I); + return (w); +} diff --git a/lib/libm/src/s_catanl.c b/lib/libm/src/s_catanl.c new file mode 100644 index 00000000000..9d3a4985d98 --- /dev/null +++ b/lib/libm/src/s_catanl.c @@ -0,0 +1,127 @@ +/* $OpenBSD: s_catanl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* catanl() + * + * Complex circular arc tangent + * + * + * + * SYNOPSIS: + * + * long double complex catanl(); + * long double complex z, w; + * + * w = catanl( z ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * 1 ( 2x ) + * Re w = - arctan(-----------) + k PI + * 2 ( 2 2) + * (1 - x - y ) + * + * ( 2 2) + * 1 (x + (y+1) ) + * Im w = - log(------------) + * 4 ( 2 2) + * (x + (y-1) ) + * + * Where k is an arbitrary integer. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 5900 1.3e-16 7.8e-18 + * IEEE -10,+10 30000 2.3e-15 8.5e-17 + * The check catan( ctan(z) ) = z, with |x| and |y| < PI/2, + * had peak relative error 1.5e-16, rms relative error + * 2.9e-17. See also clog(). + */ + +#include <complex.h> +#include <float.h> +#include <math.h> + +static long double PIL = 3.141592653589793238462643383279502884197169L; +static long double DP1 = 3.14159265358979323829596852490908531763125L; +static long double DP2 = 1.6667485837041756656403424829301998703007e-19L; +static 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 = creal(z); + y = cimag(z); + + if ((x == 0.0L) && (y > 1.0L)) + goto ovrf; + + x2 = x * x; + a = 1.0L - x2 - (y * y); + if (a == 0.0L) + goto ovrf; + + t = atan2l(2.0L * x, a) * 0.5L; + w = redupil(t); + + t = y - 1.0L; + a = x2 + (t * t); + if (a == 0.0L) + goto ovrf; + + t = y + 1.0L; + a = (x2 + (t * t)) / a; + w = w + (0.25L * logl(a)) * I; + return (w); + +ovrf: + /*mtherr( "catanl", OVERFLOW );*/ + w = LDBL_MAX + LDBL_MAX * I; + return (w); +} diff --git a/lib/libm/src/s_ccos.c b/lib/libm/src/s_ccos.c index 264143ab24c..7152bc14253 100644 --- a/lib/libm/src/s_ccos.c +++ b/lib/libm/src/s_ccos.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_ccos.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_ccos.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* ccos() * * Complex circular cosine @@ -49,7 +51,9 @@ * IEEE -10,+10 30000 3.8e-16 1.0e-16 */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> /* calculate cosh and sinh */ @@ -82,3 +86,12 @@ ccos(double complex z) w = cos(creal (z)) * ch - (sin (creal (z)) * sh) * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex ccosl(long double complex); +#else /* lint */ +__weak_alias(ccosl, ccos); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_ccosh.c b/lib/libm/src/s_ccosh.c index 0267b6f2a15..3e780e976f0 100644 --- a/lib/libm/src/s_ccosh.c +++ b/lib/libm/src/s_ccosh.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_ccosh.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_ccosh.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* ccosh * * Complex hyperbolic cosine @@ -42,7 +44,9 @@ * */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -56,3 +60,12 @@ ccosh(double complex z) w = cosh (x) * cos (y) + (sinh (x) * sin (y)) * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex ccoshl(long double complex); +#else /* lint */ +__weak_alias(ccoshl, ccosh); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_ccoshl.c b/lib/libm/src/s_ccoshl.c new file mode 100644 index 00000000000..58acdbea43b --- /dev/null +++ b/lib/libm/src/s_ccoshl.c @@ -0,0 +1,59 @@ +/* $OpenBSD: s_ccoshl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* ccoshl + * + * Complex hyperbolic cosine + * + * + * + * SYNOPSIS: + * + * long double complex ccoshl(); + * long double complex z, w; + * + * w = ccoshl (z); + * + * + * + * DESCRIPTION: + * + * ccosh(z) = cosh x cos y + i sinh x sin y . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 2.9e-16 8.1e-17 + * + */ + +#include <complex.h> +#include <math.h> + +long double complex +ccoshl(long double complex z) +{ + long double complex w; + long double x, y; + + x = creal(z); + y = cimag(z); + w = coshl(x) * cosl(y) + (sinhl(x) * sinl(y)) * I; + return (w); +} diff --git a/lib/libm/src/s_ccosl.c b/lib/libm/src/s_ccosl.c new file mode 100644 index 00000000000..02d39a2ea65 --- /dev/null +++ b/lib/libm/src/s_ccosl.c @@ -0,0 +1,82 @@ +/* $OpenBSD: s_ccosl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* ccosl() + * + * Complex circular cosine + * + * + * + * SYNOPSIS: + * + * long double complex ccosl(); + * long double complex z, w; + * + * w = ccosl( z ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * w = cos x cosh y - i sin x sinh y. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 8400 4.5e-17 1.3e-17 + * IEEE -10,+10 30000 3.8e-16 1.0e-16 + */ + +#include <complex.h> +#include <math.h> + +static void +cchshl(long double x, long double *c, long double *s) +{ + long double e, ei; + + if(fabsl(x) <= 0.5L) { + *c = coshl(x); + *s = sinhl(x); + } else { + e = expl(x); + ei = 0.5L/e; + e = 0.5L * e; + *s = e - ei; + *c = e + ei; + } +} + +long double complex +ccosl(long double complex z) +{ + long double complex w; + long double ch, sh; + + cchshl(cimag(z), &ch, &sh); + w = cosl(creal(z)) * ch + (-sinl(creal(z)) * sh) * I; + return (w); +} diff --git a/lib/libm/src/s_cexp.c b/lib/libm/src/s_cexp.c index 52241a8eb5e..14a43270f52 100644 --- a/lib/libm/src/s_cexp.c +++ b/lib/libm/src/s_cexp.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_cexp.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_cexp.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* cexp() * * Complex exponential function @@ -53,7 +55,9 @@ * */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -68,3 +72,12 @@ cexp(double complex z) w = r * cos (y) + r * sin (y) * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex cexpl(long double complex); +#else /* lint */ +__weak_alias(cexpl, cexp); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_cexpl.c b/lib/libm/src/s_cexpl.c new file mode 100644 index 00000000000..6fbd347e77d --- /dev/null +++ b/lib/libm/src/s_cexpl.c @@ -0,0 +1,70 @@ +/* $OpenBSD: s_cexpl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* cexpl() + * + * Complex exponential function + * + * + * + * SYNOPSIS: + * + * long double complex cexpl(); + * long double complex z, w; + * + * w = cexpl( z ); + * + * + * + * DESCRIPTION: + * + * Returns the exponential of the complex argument z + * into the complex result w. + * + * If + * z = x + iy, + * r = exp(x), + * + * then + * + * w = r cos y + i r sin y. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 8700 3.7e-17 1.1e-17 + * IEEE -10,+10 30000 3.0e-16 8.7e-17 + * + */ + +#include <complex.h> +#include <math.h> + +long double complex +cexpl(long double complex z) +{ + long double complex w; + long double r; + + r = expl(creal(z)); + w = r * cosl((long double)cimag(z)) + + (r * sinl((long double)cimag(z))) * I; + return (w); +} diff --git a/lib/libm/src/s_cimag.c b/lib/libm/src/s_cimag.c index 85aaa2975d1..d03a6186e76 100644 --- a/lib/libm/src/s_cimag.c +++ b/lib/libm/src/s_cimag.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_cimag.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_cimag.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org> * @@ -15,7 +15,11 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double @@ -23,3 +27,12 @@ cimag(double complex z) { return __imag__ z; } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double cimagl(long double complex); +#else /* lint */ +__weak_alias(cimagl, cimag); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_cimagl.c b/lib/libm/src/s_cimagl.c new file mode 100644 index 00000000000..b25c3402027 --- /dev/null +++ b/lib/libm/src/s_cimagl.c @@ -0,0 +1,26 @@ +/* $OpenBSD: s_cimagl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2011 Martynas Venckus <martynas@openbsd.org> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +#include <complex.h> +#include <math.h> + +long double +cimagl(long double complex z) +{ + return __imag__ z; +} diff --git a/lib/libm/src/s_clog.c b/lib/libm/src/s_clog.c index 31163700607..ce806269bfb 100644 --- a/lib/libm/src/s_clog.c +++ b/lib/libm/src/s_clog.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_clog.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_clog.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* clog.c * * Complex natural logarithm @@ -54,7 +56,9 @@ * absolute error 1.0e-16. */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -70,3 +74,12 @@ clog(double complex z) w = p + rr * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex clogl(long double complex); +#else /* lint */ +__weak_alias(clogl, clog); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_clogl.c b/lib/libm/src/s_clogl.c new file mode 100644 index 00000000000..2a70dbb42d7 --- /dev/null +++ b/lib/libm/src/s_clogl.c @@ -0,0 +1,73 @@ +/* $OpenBSD: s_clogl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* clogl.c + * + * Complex natural logarithm + * + * + * + * SYNOPSIS: + * + * long double complex clogl(); + * long double complex z, w; + * + * w = clogl( z ); + * + * + * + * DESCRIPTION: + * + * Returns complex logarithm to the base e (2.718...) of + * the complex argument x. + * + * If z = x + iy, r = sqrt( x**2 + y**2 ), + * then + * w = log(r) + i arctan(y/x). + * + * The arctangent ranges from -PI to +PI. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 7000 8.5e-17 1.9e-17 + * IEEE -10,+10 30000 5.0e-15 1.1e-16 + * + * Larger relative error can be observed for z near 1 +i0. + * In IEEE arithmetic the peak absolute error is 5.2e-16, rms + * absolute error 1.0e-16. + */ + +#include <complex.h> +#include <math.h> + +long double complex +clogl(long double complex z) +{ + long double complex w; + long double p, rr; + + /*rr = sqrt(z->r * z->r + z->i * z->i);*/ + p = cabsl(z); + p = logl(p); + rr = atan2l(cimag(z), creal(z)); + w = p + rr * I; + return (w); +} diff --git a/lib/libm/src/s_conj.c b/lib/libm/src/s_conj.c index fd38a9d034c..0ffb990a831 100644 --- a/lib/libm/src/s_conj.c +++ b/lib/libm/src/s_conj.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_conj.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_conj.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org> * @@ -15,7 +15,11 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -23,3 +27,12 @@ conj(double complex z) { return ~z; } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex conjl(long double complex); +#else /* lint */ +__weak_alias(conjl, conj); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_conjl.c b/lib/libm/src/s_conjl.c new file mode 100644 index 00000000000..8035895ef71 --- /dev/null +++ b/lib/libm/src/s_conjl.c @@ -0,0 +1,26 @@ +/* $OpenBSD: s_conjl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2011 Martynas Venckus <martynas@openbsd.org> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +#include <complex.h> +#include <math.h> + +long double complex +conjl(long double complex z) +{ + return ~z; +} diff --git a/lib/libm/src/s_cpow.c b/lib/libm/src/s_cpow.c index 6b0537393ce..a0327118e35 100644 --- a/lib/libm/src/s_cpow.c +++ b/lib/libm/src/s_cpow.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_cpow.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_cpow.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* cpow * * Complex power function @@ -44,7 +46,9 @@ * */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -57,7 +61,7 @@ cpow(double complex a, double complex z) y = cimag (z); absa = cabs (a); if (absa == 0.0) { - return (0.0 + 0.0 * I); + return (0.0 + 0.0 * I); } arga = carg (a); r = pow (absa, x); @@ -69,3 +73,12 @@ cpow(double complex a, double complex z) w = r * cos (theta) + (r * sin (theta)) * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex cpowl(long double complex, long double complex); +#else /* lint */ +__weak_alias(cpowl, cpow); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_cpowl.c b/lib/libm/src/s_cpowl.c new file mode 100644 index 00000000000..7707827fe6e --- /dev/null +++ b/lib/libm/src/s_cpowl.c @@ -0,0 +1,72 @@ +/* $OpenBSD: s_cpowl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* cpowl + * + * Complex power function + * + * + * + * SYNOPSIS: + * + * long double complex cpowl(); + * long double complex a, z, w; + * + * w = cpowl (a, z); + * + * + * + * DESCRIPTION: + * + * Raises complex A to the complex Zth power. + * Definition is per AMS55 # 4.2.8, + * analytically equivalent to cpow(a,z) = cexp(z clog(a)). + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 9.4e-15 1.5e-15 + * + */ + +#include <complex.h> +#include <math.h> + +long double complex +cpowl(long double complex a, long double complex z) +{ + long double complex w; + long double x, y, r, theta, absa, arga; + + x = creal(z); + y = cimag(z); + absa = cabsl(a); + if (absa == 0.0L) { + return (0.0L + 0.0L * I); + } + arga = cargl(a); + r = powl(absa, x); + theta = x * arga; + if (y != 0.0L) { + r = r * expl(-y * arga); + theta = theta + y * logl(absa); + } + w = r * cosl(theta) + (r * sinl(theta)) * I; + return (w); +} diff --git a/lib/libm/src/s_cproj.c b/lib/libm/src/s_cproj.c index e4fc50f8b57..543290a08c2 100644 --- a/lib/libm/src/s_cproj.c +++ b/lib/libm/src/s_cproj.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_cproj.c,v 1.2 2010/07/19 00:04:07 guenther Exp $ */ +/* $OpenBSD: s_cproj.c,v 1.3 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org> * @@ -15,7 +15,11 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -32,3 +36,12 @@ cproj(double complex z) return res; } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex cprojl(long double complex); +#else /* lint */ +__weak_alias(cprojl, cproj); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_cprojl.c b/lib/libm/src/s_cprojl.c new file mode 100644 index 00000000000..596e417a429 --- /dev/null +++ b/lib/libm/src/s_cprojl.c @@ -0,0 +1,35 @@ +/* $OpenBSD: s_cprojl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2011 Martynas Venckus <martynas@openbsd.org> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +#include <complex.h> +#include <math.h> + +long double complex +cprojl(long double complex z) +{ + long double complex res; + + if (isinf(__real__ z) || isinf(__imag__ z)) { + __real__ res = INFINITY; + __imag__ res = copysignl(0.0, __imag__ z); + } else { + res = z; + } + + return res; +} diff --git a/lib/libm/src/s_creal.c b/lib/libm/src/s_creal.c index 55c36e8d839..a45941534cc 100644 --- a/lib/libm/src/s_creal.c +++ b/lib/libm/src/s_creal.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_creal.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_creal.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org> * @@ -15,7 +15,11 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double @@ -23,3 +27,12 @@ creal(double complex z) { return __real__ z; } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double creall(long double complex); +#else /* lint */ +__weak_alias(creall, creal); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_creall.c b/lib/libm/src/s_creall.c new file mode 100644 index 00000000000..4f3f4a4194e --- /dev/null +++ b/lib/libm/src/s_creall.c @@ -0,0 +1,26 @@ +/* $OpenBSD: s_creall.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2011 Martynas Venckus <martynas@openbsd.org> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +#include <complex.h> +#include <math.h> + +long double +creall(long double complex z) +{ + return __real__ z; +} diff --git a/lib/libm/src/s_csin.c b/lib/libm/src/s_csin.c index ad249de7ac0..44f7d925148 100644 --- a/lib/libm/src/s_csin.c +++ b/lib/libm/src/s_csin.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_csin.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_csin.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* csin() * * Complex circular sine @@ -51,7 +53,9 @@ * */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> /* calculate cosh and sinh */ @@ -84,3 +88,12 @@ csin(double complex z) w = sin (creal(z)) * ch + (cos (creal(z)) * sh) * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex csinl(long double complex); +#else /* lint */ +__weak_alias(csinl, csin); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_csinh.c b/lib/libm/src/s_csinh.c index ece633cda8f..83ee14be7b8 100644 --- a/lib/libm/src/s_csinh.c +++ b/lib/libm/src/s_csinh.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_csinh.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_csinh.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* csinh * * Complex hyperbolic sine @@ -41,7 +43,9 @@ * */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -55,3 +59,12 @@ csinh(double complex z) w = sinh (x) * cos (y) + (cosh (x) * sin (y)) * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex csinhl(long double complex); +#else /* lint */ +__weak_alias(csinhl, csinh); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_csinhl.c b/lib/libm/src/s_csinhl.c new file mode 100644 index 00000000000..de5e7938d0b --- /dev/null +++ b/lib/libm/src/s_csinhl.c @@ -0,0 +1,58 @@ +/* $OpenBSD: s_csinhl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* csinhl + * + * Complex hyperbolic sine + * + * + * + * SYNOPSIS: + * + * long double complex csinhl(); + * long double complex z, w; + * + * w = csinhl (z); + * + * DESCRIPTION: + * + * csinh z = (cexp(z) - cexp(-z))/2 + * = sinh x * cos y + i cosh x * sin y . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 3.1e-16 8.2e-17 + * + */ + +#include <complex.h> +#include <math.h> + +long double complex +csinhl(long double complex z) +{ + long double complex w; + long double x, y; + + x = creal(z); + y = cimag(z); + w = sinhl(x) * cosl(y) + (coshl(x) * sinl(y)) * I; + return (w); +} diff --git a/lib/libm/src/s_csinl.c b/lib/libm/src/s_csinl.c new file mode 100644 index 00000000000..f9e20134bee --- /dev/null +++ b/lib/libm/src/s_csinl.c @@ -0,0 +1,84 @@ +/* $OpenBSD: s_csinl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* csinl() + * + * Complex circular sine + * + * + * + * SYNOPSIS: + * + * long double complex csinl(); + * long double complex z, w; + * + * w = csinl( z ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * w = sin x cosh y + i cos x sinh y. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 8400 5.3e-17 1.3e-17 + * IEEE -10,+10 30000 3.8e-16 1.0e-16 + * Also tested by csin(casin(z)) = z. + * + */ + +#include <complex.h> +#include <math.h> + +static void +cchshl(long double x, long double *c, long double *s) +{ + long double e, ei; + + if(fabsl(x) <= 0.5L) { + *c = coshl(x); + *s = sinhl(x); + } else { + e = expl(x); + ei = 0.5L/e; + e = 0.5L * e; + *s = e - ei; + *c = e + ei; + } +} + +long double complex +csinl(long double complex z) +{ + long double complex w; + long double ch, sh; + + cchshl(cimag(z), &ch, &sh); + w = sinl(creal(z)) * ch + (cosl(creal(z)) * sh) * I; + return (w); +} diff --git a/lib/libm/src/s_csqrt.c b/lib/libm/src/s_csqrt.c index 3ed4aa56076..700c17dd533 100644 --- a/lib/libm/src/s_csqrt.c +++ b/lib/libm/src/s_csqrt.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_csqrt.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_csqrt.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* csqrt() * * Complex square root @@ -58,7 +60,9 @@ * */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -129,3 +133,12 @@ csqrt(double complex z) w = t + r * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex csqrtl(long double complex); +#else /* lint */ +__weak_alias(csqrtl, csqrt); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_csqrtl.c b/lib/libm/src/s_csqrtl.c new file mode 100644 index 00000000000..86c4f380590 --- /dev/null +++ b/lib/libm/src/s_csqrtl.c @@ -0,0 +1,128 @@ +/* $OpenBSD: s_csqrtl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* csqrtl() + * + * Complex square root + * + * + * + * SYNOPSIS: + * + * long double complex csqrtl(); + * long double complex z, w; + * + * w = csqrtl( z ); + * + * + * + * DESCRIPTION: + * + * + * If z = x + iy, r = |z|, then + * + * 1/2 + * Re w = [ (r + x)/2 ] , + * + * 1/2 + * Im w = [ (r - x)/2 ] . + * + * Cancellation error in r-x or r+x is avoided by using the + * identity 2 Re w Im w = y. + * + * Note that -w is also a square root of z. The root chosen + * is always in the right half plane and Im w has the same sign as y. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 500000 1.1e-19 3.0e-20 + * + */ + +#include <complex.h> +#include <math.h> + +long double complex +csqrtl(long double complex z) +{ + long double complex w; + long double x, y, r, t, scale; + + x = creal(z); + y = cimag(z); + + if (y == 0.0L) { + if (x < 0.0L) { + w = 0.0L + sqrtl(-x) * I; + return (w); + } + else { + w = sqrtl(x) + 0.0L * I; + return (w); + } + } + + if (x == 0.0L) { + r = fabsl(y); + r = sqrtl(0.5L * r); + if (y > 0.0L) + w = r + r * I; + else + w = r - r * I; + return (w); + } + + /* Rescale to avoid internal overflow or underflow. */ + if ((fabsl(x) > 4.0L) || (fabsl(y) > 4.0L)) { + x *= 0.25L; + y *= 0.25L; + scale = 2.0L; + } + else { +#if 1 + x *= 7.3786976294838206464e19; /* 2^66 */ + y *= 7.3786976294838206464e19; + scale = 1.16415321826934814453125e-10; /* 2^-33 */ +#else + x *= 4.0L; + y *= 4.0L; + scale = 0.5L; +#endif + } + w = x + y * I; + r = cabsl(w); + if (x > 0) { + t = sqrtl(0.5L * r + 0.5L * x); + r = scale * fabsl((0.5L * y) / t); + t *= scale; + } + else { + r = sqrtl(0.5L * r - 0.5L * x); + t = scale * fabsl((0.5L * y) / r); + r *= scale; + } + if (y < 0) + w = t - r * I; + else + w = t + r * I; + return (w); +} diff --git a/lib/libm/src/s_ctan.c b/lib/libm/src/s_ctan.c index 4db6fb18e83..c98dd50b38c 100644 --- a/lib/libm/src/s_ctan.c +++ b/lib/libm/src/s_ctan.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_ctan.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_ctan.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* ctan() * * Complex circular tangent @@ -56,7 +58,9 @@ * Also tested by ctan * ccot = 1 and catan(ctan(z)) = z. */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> #define MACHEP 1.1e-16 @@ -150,3 +154,12 @@ ctan(double complex z) w = sin (2.0 * creal(z)) / d + (sinh (2.0 * cimag(z)) / d) * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex ctanl(long double complex); +#else /* lint */ +__weak_alias(ctanl, ctan); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_ctanh.c b/lib/libm/src/s_ctanh.c index 30c20155672..0b0bcb2ac02 100644 --- a/lib/libm/src/s_ctanh.c +++ b/lib/libm/src/s_ctanh.c @@ -1,4 +1,4 @@ -/* $OpenBSD: s_ctanh.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ +/* $OpenBSD: s_ctanh.c,v 1.2 2011/07/08 19:25:31 martynas Exp $ */ /* * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> * @@ -15,6 +15,8 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +/* LINTLIBRARY */ + /* ctanh * * Complex hyperbolic tangent @@ -42,7 +44,9 @@ * */ +#include <sys/cdefs.h> #include <complex.h> +#include <float.h> #include <math.h> double complex @@ -57,3 +61,12 @@ ctanh(double complex z) w = sinh (2.0 * x) / d + (sin (2.0 * y) / d) * I; return (w); } + +#if LDBL_MANT_DIG == 53 +#ifdef lint +/* PROTOLIB1 */ +long double complex ctanhl(long double complex); +#else /* lint */ +__weak_alias(ctanhl, ctanh); +#endif /* lint */ +#endif /* LDBL_MANT_DIG == 53 */ diff --git a/lib/libm/src/s_ctanhl.c b/lib/libm/src/s_ctanhl.c new file mode 100644 index 00000000000..3826dcb8cd2 --- /dev/null +++ b/lib/libm/src/s_ctanhl.c @@ -0,0 +1,60 @@ +/* $OpenBSD: s_ctanhl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* ctanhl + * + * Complex hyperbolic tangent + * + * + * + * SYNOPSIS: + * + * long double complex ctanhl(); + * long double complex z, w; + * + * w = ctanhl (z); + * + * + * + * DESCRIPTION: + * + * tanh z = (sinh 2x + i sin 2y) / (cosh 2x + cos 2y) . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.7e-14 2.4e-16 + * + */ + +#include <complex.h> +#include <math.h> + +long double complex +ctanhl(long double complex z) +{ + long double complex w; + long double x, y, d; + + x = creal(z); + y = cimag(z); + d = coshl(2.0L * x) + cosl(2.0L * y); + w = sinhl(2.0L * x) / d + (sinl(2.0L * y) / d) * I; + return (w); +} diff --git a/lib/libm/src/s_ctanl.c b/lib/libm/src/s_ctanl.c new file mode 100644 index 00000000000..1ca033f6c55 --- /dev/null +++ b/lib/libm/src/s_ctanl.c @@ -0,0 +1,157 @@ +/* $OpenBSD: s_ctanl.c,v 1.1 2011/07/08 19:25:31 martynas Exp $ */ + +/* + * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> + * + * Permission to use, copy, modify, and distribute this software for any + * purpose with or without fee is hereby granted, provided that the above + * copyright notice and this permission notice appear in all copies. + * + * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + +/* ctanl() + * + * Complex circular tangent + * + * + * + * SYNOPSIS: + * + * long double complex ctanl(); + * long double complex z, w; + * + * w = ctanl( z ); + * + * + * + * DESCRIPTION: + * + * If + * z = x + iy, + * + * then + * + * sin 2x + i sinh 2y + * w = --------------------. + * cos 2x + cosh 2y + * + * On the real axis the denominator is zero at odd multiples + * of PI/2. The denominator is evaluated by its Taylor + * series near these points. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 5200 7.1e-17 1.6e-17 + * IEEE -10,+10 30000 7.2e-16 1.2e-16 + * Also tested by ctan * ccot = 1 and catan(ctan(z)) = z. + */ + +#include <complex.h> +#include <float.h> +#include <math.h> + +#if LDBL_MANT_DIG == 64 +static long double MACHEPL= 5.42101086242752217003726400434970855712890625E-20L; +#elif LDBL_MANT_DIG == 113 +static long double MACHEPL = 9.629649721936179265279889712924636592690508e-35L; +#endif + +static long double PIL = 3.141592653589793238462643383279502884197169L; +static long double DP1 = 3.14159265358979323829596852490908531763125L; +static long double DP2 = 1.6667485837041756656403424829301998703007e-19L; +static long double DP3 = 1.8830410776607851167459095484560349402753e-39L; + +static long double +redupil(long double x) +{ + long double t; + long i; + + t = x / PIL; + if (t >= 0.0L) + t += 0.5L; + else + t -= 0.5L; + + i = t; /* the multiple */ + t = i; + t = ((x - t * DP1) - t * DP2) - t * DP3; + return (t); +} + +static long double +ctansl(long double complex z) +{ + long double f, x, x2, y, y2, rn, t; + long double d; + + x = fabsl(2.0L * creal(z)); + y = fabsl(2.0L * cimag(z)); + + x = redupil(x); + + x = x * x; + y = y * y; + x2 = 1.0L; + y2 = 1.0L; + f = 1.0L; + rn = 0.0L; + d = 0.0L; + do { + rn += 1.0L; + f *= rn; + rn += 1.0L; + f *= rn; + x2 *= x; + y2 *= y; + t = y2 + x2; + t /= f; + d += t; + + rn += 1.0L; + f *= rn; + rn += 1.0L; + f *= rn; + x2 *= x; + y2 *= y; + t = y2 - x2; + t /= f; + d += t; + } + while (fabsl(t/d) > MACHEPL); + return(d); +} + +long double complex +ctanl(long double complex z) +{ + long double complex w; + long double d, x, y; + + x = creal(z); + y = cimag(z); + d = cosl(2.0L * x) + coshl(2.0L * y); + + if (fabsl(d) < 0.25L) { + d = fabsl(d); + d = ctansl(z); + } + if (d == 0.0L) { + /*mtherr( "ctan", OVERFLOW );*/ + w = LDBL_MAX + LDBL_MAX * I; + return (w); + } + + w = sinl(2.0L * x) / d + (sinhl(2.0L * y) / d) * I; + return (w); +} |