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

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/Q/expm1l.c
          +++ new/usr/src/lib/libm/common/Q/expm1l.c
↓ open down ↓ 14 lines elided ↑ open up ↑
  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   27   * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27   28   * Use is subject to license terms.
  28   29   */
  29   30  
  30   31  #pragma weak __expm1l = expm1l
  31   32  
  32   33  #if !defined(__sparc)
  33   34  #error Unsupported architecture
  34   35  #endif
↓ open down ↓ 1 lines elided ↑ open up ↑
  36   37  /*
  37   38   * expm1l(x)
  38   39   *
  39   40   * Table driven method
  40   41   * Written by K.C. Ng, June 1995.
  41   42   * Algorithm :
  42   43   *      1. expm1(x) = x if x<2**-114
  43   44   *      2. if |x| <= 0.0625 = 1/16, use approximation
  44   45   *              expm1(x) = x + x*P/(2-P)
  45   46   * where
  46      - *      P = x - z*(P1+z*(P2+z*(P3+z*(P4+z*(P5+z*P6+z*P7))))), z = x*x;
       47 + *      P = x - z*(P1+z*(P2+z*(P3+z*(P4+z*(P5+z*P6+z*P7))))), z = x*x;
  47   48   * (this formula is derived from
  48   49   *      2-P+x = R = x*(exp(x)+1)/(exp(x)-1) ~ 2 + x*x/6 - x^4/360 + ...)
  49   50   *
  50   51   * P1 =   1.66666666666666666666666666666638500528074603030e-0001
  51   52   * P2 =  -2.77777777777777777777777759668391122822266551158e-0003
  52   53   * P3 =   6.61375661375661375657437408890138814721051293054e-0005
  53   54   * P4 =  -1.65343915343915303310185228411892601606669528828e-0006
  54   55   * P5 =   4.17535139755122945763580609663414647067443411178e-0008
  55   56   * P6 =  -1.05683795988668526689182102605260986731620026832e-0009
  56   57   * P7 =   2.67544168821852702827123344217198187229611470514e-0011
↓ open down ↓ 36 lines elided ↑ open up ↑
  93   94   * Constants:
  94   95   * Only decimal values are given. We assume that the compiler will convert
  95   96   * from decimal to binary accurately enough to produce the correct
  96   97   * hexadecimal values.
  97   98   */
  98   99  
  99  100  #include "libm.h"
 100  101  
 101  102  extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
 102  103  extern const long double _TBL_expm1lx[], _TBL_expm1l[];
 103      -
 104      -static const long double
 105      -        zero            = +0.0L,
 106      -        one             = +1.0L,
 107      -        two             = +2.0L,
 108      -        ln2_64          = +1.083042469624914545964425189778400898568e-2L,
 109      -        ovflthreshold   = +1.135652340629414394949193107797076342845e+4L,
 110      -        invln2_32       = +4.616624130844682903551758979206054839765e+1L,
 111      -        ln2_32hi        = +2.166084939249829091928849858592451515688e-2L,
 112      -        ln2_32lo        = +5.209643502595475652782654157501186731779e-27L,
 113      -        huge            = +1.0e4000L,
 114      -        tiny            = +1.0e-4000L,
      104 +static const long double zero = +0.0L,
      105 +        one = +1.0L,
      106 +        two = +2.0L,
      107 +        ln2_64 = +1.083042469624914545964425189778400898568e-2L,
      108 +        ovflthreshold = +1.135652340629414394949193107797076342845e+4L,
      109 +        invln2_32 = +4.616624130844682903551758979206054839765e+1L,
      110 +        ln2_32hi = +2.166084939249829091928849858592451515688e-2L,
      111 +        ln2_32lo = +5.209643502595475652782654157501186731779e-27L,
      112 +        huge = +1.0e4000L,
      113 +        tiny = +1.0e-4000L,
 115  114          P1 = +1.66666666666666666666666666666638500528074603030e-0001L,
 116  115          P2 = -2.77777777777777777777777759668391122822266551158e-0003L,
 117  116          P3 = +6.61375661375661375657437408890138814721051293054e-0005L,
 118  117          P4 = -1.65343915343915303310185228411892601606669528828e-0006L,
 119  118          P5 = +4.17535139755122945763580609663414647067443411178e-0008L,
 120  119          P6 = -1.05683795988668526689182102605260986731620026832e-0009L,
 121  120          P7 = +2.67544168821852702827123344217198187229611470514e-0011L,
 122  121  /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
 123  122          T1 = +1.666666666666666666666666666660876387437e-1L,
 124  123          T2 = -2.777777777777777777777707812093173478756e-3L,
 125  124          T3 = +6.613756613756613482074280932874221202424e-5L,
 126  125          T4 = -1.653439153392139954169609822742235851120e-6L,
 127  126          T5 = +4.175314851769539751387852116610973796053e-8L;
 128  127  
 129  128  long double
 130      -expm1l(long double x) {
      129 +expm1l(long double x)
      130 +{
 131  131          int hx, ix, j, k, m;
 132  132          long double t, r, s, w;
 133  133  
 134      -        hx = ((int *) &x)[HIXWORD];
      134 +        hx = ((int *)&x)[HIXWORD];
 135  135          ix = hx & ~0x80000000;
      136 +
 136  137          if (ix >= 0x7fff0000) {
 137  138                  if (x != x)
 138      -                        return (x + x); /* NaN */
      139 +                        return (x + x);                 /* NaN */
      140 +
 139  141                  if (x < zero)
 140      -                        return (-one);  /* -inf */
 141      -                return (x);     /* +inf */
      142 +                        return (-one);                  /* -inf */
      143 +
      144 +                return (x);                             /* +inf */
 142  145          }
 143      -        if (ix < 0x3fff4000) {  /* |x| < 1.25 */
 144      -                if (ix < 0x3ffb0000) {  /* |x| < 0.0625 */
      146 +
      147 +        if (ix < 0x3fff4000) {                          /* |x| < 1.25 */
      148 +                if (ix < 0x3ffb0000) {                  /* |x| < 0.0625 */
 145  149                          if (ix < 0x3f8d0000) {
 146      -                                if ((int) x == 0)
      150 +                                if ((int)x == 0)
 147  151                                          return (x);     /* |x|<2^-114 */
 148  152                          }
      153 +
 149  154                          t = x * x;
 150  155                          r = (x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t *
 151      -                                (P5 + t * (P6 + t * P7)))))));
      156 +                            (P5 + t * (P6 + t * P7)))))));
 152  157                          return (x + (x * r) / (two - r));
 153  158                  }
      159 +
 154  160                  /* compute i = [64*x] */
 155  161                  m = 0x4009 - (ix >> 16);
 156  162                  j = ((ix & 0x0000ffff) | 0x10000) >> m; /* j=4,...,67 */
      163 +
 157  164                  if (hx < 0)
 158  165                          j += 82;                        /* negative */
      166 +
 159  167                  s = x - _TBL_expm1lx[j];
 160  168                  t = s * s;
 161  169                  r = s - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))));
 162  170                  r = (s + s) / (two - r);
 163  171                  w = _TBL_expm1l[j];
 164  172                  return (w + (w + one) * r);
 165  173          }
      174 +
 166  175          if (hx > 0) {
 167  176                  if (x > ovflthreshold)
 168  177                          return (huge * huge);
 169      -                k = (int) (invln2_32 * (x + ln2_64));
      178 +
      179 +                k = (int)(invln2_32 * (x + ln2_64));
 170  180          } else {
 171  181                  if (x < -80.0)
 172  182                          return (tiny - x / x);
 173      -                k = (int) (invln2_32 * (x - ln2_64));
      183 +
      184 +                k = (int)(invln2_32 * (x - ln2_64));
 174  185          }
      186 +
 175  187          j = k & 0x1f;
 176  188          m = k >> 5;
 177      -        t = (long double) k;
      189 +        t = (long double)k;
 178  190          x = (x - t * ln2_32hi) - t * ln2_32lo;
 179  191          t = x * x;
 180  192          r = (x - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))))) - two;
 181  193          x = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
 182      -                _TBL_expl_lo[j]);
      194 +            _TBL_expl_lo[j]);
 183  195          return (scalbnl(x, m) - one);
 184  196  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX