lib/msun/ld128/s_logl.c

   1 /*-
   2  * Copyright (c) 2007-2013 Bruce D. Evans
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  * 1. Redistributions of source code must retain the above copyright
   9  *    notice unmodified, this list of conditions, and the following
  10  *    disclaimer.
  11  * 2. Redistributions in binary form must reproduce the above copyright
  12  *    notice, this list of conditions and the following disclaimer in the
  13  *    documentation and/or other materials provided with the distribution.
  14  *
  15  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  16  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
  17  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  18  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
  19  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  20  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  21  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  22  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  23  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  24  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  25  */
  26
  27 #include <sys/cdefs.h>
  28 __FBSDID("$FreeBSD$");
  29
  30 /**
  31  * Implementation of the natural logarithm of x for 128-bit format.
  32  *
  33  * First decompose x into its base 2 representation:
  34  *
  35  *    log(x) = log(X * 2**k), where X is in [1, 2)
  36  *           = log(X) + k * log(2).
  37  *
  38  * Let X = X_i + e, where X_i is the center of one of the intervals
  39  * [-1.0/256, 1.0/256), [1.0/256, 3.0/256), .... [2.0-1.0/256, 2.0+1.0/256)
  40  * and X is in this interval.  Then
  41  *
  42  *    log(X) = log(X_i + e)
  43  *           = log(X_i * (1 + e / X_i))
  44  *           = log(X_i) + log(1 + e / X_i).
  45  *
  46  * The values log(X_i) are tabulated below.  Let d = e / X_i and use
  47  *
  48  *    log(1 + d) = p(d)
  49  *
  50  * where p(d) = d - 0.5*d*d + ... is a special minimax polynomial of
  51  * suitably high degree.
  52  *
  53  * To get sufficiently small roundoff errors, k * log(2), log(X_i), and
  54  * sometimes (if |k| is not large) the first term in p(d) must be evaluated
  55  * and added up in extra precision.  Extra precision is not needed for the
  56  * rest of p(d).  In the worst case when k = 0 and log(X_i) is 0, the final
  57  * error is controlled mainly by the error in the second term in p(d).  The
  58  * error in this term itself is at most 0.5 ulps from the d*d operation in
  59  * it.  The error in this term relative to the first term is thus at most
  60  * 0.5 * |-0.5| * |d| < 1.0/1024 ulps.  We aim for an accumulated error of
  61  * at most twice this at the point of the final rounding step.  Thus the
  62  * final error should be at most 0.5 + 1.0/512 = 0.5020 ulps.  Exhaustive
  63  * testing of a float variant of this function showed a maximum final error
  64  * of 0.5008 ulps.  Non-exhaustive testing of a double variant of this
  65  * function showed a maximum final error of 0.5078 ulps (near 1+1.0/256).
  66  *
  67  * We made the maximum of |d| (and thus the total relative error and the
  68  * degree of p(d)) small by using a large number of intervals.  Using
  69  * centers of intervals instead of endpoints reduces this maximum by a
  70  * factor of 2 for a given number of intervals.  p(d) is special only
  71  * in beginning with the Taylor coefficients 0 + 1*d, which tends to happen
  72  * naturally.  The most accurate minimax polynomial of a given degree might
  73  * be different, but then we wouldn't want it since we would have to do
  74  * extra work to avoid roundoff error (especially for P0*d instead of d).
  75  */
  76
  77 #ifdef DEBUG
  78 #include <assert.h>
  79 #include <fenv.h>
  80 #endif
  81
  82 #include "fpmath.h"
  83 #include "math.h"
  84 #ifndef NO_STRUCT_RETURN
  85 #define STRUCT_RETURN
  86 #endif
  87 #include "math_private.h"
  88
  89 #if !defined(NO_UTAB) && !defined(NO_UTABL)
  90 #define USE_UTAB
  91 #endif
  92
  93 /*
  94  * Domain [-0.005280, 0.004838], range ~[-1.1577e-37, 1.1582e-37]:
  95  * |log(1 + d)/d - p(d)| < 2**-122.7
  96  */
  97 static const long double
  98 P2 = -0.5L,
  99 P3 =  3.33333333333333333333333333333233795e-1L,        /*  0x15555555555555555555555554d42.0p-114L */
 100 P4 = -2.49999999999999999999999999941139296e-1L,        /* -0x1ffffffffffffffffffffffdab14e.0p-115L */
 101 P5 =  2.00000000000000000000000085468039943e-1L,        /*  0x19999999999999999999a6d3567f4.0p-115L */
 102 P6 = -1.66666666666666666666696142372698408e-1L,        /* -0x15555555555555555567267a58e13.0p-115L */
 103 P7 =  1.42857142857142857119522943477166120e-1L,        /*  0x1249249249249248ed79a0ae434de.0p-115L */
 104 P8 = -1.24999999999999994863289015033581301e-1L;        /* -0x1fffffffffffffa13e91765e46140.0p-116L */
 105 /* Double precision gives ~ 53 + log2(P9 * max(|d|)**8) ~= 120 bits. */
 106 static const double
 107 P9 =  1.1111111111111401e-1,            /*  0x1c71c71c71c7ed.0p-56 */
 108 P10 = -1.0000000000040135e-1,           /* -0x199999999a0a92.0p-56 */
 109 P11 =  9.0909090728136258e-2,           /*  0x1745d173962111.0p-56 */
 110 P12 = -8.3333318851855284e-2,           /* -0x1555551722c7a3.0p-56 */
 111 P13 =  7.6928634666404178e-2,           /*  0x13b1985204a4ae.0p-56 */
 112 P14 = -7.1626810078462499e-2;           /* -0x12562276cdc5d0.0p-56 */
 113
 114 static volatile const double zero = 0;
 115
 116 #define INTERVALS       128
 117 #define LOG2_INTERVALS  7
 118 #define TSIZE           (INTERVALS + 1)
 119 #define G(i)            (T[(i)].G)
 120 #define F_hi(i)         (T[(i)].F_hi)
 121 #define F_lo(i)         (T[(i)].F_lo)
 122 #define ln2_hi          F_hi(TSIZE - 1)
 123 #define ln2_lo          F_lo(TSIZE - 1)
 124 #define E(i)            (U[(i)].E)
 125 #define H(i)            (U[(i)].H)
 126
 127 static const struct {
 128         float   G;                      /* 1/(1 + i/128) rounded to 8/9 bits */
 129         float   F_hi;                   /* log(1 / G_i) rounded (see below) */
 130         /* The compiler will insert 8 bytes of padding here. */
 131         long double F_lo;               /* next 113 bits for log(1 / G_i) */
 132 } T[TSIZE] = {
 133         /*
 134          * ln2_hi and each F_hi(i) are rounded to a number of bits that
 135          * makes F_hi(i) + dk*ln2_hi exact for all i and all dk.
 136          *
 137          * The last entry (for X just below 2) is used to define ln2_hi
 138          * and ln2_lo, to ensure that F_hi(i) and F_lo(i) cancel exactly
 139          * with dk*ln2_hi and dk*ln2_lo, respectively, when dk = -1.
 140          * This is needed for accuracy when x is just below 1.  (To avoid
 141          * special cases, such x are "reduced" strangely to X just below
 142          * 2 and dk = -1, and then the exact cancellation is needed
 143          * because any the error from any non-exactness would be too
 144          * large).
 145          *
 146          * The relevant range of dk is [-16445, 16383].  The maximum number
 147          * of bits in F_hi(i) that works is very dependent on i but has
 148          * a minimum of 93.  We only need about 12 bits in F_hi(i) for
 149          * it to provide enough extra precision.
 150          *
 151          * We round F_hi(i) to 24 bits so that it can have type float,
 152          * mainly to minimize the size of the table.  Using all 24 bits
 153          * in a float for it automatically satisfies the above constraints.
 154          */
 155      0x800000.0p-23,  0,               0,
 156      0xfe0000.0p-24,  0x8080ac.0p-30, -0x14ee431dae6674afa0c4bfe16e8fd.0p-144L,
 157      0xfc0000.0p-24,  0x8102b3.0p-29, -0x1db29ee2d83717be918e1119642ab.0p-144L,
 158      0xfa0000.0p-24,  0xc24929.0p-29,  0x1191957d173697cf302cc9476f561.0p-143L,
 159      0xf80000.0p-24,  0x820aec.0p-28,  0x13ce8888e02e78eba9b1113bc1c18.0p-142L,
 160      0xf60000.0p-24,  0xa33577.0p-28, -0x17a4382ce6eb7bfa509bec8da5f22.0p-142L,
 161      0xf48000.0p-24,  0xbc42cb.0p-28, -0x172a21161a107674986dcdca6709c.0p-143L,
 162      0xf30000.0p-24,  0xd57797.0p-28, -0x1e09de07cb958897a3ea46e84abb3.0p-142L,
 163      0xf10000.0p-24,  0xf7518e.0p-28,  0x1ae1eec1b036c484993c549c4bf40.0p-151L,
 164      0xef0000.0p-24,  0x8cb9df.0p-27, -0x1d7355325d560d9e9ab3d6ebab580.0p-141L,
 165      0xed8000.0p-24,  0x999ec0.0p-27, -0x1f9f02d256d5037108f4ec21e48cd.0p-142L,
 166      0xec0000.0p-24,  0xa6988b.0p-27, -0x16fc0a9d12c17a70f7a684c596b12.0p-143L,
 167      0xea0000.0p-24,  0xb80698.0p-27,  0x15d581c1e8da99ded322fb08b8462.0p-141L,
 168      0xe80000.0p-24,  0xc99af3.0p-27, -0x1535b3ba8f150ae09996d7bb4653e.0p-143L,
 169      0xe70000.0p-24,  0xd273b2.0p-27,  0x163786f5251aefe0ded34c8318f52.0p-145L,
 170      0xe50000.0p-24,  0xe442c0.0p-27,  0x1bc4b2368e32d56699c1799a244d4.0p-144L,
 171      0xe38000.0p-24,  0xf1b83f.0p-27,  0x1c6090f684e6766abceccab1d7174.0p-141L,
 172      0xe20000.0p-24,  0xff448a.0p-27, -0x1890aa69ac9f4215f93936b709efb.0p-142L,
 173      0xe08000.0p-24,  0x8673f6.0p-26,  0x1b9985194b6affd511b534b72a28e.0p-140L,
 174      0xdf0000.0p-24,  0x8d515c.0p-26, -0x1dc08d61c6ef1d9b2ef7e68680598.0p-143L,
 175      0xdd8000.0p-24,  0x943a9e.0p-26, -0x1f72a2dac729b3f46662238a9425a.0p-142L,
 176      0xdc0000.0p-24,  0x9b2fe6.0p-26, -0x1fd4dfd3a0afb9691aed4d5e3df94.0p-140L,
 177      0xda8000.0p-24,  0xa2315d.0p-26, -0x11b26121629c46c186384993e1c93.0p-142L,
 178      0xd90000.0p-24,  0xa93f2f.0p-26,  0x1286d633e8e5697dc6a402a56fce1.0p-141L,
 179      0xd78000.0p-24,  0xb05988.0p-26,  0x16128eba9367707ebfa540e45350c.0p-144L,
 180      0xd60000.0p-24,  0xb78094.0p-26,  0x16ead577390d31ef0f4c9d43f79b2.0p-140L,
 181      0xd50000.0p-24,  0xbc4c6c.0p-26,  0x151131ccf7c7b75e7d900b521c48d.0p-141L,
 182      0xd38000.0p-24,  0xc3890a.0p-26, -0x115e2cd714bd06508aeb00d2ae3e9.0p-140L,
 183      0xd20000.0p-24,  0xcad2d7.0p-26, -0x1847f406ebd3af80485c2f409633c.0p-142L,
 184      0xd10000.0p-24,  0xcfb620.0p-26,  0x1c2259904d686581799fbce0b5f19.0p-141L,
 185      0xcf8000.0p-24,  0xd71653.0p-26,  0x1ece57a8d5ae54f550444ecf8b995.0p-140L,
 186      0xce0000.0p-24,  0xde843a.0p-26, -0x1f109d4bc4595412b5d2517aaac13.0p-141L,
 187      0xcd0000.0p-24,  0xe37fde.0p-26,  0x1bc03dc271a74d3a85b5b43c0e727.0p-141L,
 188      0xcb8000.0p-24,  0xeb050c.0p-26, -0x1bf2badc0df841a71b79dd5645b46.0p-145L,
 189      0xca0000.0p-24,  0xf29878.0p-26, -0x18efededd89fbe0bcfbe6d6db9f66.0p-147L,
 190      0xc90000.0p-24,  0xf7ad6f.0p-26,  0x1373ff977baa6911c7bafcb4d84fb.0p-141L,
 191      0xc80000.0p-24,  0xfcc8e3.0p-26,  0x196766f2fb328337cc050c6d83b22.0p-140L,
 192      0xc68000.0p-24,  0x823f30.0p-25,  0x19bd076f7c434e5fcf1a212e2a91e.0p-139L,
 193      0xc58000.0p-24,  0x84d52c.0p-25, -0x1a327257af0f465e5ecab5f2a6f81.0p-139L,
 194      0xc40000.0p-24,  0x88bc74.0p-25,  0x113f23def19c5a0fe396f40f1dda9.0p-141L,
 195      0xc30000.0p-24,  0x8b5ae6.0p-25,  0x1759f6e6b37de945a049a962e66c6.0p-139L,
 196      0xc20000.0p-24,  0x8dfccb.0p-25,  0x1ad35ca6ed5147bdb6ddcaf59c425.0p-141L,
 197      0xc10000.0p-24,  0x90a22b.0p-25,  0x1a1d71a87deba46bae9827221dc98.0p-139L,
 198      0xbf8000.0p-24,  0x94a0d8.0p-25, -0x139e5210c2b730e28aba001a9b5e0.0p-140L,
 199      0xbe8000.0p-24,  0x974f16.0p-25, -0x18f6ebcff3ed72e23e13431adc4a5.0p-141L,
 200      0xbd8000.0p-24,  0x9a00f1.0p-25, -0x1aa268be39aab7148e8d80caa10b7.0p-139L,
 201      0xbc8000.0p-24,  0x9cb672.0p-25, -0x14c8815839c5663663d15faed7771.0p-139L,
 202      0xbb0000.0p-24,  0xa0cda1.0p-25,  0x1eaf46390dbb2438273918db7df5c.0p-141L,
 203      0xba0000.0p-24,  0xa38c6e.0p-25,  0x138e20d831f698298adddd7f32686.0p-141L,
 204      0xb90000.0p-24,  0xa64f05.0p-25, -0x1e8d3c41123615b147a5d47bc208f.0p-142L,
 205      0xb80000.0p-24,  0xa91570.0p-25,  0x1ce28f5f3840b263acb4351104631.0p-140L,
 206      0xb70000.0p-24,  0xabdfbb.0p-25, -0x186e5c0a42423457e22d8c650b355.0p-139L,
 207      0xb60000.0p-24,  0xaeadef.0p-25, -0x14d41a0b2a08a465dc513b13f567d.0p-143L,
 208      0xb50000.0p-24,  0xb18018.0p-25,  0x16755892770633947ffe651e7352f.0p-139L,
 209      0xb40000.0p-24,  0xb45642.0p-25, -0x16395ebe59b15228bfe8798d10ff0.0p-142L,
 210      0xb30000.0p-24,  0xb73077.0p-25,  0x1abc65c8595f088b61a335f5b688c.0p-140L,
 211      0xb20000.0p-24,  0xba0ec4.0p-25, -0x1273089d3dad88e7d353e9967d548.0p-139L,
 212      0xb10000.0p-24,  0xbcf133.0p-25,  0x10f9f67b1f4bbf45de06ecebfaf6d.0p-139L,
 213      0xb00000.0p-24,  0xbfd7d2.0p-25, -0x109fab904864092b34edda19a831e.0p-140L,
 214      0xaf0000.0p-24,  0xc2c2ac.0p-25, -0x1124680aa43333221d8a9b475a6ba.0p-139L,
 215      0xae8000.0p-24,  0xc439b3.0p-25, -0x1f360cc4710fbfe24b633f4e8d84d.0p-140L,
 216      0xad8000.0p-24,  0xc72afd.0p-25, -0x132d91f21d89c89c45003fc5d7807.0p-140L,
 217      0xac8000.0p-24,  0xca20a2.0p-25, -0x16bf9b4d1f8da8002f2449e174504.0p-139L,
 218      0xab8000.0p-24,  0xcd1aae.0p-25,  0x19deb5ce6a6a8717d5626e16acc7d.0p-141L,
 219      0xaa8000.0p-24,  0xd0192f.0p-25,  0x1a29fb48f7d3ca87dabf351aa41f4.0p-139L,
 220      0xaa0000.0p-24,  0xd19a20.0p-25,  0x1127d3c6457f9d79f51dcc73014c9.0p-141L,
 221      0xa90000.0p-24,  0xd49f6a.0p-25, -0x1ba930e486a0ac42d1bf9199188e7.0p-141L,
 222      0xa80000.0p-24,  0xd7a94b.0p-25, -0x1b6e645f31549dd1160bcc45c7e2c.0p-139L,
 223      0xa70000.0p-24,  0xdab7d0.0p-25,  0x1118a425494b610665377f15625b6.0p-140L,
 224      0xa68000.0p-24,  0xdc40d5.0p-25,  0x1966f24d29d3a2d1b2176010478be.0p-140L,
 225      0xa58000.0p-24,  0xdf566d.0p-25, -0x1d8e52eb2248f0c95dd83626d7333.0p-142L,
 226      0xa48000.0p-24,  0xe270ce.0p-25, -0x1ee370f96e6b67ccb006a5b9890ea.0p-140L,
 227      0xa40000.0p-24,  0xe3ffce.0p-25,  0x1d155324911f56db28da4d629d00a.0p-140L,
 228      0xa30000.0p-24,  0xe72179.0p-25, -0x1fe6e2f2f867d8f4d60c713346641.0p-140L,
 229      0xa20000.0p-24,  0xea4812.0p-25,  0x1b7be9add7f4d3b3d406b6cbf3ce5.0p-140L,
 230      0xa18000.0p-24,  0xebdd3d.0p-25,  0x1b3cfb3f7511dd73692609040ccc2.0p-139L,
 231      0xa08000.0p-24,  0xef0b5b.0p-25, -0x1220de1f7301901b8ad85c25afd09.0p-139L,
 232      0xa00000.0p-24,  0xf0a451.0p-25, -0x176364c9ac81cc8a4dfb804de6867.0p-140L,
 233      0x9f0000.0p-24,  0xf3da16.0p-25,  0x1eed6b9aafac8d42f78d3e65d3727.0p-141L,
 234      0x9e8000.0p-24,  0xf576e9.0p-25,  0x1d593218675af269647b783d88999.0p-139L,
 235      0x9d8000.0p-24,  0xf8b47c.0p-25, -0x13e8eb7da053e063714615f7cc91d.0p-144L,
 236      0x9d0000.0p-24,  0xfa553f.0p-25,  0x1c063259bcade02951686d5373aec.0p-139L,
 237      0x9c0000.0p-24,  0xfd9ac5.0p-25,  0x1ef491085fa3c1649349630531502.0p-139L,
 238      0x9b8000.0p-24,  0xff3f8c.0p-25,  0x1d607a7c2b8c5320619fb9433d841.0p-139L,
 239      0x9a8000.0p-24,  0x814697.0p-24, -0x12ad3817004f3f0bdff99f932b273.0p-138L,
 240      0x9a0000.0p-24,  0x821b06.0p-24, -0x189fc53117f9e54e78103a2bc1767.0p-141L,
 241      0x990000.0p-24,  0x83c5f8.0p-24,  0x14cf15a048907b7d7f47ddb45c5a3.0p-139L,
 242      0x988000.0p-24,  0x849c7d.0p-24,  0x1cbb1d35fb82873b04a9af1dd692c.0p-138L,
 243      0x978000.0p-24,  0x864ba6.0p-24,  0x1128639b814f9b9770d8cb6573540.0p-138L,
 244      0x970000.0p-24,  0x87244c.0p-24,  0x184733853300f002e836dfd47bd41.0p-139L,
 245      0x968000.0p-24,  0x87fdaa.0p-24,  0x109d23aef77dd5cd7cc94306fb3ff.0p-140L,
 246      0x958000.0p-24,  0x89b293.0p-24, -0x1a81ef367a59de2b41eeebd550702.0p-138L,
 247      0x950000.0p-24,  0x8a8e20.0p-24, -0x121ad3dbb2f45275c917a30df4ac9.0p-138L,
 248      0x948000.0p-24,  0x8b6a6a.0p-24, -0x1cfb981628af71a89df4e6df2e93b.0p-139L,
 249      0x938000.0p-24,  0x8d253a.0p-24, -0x1d21730ea76cfdec367828734cae5.0p-139L,
 250      0x930000.0p-24,  0x8e03c2.0p-24,  0x135cc00e566f76b87333891e0dec4.0p-138L,
 251      0x928000.0p-24,  0x8ee30d.0p-24, -0x10fcb5df257a263e3bf446c6e3f69.0p-140L,
 252      0x918000.0p-24,  0x90a3ee.0p-24, -0x16e171b15433d723a4c7380a448d8.0p-139L,
 253      0x910000.0p-24,  0x918587.0p-24, -0x1d050da07f3236f330972da2a7a87.0p-139L,
 254      0x908000.0p-24,  0x9267e7.0p-24,  0x1be03669a5268d21148c6002becd3.0p-139L,
 255      0x8f8000.0p-24,  0x942f04.0p-24,  0x10b28e0e26c336af90e00533323ba.0p-139L,
 256      0x8f0000.0p-24,  0x9513c3.0p-24,  0x1a1d820da57cf2f105a89060046aa.0p-138L,
 257      0x8e8000.0p-24,  0x95f950.0p-24, -0x19ef8f13ae3cf162409d8ea99d4c0.0p-139L,
 258      0x8e0000.0p-24,  0x96dfab.0p-24, -0x109e417a6e507b9dc10dac743ad7a.0p-138L,
 259      0x8d0000.0p-24,  0x98aed2.0p-24,  0x10d01a2c5b0e97c4990b23d9ac1f5.0p-139L,
 260      0x8c8000.0p-24,  0x9997a2.0p-24, -0x1d6a50d4b61ea74540bdd2aa99a42.0p-138L,
 261      0x8c0000.0p-24,  0x9a8145.0p-24,  0x1b3b190b83f9527e6aba8f2d783c1.0p-138L,
 262      0x8b8000.0p-24,  0x9b6bbf.0p-24,  0x13a69fad7e7abe7ba81c664c107e0.0p-138L,
 263      0x8b0000.0p-24,  0x9c5711.0p-24, -0x11cd12316f576aad348ae79867223.0p-138L,
 264      0x8a8000.0p-24,  0x9d433b.0p-24,  0x1c95c444b807a246726b304ccae56.0p-139L,
 265      0x898000.0p-24,  0x9f1e22.0p-24, -0x1b9c224ea698c2f9b47466d6123fe.0p-139L,
 266      0x890000.0p-24,  0xa00ce1.0p-24,  0x125ca93186cf0f38b4619a2483399.0p-141L,
 267      0x888000.0p-24,  0xa0fc80.0p-24, -0x1ee38a7bc228b3597043be78eaf49.0p-139L,
 268      0x880000.0p-24,  0xa1ed00.0p-24, -0x1a0db876613d204147dc69a07a649.0p-138L,
 269      0x878000.0p-24,  0xa2de62.0p-24,  0x193224e8516c008d3602a7b41c6e8.0p-139L,
 270      0x870000.0p-24,  0xa3d0a9.0p-24,  0x1fa28b4d2541aca7d5844606b2421.0p-139L,
 271      0x868000.0p-24,  0xa4c3d6.0p-24,  0x1c1b5760fb4571acbcfb03f16daf4.0p-138L,
 272      0x858000.0p-24,  0xa6acea.0p-24,  0x1fed5d0f65949c0a345ad743ae1ae.0p-140L,
 273      0x850000.0p-24,  0xa7a2d4.0p-24,  0x1ad270c9d749362382a7688479e24.0p-140L,
 274      0x848000.0p-24,  0xa899ab.0p-24,  0x199ff15ce532661ea9643a3a2d378.0p-139L,
 275      0x840000.0p-24,  0xa99171.0p-24,  0x1a19e15ccc45d257530a682b80490.0p-139L,
 276      0x838000.0p-24,  0xaa8a28.0p-24, -0x121a14ec532b35ba3e1f868fd0b5e.0p-140L,
 277      0x830000.0p-24,  0xab83d1.0p-24,  0x1aee319980bff3303dd481779df69.0p-139L,
 278      0x828000.0p-24,  0xac7e6f.0p-24, -0x18ffd9e3900345a85d2d86161742e.0p-140L,
 279      0x820000.0p-24,  0xad7a03.0p-24, -0x1e4db102ce29f79b026b64b42caa1.0p-140L,
 280      0x818000.0p-24,  0xae768f.0p-24,  0x17c35c55a04a82ab19f77652d977a.0p-141L,
 281      0x810000.0p-24,  0xaf7415.0p-24,  0x1448324047019b48d7b98c1cf7234.0p-138L,
 282      0x808000.0p-24,  0xb07298.0p-24, -0x1750ee3915a197e9c7359dd94152f.0p-138L,
 283      0x800000.0p-24,  0xb17218.0p-24, -0x105c610ca86c3898cff81a12a17e2.0p-141L,
 284 };
 285
 286 #ifdef USE_UTAB
 287 static const struct {
 288         float   H;                      /* 1 + i/INTERVALS (exact) */
 289         float   E;                      /* H(i) * G(i) - 1 (exact) */
 290 } U[TSIZE] = {
 291          0x800000.0p-23,  0,
 292          0x810000.0p-23, -0x800000.0p-37,
 293          0x820000.0p-23, -0x800000.0p-35,
 294          0x830000.0p-23, -0x900000.0p-34,
 295          0x840000.0p-23, -0x800000.0p-33,
 296          0x850000.0p-23, -0xc80000.0p-33,
 297          0x860000.0p-23, -0xa00000.0p-36,
 298          0x870000.0p-23,  0x940000.0p-33,
 299          0x880000.0p-23,  0x800000.0p-35,
 300          0x890000.0p-23, -0xc80000.0p-34,
 301          0x8a0000.0p-23,  0xe00000.0p-36,
 302          0x8b0000.0p-23,  0x900000.0p-33,
 303          0x8c0000.0p-23, -0x800000.0p-35,
 304          0x8d0000.0p-23, -0xe00000.0p-33,
 305          0x8e0000.0p-23,  0x880000.0p-33,
 306          0x8f0000.0p-23, -0xa80000.0p-34,
 307          0x900000.0p-23, -0x800000.0p-35,
 308          0x910000.0p-23,  0x800000.0p-37,
 309          0x920000.0p-23,  0x900000.0p-35,
 310          0x930000.0p-23,  0xd00000.0p-35,
 311          0x940000.0p-23,  0xe00000.0p-35,
 312          0x950000.0p-23,  0xc00000.0p-35,
 313          0x960000.0p-23,  0xe00000.0p-36,
 314          0x970000.0p-23, -0x800000.0p-38,
 315          0x980000.0p-23, -0xc00000.0p-35,
 316          0x990000.0p-23, -0xd00000.0p-34,
 317          0x9a0000.0p-23,  0x880000.0p-33,
 318          0x9b0000.0p-23,  0xe80000.0p-35,
 319          0x9c0000.0p-23, -0x800000.0p-35,
 320          0x9d0000.0p-23,  0xb40000.0p-33,
 321          0x9e0000.0p-23,  0x880000.0p-34,
 322          0x9f0000.0p-23, -0xe00000.0p-35,
 323          0xa00000.0p-23,  0x800000.0p-33,
 324          0xa10000.0p-23, -0x900000.0p-36,
 325          0xa20000.0p-23, -0xb00000.0p-33,
 326          0xa30000.0p-23, -0xa00000.0p-36,
 327          0xa40000.0p-23,  0x800000.0p-33,
 328          0xa50000.0p-23, -0xf80000.0p-35,
 329          0xa60000.0p-23,  0x880000.0p-34,
 330          0xa70000.0p-23, -0x900000.0p-33,
 331          0xa80000.0p-23, -0x800000.0p-35,
 332          0xa90000.0p-23,  0x900000.0p-34,
 333          0xaa0000.0p-23,  0xa80000.0p-33,
 334          0xab0000.0p-23, -0xac0000.0p-34,
 335          0xac0000.0p-23, -0x800000.0p-37,
 336          0xad0000.0p-23,  0xf80000.0p-35,
 337          0xae0000.0p-23,  0xf80000.0p-34,
 338          0xaf0000.0p-23, -0xac0000.0p-33,
 339          0xb00000.0p-23, -0x800000.0p-33,
 340          0xb10000.0p-23, -0xb80000.0p-34,
 341          0xb20000.0p-23, -0x800000.0p-34,
 342          0xb30000.0p-23, -0xb00000.0p-35,
 343          0xb40000.0p-23, -0x800000.0p-35,
 344          0xb50000.0p-23, -0xe00000.0p-36,
 345          0xb60000.0p-23, -0x800000.0p-35,
 346          0xb70000.0p-23, -0xb00000.0p-35,
 347          0xb80000.0p-23, -0x800000.0p-34,
 348          0xb90000.0p-23, -0xb80000.0p-34,
 349          0xba0000.0p-23, -0x800000.0p-33,
 350          0xbb0000.0p-23, -0xac0000.0p-33,
 351          0xbc0000.0p-23,  0x980000.0p-33,
 352          0xbd0000.0p-23,  0xbc0000.0p-34,
 353          0xbe0000.0p-23,  0xe00000.0p-36,
 354          0xbf0000.0p-23, -0xb80000.0p-35,
 355          0xc00000.0p-23, -0x800000.0p-33,
 356          0xc10000.0p-23,  0xa80000.0p-33,
 357          0xc20000.0p-23,  0x900000.0p-34,
 358          0xc30000.0p-23, -0x800000.0p-35,
 359          0xc40000.0p-23, -0x900000.0p-33,
 360          0xc50000.0p-23,  0x820000.0p-33,
 361          0xc60000.0p-23,  0x800000.0p-38,
 362          0xc70000.0p-23, -0x820000.0p-33,
 363          0xc80000.0p-23,  0x800000.0p-33,
 364          0xc90000.0p-23, -0xa00000.0p-36,
 365          0xca0000.0p-23, -0xb00000.0p-33,
 366          0xcb0000.0p-23,  0x840000.0p-34,
 367          0xcc0000.0p-23, -0xd00000.0p-34,
 368          0xcd0000.0p-23,  0x800000.0p-33,
 369          0xce0000.0p-23, -0xe00000.0p-35,
 370          0xcf0000.0p-23,  0xa60000.0p-33,
 371          0xd00000.0p-23, -0x800000.0p-35,
 372          0xd10000.0p-23,  0xb40000.0p-33,
 373          0xd20000.0p-23, -0x800000.0p-35,
 374          0xd30000.0p-23,  0xaa0000.0p-33,
 375          0xd40000.0p-23, -0xe00000.0p-35,
 376          0xd50000.0p-23,  0x880000.0p-33,
 377          0xd60000.0p-23, -0xd00000.0p-34,
 378          0xd70000.0p-23,  0x9c0000.0p-34,
 379          0xd80000.0p-23, -0xb00000.0p-33,
 380          0xd90000.0p-23, -0x800000.0p-38,
 381          0xda0000.0p-23,  0xa40000.0p-33,
 382          0xdb0000.0p-23, -0xdc0000.0p-34,
 383          0xdc0000.0p-23,  0xc00000.0p-35,
 384          0xdd0000.0p-23,  0xca0000.0p-33,
 385          0xde0000.0p-23, -0xb80000.0p-34,
 386          0xdf0000.0p-23,  0xd00000.0p-35,
 387          0xe00000.0p-23,  0xc00000.0p-33,
 388          0xe10000.0p-23, -0xf40000.0p-34,
 389          0xe20000.0p-23,  0x800000.0p-37,
 390          0xe30000.0p-23,  0x860000.0p-33,
 391          0xe40000.0p-23, -0xc80000.0p-33,
 392          0xe50000.0p-23, -0xa80000.0p-34,
 393          0xe60000.0p-23,  0xe00000.0p-36,
 394          0xe70000.0p-23,  0x880000.0p-33,
 395          0xe80000.0p-23, -0xe00000.0p-33,
 396          0xe90000.0p-23, -0xfc0000.0p-34,
 397          0xea0000.0p-23, -0x800000.0p-35,
 398          0xeb0000.0p-23,  0xe80000.0p-35,
 399          0xec0000.0p-23,  0x900000.0p-33,
 400          0xed0000.0p-23,  0xe20000.0p-33,
 401          0xee0000.0p-23, -0xac0000.0p-33,
 402          0xef0000.0p-23, -0xc80000.0p-34,
 403          0xf00000.0p-23, -0x800000.0p-35,
 404          0xf10000.0p-23,  0x800000.0p-35,
 405          0xf20000.0p-23,  0xb80000.0p-34,
 406          0xf30000.0p-23,  0x940000.0p-33,
 407          0xf40000.0p-23,  0xc80000.0p-33,
 408          0xf50000.0p-23, -0xf20000.0p-33,
 409          0xf60000.0p-23, -0xc80000.0p-33,
 410          0xf70000.0p-23, -0xa20000.0p-33,
 411          0xf80000.0p-23, -0x800000.0p-33,
 412          0xf90000.0p-23, -0xc40000.0p-34,
 413          0xfa0000.0p-23, -0x900000.0p-34,
 414          0xfb0000.0p-23, -0xc80000.0p-35,
 415          0xfc0000.0p-23, -0x800000.0p-35,
 416          0xfd0000.0p-23, -0x900000.0p-36,
 417          0xfe0000.0p-23, -0x800000.0p-37,
 418          0xff0000.0p-23, -0x800000.0p-39,
 419          0x800000.0p-22,  0,
 420 };
 421 #endif /* USE_UTAB */
 422
 423 #ifdef STRUCT_RETURN
 424 #define RETURN1(rp, v) do {     \
 425         (rp)->hi = (v);         \
 426         (rp)->lo_set = 0;       \
 427         return;                 \
 428 } while (0)
 429
 430 #define RETURN2(rp, h, l) do {  \
 431         (rp)->hi = (h);         \
 432         (rp)->lo = (l);         \
 433         (rp)->lo_set = 1;       \
 434         return;                 \
 435 } while (0)
 436
 437 struct ld {
 438         long double hi;
 439         long double lo;
 440         int     lo_set;
 441 };
 442 #else
 443 #define RETURN1(rp, v)  RETURNF(v)
 444 #define RETURN2(rp, h, l)       RETURNI((h) + (l))
 445 #endif
 446
 447 #ifdef STRUCT_RETURN
 448 static inline __always_inline void
 449 k_logl(long double x, struct ld *rp)
 450 #else
 451 long double
 452 logl(long double x)
 453 #endif
 454 {
 455         long double d, val_hi, val_lo;
 456         double dd, dk;
 457         uint64_t lx, llx;
 458         int i, k;
 459         uint16_t hx;
 460
 461         EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
 462         k = -16383;
 463 #if 0 /* Hard to do efficiently.  Don't do it until we support all modes. */
 464         if (x == 1)
 465                 RETURN1(rp, 0);         /* log(1) = +0 in all rounding modes */
 466 #endif
 467         if (hx == 0 || hx >= 0x8000) {  /* zero, negative or subnormal? */
 468                 if (((hx & 0x7fff) | lx | llx) == 0)
 469                         RETURN1(rp, -1 / zero); /* log(+-0) = -Inf */
 470                 if (hx != 0)
 471                         /* log(neg or NaN) = qNaN: */
 472                         RETURN1(rp, (x - x) / zero);
 473                 x *= 0x1.0p113;         /* subnormal; scale up x */
 474                 EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
 475                 k = -16383 - 113;
 476         } else if (hx >= 0x7fff)
 477                 RETURN1(rp, x + x);     /* log(Inf or NaN) = Inf or qNaN */
 478 #ifndef STRUCT_RETURN
 479         ENTERI();
 480 #endif
 481         k += hx;
 482         dk = k;
 483
 484         /* Scale x to be in [1, 2). */
 485         SET_LDBL_EXPSIGN(x, 0x3fff);
 486
 487         /* 0 <= i <= INTERVALS: */
 488 #define L2I     (49 - LOG2_INTERVALS)
 489         i = (lx + (1LL << (L2I - 2))) >> (L2I - 1);
 490
 491         /*
 492          * -0.005280 < d < 0.004838.  In particular, the infinite-
 493          * precision |d| is <= 2**-7.  Rounding of G(i) to 8 bits
 494          * ensures that d is representable without extra precision for
 495          * this bound on |d| (since when this calculation is expressed
 496          * as x*G(i)-1, the multiplication needs as many extra bits as
 497          * G(i) has and the subtraction cancels 8 bits).  But for
 498          * most i (107 cases out of 129), the infinite-precision |d|
 499          * is <= 2**-8.  G(i) is rounded to 9 bits for such i to give
 500          * better accuracy (this works by improving the bound on |d|,
 501          * which in turn allows rounding to 9 bits in more cases).
 502          * This is only important when the original x is near 1 -- it
 503          * lets us avoid using a special method to give the desired
 504          * accuracy for such x.
 505          */
 506         if (0)
 507                 d = x * G(i) - 1;
 508         else {
 509 #ifdef USE_UTAB
 510                 d = (x - H(i)) * G(i) + E(i);
 511 #else
 512                 long double x_hi;
 513                 double x_lo;
 514
 515                 /*
 516                  * Split x into x_hi + x_lo to calculate x*G(i)-1 exactly.
 517                  * G(i) has at most 9 bits, so the splitting point is not
 518                  * critical.
 519                  */
 520                 INSERT_LDBL128_WORDS(x_hi, 0x3fff, lx,
 521                     llx & 0xffffffffff000000ULL);
 522                 x_lo = x - x_hi;
 523                 d = x_hi * G(i) - 1 + x_lo * G(i);
 524 #endif
 525         }
 526
 527         /*
 528          * Our algorithm depends on exact cancellation of F_lo(i) and
 529          * F_hi(i) with dk*ln_2_lo and dk*ln2_hi when k is -1 and i is
 530          * at the end of the table.  This and other technical complications
 531          * make it difficult to avoid the double scaling in (dk*ln2) *
 532          * log(base) for base != e without losing more accuracy and/or
 533          * efficiency than is gained.
 534          */
 535         /*
 536          * Use double precision operations wherever possible, since long
 537          * double operations are emulated and are very slow on the only
 538          * known machines that support ld128 (sparc64).  Also, don't try
 539          * to improve parallelism by increasing the number of operations,
 540          * since any parallelism on such machines is needed for the
 541          * emulation.  Horner's method is good for this, and is also good
 542          * for accuracy.  Horner's method doesn't handle the `lo' term
 543          * well, either for efficiency or accuracy.  However, for accuracy
 544          * we evaluate d * d * P2 separately to take advantage of
 545          * by P2 being exact, and this gives a good place to sum the 'lo'
 546          * term too.
 547          */
 548         dd = (double)d;
 549         val_lo = d * d * d * (P3 +
 550             d * (P4 + d * (P5 + d * (P6 + d * (P7 + d * (P8 +
 551             dd * (P9 + dd * (P10 + dd * (P11 + dd * (P12 + dd * (P13 +
 552             dd * P14))))))))))) + (F_lo(i) + dk * ln2_lo) + d * d * P2;
 553         val_hi = d;
 554 #ifdef DEBUG
 555         if (fetestexcept(FE_UNDERFLOW))
 556                 breakpoint();
 557 #endif
 558
 559         _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
 560         RETURN2(rp, val_hi, val_lo);
 561 }
 562
 563 long double
 564 log1pl(long double x)
 565 {
 566         long double d, d_hi, f_lo, val_hi, val_lo;
 567         long double f_hi, twopminusk;
 568         double d_lo, dd, dk;
 569         uint64_t lx, llx;
 570         int i, k;
 571         int16_t ax, hx;
 572
 573         DOPRINT_START(&x);
 574         EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
 575         if (hx < 0x3fff) {              /* x < 1, or x neg NaN */
 576                 ax = hx & 0x7fff;
 577                 if (ax >= 0x3fff) {     /* x <= -1, or x neg NaN */
 578                         if (ax == 0x3fff && (lx | llx) == 0)
 579                                 RETURNP(-1 / zero);     /* log1p(-1) = -Inf */
 580                         /* log1p(x < 1, or x NaN) = qNaN: */
 581                         RETURNP((x - x) / (x - x));
 582                 }
 583                 if (ax <= 0x3f8d) {     /* |x| < 2**-113 */
 584                         if ((int)x == 0)
 585                                 RETURNP(x);     /* x with inexact if x != 0 */
 586                 }
 587                 f_hi = 1;
 588                 f_lo = x;
 589         } else if (hx >= 0x7fff) {      /* x +Inf or non-neg NaN */
 590                 RETURNP(x + x);         /* log1p(Inf or NaN) = Inf or qNaN */
 591         } else if (hx < 0x40e1) {       /* 1 <= x < 2**226 */
 592                 f_hi = x;
 593                 f_lo = 1;
 594         } else {                        /* 2**226 <= x < +Inf */
 595                 f_hi = x;
 596                 f_lo = 0;               /* avoid underflow of the P3 term */
 597         }
 598         ENTERI();
 599         x = f_hi + f_lo;
 600         f_lo = (f_hi - x) + f_lo;
 601
 602         EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
 603         k = -16383;
 604
 605         k += hx;
 606         dk = k;
 607
 608         SET_LDBL_EXPSIGN(x, 0x3fff);
 609         twopminusk = 1;
 610         SET_LDBL_EXPSIGN(twopminusk, 0x7ffe - (hx & 0x7fff));
 611         f_lo *= twopminusk;
 612
 613         i = (lx + (1LL << (L2I - 2))) >> (L2I - 1);
 614
 615         /*
 616          * x*G(i)-1 (with a reduced x) can be represented exactly, as
 617          * above, but now we need to evaluate the polynomial on d =
 618          * (x+f_lo)*G(i)-1 and extra precision is needed for that.
 619          * Since x+x_lo is a hi+lo decomposition and subtracting 1
 620          * doesn't lose too many bits, an inexact calculation for
 621          * f_lo*G(i) is good enough.
 622          */
 623         if (0)
 624                 d_hi = x * G(i) - 1;
 625         else {
 626 #ifdef USE_UTAB
 627                 d_hi = (x - H(i)) * G(i) + E(i);
 628 #else
 629                 long double x_hi;
 630                 double x_lo;
 631
 632                 INSERT_LDBL128_WORDS(x_hi, 0x3fff, lx,
 633                     llx & 0xffffffffff000000ULL);
 634                 x_lo = x - x_hi;
 635                 d_hi = x_hi * G(i) - 1 + x_lo * G(i);
 636 #endif
 637         }
 638         d_lo = f_lo * G(i);
 639
 640         /*
 641          * This is _2sumF(d_hi, d_lo) inlined.  The condition
 642          * (d_hi == 0 || |d_hi| >= |d_lo|) for using _2sumF() is not
 643          * always satisifed, so it is not clear that this works, but
 644          * it works in practice.  It works even if it gives a wrong
 645          * normalized d_lo, since |d_lo| > |d_hi| implies that i is
 646          * nonzero and d is tiny, so the F(i) term dominates d_lo.
 647          * In float precision:
 648          * (By exhaustive testing, the worst case is d_hi = 0x1.bp-25.
 649          * And if d is only a little tinier than that, we would have
 650          * another underflow problem for the P3 term; this is also ruled
 651          * out by exhaustive testing.)
 652          */
 653         d = d_hi + d_lo;
 654         d_lo = d_hi - d + d_lo;
 655         d_hi = d;
 656
 657         dd = (double)d;
 658         val_lo = d * d * d * (P3 +
 659             d * (P4 + d * (P5 + d * (P6 + d * (P7 + d * (P8 +
 660             dd * (P9 + dd * (P10 + dd * (P11 + dd * (P12 + dd * (P13 +
 661             dd * P14))))))))))) + (F_lo(i) + dk * ln2_lo + d_lo) + d * d * P2;
 662         val_hi = d_hi;
 663 #ifdef DEBUG
 664         if (fetestexcept(FE_UNDERFLOW))
 665                 breakpoint();
 666 #endif
 667
 668         _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
 669         RETURN2PI(val_hi, val_lo);
 670 }
 671
 672 #ifdef STRUCT_RETURN
 673
 674 long double
 675 logl(long double x)
 676 {
 677         struct ld r;
 678
 679         ENTERI();
 680         DOPRINT_START(&x);
 681         k_logl(x, &r);
 682         RETURNSPI(&r);
 683 }
 684
 685 /*
 686  * 29+113 bit decompositions.  The bits are distributed so that the products
 687  * of the hi terms are exact in double precision.  The types are chosen so
 688  * that the products of the hi terms are done in at least double precision,
 689  * without any explicit conversions.  More natural choices would require a
 690  * slow long double precision multiplication.
 691  */
 692 static const double
 693 invln10_hi =  4.3429448176175356e-1,            /*  0x1bcb7b15000000.0p-54 */
 694 invln2_hi =  1.4426950402557850e0;              /*  0x17154765000000.0p-52 */
 695 static const long double
 696 invln10_lo =  1.41498268538580090791605082294397000e-10L,       /*  0x137287195355baaafad33dc323ee3.0p-145L */
 697 invln2_lo =  6.33178418956604368501892137426645911e-10L;        /*  0x15c17f0bbbe87fed0691d3e88eb57.0p-143L */
 698
 699 long double
 700 log10l(long double x)
 701 {
 702         struct ld r;
 703         long double lo;
 704         float hi;
 705
 706         ENTERI();
 707         DOPRINT_START(&x);
 708         k_logl(x, &r);
 709         if (!r.lo_set)
 710                 RETURNPI(r.hi);
 711         _2sumF(r.hi, r.lo);
 712         hi = r.hi;
 713         lo = r.lo + (r.hi - hi);
 714         RETURN2PI(invln10_hi * hi,
 715             (invln10_lo + invln10_hi) * lo + invln10_lo * hi);
 716 }
 717
 718 long double
 719 log2l(long double x)
 720 {
 721         struct ld r;
 722         long double lo;
 723         float hi;
 724
 725         ENTERI();
 726         DOPRINT_START(&x);
 727         k_logl(x, &r);
 728         if (!r.lo_set)
 729                 RETURNPI(r.hi);
 730         _2sumF(r.hi, r.lo);
 731         hi = r.hi;
 732         lo = r.lo + (r.hi - hi);
 733         RETURN2PI(invln2_hi * hi,
 734             (invln2_lo + invln2_hi) * lo + invln2_lo * hi);
 735 }
 736
 737 #endif /* STRUCT_RETURN */