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  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #ifdef __LITTLE_ENDIAN
  31 #define H0(x)   *(3 + (int *) &x)
  32 #define H1(x)   *(2 + (int *) &x)
  33 #define H2(x)   *(1 + (int *) &x)
  34 #define H3(x)   *(int *) &x
  35 #else
  36 #define H0(x)   *(int *) &x
  37 #define H1(x)   *(1 + (int *) &x)
  38 #define H2(x)   *(2 + (int *) &x)
  39 #define H3(x)   *(3 + (int *) &x)
  40 #endif
  41 
  42 /*
  43  * log1pl(x)
  44  * Table look-up algorithm by modifying logl.c
  45  * By K.C. Ng, July 6, 1995
  46  *
  47  * (a). For 1+x in [31/33,33/31], using a special approximation:
  48  *      s = x/(2.0+x);  ... here |s| <= 0.03125
  49  *      z = s*s;
  50  *      return x-s*(x-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...))));
  51  *      (i.e., x is in [-2/33,2/31])
  52  *
  53  * (b). Otherwise, normalize 1+x = 2^n * 1.f.
  54  *      Here we may need a correction term for 1+x rounded.
  55  *      Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits,
  56  *      then
  57  *          log(1+x) = n*ln2 + log(1.g) + log(1.f/1.g).
  58  *      Here the leading and trailing values of log(1.g) are obtained from
  59  *      a size-64 table.
  60  *      For log(1.f/1.g), let s = (1.f-1.g)/(1.f+1.g). Note that
  61  *              1.f = 2^-n(1+x)
  62  *
  63  *      then
  64  *          log(1.f/1.g) = log((1+s)/(1-s)) = 2s + 2/3 s^3 + 2/5 s^5 +...
  65  *      Note that |s|<2**-8=0.00390625. We use an odd s-polynomial
  66  *      approximation to compute log(1.f/1.g):
  67  *              s*(A1+s^2*(A2+s^2*(A3+s^2*(A4+s^2*(A5+s^2*(A6+s^2*A7))))))
  68  *      (Precision is 2**-136.91 bits, absolute error)
  69  *
  70  *      CAUTION:
  71  *      For x>=1, compute 1+x will lost one bit (OK).
  72  *      For x in [-0.5,-1), 1+x is exact.
  73  *      For x in (-0.5,-2/33]U[2/31,1), up to 4 last bits of x will be lost
  74  *      in 1+x.  Therefore, to recover the lost bits, one need to compute
  75  *      1.f-1.g accurately.
  76  *
  77  *      Let hx = HI(x), m = (hx>>16)-0x3fff (=ilogbl(x)), note that
  78  *              -2/33 = -0.0606...= 2^-5 * 1.939...,
  79  *               2/31 =  0.09375  = 2^-4 * 1.500...,
  80  *      so for x in (-0.5,-2/33], -5<=m<=-2,  n= -1, 1+f=2*(1+x)
  81  *         for x in [2/33,1),     -4<=m<=-1,  n=  0, f=x
  82  *
  83  *      In short:
  84  *      if x>0, let g: hg= ((hx + (0x200<<(-m)))>>(10-m))<<(10-m)
  85  *      then 1.f-1.g = x-g
  86  *      if x<0, let g': hg' =((ix-(0x200)<<(-m-1))>>(9-m))<<(9-m)
  87  *      (ix=hx&0x7fffffff)
  88  *      then 1.f-1.g = 2*(g'+x),
  89  *
  90  * (c). The final result is computed by
  91  *              (n*ln2_hi+_TBL_logl_hi[j]) +
  92  *                      ( (n*ln2_lo+_TBL_logl_lo[j]) + s*(A1+...) )
  93  *
  94  * Note.
  95  *      For ln2_hi and _TBL_logl_hi[j], we force their last 32 bit to be zero
  96  *      so that n*ln2_hi + _TBL_logl_hi[j] is exact. Here
  97  *      _TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits
  98  *
  99  *
 100  * Special cases:
 101  *      log(x) is NaN with signal if x < 0 (including -INF) ;
 102  *      log(+INF) is +INF; log(0) is -INF with signal;
 103  *      log(NaN) is that NaN with no signal.
 104  *
 105  * Constants:
 106  * The hexadecimal values are the intended ones for the following constants.
 107  * The decimal values may be used, provided that the compiler will convert
 108  * from decimal to binary accurately enough to produce the hexadecimal values
 109  * shown.
 110  */
 111 
 112 #pragma weak __log1pl = log1pl
 113 
 114 #include "libm.h"
 115 
 116 extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
 117 
 118 static const long double
 119 zero    =   0.0L,
 120 one     =   1.0L,
 121 two     =   2.0L,
 122 ln2hi   =   6.931471805599453094172319547495844850203e-0001L,
 123 ln2lo   =   1.667085920830552208890449330400379754169e-0025L,
 124 A1      =   2.000000000000000000000000000000000000024e+0000L,
 125 A2      =   6.666666666666666666666666666666091393804e-0001L,
 126 A3      =   4.000000000000000000000000407167070220671e-0001L,
 127 A4      =   2.857142857142857142730077490612903681164e-0001L,
 128 A5      =   2.222222222222242577702836920812882605099e-0001L,
 129 A6      =   1.818181816435493395985912667105885828356e-0001L,
 130 A7      =   1.538537835211839751112067512805496931725e-0001L,
 131 B1      =   6.666666666666666666666666666666961498329e-0001L,
 132 B2      =   3.999999999999999999999999990037655042358e-0001L,
 133 B3      =   2.857142857142857142857273426428347457918e-0001L,
 134 B4      =   2.222222222222222221353229049747910109566e-0001L,
 135 B5      =   1.818181818181821503532559306309070138046e-0001L,
 136 B6      =   1.538461538453809210486356084587356788556e-0001L,
 137 B7      =   1.333333344463358756121456892645178795480e-0001L,
 138 B8      =   1.176460904783899064854645174603360383792e-0001L,
 139 B9      =   1.057293869956598995326368602518056990746e-0001L;
 140 
 141 long double
 142 log1pl(long double x) {
 143         long double f, s, z, qn, h, t, y, g;
 144         int i, j, ix, iy, n, hx, m;
 145 
 146         hx = H0(x);
 147         ix = hx & 0x7fffffff;
 148         if (ix < 0x3ffaf07c) {       /* |x|<2/33 */
 149                 if (ix <= 0x3f8d0000) {      /* x <= 2**-114, return x */
 150                         if ((int) x == 0)
 151                                 return (x);
 152                 }
 153                 s = x / (two + x);      /* |s|<2**-8 */
 154                 z = s * s;
 155                 return (x - s * (x - z * (B1 + z * (B2 + z * (B3 + z * (B4 +
 156                     z * (B5 + z * (B6 + z * (B7 + z * (B8 + z * B9))))))))));
 157         }
 158         if (ix >= 0x7fff0000) {      /* x is +inf or NaN */
 159                 return (x + fabsl(x));
 160         }
 161         if (hx < 0 && ix >= 0x3fff0000) {
 162                 if (ix > 0x3fff0000 || (H1(x) | H2(x) | H3(x)) != 0)
 163                         x = zero;
 164                 return (x / zero);      /* log1p(x) is NaN  if x<-1 */
 165                 /* log1p(-1) is -inf */
 166         }
 167         if (ix >= 0x7ffeffff)
 168                 y = x;          /* avoid spurious overflow */
 169         else
 170                 y = one + x;
 171         iy = H0(y);
 172         n = ((iy + 0x200) >> 16) - 0x3fff;
 173         iy = (iy & 0x0000ffff) | 0x3fff0000;        /* scale 1+x to [1,2] */
 174         H0(y) = iy;
 175         z = zero;
 176         m = (ix >> 16) - 0x3fff;
 177         /* HI(1+x) = (((hx&0xffff)|0x10000)>>(-m))|0x3fff0000 */
 178         if (n == 0) {           /* x in [2/33,1) */
 179                 g = zero;
 180                 H0(g) = ((hx + (0x200 << (-m))) >> (10 - m)) << (10 - m);
 181                 t = x - g;
 182                 i = (((((hx & 0xffff) | 0x10000) >> (-m)) | 0x3fff0000) +
 183                         0x200) >> 10;
 184                 H0(z) = i << 10;
 185 
 186         } else if ((1 + n) == 0 && (ix < 0x3ffe0000)) {      /* x in (-0.5,-2/33] */
 187                 g = zero;
 188                 H0(g) = ((ix + (0x200 << (-m - 1))) >> (9 - m)) << (9 - m);
 189                 t = g + x;
 190                 t = t + t;
 191                 /*
 192                  * HI(2*(1+x)) =
 193                  * ((0x10000-(((hx&0xffff)|0x10000)>>(-m)))<<1)|0x3fff0000
 194                  */
 195                 /*
 196                  * i =
 197                  * ((((0x10000-(((hx&0xffff)|0x10000)>>(-m)))<<1)|0x3fff0000)+
 198                  * 0x200)>>10; H0(z)=i<<10;
 199                  */
 200                 z = two * (one - g);
 201                 i = H0(z) >> 10;
 202         } else {
 203                 i = (iy + 0x200) >> 10;
 204                 H0(z) = i << 10;
 205                 t = y - z;
 206         }
 207 
 208         s = t / (y + z);
 209         j = i & 0x3f;
 210         z = s * s;
 211         qn = (long double) n;
 212         t = qn * ln2lo + _TBL_logl_lo[j];
 213         h = qn * ln2hi + _TBL_logl_hi[j];
 214         f = t + s * (A1 + z * (A2 + z * (A3 + z * (A4 + z * (A5 + z * (A6 +
 215                 z * A7))))));
 216         return (h + f);
 217 }