5261 libm should stop using synonyms.h
5298 fabs is 0-sized, confuses dis(1) and others
Reviewed by: Josef 'Jeff' Sipek <jeffpc@josefsipek.net>
Approved by: Gordon Ross <gwr@nexenta.com>

          --- old/usr/src/lib/libm/common/complex/cexp.c
          +++ new/usr/src/lib/libm/common/complex/cexp.c

   1    1  /*
   2    2   * CDDL HEADER START
   3    3   *
   4    4   * The contents of this file are subject to the terms of the
   5    5   * Common Development and Distribution License (the "License").
   6    6   * You may not use this file except in compliance with the License.
   7    7   *
   8    8   * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9    9   * or http://www.opensolaris.org/os/licensing.
  10   10   * See the License for the specific language governing permissions
  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   23   * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24   24   */
  25   25  /*
  26   26   * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27   27   * Use is subject to license terms.
  28   28   */
  29   29  
  30      -#pragma weak cexp = __cexp
       30 +#pragma weak __cexp = cexp
  31   31  
  32   32  /* INDENT OFF */
  33   33  /*
  34   34   * dcomplex cexp(dcomplex z);
  35   35   *
  36   36   *  x+iy    x
  37   37   * e     = e  (cos(y)+i*sin(y))
  38   38   *
  39   39   * Over/underflow issue
  40   40   * --------------------

  41   41   * exp(x) may be huge but cos(y) or sin(y) may be tiny. So we use
  42   42   * function __k_cexp(x,&n) to return exp(x) = __k_cexp(x,&n)*2**n.
  43   43   * Thus if exp(x+iy) = A + Bi and t = __k_cexp(x,&n), then
  44   44   *         A = t*cos(y)*2**n,   B = t*sin(y)*2**n
  45   45   *
  46   46   * Purge off all exceptional arguments:
  47   47   *      (x,0) --> (exp(x),0)         for all x, include inf and NaN
  48   48   *      (+inf, y) --> (+inf, NaN)    for inf, nan
  49   49   *      (-inf, y) --> (+-0, +-0)     for y = inf, nan
  50   50   *      (x,+-inf/NaN) --> (NaN,NaN)  for finite x
  51   51   * For all other cases, return
  52   52   *      (x,y) --> exp(x)*cos(y)+i*exp(x)*sin(y))
  53   53   *
  54   54   * Algorithm for out of range x and finite y
  55   55   *      1. compute exp(x) in factor form (t=__k_cexp(x,&n))*2**n
  56   56   *      2. compute sincos(y,&s,&c)
  57   57   *      3. compute t*s+i*(t*c), then scale back to 2**n and return.
  58   58   */
  59   59  /* INDENT ON */
  60   60  
  61   61  #include "libm.h"               /* exp/scalbn/sincos/__k_cexp */
  62   62  #include "complex_wrapper.h"
  63   63  
  64   64  static const double zero = 0.0;
  65   65  
  66   66  dcomplex
  67   67  cexp(dcomplex z) {
  68   68          dcomplex ans;
  69   69          double x, y, t, c, s;
  70   70          int n, ix, iy, hx, hy, lx, ly;
  71   71  
  72   72          x = D_RE(z);
  73   73          y = D_IM(z);
  74   74          hx = HI_WORD(x);
  75   75          lx = LO_WORD(x);
  76   76          hy = HI_WORD(y);
  77   77          ly = LO_WORD(y);
  78   78          ix = hx & 0x7fffffff;
  79   79          iy = hy & 0x7fffffff;
  80   80          if ((iy | ly) == 0) {   /* y = 0 */
  81   81                  D_RE(ans) = exp(x);
  82   82                  D_IM(ans) = y;
  83   83          } else if (ISINF(ix, lx)) {     /* x is +-inf */
  84   84                  if (hx < 0) {
  85   85                          if (iy >= 0x7ff00000) {
  86   86                                  D_RE(ans) = zero;
  87   87                                  D_IM(ans) = zero;
  88   88                          } else {
  89   89                                  sincos(y, &s, &c);
  90   90                                  D_RE(ans) = zero * c;
  91   91                                  D_IM(ans) = zero * s;
  92   92                          }
  93   93                  } else {
  94   94                          if (iy >= 0x7ff00000) {
  95   95                                  D_RE(ans) = x;
  96   96                                  D_IM(ans) = y - y;
  97   97                          } else {
  98   98                                  (void) sincos(y, &s, &c);
  99   99                                  D_RE(ans) = x * c;
 100  100                                  D_IM(ans) = x * s;
 101  101                          }
 102  102                  }
 103  103          } else {
 104  104                  (void) sincos(y, &s, &c);
 105  105                  if (ix >= 0x40862E42) { /* |x| > 709.78... ~ log(2**1024) */
 106  106                          t = __k_cexp(x, &n);
 107  107                          D_RE(ans) = scalbn(t * c, n);
 108  108                          D_IM(ans) = scalbn(t * s, n);
 109  109                  } else {
 110  110                          t = exp(x);
 111  111                          D_RE(ans) = t * c;
 112  112                          D_IM(ans) = t * s;
 113  113                  }
 114  114          }
 115  115          return (ans);
 116  116  }

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX