usr/src/lib/libm/common/Q/logl.c

   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 #pragma weak __logl = logl
  31
  32 /*
  33  * logl(x)
  34  * Table look-up algorithm
  35  * By K.C. Ng, March 6, 1989
  36  *
  37  * (a). For x in [31/33,33/31], using a special approximation:
  38  *      f = x - 1;
  39  *      s = f/(2.0+f);  ... here |s| <= 0.03125
  40  *      z = s*s;
  41  *      return f-s*(f-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...))));
  42  *
  43  * (b). Otherwise, normalize x = 2^n * 1.f.
  44  *      Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits,
  45  *      then
  46  *          log(x) = n*ln2 + log(1.g) + log(1.f/1.g).
  47  *      Here the leading and trailing values of log(1.g) are obtained from
  48  *      a size-64 table.
  49  *      For log(1.f/1.g), let s = (1.f-1.g)/(1.f+1.g), then
  50  *          log(1.f/1.g) = log((1+s)/(1-s)) = 2s + 2/3 s^3 + 2/5 s^5 +...
  51  *      Note that |s|<2**-8=0.00390625. We use an odd s-polynomial
  52  *      approximation to compute log(1.f/1.g):
  53  *              s*(A1+s^2*(A2+s^2*(A3+s^2*(A4+s^2*(A5+s^2*(A6+s^2*A7))))))
  54  *      (Precision is 2**-136.91 bits, absolute error)
  55  *
  56  * (c). The final result is computed by
  57  *              (n*ln2_hi+_TBL_logl_hi[j]) +
  58  *                      ( (n*ln2_lo+_TBL_logl_lo[j]) + s*(A1+...) )
  59  *
  60  * Note.
  61  *      For ln2_hi and _TBL_logl_hi[j], we force their last 32 bit to be zero
  62  *      so that n*ln2_hi + _TBL_logl_hi[j] is exact. Here
  63  *      _TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits
  64  *
  65  *
  66  * Special cases:
  67  *      log(x) is NaN with signal if x < 0 (including -INF) ;
  68  *      log(+INF) is +INF; log(0) is -INF with signal;
  69  *      log(NaN) is that NaN with no signal.
  70  *
  71  * Constants:
  72  * The hexadecimal values are the intended ones for the following constants.
  73  * The decimal values may be used, provided that the compiler will convert
  74  * from decimal to binary accurately enough to produce the hexadecimal values
  75  * shown.
  76  */
  77
  78 #include "libm.h"
  79
  80 extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
  81
  82 static const long double
  83         zero    =   0.0L,
  84         one     =   1.0L,
  85         two     =   2.0L,
  86         two113  =   10384593717069655257060992658440192.0L,
  87         ln2hi   =   6.931471805599453094172319547495844850203e-0001L,
  88         ln2lo   =   1.667085920830552208890449330400379754169e-0025L,
  89         A1      =   2.000000000000000000000000000000000000024e+0000L,
  90         A2      =   6.666666666666666666666666666666091393804e-0001L,
  91         A3      =   4.000000000000000000000000407167070220671e-0001L,
  92         A4      =   2.857142857142857142730077490612903681164e-0001L,
  93         A5      =   2.222222222222242577702836920812882605099e-0001L,
  94         A6      =   1.818181816435493395985912667105885828356e-0001L,
  95         A7      =   1.538537835211839751112067512805496931725e-0001L,
  96         B1      =   6.666666666666666666666666666666961498329e-0001L,
  97         B2      =   3.999999999999999999999999990037655042358e-0001L,
  98         B3      =   2.857142857142857142857273426428347457918e-0001L,
  99         B4      =   2.222222222222222221353229049747910109566e-0001L,
 100         B5      =   1.818181818181821503532559306309070138046e-0001L,
 101         B6      =   1.538461538453809210486356084587356788556e-0001L,
 102         B7      =   1.333333344463358756121456892645178795480e-0001L,
 103         B8      =   1.176460904783899064854645174603360383792e-0001L,
 104         B9      =   1.057293869956598995326368602518056990746e-0001L;
 105
 106 long double
 107 logl(long double x) {
 108         long double f, s, z, qn, h, t;
 109         int *px = (int *) &x;
 110         int *pz = (int *) &z;
 111         int i, j, ix, i0, i1, n;
 112
 113         /* get long double precision word ordering */
 114         if (*(int *) &one == 0) {
 115                 i0 = 3;
 116                 i1 = 0;
 117         } else {
 118                 i0 = 0;
 119                 i1 = 3;
 120         }
 121
 122         n = 0;
 123         ix = px[i0];
 124         if (ix > 0x3ffee0f8) {  /* if x >  31/33 */
 125                 if (ix < 0x3fff1084) {  /* if x < 33/31 */
 126                         f = x - one;
 127                         z = f * f;
 128                         if (((ix - 0x3fff0000) | px[i1] | px[2] | px[1]) == 0) {
 129                                 return (zero);  /* log(1)= +0 */
 130                         }
 131                         s = f / (two + f);      /* |s|<2**-8 */
 132                         z = s * s;
 133                         return (f - s * (f - z * (B1 + z * (B2 + z * (B3 +
 134                                 z * (B4 + z * (B5 + z * (B6 + z * (B7 +
 135                                 z * (B8 + z * B9))))))))));
 136                 }
 137                 if (ix >= 0x7fff0000)
 138                         return (x + x); /* x is +inf or NaN */
 139                 goto LARGE_N;
 140         }
 141         if (ix >= 0x00010000)
 142                 goto LARGE_N;
 143         i = ix & 0x7fffffff;
 144         if ((i | px[i1] | px[2] | px[1]) == 0) {
 145                 px[i0] |= 0x80000000;
 146                 return (one / x);       /* log(0.0) = -inf */
 147         }
 148         if (ix < 0) {
 149                 if ((unsigned) ix >= 0xffff0000)
 150                         return (x - x); /* x is -inf or NaN */
 151                 return (zero / zero);   /* log(x<0) is NaN  */
 152         }
 153         /* subnormal x */
 154         x *= two113;
 155         n = -113;
 156         ix = px[i0];
 157 LARGE_N:
 158         n += ((ix + 0x200) >> 16) - 0x3fff;
 159         ix = (ix & 0x0000ffff) | 0x3fff0000;    /* scale x to [1,2] */
 160         px[i0] = ix;
 161         i = ix + 0x200;
 162         pz[i0] = i & 0xfffffc00;
 163         pz[i1] = pz[1] = pz[2] = 0;
 164         s = (x - z) / (x + z);
 165         j = (i >> 10) & 0x3f;
 166         z = s * s;
 167         qn = (long double) n;
 168         t = qn * ln2lo + _TBL_logl_lo[j];
 169         h = qn * ln2hi + _TBL_logl_hi[j];
 170         f = t + s * (A1 + z * (A2 + z * (A3 + z * (A4 + z * (A5 +
 171                 z * (A6 + z * A7))))));
 172         return (h + f);
 173 }