1129202Scognet/*- 2129202Scognet * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl 3129202Scognet * All rights reserved. 4129202Scognet * 5129202Scognet * Redistribution and use in source and binary forms, with or without 6129202Scognet * modification, are permitted provided that the following conditions 7129202Scognet * are met: 8129202Scognet * 1. Redistributions of source code must retain the above copyright 9129202Scognet * notice unmodified, this list of conditions, and the following 10129202Scognet * disclaimer. 11129202Scognet * 2. Redistributions in binary form must reproduce the above copyright 12129202Scognet * notice, this list of conditions and the following disclaimer in the 13129202Scognet * documentation and/or other materials provided with the distribution. 14129202Scognet * 15129202Scognet * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 16129202Scognet * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 17129202Scognet * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 18129202Scognet * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 19129202Scognet * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 20129202Scognet * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 21129202Scognet * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 22129202Scognet * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23129202Scognet * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24129202Scognet * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25129202Scognet */ 26129202Scognet 27129202Scognet/* 28129202Scognet * Hyperbolic sine of a complex argument z = x + i y. 29129202Scognet * 30129202Scognet * sinh(z) = sinh(x+iy) 31129202Scognet * = sinh(x) cos(y) + i cosh(x) sin(y). 32271337Sian * 33129202Scognet * Exceptional values are noted in the comments within the source code. 34129202Scognet * These values and the return value were taken from n1124.pdf. 35129202Scognet */ 36129202Scognet 37129202Scognet#include <sys/cdefs.h> 38129202Scognet__FBSDID("$FreeBSD$"); 39129202Scognet 40129202Scognet#include <complex.h> 41129202Scognet#include <math.h> 42129202Scognet 43129202Scognet#include "math_private.h" 44129202Scognet 45129202Scognetstatic const double huge = 0x1p1023; 46129202Scognet 47129202Scognetdouble complex 48129202Scognetcsinh(double complex z) 49129202Scognet{ 50129202Scognet double x, y, h; 51137464Scognet int32_t hx, hy, ix, iy, lx, ly; 52271337Sian 53129202Scognet x = creal(z); 54129202Scognet y = cimag(z); 55129202Scognet 56129202Scognet EXTRACT_WORDS(hx, lx, x); 57129202Scognet EXTRACT_WORDS(hy, ly, y); 58129202Scognet 59129202Scognet ix = 0x7fffffff & hx; 60129202Scognet iy = 0x7fffffff & hy; 61129202Scognet 62129202Scognet /* Handle the nearly-non-exceptional cases where x and y are finite. */ 63129202Scognet if (ix < 0x7ff00000 && iy < 0x7ff00000) { 64129202Scognet if ((iy | ly) == 0) 65129202Scognet return (cpack(sinh(x), y)); 66129202Scognet if (ix < 0x40360000) /* small x: normal case */ 67129202Scognet return (cpack(sinh(x) * cos(y), cosh(x) * sin(y))); 68129202Scognet 69129202Scognet /* |x| >= 22, so cosh(x) ~= exp(|x|) */ 70129202Scognet if (ix < 0x40862e42) { 71129202Scognet /* x < 710: exp(|x|) won't overflow */ 72129202Scognet h = exp(fabs(x)) * 0.5; 73129202Scognet return (cpack(copysign(h, x) * cos(y), h * sin(y))); 74137464Scognet } else if (ix < 0x4096bbaa) { 75271337Sian /* x < 1455: scale to avoid overflow */ 76129202Scognet z = __ldexp_cexp(cpack(fabs(x), y), -1); 77129202Scognet return (cpack(creal(z) * copysign(1, x), cimag(z))); 78129202Scognet } else { 79129202Scognet /* x >= 1455: the result always overflows */ 80129202Scognet h = huge * x; 81129202Scognet return (cpack(h * cos(y), h * h * sin(y))); 82129202Scognet } 83129202Scognet } 84129202Scognet 85129202Scognet /* 86129202Scognet * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN. 87129202Scognet * The sign of 0 in the result is unspecified. Choice = normally 88129202Scognet * the same as dNaN. Raise the invalid floating-point exception. 89129202Scognet * 90129202Scognet * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN). 91129202Scognet * The sign of 0 in the result is unspecified. Choice = normally 92129202Scognet * the same as d(NaN). 93129202Scognet */ 94129202Scognet if ((ix | lx) == 0 && iy >= 0x7ff00000) 95129202Scognet return (cpack(copysign(0, x * (y - y)), y - y)); 96129202Scognet 97129202Scognet /* 98129202Scognet * sinh(+-Inf +- I 0) = +-Inf + I +-0. 99129202Scognet * 100129202Scognet * sinh(NaN +- I 0) = d(NaN) + I +-0. 101129202Scognet */ 102129202Scognet if ((iy | ly) == 0 && ix >= 0x7ff00000) { 103129202Scognet if (((hx & 0xfffff) | lx) == 0) 104129202Scognet return (cpack(x, y)); 105129202Scognet return (cpack(x, copysign(0, y))); 106129202Scognet } 107129202Scognet 108129202Scognet /* 109129202Scognet * sinh(x +- I Inf) = dNaN + I dNaN. 110129202Scognet * Raise the invalid floating-point exception for finite nonzero x. 111129202Scognet * 112129202Scognet * sinh(x + I NaN) = d(NaN) + I d(NaN). 113129202Scognet * Optionally raises the invalid floating-point exception for finite 114129202Scognet * nonzero x. Choice = don't raise (except for signaling NaNs). 115129202Scognet */ 116129202Scognet if (ix < 0x7ff00000 && iy >= 0x7ff00000) 117129202Scognet return (cpack(y - y, x * (y - y))); 118129202Scognet 119129202Scognet /* 120129202Scognet * sinh(+-Inf + I NaN) = +-Inf + I d(NaN). 121129202Scognet * The sign of Inf in the result is unspecified. Choice = normally 122129202Scognet * the same as d(NaN). 123129202Scognet * 124129202Scognet * sinh(+-Inf +- I Inf) = +Inf + I dNaN. 125129202Scognet * The sign of Inf in the result is unspecified. Choice = always +. 126129202Scognet * Raise the invalid floating-point exception. 127129202Scognet * 128129202Scognet * sinh(+-Inf + I y) = +-Inf cos(y) + I Inf sin(y) 129129202Scognet */ 130129202Scognet if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { 131129202Scognet if (iy >= 0x7ff00000) 132129202Scognet return (cpack(x * x, x * (y - y))); 133129202Scognet return (cpack(x * cos(y), INFINITY * sin(y))); 134129202Scognet } 135129202Scognet 136129202Scognet /* 137129202Scognet * sinh(NaN + I NaN) = d(NaN) + I d(NaN). 138129202Scognet * 139129202Scognet * sinh(NaN +- I Inf) = d(NaN) + I d(NaN). 140129202Scognet * Optionally raises the invalid floating-point exception. 141129202Scognet * Choice = raise. 142129202Scognet * 143129202Scognet * sinh(NaN + I y) = d(NaN) + I d(NaN). 144129202Scognet * Optionally raises the invalid floating-point exception for finite 145129202Scognet * nonzero y. Choice = don't raise (except for signaling NaNs). 146129202Scognet */ 147129202Scognet return (cpack((x * x) * (y - y), (x + x) * (y - y))); 148129202Scognet} 149129202Scognet 150129202Scognetdouble complex 151129202Scognetcsin(double complex z) 152129202Scognet{ 153129202Scognet 154129202Scognet /* csin(z) = -I * csinh(I * z) */ 155129202Scognet z = csinh(cpack(-cimag(z), creal(z))); 156129202Scognet return (cpack(cimag(z), -creal(z))); 157129202Scognet} 158129202Scognet