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 __csqrt = csqrt
  31 
  32 /* INDENT OFF */
  33 /*
  34  * dcomplex csqrt(dcomplex z);
  35  *
  36  *                                         2    2    2
  37  * Let w=r+i*s = sqrt(x+iy). Then (r + i s)  = r  - s  + i 2sr = x + i y.
  38  *
  39  * Hence x = r*r-s*s, y = 2sr.
  40  *
  41  * Note that x*x+y*y = (s*s+r*r)**2. Thus, we have
  42  *                        ________
  43  *            2    2     / 2    2
  44  *      (1) r  + s  = \/ x  + y  ,
  45  *
  46  *            2    2
  47  *       (2) r  - s  = x
  48  *
  49  *      (3) 2sr = y.
  50  *
  51  * Perform (1)-(2) and (1)+(2), we obtain
  52  *


  83  *             2     \/  x  + y    -  x                y
  84  *            s =   ---------------------,      r = -------;    (7)
  85  *                             2                      2 s
  86  *
  87  * EXCEPTION:
  88  *
  89  * One may use the polar coordinate of a complex number to justify the
  90  * following exception cases:
  91  *
  92  * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
  93  *    csqrt(+-0+ i 0   ) =  0    + i 0
  94  *    csqrt( x + i inf ) =  inf  + i inf for all x (including NaN)
  95  *    csqrt( x + i NaN ) =  NaN  + i NaN with invalid for finite x
  96  *    csqrt(-inf+ iy   ) =  0    + i inf for finite positive-signed y
  97  *    csqrt(+inf+ iy   ) =  inf  + i 0   for finite positive-signed y
  98  *    csqrt(-inf+ i NaN) =  NaN  +-i inf
  99  *    csqrt(+inf+ i NaN) =  inf  + i NaN
 100  *    csqrt(NaN + i y  ) =  NaN  + i NaN for finite y
 101  *    csqrt(NaN + i NaN) =  NaN  + i NaN
 102  */
 103 /* INDENT ON */
 104 
 105 #include "libm.h"               /* fabs/sqrt */
 106 #include "complex_wrapper.h"
 107 
 108 /* INDENT OFF */
 109 static const double
 110         two300 = 2.03703597633448608627e+90,
 111         twom300 = 4.90909346529772655310e-91,
 112         two599 = 2.07475778444049647926e+180,
 113         twom601 = 1.20495993255144205887e-181,
 114         two = 2.0,
 115         zero = 0.0,
 116         half = 0.5;
 117 /* INDENT ON */
 118 
 119 dcomplex
 120 csqrt(dcomplex z) {

 121         dcomplex ans;
 122         double x, y, t, ax, ay;
 123         int n, ix, iy, hx, hy, lx, ly;
 124 
 125         x = D_RE(z);
 126         y = D_IM(z);
 127         hx = HI_WORD(x);
 128         lx = LO_WORD(x);
 129         hy = HI_WORD(y);
 130         ly = LO_WORD(y);
 131         ix = hx & 0x7fffffff;
 132         iy = hy & 0x7fffffff;
 133         ay = fabs(y);
 134         ax = fabs(x);

 135         if (ix >= 0x7ff00000 || iy >= 0x7ff00000) {
 136                 /* x or y is Inf or NaN */
 137                 if (ISINF(iy, ly))
 138                         D_IM(ans) = D_RE(ans) = ay;
 139                 else if (ISINF(ix, lx)) {
 140                         if (hx > 0) {
 141                                 D_RE(ans) = ax;
 142                                 D_IM(ans) = ay * zero;
 143                         } else {
 144                                 D_RE(ans) = ay * zero;
 145                                 D_IM(ans) = ax;
 146                         }
 147                 } else
 148                         D_IM(ans) = D_RE(ans) = ax + ay;

 149         } else if ((iy | ly) == 0) {    /* y = 0 */
 150                 if (hx >= 0) {
 151                         D_RE(ans) = sqrt(ax);
 152                         D_IM(ans) = zero;
 153                 } else {
 154                         D_IM(ans) = sqrt(ax);
 155                         D_RE(ans) = zero;
 156                 }
 157         } else if (ix >= iy) {
 158                 n = (ix - iy) >> 20;

 159                 if (n >= 30) {       /* x >> y or y=0 */
 160                         t = sqrt(ax);
 161                 } else if (ix >= 0x5f300000) {       /* x > 2**500 */
 162                         ax *= twom601;
 163                         y *= twom601;
 164                         t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
 165                 } else if (iy < 0x20b00000) {        /* y < 2**-500 */
 166                         ax *= two599;
 167                         y *= two599;
 168                         t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
 169                 } else
 170                         t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));


 171                 if (hx >= 0) {
 172                         D_RE(ans) = t;
 173                         D_IM(ans) = ay / (t + t);
 174                 } else {
 175                         D_IM(ans) = t;
 176                         D_RE(ans) = ay / (t + t);
 177                 }
 178         } else {
 179                 n = (iy - ix) >> 20;

 180                 if (n >= 30) {       /* y >> x */
 181                         if (n >= 60)
 182                                 t = sqrt(half * ay);
 183                         else if (iy >= 0x7fe00000)
 184                                 t = sqrt(half * ay + half * ax);
 185                         else if (ix <= 0x00100000)
 186                                 t = half * sqrt(two * (ay + ax));
 187                         else
 188                                 t = sqrt(half * (ay + ax));
 189                 } else if (iy >= 0x5f300000) {       /* y > 2**500 */
 190                         ax *= twom601;
 191                         y *= twom601;
 192                         t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
 193                 } else if (ix < 0x20b00000) {        /* x < 2**-500 */
 194                         ax *= two599;
 195                         y *= two599;
 196                         t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
 197                 } else
 198                         t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));


 199                 if (hx >= 0) {
 200                         D_RE(ans) = t;
 201                         D_IM(ans) = ay / (t + t);
 202                 } else {
 203                         D_IM(ans) = t;
 204                         D_RE(ans) = ay / (t + t);
 205                 }
 206         }

 207         if (hy < 0)
 208                 D_IM(ans) = -D_IM(ans);

 209         return (ans);
 210 }


   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 __csqrt = csqrt
  32 
  33 
  34 /*
  35  * dcomplex csqrt(dcomplex z);
  36  *
  37  *                                         2    2    2
  38  * Let w=r+i*s = sqrt(x+iy). Then (r + i s)  = r  - s  + i 2sr = x + i y.
  39  *
  40  * Hence x = r*r-s*s, y = 2sr.
  41  *
  42  * Note that x*x+y*y = (s*s+r*r)**2. Thus, we have
  43  *                        ________
  44  *            2    2     / 2    2
  45  *      (1) r  + s  = \/ x  + y  ,
  46  *
  47  *            2    2
  48  *       (2) r  - s  = x
  49  *
  50  *      (3) 2sr = y.
  51  *
  52  * Perform (1)-(2) and (1)+(2), we obtain
  53  *


  84  *             2     \/  x  + y    -  x                y
  85  *            s =   ---------------------,      r = -------;    (7)
  86  *                             2                      2 s
  87  *
  88  * EXCEPTION:
  89  *
  90  * One may use the polar coordinate of a complex number to justify the
  91  * following exception cases:
  92  *
  93  * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
  94  *    csqrt(+-0+ i 0   ) =  0    + i 0
  95  *    csqrt( x + i inf ) =  inf  + i inf for all x (including NaN)
  96  *    csqrt( x + i NaN ) =  NaN  + i NaN with invalid for finite x
  97  *    csqrt(-inf+ iy   ) =  0    + i inf for finite positive-signed y
  98  *    csqrt(+inf+ iy   ) =  inf  + i 0   for finite positive-signed y
  99  *    csqrt(-inf+ i NaN) =  NaN  +-i inf
 100  *    csqrt(+inf+ i NaN) =  inf  + i NaN
 101  *    csqrt(NaN + i y  ) =  NaN  + i NaN for finite y
 102  *    csqrt(NaN + i NaN) =  NaN  + i NaN
 103  */

 104 
 105 #include "libm.h"                       /* fabs/sqrt */
 106 #include "complex_wrapper.h"
 107 
 108 static const double two300 = 2.03703597633448608627e+90,


 109         twom300 = 4.90909346529772655310e-91,
 110         two599 = 2.07475778444049647926e+180,
 111         twom601 = 1.20495993255144205887e-181,
 112         two = 2.0,
 113         zero = 0.0,
 114         half = 0.5;
 115 
 116 
 117 dcomplex
 118 csqrt(dcomplex z)
 119 {
 120         dcomplex ans;
 121         double x, y, t, ax, ay;
 122         int n, ix, iy, hx, hy, lx, ly;
 123 
 124         x = D_RE(z);
 125         y = D_IM(z);
 126         hx = HI_WORD(x);
 127         lx = LO_WORD(x);
 128         hy = HI_WORD(y);
 129         ly = LO_WORD(y);
 130         ix = hx & 0x7fffffff;
 131         iy = hy & 0x7fffffff;
 132         ay = fabs(y);
 133         ax = fabs(x);
 134 
 135         if (ix >= 0x7ff00000 || iy >= 0x7ff00000) {
 136                 /* x or y is Inf or NaN */
 137                 if (ISINF(iy, ly)) {
 138                         D_IM(ans) = D_RE(ans) = ay;
 139                 } else if (ISINF(ix, lx)) {
 140                         if (hx > 0) {
 141                                 D_RE(ans) = ax;
 142                                 D_IM(ans) = ay * zero;
 143                         } else {
 144                                 D_RE(ans) = ay * zero;
 145                                 D_IM(ans) = ax;
 146                         }
 147                 } else {
 148                         D_IM(ans) = D_RE(ans) = ax + ay;
 149                 }
 150         } else if ((iy | ly) == 0) {    /* y = 0 */
 151                 if (hx >= 0) {
 152                         D_RE(ans) = sqrt(ax);
 153                         D_IM(ans) = zero;
 154                 } else {
 155                         D_IM(ans) = sqrt(ax);
 156                         D_RE(ans) = zero;
 157                 }
 158         } else if (ix >= iy) {
 159                 n = (ix - iy) >> 20;
 160 
 161                 if (n >= 30) {                       /* x >> y or y=0 */
 162                         t = sqrt(ax);
 163                 } else if (ix >= 0x5f300000) {       /* x > 2**500 */
 164                         ax *= twom601;
 165                         y *= twom601;
 166                         t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
 167                 } else if (iy < 0x20b00000) {        /* y < 2**-500 */
 168                         ax *= two599;
 169                         y *= two599;
 170                         t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
 171                 } else {
 172                         t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
 173                 }
 174 
 175                 if (hx >= 0) {
 176                         D_RE(ans) = t;
 177                         D_IM(ans) = ay / (t + t);
 178                 } else {
 179                         D_IM(ans) = t;
 180                         D_RE(ans) = ay / (t + t);
 181                 }
 182         } else {
 183                 n = (iy - ix) >> 20;
 184 
 185                 if (n >= 30) {               /* y >> x */
 186                         if (n >= 60)
 187                                 t = sqrt(half * ay);
 188                         else if (iy >= 0x7fe00000)
 189                                 t = sqrt(half * ay + half * ax);
 190                         else if (ix <= 0x00100000)
 191                                 t = half * sqrt(two * (ay + ax));
 192                         else
 193                                 t = sqrt(half * (ay + ax));
 194                 } else if (iy >= 0x5f300000) {       /* y > 2**500 */
 195                         ax *= twom601;
 196                         y *= twom601;
 197                         t = two300 * sqrt(ax + sqrt(ax * ax + y * y));
 198                 } else if (ix < 0x20b00000) {        /* x < 2**-500 */
 199                         ax *= two599;
 200                         y *= two599;
 201                         t = twom300 * sqrt(ax + sqrt(ax * ax + y * y));
 202                 } else {
 203                         t = sqrt(half * (ax + sqrt(ax * ax + ay * ay)));
 204                 }
 205 
 206                 if (hx >= 0) {
 207                         D_RE(ans) = t;
 208                         D_IM(ans) = ay / (t + t);
 209                 } else {
 210                         D_IM(ans) = t;
 211                         D_RE(ans) = ay / (t + t);
 212                 }
 213         }
 214 
 215         if (hy < 0)
 216                 D_IM(ans) = -D_IM(ans);
 217 
 218         return (ans);
 219 }