usr/src/lib/libm/common/complex/csqrtl.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 __csqrtl = csqrtl
  31
  32 #include "libm.h"               /* fabsl/isinfl/sqrtl */
  33 #include "complex_wrapper.h"
  34 #include "longdouble.h"
  35
  36 /* INDENT OFF */
  37 static const long double
  38         twom9001 = 2.6854002716003034957421765100615693043656e-2710L,
  39         twom4500 = 2.3174987687592429423263242862381544149252e-1355L,
  40         two8999 = 9.3095991180122343502582347372163290310934e+2708L,
  41         two4500 = 4.3149968987270974283777803545571722250806e+1354L,
  42         zero = 0.0L,
  43         half = 0.5L,
  44         two = 2.0L;
  45 /* INDENT ON */
  46
  47 ldcomplex
  48 csqrtl(ldcomplex z) {
  49         ldcomplex ans;
  50         long double x, y, t, ax, ay;
  51         int n, ix, iy, hx, hy;
  52
  53         x = LD_RE(z);
  54         y = LD_IM(z);
  55         hx = HI_XWORD(x);
  56         hy = HI_XWORD(y);
  57         ix = hx & 0x7fffffff;
  58         iy = hy & 0x7fffffff;
  59         ay = fabsl(y);
  60         ax = fabsl(x);
  61         if (ix >= 0x7fff0000 || iy >= 0x7fff0000) {
  62                 /* x or y is Inf or NaN */
  63                 if (isinfl(y))
  64                         LD_IM(ans) = LD_RE(ans) = ay;
  65                 else if (isinfl(x)) {
  66                         if (hx > 0) {
  67                                 LD_RE(ans) = ax;
  68                                 LD_IM(ans) = ay * zero;
  69                         } else {
  70                                 LD_RE(ans) = ay * zero;
  71                                 LD_IM(ans) = ax;
  72                         }
  73                 } else
  74                         LD_IM(ans) = LD_RE(ans) = ax + ay;
  75         } else if (y == zero) {
  76                 if (hx >= 0) {
  77                         LD_RE(ans) = sqrtl(ax);
  78                         LD_IM(ans) = zero;
  79                 } else {
  80                         LD_IM(ans) = sqrtl(ax);
  81                         LD_RE(ans) = zero;
  82                 }
  83         } else if (ix >= iy) {
  84                 n = (ix - iy) >> 16;
  85 #if defined(__x86)              /* 64 significant bits */
  86                 if (n >= 35)
  87 #else                           /* 113 significant bits  */
  88                 if (n >= 60)
  89 #endif
  90                         t = sqrtl(ax);
  91                 else if (ix >= 0x5f3f0000) {    /* x > 2**8000 */
  92                         ax *= twom9001;
  93                         y *= twom9001;
  94                         t = two4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
  95                 } else if (iy <= 0x20bf0000) {  /* y < 2**-8000 */
  96                         ax *= two8999;
  97                         y *= two8999;
  98                         t = twom4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
  99                 } else
 100                         t = sqrtl(half * (ax + sqrtl(ax * ax + y * y)));
 101
 102                 if (hx >= 0) {
 103                         LD_RE(ans) = t;
 104                         LD_IM(ans) = ay / (t + t);
 105                 } else {
 106                         LD_IM(ans) = t;
 107                         LD_RE(ans) = ay / (t + t);
 108                 }
 109         } else {
 110                 n = (iy - ix) >> 16;
 111 #if defined(__x86)              /* 64 significant bits */
 112                 if (n >= 35) {  /* } */
 113 #else                           /* 113 significant bits  */
 114                 if (n >= 60) {
 115 #endif
 116                         if (n >= 120)
 117                                 t = sqrtl(half * ay);
 118                         else if (iy >= 0x7ffe0000)
 119                                 t = sqrtl(half * ay + half * ax);
 120                         else if (ix <= 0x00010000)
 121                                 t = half * (sqrtl(two * (ax + ay)));
 122                         else
 123                                 t = sqrtl(half * (ax + ay));
 124                 } else if (iy >= 0x5f3f0000) {  /* y > 2**8000 */
 125                         ax *= twom9001;
 126                         y *= twom9001;
 127                         t = two4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
 128                 } else if (ix <= 0x20bf0000) {
 129                         ax *= two8999;
 130                         y *= two8999;
 131                         t = twom4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
 132                 } else
 133                         t = sqrtl(half * (ax + sqrtl(ax * ax + y * y)));
 134
 135                 if (hx >= 0) {
 136                         LD_RE(ans) = t;
 137                         LD_IM(ans) = ay / (t + t);
 138                 } else {
 139                         LD_IM(ans) = t;
 140                         LD_RE(ans) = ay / (t + t);
 141                 }
 142         }
 143         if (hy < 0)
 144                 LD_IM(ans) = -LD_IM(ans);
 145         return (ans);
 146 }