Print this page
11210 libm should be cstyle(1ONBLD) clean
   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  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23  */

  24 /*
  25  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26  * Use is subject to license terms.
  27  */
  28 
  29 #pragma weak __cpowf = cpowf
  30 
  31 #include "libm.h"
  32 #include "complex_wrapper.h"
  33 
  34 extern void sincospi(double, double *, double *);
  35 extern void sincospif(float, float *, float *);
  36 extern double atan2pi(double, double);
  37 extern float atan2pif(float, float);
  38 
  39 #if defined(__i386) && !defined(__amd64)
  40 extern int __swapRP(int);
  41 #endif
  42 
  43 static const double
  44         dpi = 3.1415926535897931160E0,  /* Hex 2^ 1 * 1.921FB54442D18 */
  45         dhalf = 0.5,

  46         dsqrt2 = 1.41421356237309514547,        /* 3FF6A09E 667F3BCD */
  47         dinvpi = 0.3183098861837906715377675;
  48 
  49 static const float one = 1.0F, zero = 0.0F;
  50 
  51 #define hiinf   0x7f800000
  52 
  53 fcomplex
  54 cpowf(fcomplex z, fcomplex w) {

  55         fcomplex        ans;
  56         float           x, y, u, v, t, c, s;
  57         double          dx, dy, du, dv, dt, dc, ds, dp, dq, dr;
  58         int             ix, iy, hx, hy, hv, hu, iu, iv, j;
  59 
  60         x = F_RE(z);
  61         y = F_IM(z);
  62         u = F_RE(w);
  63         v = F_IM(w);
  64         hx = THE_WORD(x);
  65         hy = THE_WORD(y);
  66         hu = THE_WORD(u);
  67         hv = THE_WORD(v);
  68         ix = hx & 0x7fffffff;
  69         iy = hy & 0x7fffffff;
  70         iu = hu & 0x7fffffff;
  71         iv = hv & 0x7fffffff;
  72 
  73         j = 0;

  74         if (iv == 0) {          /* z**(real) */
  75                 if (hu == 0x3f800000) { /* (anything) ** 1  is itself */
  76                         F_RE(ans) = x;
  77                         F_IM(ans) = y;
  78                 } else if (iu == 0) {   /* (anything) ** 0  is 1 */
  79                         F_RE(ans) = one;
  80                         F_IM(ans) = zero;
  81                 } else if (iy == 0) {   /* (real)**(real) */
  82                         F_IM(ans) = zero;

  83                         if (hx < 0 && ix < hiinf && iu < hiinf) {
  84                                 /* -x ** u  is exp(i*pi*u)*pow(x,u) */
  85                                 t = powf(-x, u);
  86                                 sincospif(u, &s, &c);
  87                                 F_RE(ans) = (c == zero)? c: c * t;
  88                                 F_IM(ans) = (s == zero)? s: s * t;
  89                         } else {
  90                                 F_RE(ans) = powf(x, u);
  91                         }
  92                 } else if (ix == 0 || ix >= hiinf || iy >= hiinf) {
  93                         if (ix > hiinf || iy > hiinf || iu > hiinf) {
  94                                 F_RE(ans) = F_IM(ans) = x + y + u;
  95                         } else {
  96                                 v = fabsf(y);

  97                                 if (ix != 0)
  98                                         v += fabsf(x);

  99                                 t = atan2pif(y, x);
 100                                 sincospif(t * u, &s, &c);
 101                                 F_RE(ans) = (c == zero)? c: c * v;
 102                                 F_IM(ans) = (s == zero)? s: s * v;
 103                         }
 104                 } else if (ix == iy) {  /* if |x| == |y| */
 105 #if defined(__i386) && !defined(__amd64)
 106                         int     rp = __swapRP(fp_extended);
 107 #endif
 108                         dx = (double)x;
 109                         du = (double)u;
 110                         dt = (hx >= 0)? 0.25 : 0.75;

 111                         if (hy < 0)
 112                                 dt = -dt;

 113                         dr = pow(dsqrt2 * dx, du);
 114                         sincospi(dt * du, &ds, &dc);
 115                         F_RE(ans) = (float)(dr * dc);
 116                         F_IM(ans) = (float)(dr * ds);
 117 #if defined(__i386) && !defined(__amd64)
 118                         if (rp != fp_extended)
 119                                 (void) __swapRP(rp);
 120 #endif
 121                 } else {
 122                         j = 1;
 123                 }

 124                 if (j == 0)
 125                         return (ans);
 126         }

 127         if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) {
 128                 /*
 129                  * non-zero imaginery part(s) with inf component(s) yields NaN
 130                  */
 131                 t = fabsf(x) + fabsf(y) + fabsf(u) + fabsf(v);
 132                 F_RE(ans) = F_IM(ans) = t - t;
 133         } else {
 134 #if defined(__i386) && !defined(__amd64)
 135                 int     rp = __swapRP(fp_extended);
 136 #endif
 137                 /* INDENT OFF */
 138                 /*
 139                  * r = u*log(hypot(x,y))-v*atan2(y,x),
 140                  * q = u*atan2(y,x)+v*log(hypot(x,y))
 141                  * or
 142                  * r = u*log(hypot(x,y))-v*pi*atan2pi(y,x),
 143                  * q/pi = u*atan2pi(y,x)+v*log(hypot(x,y))/pi
 144                  * ans = exp(r)*(cospi(q/pi)  + i sinpi(q/pi))
 145                  */
 146                 /* INDENT ON */
 147                 dx = (double)x;
 148                 dy = (double)y;
 149                 du = (double)u;
 150                 dv = (double)v;

 151                 if (ix > 0x3f000000 && ix < 0x40000000)   /* .5 < |x| < 2 */
 152                         dt = dhalf * log1p((dx - 1.0) * (dx + 1.0) + dy * dy);
 153                 else if (iy > 0x3f000000 && iy < 0x40000000) /* .5 < |y| < 2 */
 154                         dt = dhalf * log1p((dy - 1.0) * (dy + 1.0) + dx * dx);
 155                 else
 156                         dt = dhalf * log(dx * dx + dy * dy);

 157                 dp = atan2pi(dy, dx);

 158                 if (iv == 0) {  /* dv = 0 */
 159                         dr = exp(du * dt);
 160                         dq = du * dp;
 161                 } else {
 162                         dr = exp(du * dt - dv * dp * dpi);
 163                         dq = du * dp + dv * dt * dinvpi;
 164                 }

 165                 sincospi(dq, &ds, &dc);
 166                 F_RE(ans) = (float)(dr * dc);
 167                 F_IM(ans) = (float)(dr * ds);
 168 #if defined(__i386) && !defined(__amd64)
 169                 if (rp != fp_extended)
 170                         (void) __swapRP(rp);
 171 #endif
 172         }

 173         return (ans);
 174 }
   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 /*
  27  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __cpowf = cpowf
  32 
  33 #include "libm.h"
  34 #include "complex_wrapper.h"
  35 
  36 extern void sincospi(double, double *, double *);
  37 extern void sincospif(float, float *, float *);
  38 extern double atan2pi(double, double);
  39 extern float atan2pif(float, float);
  40 
  41 #if defined(__i386) && !defined(__amd64)
  42 extern int __swapRP(int);
  43 #endif
  44 
  45 /* Hex 2^ 1 * 1.921FB54442D18 */
  46 static const double dpi = 3.1415926535897931160E0;
  47 
  48 static const double dhalf = 0.5,
  49         dsqrt2 = 1.41421356237309514547,        /* 3FF6A09E 667F3BCD */
  50         dinvpi = 0.3183098861837906715377675;
  51 
  52 static const float one = 1.0F, zero = 0.0F;
  53 
  54 #define hiinf   0x7f800000
  55 
  56 fcomplex
  57 cpowf(fcomplex z, fcomplex w)
  58 {
  59         fcomplex ans;
  60         float x, y, u, v, t, c, s;
  61         double dx, dy, du, dv, dt, dc, ds, dp, dq, dr;
  62         int ix, iy, hx, hy, hv, hu, iu, iv, j;
  63 
  64         x = F_RE(z);
  65         y = F_IM(z);
  66         u = F_RE(w);
  67         v = F_IM(w);
  68         hx = THE_WORD(x);
  69         hy = THE_WORD(y);
  70         hu = THE_WORD(u);
  71         hv = THE_WORD(v);
  72         ix = hx & 0x7fffffff;
  73         iy = hy & 0x7fffffff;
  74         iu = hu & 0x7fffffff;
  75         iv = hv & 0x7fffffff;
  76 
  77         j = 0;
  78 
  79         if (iv == 0) {                  /* z**(real) */
  80                 if (hu == 0x3f800000) { /* (anything) ** 1  is itself */
  81                         F_RE(ans) = x;
  82                         F_IM(ans) = y;
  83                 } else if (iu == 0) {   /* (anything) ** 0  is 1 */
  84                         F_RE(ans) = one;
  85                         F_IM(ans) = zero;
  86                 } else if (iy == 0) {   /* (real)**(real) */
  87                         F_IM(ans) = zero;
  88 
  89                         if (hx < 0 && ix < hiinf && iu < hiinf) {
  90                                 /* -x ** u  is exp(i*pi*u)*pow(x,u) */
  91                                 t = powf(-x, u);
  92                                 sincospif(u, &s, &c);
  93                                 F_RE(ans) = (c == zero) ? c : c *t;
  94                                 F_IM(ans) = (s == zero) ? s : s *t;
  95                         } else {
  96                                 F_RE(ans) = powf(x, u);
  97                         }
  98                 } else if (ix == 0 || ix >= hiinf || iy >= hiinf) {
  99                         if (ix > hiinf || iy > hiinf || iu > hiinf) {
 100                                 F_RE(ans) = F_IM(ans) = x + y + u;
 101                         } else {
 102                                 v = fabsf(y);
 103 
 104                                 if (ix != 0)
 105                                         v += fabsf(x);
 106 
 107                                 t = atan2pif(y, x);
 108                                 sincospif(t * u, &s, &c);
 109                                 F_RE(ans) = (c == zero) ? c : c *v;
 110                                 F_IM(ans) = (s == zero) ? s : s *v;
 111                         }
 112                 } else if (ix == iy) {  /* if |x| == |y| */
 113 #if defined(__i386) && !defined(__amd64)
 114                         int rp = __swapRP(fp_extended);
 115 #endif
 116                         dx = (double)x;
 117                         du = (double)u;
 118                         dt = (hx >= 0) ? 0.25 : 0.75;
 119 
 120                         if (hy < 0)
 121                                 dt = -dt;
 122 
 123                         dr = pow(dsqrt2 * dx, du);
 124                         sincospi(dt * du, &ds, &dc);
 125                         F_RE(ans) = (float)(dr * dc);
 126                         F_IM(ans) = (float)(dr * ds);
 127 #if defined(__i386) && !defined(__amd64)
 128                         if (rp != fp_extended)
 129                                 (void) __swapRP(rp);
 130 #endif
 131                 } else {
 132                         j = 1;
 133                 }
 134 
 135                 if (j == 0)
 136                         return (ans);
 137         }
 138 
 139         if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) {
 140                 /*
 141                  * non-zero imaginery part(s) with inf component(s) yields NaN
 142                  */
 143                 t = fabsf(x) + fabsf(y) + fabsf(u) + fabsf(v);
 144                 F_RE(ans) = F_IM(ans) = t - t;
 145         } else {
 146 #if defined(__i386) && !defined(__amd64)
 147                 int rp = __swapRP(fp_extended);
 148 #endif
 149 
 150                 /*
 151                  * r = u*log(hypot(x,y))-v*atan2(y,x),
 152                  * q = u*atan2(y,x)+v*log(hypot(x,y))
 153                  * or
 154                  * r = u*log(hypot(x,y))-v*pi*atan2pi(y,x),
 155                  * q/pi = u*atan2pi(y,x)+v*log(hypot(x,y))/pi
 156                  * ans = exp(r)*(cospi(q/pi)  + i sinpi(q/pi))
 157                  */

 158                 dx = (double)x;
 159                 dy = (double)y;
 160                 du = (double)u;
 161                 dv = (double)v;
 162 
 163                 if (ix > 0x3f000000 && ix < 0x40000000)   /* .5 < |x| < 2 */
 164                         dt = dhalf * log1p((dx - 1.0) * (dx + 1.0) + dy * dy);
 165                 else if (iy > 0x3f000000 && iy < 0x40000000) /* .5 < |y| < 2 */
 166                         dt = dhalf * log1p((dy - 1.0) * (dy + 1.0) + dx * dx);
 167                 else
 168                         dt = dhalf * log(dx * dx + dy * dy);
 169 
 170                 dp = atan2pi(dy, dx);
 171 
 172                 if (iv == 0) {          /* dv = 0 */
 173                         dr = exp(du * dt);
 174                         dq = du * dp;
 175                 } else {
 176                         dr = exp(du * dt - dv * dp * dpi);
 177                         dq = du * dp + dv * dt * dinvpi;
 178                 }
 179 
 180                 sincospi(dq, &ds, &dc);
 181                 F_RE(ans) = (float)(dr * dc);
 182                 F_IM(ans) = (float)(dr * ds);
 183 #if defined(__i386) && !defined(__amd64)
 184                 if (rp != fp_extended)
 185                         (void) __swapRP(rp);
 186 #endif
 187         }
 188 
 189         return (ans);
 190 }