usr/src/lib/libm/common/complex/casinl.c

   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29
  30 #pragma weak __casinl = casinl
  31
  32 #include "libm.h"               /* asinl/atanl/fabsl/isinfl/log1pl/logl/sqrtl */
  33 #include "complex_wrapper.h"
  34 #include "longdouble.h"
  35
  36 /* INDENT OFF */
  37 static const long double
  38 zero = 0.0L,
  39 one = 1.0L,
  40 Acrossover = 1.5L,
  41 Bcrossover = 0.6417L,
  42 half = 0.5L,
  43 ln2 = 6.931471805599453094172321214581765680755e-0001L,
  44 Foursqrtu = 7.3344154702193886624856495681939326638255e-2466L,  /* 2**-8189 */
  45 #if defined(__x86)
  46 E = 5.4210108624275221700372640043497085571289e-20L,    /* 2**-64 */
  47 pi_4 = 0.7853981633974483095739921312272713294078130L,
  48 pi_4_l = 4.1668714592604391641479322342670193036704898e-20L,
  49 pi_2 = 1.5707963267948966191479842624545426588156260L,
  50 pi_2_l = 8.3337429185208783282958644685340386073409796e-20L;
  51
  52 #else
  53 E = 9.6296497219361792652798897129246365926905e-35L,    /* 2**-113 */
  54 pi_4 = 0.7853981633974483096156608458198756993697670L,
  55 pi_4_l = 2.1679525325309452561992610065108379921905808e-35L,
  56 pi_2 = 1.5707963267948966192313216916397513987395340L,
  57 pi_2_l = 4.3359050650618905123985220130216759843811616e-35L;
  58
  59 #endif
  60 /* INDENT ON */
  61
  62 #if defined(__x86)
  63 static const int ip1 = 0x40400000;      /* 2**65 */
  64 #else
  65 static const int ip1 = 0x40710000;      /* 2**114 */
  66 #endif
  67
  68 ldcomplex
  69 casinl(ldcomplex z) {
  70         long double x, y, t, R, S, A, Am1, B, y2, xm1, xp1, Apx;
  71         int ix, iy, hx, hy;
  72         ldcomplex ans;
  73
  74         x = LD_RE(z);
  75         y = LD_IM(z);
  76         hx = HI_XWORD(x);
  77         hy = HI_XWORD(y);
  78         ix = hx & 0x7fffffff;
  79         iy = hy & 0x7fffffff;
  80         x = fabsl(x);
  81         y = fabsl(y);
  82
  83         /* special cases */
  84
  85         /* x is inf or NaN */
  86         if (ix >= 0x7fff0000) { /* x is inf or NaN */
  87                 if (isinfl(x)) {        /* x is INF */
  88                         LD_IM(ans) = x;
  89                         if (iy >= 0x7fff0000) {
  90                                 if (isinfl(y))
  91                                         /* casin(inf + i inf) = pi/4 + i inf */
  92                                         LD_RE(ans) = pi_4 + pi_4_l;
  93                                 else    /* casin(inf + i NaN) = NaN + i inf */
  94                                         LD_RE(ans) = y + y;
  95                         } else  /* casin(inf + iy) = pi/2 + i inf */
  96                                 LD_RE(ans) = pi_2 + pi_2_l;
  97                 } else {                /* x is NaN */
  98                         if (iy >= 0x7fff0000) {
  99                                 /* INDENT OFF */
 100                                 /*
 101                                  * casin(NaN + i inf) = NaN  + i inf
 102                                  * casin(NaN + i NaN) = NaN  + i NaN
 103                                  */
 104                                 /* INDENT ON */
 105                                 LD_IM(ans) = y + y;
 106                                 LD_RE(ans) = x + x;
 107                         } else {
 108                                 /* INDENT OFF */
 109                                 /* casin(NaN + i y ) = NaN  + i NaN */
 110                                 /* INDENT ON */
 111                                 LD_IM(ans) = LD_RE(ans) = x + y;
 112                         }
 113                 }
 114                 if (hx < 0)
 115                         LD_RE(ans) = -LD_RE(ans);
 116                 if (hy < 0)
 117                         LD_IM(ans) = -LD_IM(ans);
 118                 return (ans);
 119         }
 120
 121         /* casin(+0 + i 0) = 0 + i 0. */
 122         if (x == zero && y == zero)
 123                 return (z);
 124
 125         if (iy >= 0x7fff0000) { /* y is inf or NaN */
 126                 if (isinfl(y)) {        /* casin(x + i inf) = 0 + i inf */
 127                         LD_IM(ans) = y;
 128                         LD_RE(ans) = zero;
 129                 } else {                /* casin(x + i NaN) = NaN + i NaN */
 130                         LD_IM(ans) = x + y;
 131                         if (x == zero)
 132                                 LD_RE(ans) = x;
 133                         else
 134                                 LD_RE(ans) = y;
 135                 }
 136                 if (hx < 0)
 137                         LD_RE(ans) = -LD_RE(ans);
 138                 if (hy < 0)
 139                         LD_IM(ans) = -LD_IM(ans);
 140                 return (ans);
 141         }
 142
 143         if (y == zero) {        /* region 1: y=0 */
 144                 if (ix < 0x3fff0000) {  /* |x| < 1 */
 145                         LD_RE(ans) = asinl(x);
 146                         LD_IM(ans) = zero;
 147                 } else {
 148                         LD_RE(ans) = pi_2 + pi_2_l;
 149                         if (ix >= ip1)  /* |x| >= i386 ? 2**65 : 2**114 */
 150                                 LD_IM(ans) = ln2 + logl(x);
 151                         else if (ix >= 0x3fff8000)      /* x > Acrossover */
 152                                 LD_IM(ans) = logl(x + sqrtl((x - one) * (x +
 153                                         one)));
 154                         else {
 155                                 xm1 = x - one;
 156                                 LD_IM(ans) = log1pl(xm1 + sqrtl(xm1 * (x +
 157                                         one)));
 158                         }
 159                 }
 160         } else if (y <= E * fabsl(x - one)) {   /* region 2: y < tiny*|x-1| */
 161                 if (ix < 0x3fff0000) {  /* x < 1 */
 162                         LD_RE(ans) = asinl(x);
 163                         LD_IM(ans) = y / sqrtl((one + x) * (one - x));
 164                 } else {
 165                         LD_RE(ans) = pi_2 + pi_2_l;
 166                         if (ix >= ip1)  /* i386 ? 2**65 : 2**114 */
 167                                 LD_IM(ans) = ln2 + logl(x);
 168                         else if (ix >= 0x3fff8000)      /* x > Acrossover */
 169                                 LD_IM(ans) = logl(x + sqrtl((x - one) * (x +
 170                                         one)));
 171                         else
 172                                 LD_IM(ans) = log1pl((x - one) + sqrtl((x -
 173                                         one) * (x + one)));
 174                 }
 175         } else if (y < Foursqrtu) {     /* region 3 */
 176                 t = sqrtl(y);
 177                 LD_RE(ans) = pi_2 - (t - pi_2_l);
 178                 LD_IM(ans) = t;
 179         } else if (E * y - one >= x) {  /* region 4 */
 180                 LD_RE(ans) = x / y;     /* need to fix underflow cases */
 181                 LD_IM(ans) = ln2 + logl(y);
 182         } else if (ix >= 0x5ffb0000 || iy >= 0x5ffb0000) {
 183                 /* region 5: x+1 and y are both (>= sqrt(max)/8) i.e. 2**8188 */
 184                 t = x / y;
 185                 LD_RE(ans) = atanl(t);
 186                 LD_IM(ans) = ln2 + logl(y) + half * log1pl(t * t);
 187         } else if (x < Foursqrtu) {
 188                 /* region 6: x is very small, < 4sqrt(min) */
 189                 A = sqrtl(one + y * y);
 190                 LD_RE(ans) = x / A;     /* may underflow */
 191                 if (iy >= 0x3fff8000)   /* if y > Acrossover */
 192                         LD_IM(ans) = logl(y + A);
 193                 else
 194                         LD_IM(ans) = half * log1pl((y + y) * (y + A));
 195         } else {        /* safe region */
 196                 y2 = y * y;
 197                 xp1 = x + one;
 198                 xm1 = x - one;
 199                 R = sqrtl(xp1 * xp1 + y2);
 200                 S = sqrtl(xm1 * xm1 + y2);
 201                 A = half * (R + S);
 202                 B = x / A;
 203                 if (B <= Bcrossover)
 204                         LD_RE(ans) = asinl(B);
 205                 else {          /* use atan and an accurate approx to a-x */
 206                         Apx = A + x;
 207                         if (x <= one)
 208                                 LD_RE(ans) = atanl(x / sqrtl(half * Apx * (y2 /
 209                                         (R + xp1) + (S - xm1))));
 210                         else
 211                                 LD_RE(ans) = atanl(x / (y * sqrtl(half * (Apx /
 212                                         (R + xp1) + Apx / (S + xm1)))));
 213                 }
 214                 if (A <= Acrossover) {
 215                         /* use log1p and an accurate approx to A-1 */
 216                         if (x < one)
 217                                 Am1 = half * (y2 / (R + xp1) + y2 / (S - xm1));
 218                         else
 219                                 Am1 = half * (y2 / (R + xp1) + (S + xm1));
 220                         LD_IM(ans) = log1pl(Am1 + sqrtl(Am1 * (A + one)));
 221                 } else {
 222                         LD_IM(ans) = logl(A + sqrtl(A * A - one));
 223                 }
 224         }
 225
 226         if (hx < 0)
 227                 LD_RE(ans) = -LD_RE(ans);
 228         if (hy < 0)
 229                 LD_IM(ans) = -LD_IM(ans);
 230
 231         return (ans);
 232 }