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

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/complex/cpowf.c
          +++ new/usr/src/lib/libm/common/complex/cpowf.c
↓ open down ↓ 10 lines elided ↑ open up ↑
  11   11   * and limitations under the License.
  12   12   *
  13   13   * When distributing Covered Code, include this CDDL HEADER in each
  14   14   * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15   15   * If applicable, add the following below this CDDL HEADER, with the
  16   16   * fields enclosed by brackets "[]" replaced with your own identifying
  17   17   * information: Portions Copyright [yyyy] [name of copyright owner]
  18   18   *
  19   19   * CDDL HEADER END
  20   20   */
       21 +
  21   22  /*
  22   23   * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23   24   */
       25 +
  24   26  /*
  25   27   * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26   28   * Use is subject to license terms.
  27   29   */
  28   30  
  29   31  #pragma weak __cpowf = cpowf
  30   32  
  31   33  #include "libm.h"
  32   34  #include "complex_wrapper.h"
  33   35  
  34   36  extern void sincospi(double, double *, double *);
  35   37  extern void sincospif(float, float *, float *);
  36   38  extern double atan2pi(double, double);
  37   39  extern float atan2pif(float, float);
  38   40  
  39   41  #if defined(__i386) && !defined(__amd64)
  40   42  extern int __swapRP(int);
  41   43  #endif
  42   44  
  43      -static const double
  44      -        dpi = 3.1415926535897931160E0,  /* Hex 2^ 1 * 1.921FB54442D18 */
  45      -        dhalf = 0.5,
       45 +/* Hex 2^ 1 * 1.921FB54442D18 */
       46 +static const double dpi = 3.1415926535897931160E0;
       47 +
       48 +static const double dhalf = 0.5,
  46   49          dsqrt2 = 1.41421356237309514547,        /* 3FF6A09E 667F3BCD */
  47   50          dinvpi = 0.3183098861837906715377675;
  48   51  
  49   52  static const float one = 1.0F, zero = 0.0F;
  50   53  
  51   54  #define hiinf   0x7f800000
  52   55  
  53   56  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;
       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;
  59   63  
  60   64          x = F_RE(z);
  61   65          y = F_IM(z);
  62   66          u = F_RE(w);
  63   67          v = F_IM(w);
  64   68          hx = THE_WORD(x);
  65   69          hy = THE_WORD(y);
  66   70          hu = THE_WORD(u);
  67   71          hv = THE_WORD(v);
  68   72          ix = hx & 0x7fffffff;
  69   73          iy = hy & 0x7fffffff;
  70   74          iu = hu & 0x7fffffff;
  71   75          iv = hv & 0x7fffffff;
  72   76  
  73   77          j = 0;
  74      -        if (iv == 0) {          /* z**(real) */
       78 +
       79 +        if (iv == 0) {                  /* z**(real) */
  75   80                  if (hu == 0x3f800000) { /* (anything) ** 1  is itself */
  76   81                          F_RE(ans) = x;
  77   82                          F_IM(ans) = y;
  78   83                  } else if (iu == 0) {   /* (anything) ** 0  is 1 */
  79   84                          F_RE(ans) = one;
  80   85                          F_IM(ans) = zero;
  81   86                  } else if (iy == 0) {   /* (real)**(real) */
  82   87                          F_IM(ans) = zero;
       88 +
  83   89                          if (hx < 0 && ix < hiinf && iu < hiinf) {
  84   90                                  /* -x ** u  is exp(i*pi*u)*pow(x,u) */
  85   91                                  t = powf(-x, u);
  86   92                                  sincospif(u, &s, &c);
  87      -                                F_RE(ans) = (c == zero)? c: c * t;
  88      -                                F_IM(ans) = (s == zero)? s: s * t;
       93 +                                F_RE(ans) = (c == zero) ? c : c *t;
       94 +                                F_IM(ans) = (s == zero) ? s : s *t;
  89   95                          } else {
  90   96                                  F_RE(ans) = powf(x, u);
  91   97                          }
  92   98                  } else if (ix == 0 || ix >= hiinf || iy >= hiinf) {
  93   99                          if (ix > hiinf || iy > hiinf || iu > hiinf) {
  94  100                                  F_RE(ans) = F_IM(ans) = x + y + u;
  95  101                          } else {
  96  102                                  v = fabsf(y);
      103 +
  97  104                                  if (ix != 0)
  98  105                                          v += fabsf(x);
      106 +
  99  107                                  t = atan2pif(y, x);
 100  108                                  sincospif(t * u, &s, &c);
 101      -                                F_RE(ans) = (c == zero)? c: c * v;
 102      -                                F_IM(ans) = (s == zero)? s: s * v;
      109 +                                F_RE(ans) = (c == zero) ? c : c *v;
      110 +                                F_IM(ans) = (s == zero) ? s : s *v;
 103  111                          }
 104  112                  } else if (ix == iy) {  /* if |x| == |y| */
 105  113  #if defined(__i386) && !defined(__amd64)
 106      -                        int     rp = __swapRP(fp_extended);
      114 +                        int rp = __swapRP(fp_extended);
 107  115  #endif
 108  116                          dx = (double)x;
 109  117                          du = (double)u;
 110      -                        dt = (hx >= 0)? 0.25 : 0.75;
      118 +                        dt = (hx >= 0) ? 0.25 : 0.75;
      119 +
 111  120                          if (hy < 0)
 112  121                                  dt = -dt;
      122 +
 113  123                          dr = pow(dsqrt2 * dx, du);
 114  124                          sincospi(dt * du, &ds, &dc);
 115  125                          F_RE(ans) = (float)(dr * dc);
 116  126                          F_IM(ans) = (float)(dr * ds);
 117  127  #if defined(__i386) && !defined(__amd64)
 118  128                          if (rp != fp_extended)
 119  129                                  (void) __swapRP(rp);
 120  130  #endif
 121  131                  } else {
 122  132                          j = 1;
 123  133                  }
      134 +
 124  135                  if (j == 0)
 125  136                          return (ans);
 126  137          }
      138 +
 127  139          if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) {
 128  140                  /*
 129  141                   * non-zero imaginery part(s) with inf component(s) yields NaN
 130  142                   */
 131  143                  t = fabsf(x) + fabsf(y) + fabsf(u) + fabsf(v);
 132  144                  F_RE(ans) = F_IM(ans) = t - t;
 133  145          } else {
 134  146  #if defined(__i386) && !defined(__amd64)
 135      -                int     rp = __swapRP(fp_extended);
      147 +                int rp = __swapRP(fp_extended);
 136  148  #endif
 137      -                /* INDENT OFF */
      149 +
 138  150                  /*
 139  151                   * r = u*log(hypot(x,y))-v*atan2(y,x),
 140  152                   * q = u*atan2(y,x)+v*log(hypot(x,y))
 141  153                   * or
 142  154                   * r = u*log(hypot(x,y))-v*pi*atan2pi(y,x),
 143  155                   * q/pi = u*atan2pi(y,x)+v*log(hypot(x,y))/pi
 144  156                   * ans = exp(r)*(cospi(q/pi)  + i sinpi(q/pi))
 145  157                   */
 146      -                /* INDENT ON */
 147  158                  dx = (double)x;
 148  159                  dy = (double)y;
 149  160                  du = (double)u;
 150  161                  dv = (double)v;
      162 +
 151  163                  if (ix > 0x3f000000 && ix < 0x40000000) /* .5 < |x| < 2 */
 152  164                          dt = dhalf * log1p((dx - 1.0) * (dx + 1.0) + dy * dy);
 153  165                  else if (iy > 0x3f000000 && iy < 0x40000000) /* .5 < |y| < 2 */
 154  166                          dt = dhalf * log1p((dy - 1.0) * (dy + 1.0) + dx * dx);
 155  167                  else
 156  168                          dt = dhalf * log(dx * dx + dy * dy);
      169 +
 157  170                  dp = atan2pi(dy, dx);
 158      -                if (iv == 0) {  /* dv = 0 */
      171 +
      172 +                if (iv == 0) {          /* dv = 0 */
 159  173                          dr = exp(du * dt);
 160  174                          dq = du * dp;
 161  175                  } else {
 162  176                          dr = exp(du * dt - dv * dp * dpi);
 163  177                          dq = du * dp + dv * dt * dinvpi;
 164  178                  }
      179 +
 165  180                  sincospi(dq, &ds, &dc);
 166  181                  F_RE(ans) = (float)(dr * dc);
 167  182                  F_IM(ans) = (float)(dr * ds);
 168  183  #if defined(__i386) && !defined(__amd64)
 169  184                  if (rp != fp_extended)
 170  185                          (void) __swapRP(rp);
 171  186  #endif
 172  187          }
      188 +
 173  189          return (ans);
 174  190  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX