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 /*
  31  * expl(x)
  32  * Table driven method
  33  * Written by K.C. Ng, November 1988.
  34  * Algorithm :
  35  *      1. Argument Reduction: given the input x, find r and integer k
  36  *         and j such that
  37  *                   x = (32k+j)*ln2 + r,  |r| <= (1/64)*ln2 .
  38  *
  39  *      2. expl(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r))
  40  *         Note:
  41  *         a. expm1(r) = (2r)/(2-R), R = r - r^2*(t1 + t2*r^2)
  42  *         b. 2^(j/32) is represented as
  43  *                      _TBL_expl_hi[j]+_TBL_expl_lo[j]
  44  *         where


  55  *      an ulp (unit in the last place).
  56  *
  57  * Misc. info.
  58  *      For 113 bit long double
  59  *              if x >  1.135652340629414394949193107797076342845e+4
  60  *      then expl(x) overflow;
  61  *              if x < -1.143346274333629787883724384345262150341e+4
  62  *      then expl(x) underflow
  63  *
  64  * Constants:
  65  * Only decimal values are given. We assume that the compiler will convert
  66  * from decimal to binary accurately enough to produce the correct
  67  * hexadecimal values.
  68  */
  69 
  70 #pragma weak __expl = expl
  71 
  72 #include "libm.h"
  73 
  74 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
  75 
  76 static const long double
  77 one             =  1.0L,
  78 two             =  2.0L,
  79 ln2_64          =  1.083042469624914545964425189778400898568e-2L,
  80 ovflthreshold   =  1.135652340629414394949193107797076342845e+4L,
  81 unflthreshold   = -1.143346274333629787883724384345262150341e+4L,
  82 invln2_32       =  4.616624130844682903551758979206054839765e+1L,
  83 ln2_32hi        =  2.166084939249829091928849858592451515688e-2L,
  84 ln2_32lo        =  5.209643502595475652782654157501186731779e-27L;
  85 
  86 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
  87 static const long double
  88 t1 =   1.666666666666666666666666666660876387437e-1L,
  89 t2 =  -2.777777777777777777777707812093173478756e-3L,
  90 t3 =   6.613756613756613482074280932874221202424e-5L,
  91 t4 =  -1.653439153392139954169609822742235851120e-6L,
  92 t5 =   4.175314851769539751387852116610973796053e-8L;
  93 
  94 long double
  95 expl(long double x) {
  96         int *px = (int *) &x, ix, j, k, m;

  97         long double t, r;
  98 
  99         ix = px[0];                             /* high word of x */

 100         if (ix >= 0x7fff0000)
 101                 return (x + x);                 /* NaN of +inf */
 102         if (((unsigned) ix) >= 0xffff0000)

 103                 return (-one / x);              /* NaN or -inf */

 104         if ((ix & 0x7fffffff) < 0x3fc30000) {
 105                 if ((int) x < 1)
 106                         return (one + x);       /* |x|<2^-60 */
 107         }

 108         if (ix > 0) {
 109                 if (x > ovflthreshold)
 110                         return (scalbnl(x, 20000));
 111                 k = (int) (invln2_32 * (x + ln2_64));

 112         } else {
 113                 if (x < unflthreshold)
 114                         return (scalbnl(-x, -40000));
 115                 k = (int) (invln2_32 * (x - ln2_64));

 116         }
 117         j  = k&0x1f;
 118         m  = k>>5;
 119         t  = (long double) k;

 120         x  = (x - t * ln2_32hi) - t * ln2_32lo;
 121         t  = x * x;
 122         r  = (x - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - two;
 123         x  = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
 124                 _TBL_expl_lo[j]);
 125         return (scalbnl(x, m));
 126 }


   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 /*
  32  * expl(x)
  33  * Table driven method
  34  * Written by K.C. Ng, November 1988.
  35  * Algorithm :
  36  *      1. Argument Reduction: given the input x, find r and integer k
  37  *         and j such that
  38  *                   x = (32k+j)*ln2 + r,  |r| <= (1/64)*ln2 .
  39  *
  40  *      2. expl(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r))
  41  *         Note:
  42  *         a. expm1(r) = (2r)/(2-R), R = r - r^2*(t1 + t2*r^2)
  43  *         b. 2^(j/32) is represented as
  44  *                      _TBL_expl_hi[j]+_TBL_expl_lo[j]
  45  *         where


  56  *      an ulp (unit in the last place).
  57  *
  58  * Misc. info.
  59  *      For 113 bit long double
  60  *              if x >  1.135652340629414394949193107797076342845e+4
  61  *      then expl(x) overflow;
  62  *              if x < -1.143346274333629787883724384345262150341e+4
  63  *      then expl(x) underflow
  64  *
  65  * Constants:
  66  * Only decimal values are given. We assume that the compiler will convert
  67  * from decimal to binary accurately enough to produce the correct
  68  * hexadecimal values.
  69  */
  70 
  71 #pragma weak __expl = expl
  72 
  73 #include "libm.h"
  74 
  75 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
  76 static const long double one = 1.0L,
  77         two = 2.0L,
  78         ln2_64 = 1.083042469624914545964425189778400898568e-2L,
  79         ovflthreshold = 1.135652340629414394949193107797076342845e+4L,
  80         unflthreshold = -1.143346274333629787883724384345262150341e+4L,
  81         invln2_32 = 4.616624130844682903551758979206054839765e+1L,
  82         ln2_32hi = 2.166084939249829091928849858592451515688e-2L,
  83         ln2_32lo = 5.209643502595475652782654157501186731779e-27L;


  84 
  85 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
  86 static const long double t1 = 1.666666666666666666666666666660876387437e-1L,
  87         t2 = -2.777777777777777777777707812093173478756e-3L,
  88         t3 = 6.613756613756613482074280932874221202424e-5L,
  89         t4 = -1.653439153392139954169609822742235851120e-6L,
  90         t5 = 4.175314851769539751387852116610973796053e-8L;

  91 
  92 long double
  93 expl(long double x)
  94 {
  95         int *px = (int *)&x, ix, j, k, m;
  96         long double t, r;
  97 
  98         ix = px[0];                     /* high word of x */
  99 
 100         if (ix >= 0x7fff0000)
 101                 return (x + x);         /* NaN of +inf */
 102 
 103         if (((unsigned)ix) >= 0xffff0000)
 104                 return (-one / x);      /* NaN or -inf */
 105 
 106         if ((ix & 0x7fffffff) < 0x3fc30000) {
 107                 if ((int)x < 1)
 108                         return (one + x);       /* |x|<2^-60 */
 109         }
 110 
 111         if (ix > 0) {
 112                 if (x > ovflthreshold)
 113                         return (scalbnl(x, 20000));
 114 
 115                 k = (int)(invln2_32 * (x + ln2_64));
 116         } else {
 117                 if (x < unflthreshold)
 118                         return (scalbnl(-x, -40000));
 119 
 120                 k = (int)(invln2_32 * (x - ln2_64));
 121         }
 122 
 123         j = k & 0x1f;
 124         m = k >> 5;
 125         t = (long double)k;
 126         x = (x - t * ln2_32hi) - t * ln2_32lo;
 127         t = x * x;
 128         r = (x - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - two;
 129         x = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
 130             _TBL_expl_lo[j]);
 131         return (scalbnl(x, m));
 132 }