1251292Sdas/*- 2251292Sdas * Copyright (c) 2007-2013 Bruce D. Evans 3251292Sdas * All rights reserved. 4251292Sdas * 5251292Sdas * Redistribution and use in source and binary forms, with or without 6251292Sdas * modification, are permitted provided that the following conditions 7251292Sdas * are met: 8251292Sdas * 1. Redistributions of source code must retain the above copyright 9251292Sdas * notice unmodified, this list of conditions, and the following 10251292Sdas * disclaimer. 11251292Sdas * 2. Redistributions in binary form must reproduce the above copyright 12251292Sdas * notice, this list of conditions and the following disclaimer in the 13251292Sdas * documentation and/or other materials provided with the distribution. 14251292Sdas * 15251292Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 16251292Sdas * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 17251292Sdas * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 18251292Sdas * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 19251292Sdas * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 20251292Sdas * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 21251292Sdas * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 22251292Sdas * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23251292Sdas * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24251292Sdas * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25251292Sdas */ 26251292Sdas 27251292Sdas#include <sys/cdefs.h> 28251292Sdas__FBSDID("$FreeBSD$"); 29251292Sdas 30251292Sdas/** 31251292Sdas * Implementation of the natural logarithm of x for 128-bit format. 32251292Sdas * 33251292Sdas * First decompose x into its base 2 representation: 34251292Sdas * 35251292Sdas * log(x) = log(X * 2**k), where X is in [1, 2) 36251292Sdas * = log(X) + k * log(2). 37251292Sdas * 38251292Sdas * Let X = X_i + e, where X_i is the center of one of the intervals 39251292Sdas * [-1.0/256, 1.0/256), [1.0/256, 3.0/256), .... [2.0-1.0/256, 2.0+1.0/256) 40251292Sdas * and X is in this interval. Then 41251292Sdas * 42251292Sdas * log(X) = log(X_i + e) 43251292Sdas * = log(X_i * (1 + e / X_i)) 44251292Sdas * = log(X_i) + log(1 + e / X_i). 45251292Sdas * 46251292Sdas * The values log(X_i) are tabulated below. Let d = e / X_i and use 47251292Sdas * 48251292Sdas * log(1 + d) = p(d) 49251292Sdas * 50251292Sdas * where p(d) = d - 0.5*d*d + ... is a special minimax polynomial of 51251292Sdas * suitably high degree. 52251292Sdas * 53251292Sdas * To get sufficiently small roundoff errors, k * log(2), log(X_i), and 54251292Sdas * sometimes (if |k| is not large) the first term in p(d) must be evaluated 55251292Sdas * and added up in extra precision. Extra precision is not needed for the 56251292Sdas * rest of p(d). In the worst case when k = 0 and log(X_i) is 0, the final 57251292Sdas * error is controlled mainly by the error in the second term in p(d). The 58251292Sdas * error in this term itself is at most 0.5 ulps from the d*d operation in 59251292Sdas * it. The error in this term relative to the first term is thus at most 60251292Sdas * 0.5 * |-0.5| * |d| < 1.0/1024 ulps. We aim for an accumulated error of 61251292Sdas * at most twice this at the point of the final rounding step. Thus the 62251292Sdas * final error should be at most 0.5 + 1.0/512 = 0.5020 ulps. Exhaustive 63251292Sdas * testing of a float variant of this function showed a maximum final error 64251292Sdas * of 0.5008 ulps. Non-exhaustive testing of a double variant of this 65251292Sdas * function showed a maximum final error of 0.5078 ulps (near 1+1.0/256). 66251292Sdas * 67251292Sdas * We made the maximum of |d| (and thus the total relative error and the 68251292Sdas * degree of p(d)) small by using a large number of intervals. Using 69251292Sdas * centers of intervals instead of endpoints reduces this maximum by a 70251292Sdas * factor of 2 for a given number of intervals. p(d) is special only 71251292Sdas * in beginning with the Taylor coefficients 0 + 1*d, which tends to happen 72251292Sdas * naturally. The most accurate minimax polynomial of a given degree might 73251292Sdas * be different, but then we wouldn't want it since we would have to do 74251292Sdas * extra work to avoid roundoff error (especially for P0*d instead of d). 75251292Sdas */ 76251292Sdas 77251292Sdas#ifdef DEBUG 78251292Sdas#include <assert.h> 79251292Sdas#include <fenv.h> 80251292Sdas#endif 81251292Sdas 82251292Sdas#include "fpmath.h" 83251292Sdas#include "math.h" 84251292Sdas#ifndef NO_STRUCT_RETURN 85251292Sdas#define STRUCT_RETURN 86251292Sdas#endif 87251292Sdas#include "math_private.h" 88251292Sdas 89251292Sdas#if !defined(NO_UTAB) && !defined(NO_UTABL) 90251292Sdas#define USE_UTAB 91251292Sdas#endif 92251292Sdas 93251292Sdas/* 94251292Sdas * Domain [-0.005280, 0.004838], range ~[-1.1577e-37, 1.1582e-37]: 95251292Sdas * |log(1 + d)/d - p(d)| < 2**-122.7 96251292Sdas */ 97251292Sdasstatic const long double 98251292SdasP2 = -0.5L, 99251292SdasP3 = 3.33333333333333333333333333333233795e-1L, /* 0x15555555555555555555555554d42.0p-114L */ 100251292SdasP4 = -2.49999999999999999999999999941139296e-1L, /* -0x1ffffffffffffffffffffffdab14e.0p-115L */ 101251292SdasP5 = 2.00000000000000000000000085468039943e-1L, /* 0x19999999999999999999a6d3567f4.0p-115L */ 102251292SdasP6 = -1.66666666666666666666696142372698408e-1L, /* -0x15555555555555555567267a58e13.0p-115L */ 103251292SdasP7 = 1.42857142857142857119522943477166120e-1L, /* 0x1249249249249248ed79a0ae434de.0p-115L */ 104251292SdasP8 = -1.24999999999999994863289015033581301e-1L; /* -0x1fffffffffffffa13e91765e46140.0p-116L */ 105251292Sdas/* Double precision gives ~ 53 + log2(P9 * max(|d|)**8) ~= 120 bits. */ 106251292Sdasstatic const double 107251292SdasP9 = 1.1111111111111401e-1, /* 0x1c71c71c71c7ed.0p-56 */ 108251292SdasP10 = -1.0000000000040135e-1, /* -0x199999999a0a92.0p-56 */ 109251292SdasP11 = 9.0909090728136258e-2, /* 0x1745d173962111.0p-56 */ 110251292SdasP12 = -8.3333318851855284e-2, /* -0x1555551722c7a3.0p-56 */ 111251292SdasP13 = 7.6928634666404178e-2, /* 0x13b1985204a4ae.0p-56 */ 112251292SdasP14 = -7.1626810078462499e-2; /* -0x12562276cdc5d0.0p-56 */ 113251292Sdas 114251292Sdasstatic volatile const double zero = 0; 115251292Sdas 116251292Sdas#define INTERVALS 128 117251292Sdas#define LOG2_INTERVALS 7 118251292Sdas#define TSIZE (INTERVALS + 1) 119251292Sdas#define G(i) (T[(i)].G) 120251292Sdas#define F_hi(i) (T[(i)].F_hi) 121251292Sdas#define F_lo(i) (T[(i)].F_lo) 122251292Sdas#define ln2_hi F_hi(TSIZE - 1) 123251292Sdas#define ln2_lo F_lo(TSIZE - 1) 124251292Sdas#define E(i) (U[(i)].E) 125251292Sdas#define H(i) (U[(i)].H) 126251292Sdas 127251292Sdasstatic const struct { 128251292Sdas float G; /* 1/(1 + i/128) rounded to 8/9 bits */ 129251292Sdas float F_hi; /* log(1 / G_i) rounded (see below) */ 130251292Sdas /* The compiler will insert 8 bytes of padding here. */ 131251292Sdas long double F_lo; /* next 113 bits for log(1 / G_i) */ 132251292Sdas} T[TSIZE] = { 133251292Sdas /* 134251292Sdas * ln2_hi and each F_hi(i) are rounded to a number of bits that 135251292Sdas * makes F_hi(i) + dk*ln2_hi exact for all i and all dk. 136251292Sdas * 137251292Sdas * The last entry (for X just below 2) is used to define ln2_hi 138251292Sdas * and ln2_lo, to ensure that F_hi(i) and F_lo(i) cancel exactly 139251292Sdas * with dk*ln2_hi and dk*ln2_lo, respectively, when dk = -1. 140251292Sdas * This is needed for accuracy when x is just below 1. (To avoid 141251292Sdas * special cases, such x are "reduced" strangely to X just below 142251292Sdas * 2 and dk = -1, and then the exact cancellation is needed 143251292Sdas * because any the error from any non-exactness would be too 144251292Sdas * large). 145251292Sdas * 146251292Sdas * The relevant range of dk is [-16445, 16383]. The maximum number 147251292Sdas * of bits in F_hi(i) that works is very dependent on i but has 148251292Sdas * a minimum of 93. We only need about 12 bits in F_hi(i) for 149251292Sdas * it to provide enough extra precision. 150251292Sdas * 151251292Sdas * We round F_hi(i) to 24 bits so that it can have type float, 152251292Sdas * mainly to minimize the size of the table. Using all 24 bits 153251292Sdas * in a float for it automatically satisfies the above constraints. 154251292Sdas */ 155251292Sdas 0x800000.0p-23, 0, 0, 156251292Sdas 0xfe0000.0p-24, 0x8080ac.0p-30, -0x14ee431dae6674afa0c4bfe16e8fd.0p-144L, 157251292Sdas 0xfc0000.0p-24, 0x8102b3.0p-29, -0x1db29ee2d83717be918e1119642ab.0p-144L, 158251292Sdas 0xfa0000.0p-24, 0xc24929.0p-29, 0x1191957d173697cf302cc9476f561.0p-143L, 159251292Sdas 0xf80000.0p-24, 0x820aec.0p-28, 0x13ce8888e02e78eba9b1113bc1c18.0p-142L, 160251292Sdas 0xf60000.0p-24, 0xa33577.0p-28, -0x17a4382ce6eb7bfa509bec8da5f22.0p-142L, 161251292Sdas 0xf48000.0p-24, 0xbc42cb.0p-28, -0x172a21161a107674986dcdca6709c.0p-143L, 162251292Sdas 0xf30000.0p-24, 0xd57797.0p-28, -0x1e09de07cb958897a3ea46e84abb3.0p-142L, 163251292Sdas 0xf10000.0p-24, 0xf7518e.0p-28, 0x1ae1eec1b036c484993c549c4bf40.0p-151L, 164251292Sdas 0xef0000.0p-24, 0x8cb9df.0p-27, -0x1d7355325d560d9e9ab3d6ebab580.0p-141L, 165251292Sdas 0xed8000.0p-24, 0x999ec0.0p-27, -0x1f9f02d256d5037108f4ec21e48cd.0p-142L, 166251292Sdas 0xec0000.0p-24, 0xa6988b.0p-27, -0x16fc0a9d12c17a70f7a684c596b12.0p-143L, 167251292Sdas 0xea0000.0p-24, 0xb80698.0p-27, 0x15d581c1e8da99ded322fb08b8462.0p-141L, 168251292Sdas 0xe80000.0p-24, 0xc99af3.0p-27, -0x1535b3ba8f150ae09996d7bb4653e.0p-143L, 169251292Sdas 0xe70000.0p-24, 0xd273b2.0p-27, 0x163786f5251aefe0ded34c8318f52.0p-145L, 170251292Sdas 0xe50000.0p-24, 0xe442c0.0p-27, 0x1bc4b2368e32d56699c1799a244d4.0p-144L, 171251292Sdas 0xe38000.0p-24, 0xf1b83f.0p-27, 0x1c6090f684e6766abceccab1d7174.0p-141L, 172251292Sdas 0xe20000.0p-24, 0xff448a.0p-27, -0x1890aa69ac9f4215f93936b709efb.0p-142L, 173251292Sdas 0xe08000.0p-24, 0x8673f6.0p-26, 0x1b9985194b6affd511b534b72a28e.0p-140L, 174251292Sdas 0xdf0000.0p-24, 0x8d515c.0p-26, -0x1dc08d61c6ef1d9b2ef7e68680598.0p-143L, 175251292Sdas 0xdd8000.0p-24, 0x943a9e.0p-26, -0x1f72a2dac729b3f46662238a9425a.0p-142L, 176251292Sdas 0xdc0000.0p-24, 0x9b2fe6.0p-26, -0x1fd4dfd3a0afb9691aed4d5e3df94.0p-140L, 177251292Sdas 0xda8000.0p-24, 0xa2315d.0p-26, -0x11b26121629c46c186384993e1c93.0p-142L, 178251292Sdas 0xd90000.0p-24, 0xa93f2f.0p-26, 0x1286d633e8e5697dc6a402a56fce1.0p-141L, 179251292Sdas 0xd78000.0p-24, 0xb05988.0p-26, 0x16128eba9367707ebfa540e45350c.0p-144L, 180251292Sdas 0xd60000.0p-24, 0xb78094.0p-26, 0x16ead577390d31ef0f4c9d43f79b2.0p-140L, 181251292Sdas 0xd50000.0p-24, 0xbc4c6c.0p-26, 0x151131ccf7c7b75e7d900b521c48d.0p-141L, 182251292Sdas 0xd38000.0p-24, 0xc3890a.0p-26, -0x115e2cd714bd06508aeb00d2ae3e9.0p-140L, 183251292Sdas 0xd20000.0p-24, 0xcad2d7.0p-26, -0x1847f406ebd3af80485c2f409633c.0p-142L, 184251292Sdas 0xd10000.0p-24, 0xcfb620.0p-26, 0x1c2259904d686581799fbce0b5f19.0p-141L, 185251292Sdas 0xcf8000.0p-24, 0xd71653.0p-26, 0x1ece57a8d5ae54f550444ecf8b995.0p-140L, 186251292Sdas 0xce0000.0p-24, 0xde843a.0p-26, -0x1f109d4bc4595412b5d2517aaac13.0p-141L, 187251292Sdas 0xcd0000.0p-24, 0xe37fde.0p-26, 0x1bc03dc271a74d3a85b5b43c0e727.0p-141L, 188251292Sdas 0xcb8000.0p-24, 0xeb050c.0p-26, -0x1bf2badc0df841a71b79dd5645b46.0p-145L, 189251292Sdas 0xca0000.0p-24, 0xf29878.0p-26, -0x18efededd89fbe0bcfbe6d6db9f66.0p-147L, 190251292Sdas 0xc90000.0p-24, 0xf7ad6f.0p-26, 0x1373ff977baa6911c7bafcb4d84fb.0p-141L, 191251292Sdas 0xc80000.0p-24, 0xfcc8e3.0p-26, 0x196766f2fb328337cc050c6d83b22.0p-140L, 192251292Sdas 0xc68000.0p-24, 0x823f30.0p-25, 0x19bd076f7c434e5fcf1a212e2a91e.0p-139L, 193251292Sdas 0xc58000.0p-24, 0x84d52c.0p-25, -0x1a327257af0f465e5ecab5f2a6f81.0p-139L, 194251292Sdas 0xc40000.0p-24, 0x88bc74.0p-25, 0x113f23def19c5a0fe396f40f1dda9.0p-141L, 195251292Sdas 0xc30000.0p-24, 0x8b5ae6.0p-25, 0x1759f6e6b37de945a049a962e66c6.0p-139L, 196251292Sdas 0xc20000.0p-24, 0x8dfccb.0p-25, 0x1ad35ca6ed5147bdb6ddcaf59c425.0p-141L, 197251292Sdas 0xc10000.0p-24, 0x90a22b.0p-25, 0x1a1d71a87deba46bae9827221dc98.0p-139L, 198251292Sdas 0xbf8000.0p-24, 0x94a0d8.0p-25, -0x139e5210c2b730e28aba001a9b5e0.0p-140L, 199251292Sdas 0xbe8000.0p-24, 0x974f16.0p-25, -0x18f6ebcff3ed72e23e13431adc4a5.0p-141L, 200251292Sdas 0xbd8000.0p-24, 0x9a00f1.0p-25, -0x1aa268be39aab7148e8d80caa10b7.0p-139L, 201251292Sdas 0xbc8000.0p-24, 0x9cb672.0p-25, -0x14c8815839c5663663d15faed7771.0p-139L, 202251292Sdas 0xbb0000.0p-24, 0xa0cda1.0p-25, 0x1eaf46390dbb2438273918db7df5c.0p-141L, 203251292Sdas 0xba0000.0p-24, 0xa38c6e.0p-25, 0x138e20d831f698298adddd7f32686.0p-141L, 204251292Sdas 0xb90000.0p-24, 0xa64f05.0p-25, -0x1e8d3c41123615b147a5d47bc208f.0p-142L, 205251292Sdas 0xb80000.0p-24, 0xa91570.0p-25, 0x1ce28f5f3840b263acb4351104631.0p-140L, 206251292Sdas 0xb70000.0p-24, 0xabdfbb.0p-25, -0x186e5c0a42423457e22d8c650b355.0p-139L, 207251292Sdas 0xb60000.0p-24, 0xaeadef.0p-25, -0x14d41a0b2a08a465dc513b13f567d.0p-143L, 208251292Sdas 0xb50000.0p-24, 0xb18018.0p-25, 0x16755892770633947ffe651e7352f.0p-139L, 209251292Sdas 0xb40000.0p-24, 0xb45642.0p-25, -0x16395ebe59b15228bfe8798d10ff0.0p-142L, 210251292Sdas 0xb30000.0p-24, 0xb73077.0p-25, 0x1abc65c8595f088b61a335f5b688c.0p-140L, 211251292Sdas 0xb20000.0p-24, 0xba0ec4.0p-25, -0x1273089d3dad88e7d353e9967d548.0p-139L, 212251292Sdas 0xb10000.0p-24, 0xbcf133.0p-25, 0x10f9f67b1f4bbf45de06ecebfaf6d.0p-139L, 213251292Sdas 0xb00000.0p-24, 0xbfd7d2.0p-25, -0x109fab904864092b34edda19a831e.0p-140L, 214251292Sdas 0xaf0000.0p-24, 0xc2c2ac.0p-25, -0x1124680aa43333221d8a9b475a6ba.0p-139L, 215251292Sdas 0xae8000.0p-24, 0xc439b3.0p-25, -0x1f360cc4710fbfe24b633f4e8d84d.0p-140L, 216251292Sdas 0xad8000.0p-24, 0xc72afd.0p-25, -0x132d91f21d89c89c45003fc5d7807.0p-140L, 217251292Sdas 0xac8000.0p-24, 0xca20a2.0p-25, -0x16bf9b4d1f8da8002f2449e174504.0p-139L, 218251292Sdas 0xab8000.0p-24, 0xcd1aae.0p-25, 0x19deb5ce6a6a8717d5626e16acc7d.0p-141L, 219251292Sdas 0xaa8000.0p-24, 0xd0192f.0p-25, 0x1a29fb48f7d3ca87dabf351aa41f4.0p-139L, 220251292Sdas 0xaa0000.0p-24, 0xd19a20.0p-25, 0x1127d3c6457f9d79f51dcc73014c9.0p-141L, 221251292Sdas 0xa90000.0p-24, 0xd49f6a.0p-25, -0x1ba930e486a0ac42d1bf9199188e7.0p-141L, 222251292Sdas 0xa80000.0p-24, 0xd7a94b.0p-25, -0x1b6e645f31549dd1160bcc45c7e2c.0p-139L, 223251292Sdas 0xa70000.0p-24, 0xdab7d0.0p-25, 0x1118a425494b610665377f15625b6.0p-140L, 224251292Sdas 0xa68000.0p-24, 0xdc40d5.0p-25, 0x1966f24d29d3a2d1b2176010478be.0p-140L, 225251292Sdas 0xa58000.0p-24, 0xdf566d.0p-25, -0x1d8e52eb2248f0c95dd83626d7333.0p-142L, 226251292Sdas 0xa48000.0p-24, 0xe270ce.0p-25, -0x1ee370f96e6b67ccb006a5b9890ea.0p-140L, 227251292Sdas 0xa40000.0p-24, 0xe3ffce.0p-25, 0x1d155324911f56db28da4d629d00a.0p-140L, 228251292Sdas 0xa30000.0p-24, 0xe72179.0p-25, -0x1fe6e2f2f867d8f4d60c713346641.0p-140L, 229251292Sdas 0xa20000.0p-24, 0xea4812.0p-25, 0x1b7be9add7f4d3b3d406b6cbf3ce5.0p-140L, 230251292Sdas 0xa18000.0p-24, 0xebdd3d.0p-25, 0x1b3cfb3f7511dd73692609040ccc2.0p-139L, 231251292Sdas 0xa08000.0p-24, 0xef0b5b.0p-25, -0x1220de1f7301901b8ad85c25afd09.0p-139L, 232251292Sdas 0xa00000.0p-24, 0xf0a451.0p-25, -0x176364c9ac81cc8a4dfb804de6867.0p-140L, 233251292Sdas 0x9f0000.0p-24, 0xf3da16.0p-25, 0x1eed6b9aafac8d42f78d3e65d3727.0p-141L, 234251292Sdas 0x9e8000.0p-24, 0xf576e9.0p-25, 0x1d593218675af269647b783d88999.0p-139L, 235251292Sdas 0x9d8000.0p-24, 0xf8b47c.0p-25, -0x13e8eb7da053e063714615f7cc91d.0p-144L, 236251292Sdas 0x9d0000.0p-24, 0xfa553f.0p-25, 0x1c063259bcade02951686d5373aec.0p-139L, 237251292Sdas 0x9c0000.0p-24, 0xfd9ac5.0p-25, 0x1ef491085fa3c1649349630531502.0p-139L, 238251292Sdas 0x9b8000.0p-24, 0xff3f8c.0p-25, 0x1d607a7c2b8c5320619fb9433d841.0p-139L, 239251292Sdas 0x9a8000.0p-24, 0x814697.0p-24, -0x12ad3817004f3f0bdff99f932b273.0p-138L, 240251292Sdas 0x9a0000.0p-24, 0x821b06.0p-24, -0x189fc53117f9e54e78103a2bc1767.0p-141L, 241251292Sdas 0x990000.0p-24, 0x83c5f8.0p-24, 0x14cf15a048907b7d7f47ddb45c5a3.0p-139L, 242251292Sdas 0x988000.0p-24, 0x849c7d.0p-24, 0x1cbb1d35fb82873b04a9af1dd692c.0p-138L, 243251292Sdas 0x978000.0p-24, 0x864ba6.0p-24, 0x1128639b814f9b9770d8cb6573540.0p-138L, 244251292Sdas 0x970000.0p-24, 0x87244c.0p-24, 0x184733853300f002e836dfd47bd41.0p-139L, 245251292Sdas 0x968000.0p-24, 0x87fdaa.0p-24, 0x109d23aef77dd5cd7cc94306fb3ff.0p-140L, 246251292Sdas 0x958000.0p-24, 0x89b293.0p-24, -0x1a81ef367a59de2b41eeebd550702.0p-138L, 247251292Sdas 0x950000.0p-24, 0x8a8e20.0p-24, -0x121ad3dbb2f45275c917a30df4ac9.0p-138L, 248251292Sdas 0x948000.0p-24, 0x8b6a6a.0p-24, -0x1cfb981628af71a89df4e6df2e93b.0p-139L, 249251292Sdas 0x938000.0p-24, 0x8d253a.0p-24, -0x1d21730ea76cfdec367828734cae5.0p-139L, 250251292Sdas 0x930000.0p-24, 0x8e03c2.0p-24, 0x135cc00e566f76b87333891e0dec4.0p-138L, 251251292Sdas 0x928000.0p-24, 0x8ee30d.0p-24, -0x10fcb5df257a263e3bf446c6e3f69.0p-140L, 252251292Sdas 0x918000.0p-24, 0x90a3ee.0p-24, -0x16e171b15433d723a4c7380a448d8.0p-139L, 253251292Sdas 0x910000.0p-24, 0x918587.0p-24, -0x1d050da07f3236f330972da2a7a87.0p-139L, 254251292Sdas 0x908000.0p-24, 0x9267e7.0p-24, 0x1be03669a5268d21148c6002becd3.0p-139L, 255251292Sdas 0x8f8000.0p-24, 0x942f04.0p-24, 0x10b28e0e26c336af90e00533323ba.0p-139L, 256251292Sdas 0x8f0000.0p-24, 0x9513c3.0p-24, 0x1a1d820da57cf2f105a89060046aa.0p-138L, 257251292Sdas 0x8e8000.0p-24, 0x95f950.0p-24, -0x19ef8f13ae3cf162409d8ea99d4c0.0p-139L, 258251292Sdas 0x8e0000.0p-24, 0x96dfab.0p-24, -0x109e417a6e507b9dc10dac743ad7a.0p-138L, 259251292Sdas 0x8d0000.0p-24, 0x98aed2.0p-24, 0x10d01a2c5b0e97c4990b23d9ac1f5.0p-139L, 260251292Sdas 0x8c8000.0p-24, 0x9997a2.0p-24, -0x1d6a50d4b61ea74540bdd2aa99a42.0p-138L, 261251292Sdas 0x8c0000.0p-24, 0x9a8145.0p-24, 0x1b3b190b83f9527e6aba8f2d783c1.0p-138L, 262251292Sdas 0x8b8000.0p-24, 0x9b6bbf.0p-24, 0x13a69fad7e7abe7ba81c664c107e0.0p-138L, 263251292Sdas 0x8b0000.0p-24, 0x9c5711.0p-24, -0x11cd12316f576aad348ae79867223.0p-138L, 264251292Sdas 0x8a8000.0p-24, 0x9d433b.0p-24, 0x1c95c444b807a246726b304ccae56.0p-139L, 265251292Sdas 0x898000.0p-24, 0x9f1e22.0p-24, -0x1b9c224ea698c2f9b47466d6123fe.0p-139L, 266251292Sdas 0x890000.0p-24, 0xa00ce1.0p-24, 0x125ca93186cf0f38b4619a2483399.0p-141L, 267251292Sdas 0x888000.0p-24, 0xa0fc80.0p-24, -0x1ee38a7bc228b3597043be78eaf49.0p-139L, 268251292Sdas 0x880000.0p-24, 0xa1ed00.0p-24, -0x1a0db876613d204147dc69a07a649.0p-138L, 269251292Sdas 0x878000.0p-24, 0xa2de62.0p-24, 0x193224e8516c008d3602a7b41c6e8.0p-139L, 270251292Sdas 0x870000.0p-24, 0xa3d0a9.0p-24, 0x1fa28b4d2541aca7d5844606b2421.0p-139L, 271251292Sdas 0x868000.0p-24, 0xa4c3d6.0p-24, 0x1c1b5760fb4571acbcfb03f16daf4.0p-138L, 272251292Sdas 0x858000.0p-24, 0xa6acea.0p-24, 0x1fed5d0f65949c0a345ad743ae1ae.0p-140L, 273251292Sdas 0x850000.0p-24, 0xa7a2d4.0p-24, 0x1ad270c9d749362382a7688479e24.0p-140L, 274251292Sdas 0x848000.0p-24, 0xa899ab.0p-24, 0x199ff15ce532661ea9643a3a2d378.0p-139L, 275251292Sdas 0x840000.0p-24, 0xa99171.0p-24, 0x1a19e15ccc45d257530a682b80490.0p-139L, 276251292Sdas 0x838000.0p-24, 0xaa8a28.0p-24, -0x121a14ec532b35ba3e1f868fd0b5e.0p-140L, 277251292Sdas 0x830000.0p-24, 0xab83d1.0p-24, 0x1aee319980bff3303dd481779df69.0p-139L, 278251292Sdas 0x828000.0p-24, 0xac7e6f.0p-24, -0x18ffd9e3900345a85d2d86161742e.0p-140L, 279251292Sdas 0x820000.0p-24, 0xad7a03.0p-24, -0x1e4db102ce29f79b026b64b42caa1.0p-140L, 280251292Sdas 0x818000.0p-24, 0xae768f.0p-24, 0x17c35c55a04a82ab19f77652d977a.0p-141L, 281251292Sdas 0x810000.0p-24, 0xaf7415.0p-24, 0x1448324047019b48d7b98c1cf7234.0p-138L, 282251292Sdas 0x808000.0p-24, 0xb07298.0p-24, -0x1750ee3915a197e9c7359dd94152f.0p-138L, 283251292Sdas 0x800000.0p-24, 0xb17218.0p-24, -0x105c610ca86c3898cff81a12a17e2.0p-141L, 284251292Sdas}; 285251292Sdas 286251292Sdas#ifdef USE_UTAB 287251292Sdasstatic const struct { 288251292Sdas float H; /* 1 + i/INTERVALS (exact) */ 289251292Sdas float E; /* H(i) * G(i) - 1 (exact) */ 290251292Sdas} U[TSIZE] = { 291251292Sdas 0x800000.0p-23, 0, 292251292Sdas 0x810000.0p-23, -0x800000.0p-37, 293251292Sdas 0x820000.0p-23, -0x800000.0p-35, 294251292Sdas 0x830000.0p-23, -0x900000.0p-34, 295251292Sdas 0x840000.0p-23, -0x800000.0p-33, 296251292Sdas 0x850000.0p-23, -0xc80000.0p-33, 297251292Sdas 0x860000.0p-23, -0xa00000.0p-36, 298251292Sdas 0x870000.0p-23, 0x940000.0p-33, 299251292Sdas 0x880000.0p-23, 0x800000.0p-35, 300251292Sdas 0x890000.0p-23, -0xc80000.0p-34, 301251292Sdas 0x8a0000.0p-23, 0xe00000.0p-36, 302251292Sdas 0x8b0000.0p-23, 0x900000.0p-33, 303251292Sdas 0x8c0000.0p-23, -0x800000.0p-35, 304251292Sdas 0x8d0000.0p-23, -0xe00000.0p-33, 305251292Sdas 0x8e0000.0p-23, 0x880000.0p-33, 306251292Sdas 0x8f0000.0p-23, -0xa80000.0p-34, 307251292Sdas 0x900000.0p-23, -0x800000.0p-35, 308251292Sdas 0x910000.0p-23, 0x800000.0p-37, 309251292Sdas 0x920000.0p-23, 0x900000.0p-35, 310251292Sdas 0x930000.0p-23, 0xd00000.0p-35, 311251292Sdas 0x940000.0p-23, 0xe00000.0p-35, 312251292Sdas 0x950000.0p-23, 0xc00000.0p-35, 313251292Sdas 0x960000.0p-23, 0xe00000.0p-36, 314251292Sdas 0x970000.0p-23, -0x800000.0p-38, 315251292Sdas 0x980000.0p-23, -0xc00000.0p-35, 316251292Sdas 0x990000.0p-23, -0xd00000.0p-34, 317251292Sdas 0x9a0000.0p-23, 0x880000.0p-33, 318251292Sdas 0x9b0000.0p-23, 0xe80000.0p-35, 319251292Sdas 0x9c0000.0p-23, -0x800000.0p-35, 320251292Sdas 0x9d0000.0p-23, 0xb40000.0p-33, 321251292Sdas 0x9e0000.0p-23, 0x880000.0p-34, 322251292Sdas 0x9f0000.0p-23, -0xe00000.0p-35, 323251292Sdas 0xa00000.0p-23, 0x800000.0p-33, 324251292Sdas 0xa10000.0p-23, -0x900000.0p-36, 325251292Sdas 0xa20000.0p-23, -0xb00000.0p-33, 326251292Sdas 0xa30000.0p-23, -0xa00000.0p-36, 327251292Sdas 0xa40000.0p-23, 0x800000.0p-33, 328251292Sdas 0xa50000.0p-23, -0xf80000.0p-35, 329251292Sdas 0xa60000.0p-23, 0x880000.0p-34, 330251292Sdas 0xa70000.0p-23, -0x900000.0p-33, 331251292Sdas 0xa80000.0p-23, -0x800000.0p-35, 332251292Sdas 0xa90000.0p-23, 0x900000.0p-34, 333251292Sdas 0xaa0000.0p-23, 0xa80000.0p-33, 334251292Sdas 0xab0000.0p-23, -0xac0000.0p-34, 335251292Sdas 0xac0000.0p-23, -0x800000.0p-37, 336251292Sdas 0xad0000.0p-23, 0xf80000.0p-35, 337251292Sdas 0xae0000.0p-23, 0xf80000.0p-34, 338251292Sdas 0xaf0000.0p-23, -0xac0000.0p-33, 339251292Sdas 0xb00000.0p-23, -0x800000.0p-33, 340251292Sdas 0xb10000.0p-23, -0xb80000.0p-34, 341251292Sdas 0xb20000.0p-23, -0x800000.0p-34, 342251292Sdas 0xb30000.0p-23, -0xb00000.0p-35, 343251292Sdas 0xb40000.0p-23, -0x800000.0p-35, 344251292Sdas 0xb50000.0p-23, -0xe00000.0p-36, 345251292Sdas 0xb60000.0p-23, -0x800000.0p-35, 346251292Sdas 0xb70000.0p-23, -0xb00000.0p-35, 347251292Sdas 0xb80000.0p-23, -0x800000.0p-34, 348251292Sdas 0xb90000.0p-23, -0xb80000.0p-34, 349251292Sdas 0xba0000.0p-23, -0x800000.0p-33, 350251292Sdas 0xbb0000.0p-23, -0xac0000.0p-33, 351251292Sdas 0xbc0000.0p-23, 0x980000.0p-33, 352251292Sdas 0xbd0000.0p-23, 0xbc0000.0p-34, 353251292Sdas 0xbe0000.0p-23, 0xe00000.0p-36, 354251292Sdas 0xbf0000.0p-23, -0xb80000.0p-35, 355251292Sdas 0xc00000.0p-23, -0x800000.0p-33, 356251292Sdas 0xc10000.0p-23, 0xa80000.0p-33, 357251292Sdas 0xc20000.0p-23, 0x900000.0p-34, 358251292Sdas 0xc30000.0p-23, -0x800000.0p-35, 359251292Sdas 0xc40000.0p-23, -0x900000.0p-33, 360251292Sdas 0xc50000.0p-23, 0x820000.0p-33, 361251292Sdas 0xc60000.0p-23, 0x800000.0p-38, 362251292Sdas 0xc70000.0p-23, -0x820000.0p-33, 363251292Sdas 0xc80000.0p-23, 0x800000.0p-33, 364251292Sdas 0xc90000.0p-23, -0xa00000.0p-36, 365251292Sdas 0xca0000.0p-23, -0xb00000.0p-33, 366251292Sdas 0xcb0000.0p-23, 0x840000.0p-34, 367251292Sdas 0xcc0000.0p-23, -0xd00000.0p-34, 368251292Sdas 0xcd0000.0p-23, 0x800000.0p-33, 369251292Sdas 0xce0000.0p-23, -0xe00000.0p-35, 370251292Sdas 0xcf0000.0p-23, 0xa60000.0p-33, 371251292Sdas 0xd00000.0p-23, -0x800000.0p-35, 372251292Sdas 0xd10000.0p-23, 0xb40000.0p-33, 373251292Sdas 0xd20000.0p-23, -0x800000.0p-35, 374251292Sdas 0xd30000.0p-23, 0xaa0000.0p-33, 375251292Sdas 0xd40000.0p-23, -0xe00000.0p-35, 376251292Sdas 0xd50000.0p-23, 0x880000.0p-33, 377251292Sdas 0xd60000.0p-23, -0xd00000.0p-34, 378251292Sdas 0xd70000.0p-23, 0x9c0000.0p-34, 379251292Sdas 0xd80000.0p-23, -0xb00000.0p-33, 380251292Sdas 0xd90000.0p-23, -0x800000.0p-38, 381251292Sdas 0xda0000.0p-23, 0xa40000.0p-33, 382251292Sdas 0xdb0000.0p-23, -0xdc0000.0p-34, 383251292Sdas 0xdc0000.0p-23, 0xc00000.0p-35, 384251292Sdas 0xdd0000.0p-23, 0xca0000.0p-33, 385251292Sdas 0xde0000.0p-23, -0xb80000.0p-34, 386251292Sdas 0xdf0000.0p-23, 0xd00000.0p-35, 387251292Sdas 0xe00000.0p-23, 0xc00000.0p-33, 388251292Sdas 0xe10000.0p-23, -0xf40000.0p-34, 389251292Sdas 0xe20000.0p-23, 0x800000.0p-37, 390251292Sdas 0xe30000.0p-23, 0x860000.0p-33, 391251292Sdas 0xe40000.0p-23, -0xc80000.0p-33, 392251292Sdas 0xe50000.0p-23, -0xa80000.0p-34, 393251292Sdas 0xe60000.0p-23, 0xe00000.0p-36, 394251292Sdas 0xe70000.0p-23, 0x880000.0p-33, 395251292Sdas 0xe80000.0p-23, -0xe00000.0p-33, 396251292Sdas 0xe90000.0p-23, -0xfc0000.0p-34, 397251292Sdas 0xea0000.0p-23, -0x800000.0p-35, 398251292Sdas 0xeb0000.0p-23, 0xe80000.0p-35, 399251292Sdas 0xec0000.0p-23, 0x900000.0p-33, 400251292Sdas 0xed0000.0p-23, 0xe20000.0p-33, 401251292Sdas 0xee0000.0p-23, -0xac0000.0p-33, 402251292Sdas 0xef0000.0p-23, -0xc80000.0p-34, 403251292Sdas 0xf00000.0p-23, -0x800000.0p-35, 404251292Sdas 0xf10000.0p-23, 0x800000.0p-35, 405251292Sdas 0xf20000.0p-23, 0xb80000.0p-34, 406251292Sdas 0xf30000.0p-23, 0x940000.0p-33, 407251292Sdas 0xf40000.0p-23, 0xc80000.0p-33, 408251292Sdas 0xf50000.0p-23, -0xf20000.0p-33, 409251292Sdas 0xf60000.0p-23, -0xc80000.0p-33, 410251292Sdas 0xf70000.0p-23, -0xa20000.0p-33, 411251292Sdas 0xf80000.0p-23, -0x800000.0p-33, 412251292Sdas 0xf90000.0p-23, -0xc40000.0p-34, 413251292Sdas 0xfa0000.0p-23, -0x900000.0p-34, 414251292Sdas 0xfb0000.0p-23, -0xc80000.0p-35, 415251292Sdas 0xfc0000.0p-23, -0x800000.0p-35, 416251292Sdas 0xfd0000.0p-23, -0x900000.0p-36, 417251292Sdas 0xfe0000.0p-23, -0x800000.0p-37, 418251292Sdas 0xff0000.0p-23, -0x800000.0p-39, 419251292Sdas 0x800000.0p-22, 0, 420251292Sdas}; 421251292Sdas#endif /* USE_UTAB */ 422251292Sdas 423251292Sdas#ifdef STRUCT_RETURN 424251292Sdas#define RETURN1(rp, v) do { \ 425251292Sdas (rp)->hi = (v); \ 426251292Sdas (rp)->lo_set = 0; \ 427251292Sdas return; \ 428251292Sdas} while (0) 429251292Sdas 430251292Sdas#define RETURN2(rp, h, l) do { \ 431251292Sdas (rp)->hi = (h); \ 432251292Sdas (rp)->lo = (l); \ 433251292Sdas (rp)->lo_set = 1; \ 434251292Sdas return; \ 435251292Sdas} while (0) 436251292Sdas 437251292Sdasstruct ld { 438251292Sdas long double hi; 439251292Sdas long double lo; 440251292Sdas int lo_set; 441251292Sdas}; 442251292Sdas#else 443251292Sdas#define RETURN1(rp, v) RETURNF(v) 444251292Sdas#define RETURN2(rp, h, l) RETURNI((h) + (l)) 445251292Sdas#endif 446251292Sdas 447251292Sdas#ifdef STRUCT_RETURN 448251292Sdasstatic inline __always_inline void 449251292Sdask_logl(long double x, struct ld *rp) 450251292Sdas#else 451251292Sdaslong double 452251292Sdaslogl(long double x) 453251292Sdas#endif 454251292Sdas{ 455251292Sdas long double d, val_hi, val_lo; 456251292Sdas double dd, dk; 457251292Sdas uint64_t lx, llx; 458251292Sdas int i, k; 459251292Sdas uint16_t hx; 460251292Sdas 461251292Sdas EXTRACT_LDBL128_WORDS(hx, lx, llx, x); 462251292Sdas k = -16383; 463251292Sdas#if 0 /* Hard to do efficiently. Don't do it until we support all modes. */ 464251292Sdas if (x == 1) 465251292Sdas RETURN1(rp, 0); /* log(1) = +0 in all rounding modes */ 466251292Sdas#endif 467251292Sdas if (hx == 0 || hx >= 0x8000) { /* zero, negative or subnormal? */ 468251292Sdas if (((hx & 0x7fff) | lx | llx) == 0) 469251292Sdas RETURN1(rp, -1 / zero); /* log(+-0) = -Inf */ 470251292Sdas if (hx != 0) 471251292Sdas /* log(neg or NaN) = qNaN: */ 472251292Sdas RETURN1(rp, (x - x) / zero); 473251292Sdas x *= 0x1.0p113; /* subnormal; scale up x */ 474251292Sdas EXTRACT_LDBL128_WORDS(hx, lx, llx, x); 475251292Sdas k = -16383 - 113; 476251292Sdas } else if (hx >= 0x7fff) 477251292Sdas RETURN1(rp, x + x); /* log(Inf or NaN) = Inf or qNaN */ 478251292Sdas#ifndef STRUCT_RETURN 479251292Sdas ENTERI(); 480251292Sdas#endif 481251292Sdas k += hx; 482251292Sdas dk = k; 483251292Sdas 484251292Sdas /* Scale x to be in [1, 2). */ 485251292Sdas SET_LDBL_EXPSIGN(x, 0x3fff); 486251292Sdas 487251292Sdas /* 0 <= i <= INTERVALS: */ 488251292Sdas#define L2I (49 - LOG2_INTERVALS) 489251292Sdas i = (lx + (1LL << (L2I - 2))) >> (L2I - 1); 490251292Sdas 491251292Sdas /* 492251292Sdas * -0.005280 < d < 0.004838. In particular, the infinite- 493251292Sdas * precision |d| is <= 2**-7. Rounding of G(i) to 8 bits 494251292Sdas * ensures that d is representable without extra precision for 495251292Sdas * this bound on |d| (since when this calculation is expressed 496251292Sdas * as x*G(i)-1, the multiplication needs as many extra bits as 497251292Sdas * G(i) has and the subtraction cancels 8 bits). But for 498251292Sdas * most i (107 cases out of 129), the infinite-precision |d| 499251292Sdas * is <= 2**-8. G(i) is rounded to 9 bits for such i to give 500251292Sdas * better accuracy (this works by improving the bound on |d|, 501251292Sdas * which in turn allows rounding to 9 bits in more cases). 502251292Sdas * This is only important when the original x is near 1 -- it 503251292Sdas * lets us avoid using a special method to give the desired 504251292Sdas * accuracy for such x. 505251292Sdas */ 506251292Sdas if (0) 507251292Sdas d = x * G(i) - 1; 508251292Sdas else { 509251292Sdas#ifdef USE_UTAB 510251292Sdas d = (x - H(i)) * G(i) + E(i); 511251292Sdas#else 512251292Sdas long double x_hi; 513251292Sdas double x_lo; 514251292Sdas 515251292Sdas /* 516251292Sdas * Split x into x_hi + x_lo to calculate x*G(i)-1 exactly. 517251292Sdas * G(i) has at most 9 bits, so the splitting point is not 518251292Sdas * critical. 519251292Sdas */ 520251292Sdas INSERT_LDBL128_WORDS(x_hi, 0x3fff, lx, 521251292Sdas llx & 0xffffffffff000000ULL); 522251292Sdas x_lo = x - x_hi; 523251292Sdas d = x_hi * G(i) - 1 + x_lo * G(i); 524251292Sdas#endif 525251292Sdas } 526251292Sdas 527251292Sdas /* 528251292Sdas * Our algorithm depends on exact cancellation of F_lo(i) and 529251292Sdas * F_hi(i) with dk*ln_2_lo and dk*ln2_hi when k is -1 and i is 530251292Sdas * at the end of the table. This and other technical complications 531251292Sdas * make it difficult to avoid the double scaling in (dk*ln2) * 532251292Sdas * log(base) for base != e without losing more accuracy and/or 533251292Sdas * efficiency than is gained. 534251292Sdas */ 535251292Sdas /* 536251292Sdas * Use double precision operations wherever possible, since long 537251292Sdas * double operations are emulated and are very slow on the only 538251292Sdas * known machines that support ld128 (sparc64). Also, don't try 539251292Sdas * to improve parallelism by increasing the number of operations, 540251292Sdas * since any parallelism on such machines is needed for the 541251292Sdas * emulation. Horner's method is good for this, and is also good 542251292Sdas * for accuracy. Horner's method doesn't handle the `lo' term 543251292Sdas * well, either for efficiency or accuracy. However, for accuracy 544251292Sdas * we evaluate d * d * P2 separately to take advantage of 545251292Sdas * by P2 being exact, and this gives a good place to sum the 'lo' 546251292Sdas * term too. 547251292Sdas */ 548251292Sdas dd = (double)d; 549251292Sdas val_lo = d * d * d * (P3 + 550251292Sdas d * (P4 + d * (P5 + d * (P6 + d * (P7 + d * (P8 + 551251292Sdas dd * (P9 + dd * (P10 + dd * (P11 + dd * (P12 + dd * (P13 + 552251292Sdas dd * P14))))))))))) + (F_lo(i) + dk * ln2_lo) + d * d * P2; 553251292Sdas val_hi = d; 554251292Sdas#ifdef DEBUG 555251292Sdas if (fetestexcept(FE_UNDERFLOW)) 556251292Sdas breakpoint(); 557251292Sdas#endif 558251292Sdas 559251292Sdas _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi); 560251292Sdas RETURN2(rp, val_hi, val_lo); 561251292Sdas} 562251292Sdas 563251292Sdaslong double 564251292Sdaslog1pl(long double x) 565251292Sdas{ 566251292Sdas long double d, d_hi, f_lo, val_hi, val_lo; 567251292Sdas long double f_hi, twopminusk; 568251292Sdas double d_lo, dd, dk; 569251292Sdas uint64_t lx, llx; 570251292Sdas int i, k; 571251292Sdas int16_t ax, hx; 572251292Sdas 573251292Sdas DOPRINT_START(&x); 574251292Sdas EXTRACT_LDBL128_WORDS(hx, lx, llx, x); 575251292Sdas if (hx < 0x3fff) { /* x < 1, or x neg NaN */ 576251292Sdas ax = hx & 0x7fff; 577251292Sdas if (ax >= 0x3fff) { /* x <= -1, or x neg NaN */ 578251292Sdas if (ax == 0x3fff && (lx | llx) == 0) 579251292Sdas RETURNP(-1 / zero); /* log1p(-1) = -Inf */ 580251292Sdas /* log1p(x < 1, or x NaN) = qNaN: */ 581251292Sdas RETURNP((x - x) / (x - x)); 582251292Sdas } 583251292Sdas if (ax <= 0x3f8d) { /* |x| < 2**-113 */ 584251292Sdas if ((int)x == 0) 585251292Sdas RETURNP(x); /* x with inexact if x != 0 */ 586251292Sdas } 587251292Sdas f_hi = 1; 588251292Sdas f_lo = x; 589251292Sdas } else if (hx >= 0x7fff) { /* x +Inf or non-neg NaN */ 590251292Sdas RETURNP(x + x); /* log1p(Inf or NaN) = Inf or qNaN */ 591251292Sdas } else if (hx < 0x40e1) { /* 1 <= x < 2**226 */ 592251292Sdas f_hi = x; 593251292Sdas f_lo = 1; 594251292Sdas } else { /* 2**226 <= x < +Inf */ 595251292Sdas f_hi = x; 596251292Sdas f_lo = 0; /* avoid underflow of the P3 term */ 597251292Sdas } 598251292Sdas ENTERI(); 599251292Sdas x = f_hi + f_lo; 600251292Sdas f_lo = (f_hi - x) + f_lo; 601251292Sdas 602251292Sdas EXTRACT_LDBL128_WORDS(hx, lx, llx, x); 603251292Sdas k = -16383; 604251292Sdas 605251292Sdas k += hx; 606251292Sdas dk = k; 607251292Sdas 608251292Sdas SET_LDBL_EXPSIGN(x, 0x3fff); 609251292Sdas twopminusk = 1; 610251292Sdas SET_LDBL_EXPSIGN(twopminusk, 0x7ffe - (hx & 0x7fff)); 611251292Sdas f_lo *= twopminusk; 612251292Sdas 613251292Sdas i = (lx + (1LL << (L2I - 2))) >> (L2I - 1); 614251292Sdas 615251292Sdas /* 616251292Sdas * x*G(i)-1 (with a reduced x) can be represented exactly, as 617251292Sdas * above, but now we need to evaluate the polynomial on d = 618251292Sdas * (x+f_lo)*G(i)-1 and extra precision is needed for that. 619251292Sdas * Since x+x_lo is a hi+lo decomposition and subtracting 1 620251292Sdas * doesn't lose too many bits, an inexact calculation for 621251292Sdas * f_lo*G(i) is good enough. 622251292Sdas */ 623251292Sdas if (0) 624251292Sdas d_hi = x * G(i) - 1; 625251292Sdas else { 626251292Sdas#ifdef USE_UTAB 627251292Sdas d_hi = (x - H(i)) * G(i) + E(i); 628251292Sdas#else 629251292Sdas long double x_hi; 630251292Sdas double x_lo; 631251292Sdas 632251292Sdas INSERT_LDBL128_WORDS(x_hi, 0x3fff, lx, 633251292Sdas llx & 0xffffffffff000000ULL); 634251292Sdas x_lo = x - x_hi; 635251292Sdas d_hi = x_hi * G(i) - 1 + x_lo * G(i); 636251292Sdas#endif 637251292Sdas } 638251292Sdas d_lo = f_lo * G(i); 639251292Sdas 640251292Sdas /* 641251292Sdas * This is _2sumF(d_hi, d_lo) inlined. The condition 642251292Sdas * (d_hi == 0 || |d_hi| >= |d_lo|) for using _2sumF() is not 643251292Sdas * always satisifed, so it is not clear that this works, but 644251292Sdas * it works in practice. It works even if it gives a wrong 645251292Sdas * normalized d_lo, since |d_lo| > |d_hi| implies that i is 646251292Sdas * nonzero and d is tiny, so the F(i) term dominates d_lo. 647251292Sdas * In float precision: 648251292Sdas * (By exhaustive testing, the worst case is d_hi = 0x1.bp-25. 649251292Sdas * And if d is only a little tinier than that, we would have 650251292Sdas * another underflow problem for the P3 term; this is also ruled 651251292Sdas * out by exhaustive testing.) 652251292Sdas */ 653251292Sdas d = d_hi + d_lo; 654251292Sdas d_lo = d_hi - d + d_lo; 655251292Sdas d_hi = d; 656251292Sdas 657251292Sdas dd = (double)d; 658251292Sdas val_lo = d * d * d * (P3 + 659251292Sdas d * (P4 + d * (P5 + d * (P6 + d * (P7 + d * (P8 + 660251292Sdas dd * (P9 + dd * (P10 + dd * (P11 + dd * (P12 + dd * (P13 + 661251292Sdas dd * P14))))))))))) + (F_lo(i) + dk * ln2_lo + d_lo) + d * d * P2; 662251292Sdas val_hi = d_hi; 663251292Sdas#ifdef DEBUG 664251292Sdas if (fetestexcept(FE_UNDERFLOW)) 665251292Sdas breakpoint(); 666251292Sdas#endif 667251292Sdas 668251292Sdas _3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi); 669251292Sdas RETURN2PI(val_hi, val_lo); 670251292Sdas} 671251292Sdas 672251292Sdas#ifdef STRUCT_RETURN 673251292Sdas 674251292Sdaslong double 675251292Sdaslogl(long double x) 676251292Sdas{ 677251292Sdas struct ld r; 678251292Sdas 679251292Sdas ENTERI(); 680251292Sdas DOPRINT_START(&x); 681251292Sdas k_logl(x, &r); 682251292Sdas RETURNSPI(&r); 683251292Sdas} 684251292Sdas 685251292Sdas/* 686251292Sdas * 29+113 bit decompositions. The bits are distributed so that the products 687251292Sdas * of the hi terms are exact in double precision. The types are chosen so 688251292Sdas * that the products of the hi terms are done in at least double precision, 689251292Sdas * without any explicit conversions. More natural choices would require a 690251292Sdas * slow long double precision multiplication. 691251292Sdas */ 692251292Sdasstatic const double 693251292Sdasinvln10_hi = 4.3429448176175356e-1, /* 0x1bcb7b15000000.0p-54 */ 694251292Sdasinvln2_hi = 1.4426950402557850e0; /* 0x17154765000000.0p-52 */ 695251292Sdasstatic const long double 696251292Sdasinvln10_lo = 1.41498268538580090791605082294397000e-10L, /* 0x137287195355baaafad33dc323ee3.0p-145L */ 697251292Sdasinvln2_lo = 6.33178418956604368501892137426645911e-10L; /* 0x15c17f0bbbe87fed0691d3e88eb57.0p-143L */ 698251292Sdas 699251292Sdaslong double 700251292Sdaslog10l(long double x) 701251292Sdas{ 702251292Sdas struct ld r; 703251292Sdas long double lo; 704251292Sdas float hi; 705251292Sdas 706251292Sdas ENTERI(); 707251292Sdas DOPRINT_START(&x); 708251292Sdas k_logl(x, &r); 709251292Sdas if (!r.lo_set) 710251292Sdas RETURNPI(r.hi); 711251292Sdas _2sumF(r.hi, r.lo); 712251292Sdas hi = r.hi; 713251292Sdas lo = r.lo + (r.hi - hi); 714251292Sdas RETURN2PI(invln10_hi * hi, 715251292Sdas (invln10_lo + invln10_hi) * lo + invln10_lo * hi); 716251292Sdas} 717251292Sdas 718251292Sdaslong double 719251292Sdaslog2l(long double x) 720251292Sdas{ 721251292Sdas struct ld r; 722251292Sdas long double lo; 723251292Sdas float hi; 724251292Sdas 725251292Sdas ENTERI(); 726251292Sdas DOPRINT_START(&x); 727251292Sdas k_logl(x, &r); 728251292Sdas if (!r.lo_set) 729251292Sdas RETURNPI(r.hi); 730251292Sdas _2sumF(r.hi, r.lo); 731251292Sdas hi = r.hi; 732251292Sdas lo = r.lo + (r.hi - hi); 733251292Sdas RETURN2PI(invln2_hi * hi, 734251292Sdas (invln2_lo + invln2_hi) * lo + invln2_lo * hi); 735251292Sdas} 736251292Sdas 737251292Sdas#endif /* STRUCT_RETURN */ 738