/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8 -*-│ │vi: set et ft=c ts=8 tw=8 fenc=utf-8 :vi│ ╚──────────────────────────────────────────────────────────────────────────────╝ │ │ │ Musl Libc │ │ Copyright © 2005-2014 Rich Felker, et al. │ │ │ │ Permission is hereby granted, free of charge, to any person obtaining │ │ a copy of this software and associated documentation files (the │ │ "Software"), to deal in the Software without restriction, including │ │ without limitation the rights to use, copy, modify, merge, publish, │ │ distribute, sublicense, and/or sell copies of the Software, and to │ │ permit persons to whom the Software is furnished to do so, subject to │ │ the following conditions: │ │ │ │ The above copyright notice and this permission notice shall be │ │ included in all copies or substantial portions of the Software. │ │ │ │ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, │ │ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF │ │ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. │ │ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY │ │ CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, │ │ TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE │ │ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. │ │ │ ╚─────────────────────────────────────────────────────────────────────────────*/ #include "libc/complex.h" #include "libc/math.h" #include "libc/tinymath/complex.internal.h" asm(".ident\t\"\\n\\n\ FreeBSD libm (BSD-2 License)\\n\ Copyright (c) 2005-2011, Bruce D. Evans, Steven G. Kargl, David Schultz.\""); asm(".ident\t\"\\n\\n\ Musl libc (MIT License)\\n\ Copyright 2005-2014 Rich Felker, et. al.\""); asm(".include \"libc/disclaimer.inc\""); // clang-format off /* origin: FreeBSD /usr/src/lib/msun/src/s_ccosh.c */ /*- * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Hyperbolic cosine of a complex argument z = x + i y. * * cosh(z) = cosh(x+iy) * = cosh(x) cos(y) + i sinh(x) sin(y). * * Exceptional values are noted in the comments within the source code. * These values and the return value were taken from n1124.pdf. */ static const double huge = 0x1p1023; double complex ccosh(double complex z) { double x, y, h; int32_t hx, hy, ix, iy, lx, ly; x = creal(z); y = cimag(z); EXTRACT_WORDS(hx, lx, x); EXTRACT_WORDS(hy, ly, y); ix = 0x7fffffff & hx; iy = 0x7fffffff & hy; /* Handle the nearly-non-exceptional cases where x and y are finite. */ if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) return CMPLX(cosh(x), x * y); if (ix < 0x40360000) /* small x: normal case */ return CMPLX(cosh(x) * cos(y), sinh(x) * sin(y)); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ if (ix < 0x40862e42) { /* x < 710: exp(|x|) won't overflow */ h = exp(fabs(x)) * 0.5; return CMPLX(h * cos(y), copysign(h, x) * sin(y)); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ z = __ldexp_cexp(CMPLX(fabs(x), y), -1); return CMPLX(creal(z), cimag(z) * copysign(1, x)); } else { /* x >= 1455: the result always overflows */ h = huge * x; return CMPLX(h * h * cos(y), h * sin(y)); } } /* * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0. * The sign of 0 in the result is unspecified. Choice = normally * the same as dNaN. Raise the invalid floating-point exception. * * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0. * The sign of 0 in the result is unspecified. Choice = normally * the same as d(NaN). */ if ((ix | lx) == 0 && iy >= 0x7ff00000) return CMPLX(y - y, copysign(0, x * (y - y))); /* * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0. * * cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0. * The sign of 0 in the result is unspecified. */ if ((iy | ly) == 0 && ix >= 0x7ff00000) { if (((hx & 0xfffff) | lx) == 0) return CMPLX(x * x, copysign(0, x) * y); return CMPLX(x * x, copysign(0, (x + x) * y)); } /* * cosh(x +- I Inf) = dNaN + I dNaN. * Raise the invalid floating-point exception for finite nonzero x. * * cosh(x + I NaN) = d(NaN) + I d(NaN). * Optionally raises the invalid floating-point exception for finite * nonzero x. Choice = don't raise (except for signaling NaNs). */ if (ix < 0x7ff00000 && iy >= 0x7ff00000) return CMPLX(y - y, x * (y - y)); /* * cosh(+-Inf + I NaN) = +Inf + I d(NaN). * * cosh(+-Inf +- I Inf) = +Inf + I dNaN. * The sign of Inf in the result is unspecified. Choice = always +. * Raise the invalid floating-point exception. * * cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y) */ if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { if (iy >= 0x7ff00000) return CMPLX(x * x, x * (y - y)); return CMPLX((x * x) * cos(y), x * sin(y)); } /* * cosh(NaN + I NaN) = d(NaN) + I d(NaN). * * cosh(NaN +- I Inf) = d(NaN) + I d(NaN). * Optionally raises the invalid floating-point exception. * Choice = raise. * * cosh(NaN + I y) = d(NaN) + I d(NaN). * Optionally raises the invalid floating-point exception for finite * nonzero y. Choice = don't raise (except for signaling NaNs). */ return CMPLX((x * x) * (y - y), (x + x) * (y - y)); }