Print this page
11210 libm should be cstyle(1ONBLD) clean


   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 }


   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 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __csqrtl = csqrtl
  32 
  33 #include "libm.h"                       /* fabsl/isinfl/sqrtl */
  34 #include "complex_wrapper.h"
  35 #include "longdouble.h"
  36 

  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 
  46 
  47 ldcomplex
  48 csqrtl(ldcomplex z)
  49 {
  50         ldcomplex ans;
  51         long double x, y, t, ax, ay;
  52         int n, ix, iy, hx, hy;
  53 
  54         x = LD_RE(z);
  55         y = LD_IM(z);
  56         hx = HI_XWORD(x);
  57         hy = HI_XWORD(y);
  58         ix = hx & 0x7fffffff;
  59         iy = hy & 0x7fffffff;
  60         ay = fabsl(y);
  61         ax = fabsl(x);
  62 
  63         if (ix >= 0x7fff0000 || iy >= 0x7fff0000) {
  64                 /* x or y is Inf or NaN */
  65                 if (isinfl(y)) {
  66                         LD_IM(ans) = LD_RE(ans) = ay;
  67                 } else if (isinfl(x)) {
  68                         if (hx > 0) {
  69                                 LD_RE(ans) = ax;
  70                                 LD_IM(ans) = ay * zero;
  71                         } else {
  72                                 LD_RE(ans) = ay * zero;
  73                                 LD_IM(ans) = ax;
  74                         }
  75                 } else {
  76                         LD_IM(ans) = LD_RE(ans) = ax + ay;
  77                 }
  78         } else if (y == zero) {
  79                 if (hx >= 0) {
  80                         LD_RE(ans) = sqrtl(ax);
  81                         LD_IM(ans) = zero;
  82                 } else {
  83                         LD_IM(ans) = sqrtl(ax);
  84                         LD_RE(ans) = zero;
  85                 }
  86         } else if (ix >= iy) {
  87                 n = (ix - iy) >> 16;
  88 #if defined(__x86)                      /* 64 significant bits */
  89                 if (n >= 35)
  90 #else /* 113 significant bits  */
  91                 if (n >= 60)
  92 #endif
  93                 {
  94                         t = sqrtl(ax);
  95                 } else if (ix >= 0x5f3f0000) {       /* x > 2**8000 */
  96                         ax *= twom9001;
  97                         y *= twom9001;
  98                         t = two4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
  99                 } else if (iy <= 0x20bf0000) {       /* y < 2**-8000 */
 100                         ax *= two8999;
 101                         y *= two8999;
 102                         t = twom4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
 103                 } else {
 104                         t = sqrtl(half * (ax + sqrtl(ax * ax + y * y)));
 105                 }
 106 
 107                 if (hx >= 0) {
 108                         LD_RE(ans) = t;
 109                         LD_IM(ans) = ay / (t + t);
 110                 } else {
 111                         LD_IM(ans) = t;
 112                         LD_RE(ans) = ay / (t + t);
 113                 }
 114         } else {
 115                 n = (iy - ix) >> 16;
 116 #if defined(__x86)                      /* 64 significant bits */
 117                 if (n >= 35) {               /* } */
 118 #else /* 113 significant bits  */
 119                 if (n >= 60) {
 120 #endif
 121 
 122                         if (n >= 120)
 123                                 t = sqrtl(half * ay);
 124                         else if (iy >= 0x7ffe0000)
 125                                 t = sqrtl(half * ay + half * ax);
 126                         else if (ix <= 0x00010000)
 127                                 t = half * (sqrtl(two * (ax + ay)));
 128                         else
 129                                 t = sqrtl(half * (ax + ay));
 130                 } else if (iy >= 0x5f3f0000) {       /* y > 2**8000 */
 131                         ax *= twom9001;
 132                         y *= twom9001;
 133                         t = two4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
 134                 } else if (ix <= 0x20bf0000) {
 135                         ax *= two8999;
 136                         y *= two8999;
 137                         t = twom4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
 138                 } else {
 139                         t = sqrtl(half * (ax + sqrtl(ax * ax + y * y)));
 140                 }
 141 
 142                 if (hx >= 0) {
 143                         LD_RE(ans) = t;
 144                         LD_IM(ans) = ay / (t + t);
 145                 } else {
 146                         LD_IM(ans) = t;
 147                         LD_RE(ans) = ay / (t + t);
 148                 }
 149         }
 150 
 151         if (hy < 0)
 152                 LD_IM(ans) = -LD_IM(ans);
 153 
 154         return (ans);
 155 }