Home | History | Annotate | Line # | Download | only in ld80
      1  1.1  christos /*-
      2  1.1  christos  * SPDX-License-Identifier: BSD-2-Clause
      3  1.1  christos  *
      4  1.1  christos  * Copyright (c) 2007-2013 Bruce D. Evans
      5  1.1  christos  * All rights reserved.
      6  1.1  christos  *
      7  1.1  christos  * Redistribution and use in source and binary forms, with or without
      8  1.1  christos  * modification, are permitted provided that the following conditions
      9  1.1  christos  * are met:
     10  1.1  christos  * 1. Redistributions of source code must retain the above copyright
     11  1.1  christos  *    notice unmodified, this list of conditions, and the following
     12  1.1  christos  *    disclaimer.
     13  1.1  christos  * 2. Redistributions in binary form must reproduce the above copyright
     14  1.1  christos  *    notice, this list of conditions and the following disclaimer in the
     15  1.1  christos  *    documentation and/or other materials provided with the distribution.
     16  1.1  christos  *
     17  1.1  christos  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
     18  1.1  christos  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
     19  1.1  christos  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
     20  1.1  christos  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
     21  1.1  christos  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
     22  1.1  christos  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
     23  1.1  christos  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
     24  1.1  christos  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
     25  1.1  christos  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
     26  1.1  christos  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
     27  1.1  christos  */
     28  1.1  christos 
     29  1.1  christos #include <sys/cdefs.h>
     30  1.1  christos /**
     31  1.1  christos  * Implementation of the natural logarithm of x for Intel 80-bit format.
     32  1.1  christos  *
     33  1.1  christos  * First decompose x into its base 2 representation:
     34  1.1  christos  *
     35  1.1  christos  *    log(x) = log(X * 2**k), where X is in [1, 2)
     36  1.1  christos  *           = log(X) + k * log(2).
     37  1.1  christos  *
     38  1.1  christos  * Let X = X_i + e, where X_i is the center of one of the intervals
     39  1.1  christos  * [-1.0/256, 1.0/256), [1.0/256, 3.0/256), .... [2.0-1.0/256, 2.0+1.0/256)
     40  1.1  christos  * and X is in this interval.  Then
     41  1.1  christos  *
     42  1.1  christos  *    log(X) = log(X_i + e)
     43  1.1  christos  *           = log(X_i * (1 + e / X_i))
     44  1.1  christos  *           = log(X_i) + log(1 + e / X_i).
     45  1.1  christos  *
     46  1.1  christos  * The values log(X_i) are tabulated below.  Let d = e / X_i and use
     47  1.1  christos  *
     48  1.1  christos  *    log(1 + d) = p(d)
     49  1.1  christos  *
     50  1.1  christos  * where p(d) = d - 0.5*d*d + ... is a special minimax polynomial of
     51  1.1  christos  * suitably high degree.
     52  1.1  christos  *
     53  1.1  christos  * To get sufficiently small roundoff errors, k * log(2), log(X_i), and
     54  1.1  christos  * sometimes (if |k| is not large) the first term in p(d) must be evaluated
     55  1.1  christos  * and added up in extra precision.  Extra precision is not needed for the
     56  1.1  christos  * rest of p(d).  In the worst case when k = 0 and log(X_i) is 0, the final
     57  1.1  christos  * error is controlled mainly by the error in the second term in p(d).  The
     58  1.1  christos  * error in this term itself is at most 0.5 ulps from the d*d operation in
     59  1.1  christos  * it.  The error in this term relative to the first term is thus at most
     60  1.1  christos  * 0.5 * |-0.5| * |d| < 1.0/1024 ulps.  We aim for an accumulated error of
     61  1.1  christos  * at most twice this at the point of the final rounding step.  Thus the
     62  1.1  christos  * final error should be at most 0.5 + 1.0/512 = 0.5020 ulps.  Exhaustive
     63  1.1  christos  * testing of a float variant of this function showed a maximum final error
     64  1.1  christos  * of 0.5008 ulps.  Non-exhaustive testing of a double variant of this
     65  1.1  christos  * function showed a maximum final error of 0.5078 ulps (near 1+1.0/256).
     66  1.1  christos  *
     67  1.1  christos  * We made the maximum of |d| (and thus the total relative error and the
     68  1.1  christos  * degree of p(d)) small by using a large number of intervals.  Using
     69  1.1  christos  * centers of intervals instead of endpoints reduces this maximum by a
     70  1.1  christos  * factor of 2 for a given number of intervals.  p(d) is special only
     71  1.1  christos  * in beginning with the Taylor coefficients 0 + 1*d, which tends to happen
     72  1.1  christos  * naturally.  The most accurate minimax polynomial of a given degree might
     73  1.1  christos  * be different, but then we wouldn't want it since we would have to do
     74  1.1  christos  * extra work to avoid roundoff error (especially for P0*d instead of d).
     75  1.1  christos  */
     76  1.1  christos 
     77  1.1  christos #ifdef DEBUG
     78  1.1  christos #include <fenv.h>
     79  1.1  christos #endif
     80  1.1  christos 
     81  1.1  christos #ifdef __FreeBSD__
     82  1.1  christos #include "fpmath.h"
     83  1.1  christos #endif
     84  1.1  christos #include "math.h"
     85  1.1  christos #define	i386_SSE_GOOD
     86  1.1  christos #ifndef NO_STRUCT_RETURN
     87  1.1  christos #define	STRUCT_RETURN
     88  1.1  christos #endif
     89  1.1  christos #include "math_private.h"
     90  1.1  christos 
     91  1.1  christos #if !defined(NO_UTAB) && !defined(NO_UTABL)
     92  1.1  christos #define	USE_UTAB
     93  1.1  christos #endif
     94  1.1  christos 
     95  1.1  christos /*
     96  1.1  christos  * Domain [-0.005280, 0.004838], range ~[-5.1736e-22, 5.1738e-22]:
     97  1.1  christos  * |log(1 + d)/d - p(d)| < 2**-70.7
     98  1.1  christos  */
     99  1.1  christos static const double
    100  1.1  christos P2 = -0.5,
    101  1.1  christos P3 =  3.3333333333333359e-1,		/*  0x1555555555555a.0p-54 */
    102  1.1  christos P4 = -2.5000000000004424e-1,		/* -0x1000000000031d.0p-54 */
    103  1.1  christos P5 =  1.9999999992970016e-1,		/*  0x1999999972f3c7.0p-55 */
    104  1.1  christos P6 = -1.6666666072191585e-1,		/* -0x15555548912c09.0p-55 */
    105  1.1  christos P7 =  1.4286227413310518e-1,		/*  0x12494f9d9def91.0p-55 */
    106  1.1  christos P8 = -1.2518388626763144e-1;		/* -0x1006068cc0b97c.0p-55 */
    107  1.1  christos 
    108  1.1  christos static volatile const double zero = 0;
    109  1.1  christos 
    110  1.1  christos #define	INTERVALS	128
    111  1.1  christos #define	LOG2_INTERVALS	7
    112  1.1  christos #define	TSIZE		(INTERVALS + 1)
    113  1.1  christos #define	G(i)		(T[(i)].G)
    114  1.1  christos #define	F_hi(i)		(T[(i)].F_hi)
    115  1.1  christos #define	F_lo(i)		(T[(i)].F_lo)
    116  1.1  christos #define	ln2_hi		F_hi(TSIZE - 1)
    117  1.1  christos #define	ln2_lo		F_lo(TSIZE - 1)
    118  1.1  christos #define	E(i)		(U[(i)].E)
    119  1.1  christos #define	H(i)		(U[(i)].H)
    120  1.1  christos 
    121  1.1  christos static const struct {
    122  1.1  christos 	float	G;			/* 1/(1 + i/128) rounded to 8/9 bits */
    123  1.1  christos 	float	F_hi;			/* log(1 / G_i) rounded (see below) */
    124  1.1  christos 	double	F_lo;			/* next 53 bits for log(1 / G_i) */
    125  1.1  christos } T[TSIZE] = {
    126  1.1  christos 	/*
    127  1.1  christos 	 * ln2_hi and each F_hi(i) are rounded to a number of bits that
    128  1.1  christos 	 * makes F_hi(i) + dk*ln2_hi exact for all i and all dk.
    129  1.1  christos 	 *
    130  1.1  christos 	 * The last entry (for X just below 2) is used to define ln2_hi
    131  1.1  christos 	 * and ln2_lo, to ensure that F_hi(i) and F_lo(i) cancel exactly
    132  1.1  christos 	 * with dk*ln2_hi and dk*ln2_lo, respectively, when dk = -1.
    133  1.1  christos 	 * This is needed for accuracy when x is just below 1.  (To avoid
    134  1.1  christos 	 * special cases, such x are "reduced" strangely to X just below
    135  1.1  christos 	 * 2 and dk = -1, and then the exact cancellation is needed
    136  1.1  christos 	 * because any the error from any non-exactness would be too
    137  1.1  christos 	 * large).
    138  1.1  christos 	 *
    139  1.1  christos 	 * We want to share this table between double precision and ld80,
    140  1.1  christos 	 * so the relevant range of dk is the larger one of ld80
    141  1.1  christos 	 * ([-16445, 16383]) and the relevant exactness requirement is
    142  1.1  christos 	 * the stricter one of double precision.  The maximum number of
    143  1.1  christos 	 * bits in F_hi(i) that works is very dependent on i but has
    144  1.1  christos 	 * a minimum of 33.  We only need about 12 bits in F_hi(i) for
    145  1.1  christos 	 * it to provide enough extra precision in double precision (11
    146  1.1  christos 	 * more than that are required for ld80).
    147  1.1  christos 	 *
    148  1.1  christos 	 * We round F_hi(i) to 24 bits so that it can have type float,
    149  1.1  christos 	 * mainly to minimize the size of the table.  Using all 24 bits
    150  1.1  christos 	 * in a float for it automatically satisfies the above constraints.
    151  1.1  christos 	 */
    152  1.1  christos 	 { 0x800000.0p-23,  0,               0 },
    153  1.1  christos 	 { 0xfe0000.0p-24,  0x8080ac.0p-30, -0x14ee431dae6675.0p-84 },
    154  1.1  christos 	 { 0xfc0000.0p-24,  0x8102b3.0p-29, -0x1db29ee2d83718.0p-84 },
    155  1.1  christos 	 { 0xfa0000.0p-24,  0xc24929.0p-29,  0x1191957d173698.0p-83 },
    156  1.1  christos 	 { 0xf80000.0p-24,  0x820aec.0p-28,  0x13ce8888e02e79.0p-82 },
    157  1.1  christos 	 { 0xf60000.0p-24,  0xa33577.0p-28, -0x17a4382ce6eb7c.0p-82 },
    158  1.1  christos 	 { 0xf48000.0p-24,  0xbc42cb.0p-28, -0x172a21161a1076.0p-83 },
    159  1.1  christos 	 { 0xf30000.0p-24,  0xd57797.0p-28, -0x1e09de07cb9589.0p-82 },
    160  1.1  christos 	 { 0xf10000.0p-24,  0xf7518e.0p-28,  0x1ae1eec1b036c5.0p-91 },
    161  1.1  christos 	 { 0xef0000.0p-24,  0x8cb9df.0p-27, -0x1d7355325d560e.0p-81 },
    162  1.1  christos 	 { 0xed8000.0p-24,  0x999ec0.0p-27, -0x1f9f02d256d503.0p-82 },
    163  1.1  christos 	 { 0xec0000.0p-24,  0xa6988b.0p-27, -0x16fc0a9d12c17a.0p-83 },
    164  1.1  christos 	 { 0xea0000.0p-24,  0xb80698.0p-27,  0x15d581c1e8da9a.0p-81 },
    165  1.1  christos 	 { 0xe80000.0p-24,  0xc99af3.0p-27, -0x1535b3ba8f150b.0p-83 },
    166  1.1  christos 	 { 0xe70000.0p-24,  0xd273b2.0p-27,  0x163786f5251af0.0p-85 },
    167  1.1  christos 	 { 0xe50000.0p-24,  0xe442c0.0p-27,  0x1bc4b2368e32d5.0p-84 },
    168  1.1  christos 	 { 0xe38000.0p-24,  0xf1b83f.0p-27,  0x1c6090f684e676.0p-81 },
    169  1.1  christos 	 { 0xe20000.0p-24,  0xff448a.0p-27, -0x1890aa69ac9f42.0p-82 },
    170  1.1  christos 	 { 0xe08000.0p-24,  0x8673f6.0p-26,  0x1b9985194b6b00.0p-80 },
    171  1.1  christos 	 { 0xdf0000.0p-24,  0x8d515c.0p-26, -0x1dc08d61c6ef1e.0p-83 },
    172  1.1  christos 	 { 0xdd8000.0p-24,  0x943a9e.0p-26, -0x1f72a2dac729b4.0p-82 },
    173  1.1  christos 	 { 0xdc0000.0p-24,  0x9b2fe6.0p-26, -0x1fd4dfd3a0afb9.0p-80 },
    174  1.1  christos 	 { 0xda8000.0p-24,  0xa2315d.0p-26, -0x11b26121629c47.0p-82 },
    175  1.1  christos 	 { 0xd90000.0p-24,  0xa93f2f.0p-26,  0x1286d633e8e569.0p-81 },
    176  1.1  christos 	 { 0xd78000.0p-24,  0xb05988.0p-26,  0x16128eba936770.0p-84 },
    177  1.1  christos 	 { 0xd60000.0p-24,  0xb78094.0p-26,  0x16ead577390d32.0p-80 },
    178  1.1  christos 	 { 0xd50000.0p-24,  0xbc4c6c.0p-26,  0x151131ccf7c7b7.0p-81 },
    179  1.1  christos 	 { 0xd38000.0p-24,  0xc3890a.0p-26, -0x115e2cd714bd06.0p-80 },
    180  1.1  christos 	 { 0xd20000.0p-24,  0xcad2d7.0p-26, -0x1847f406ebd3b0.0p-82 },
    181  1.1  christos 	 { 0xd10000.0p-24,  0xcfb620.0p-26,  0x1c2259904d6866.0p-81 },
    182  1.1  christos 	 { 0xcf8000.0p-24,  0xd71653.0p-26,  0x1ece57a8d5ae55.0p-80 },
    183  1.1  christos 	 { 0xce0000.0p-24,  0xde843a.0p-26, -0x1f109d4bc45954.0p-81 },
    184  1.1  christos 	 { 0xcd0000.0p-24,  0xe37fde.0p-26,  0x1bc03dc271a74d.0p-81 },
    185  1.1  christos 	 { 0xcb8000.0p-24,  0xeb050c.0p-26, -0x1bf2badc0df842.0p-85 },
    186  1.1  christos 	 { 0xca0000.0p-24,  0xf29878.0p-26, -0x18efededd89fbe.0p-87 },
    187  1.1  christos 	 { 0xc90000.0p-24,  0xf7ad6f.0p-26,  0x1373ff977baa69.0p-81 },
    188  1.1  christos 	 { 0xc80000.0p-24,  0xfcc8e3.0p-26,  0x196766f2fb3283.0p-80 },
    189  1.1  christos 	 { 0xc68000.0p-24,  0x823f30.0p-25,  0x19bd076f7c434e.0p-79 },
    190  1.1  christos 	 { 0xc58000.0p-24,  0x84d52c.0p-25, -0x1a327257af0f46.0p-79 },
    191  1.1  christos 	 { 0xc40000.0p-24,  0x88bc74.0p-25,  0x113f23def19c5a.0p-81 },
    192  1.1  christos 	 { 0xc30000.0p-24,  0x8b5ae6.0p-25,  0x1759f6e6b37de9.0p-79 },
    193  1.1  christos 	 { 0xc20000.0p-24,  0x8dfccb.0p-25,  0x1ad35ca6ed5148.0p-81 },
    194  1.1  christos 	 { 0xc10000.0p-24,  0x90a22b.0p-25,  0x1a1d71a87deba4.0p-79 },
    195  1.1  christos 	 { 0xbf8000.0p-24,  0x94a0d8.0p-25, -0x139e5210c2b731.0p-80 },
    196  1.1  christos 	 { 0xbe8000.0p-24,  0x974f16.0p-25, -0x18f6ebcff3ed73.0p-81 },
    197  1.1  christos 	 { 0xbd8000.0p-24,  0x9a00f1.0p-25, -0x1aa268be39aab7.0p-79 },
    198  1.1  christos 	 { 0xbc8000.0p-24,  0x9cb672.0p-25, -0x14c8815839c566.0p-79 },
    199  1.1  christos 	 { 0xbb0000.0p-24,  0xa0cda1.0p-25,  0x1eaf46390dbb24.0p-81 },
    200  1.1  christos 	 { 0xba0000.0p-24,  0xa38c6e.0p-25,  0x138e20d831f698.0p-81 },
    201  1.1  christos 	 { 0xb90000.0p-24,  0xa64f05.0p-25, -0x1e8d3c41123616.0p-82 },
    202  1.1  christos 	 { 0xb80000.0p-24,  0xa91570.0p-25,  0x1ce28f5f3840b2.0p-80 },
    203  1.1  christos 	 { 0xb70000.0p-24,  0xabdfbb.0p-25, -0x186e5c0a424234.0p-79 },
    204  1.1  christos 	 { 0xb60000.0p-24,  0xaeadef.0p-25, -0x14d41a0b2a08a4.0p-83 },
    205  1.1  christos 	 { 0xb50000.0p-24,  0xb18018.0p-25,  0x16755892770634.0p-79 },
    206  1.1  christos 	 { 0xb40000.0p-24,  0xb45642.0p-25, -0x16395ebe59b152.0p-82 },
    207  1.1  christos 	 { 0xb30000.0p-24,  0xb73077.0p-25,  0x1abc65c8595f09.0p-80 },
    208  1.1  christos 	 { 0xb20000.0p-24,  0xba0ec4.0p-25, -0x1273089d3dad89.0p-79 },
    209  1.1  christos 	 { 0xb10000.0p-24,  0xbcf133.0p-25,  0x10f9f67b1f4bbf.0p-79 },
    210  1.1  christos 	 { 0xb00000.0p-24,  0xbfd7d2.0p-25, -0x109fab90486409.0p-80 },
    211  1.1  christos 	 { 0xaf0000.0p-24,  0xc2c2ac.0p-25, -0x1124680aa43333.0p-79 },
    212  1.1  christos 	 { 0xae8000.0p-24,  0xc439b3.0p-25, -0x1f360cc4710fc0.0p-80 },
    213  1.1  christos 	 { 0xad8000.0p-24,  0xc72afd.0p-25, -0x132d91f21d89c9.0p-80 },
    214  1.1  christos 	 { 0xac8000.0p-24,  0xca20a2.0p-25, -0x16bf9b4d1f8da8.0p-79 },
    215  1.1  christos 	 { 0xab8000.0p-24,  0xcd1aae.0p-25,  0x19deb5ce6a6a87.0p-81 },
    216  1.1  christos 	 { 0xaa8000.0p-24,  0xd0192f.0p-25,  0x1a29fb48f7d3cb.0p-79 },
    217  1.1  christos 	 { 0xaa0000.0p-24,  0xd19a20.0p-25,  0x1127d3c6457f9d.0p-81 },
    218  1.1  christos 	 { 0xa90000.0p-24,  0xd49f6a.0p-25, -0x1ba930e486a0ac.0p-81 },
    219  1.1  christos 	 { 0xa80000.0p-24,  0xd7a94b.0p-25, -0x1b6e645f31549e.0p-79 },
    220  1.1  christos 	 { 0xa70000.0p-24,  0xdab7d0.0p-25,  0x1118a425494b61.0p-80 },
    221  1.1  christos 	 { 0xa68000.0p-24,  0xdc40d5.0p-25,  0x1966f24d29d3a3.0p-80 },
    222  1.1  christos 	 { 0xa58000.0p-24,  0xdf566d.0p-25, -0x1d8e52eb2248f1.0p-82 },
    223  1.1  christos 	 { 0xa48000.0p-24,  0xe270ce.0p-25, -0x1ee370f96e6b68.0p-80 },
    224  1.1  christos 	 { 0xa40000.0p-24,  0xe3ffce.0p-25,  0x1d155324911f57.0p-80 },
    225  1.1  christos 	 { 0xa30000.0p-24,  0xe72179.0p-25, -0x1fe6e2f2f867d9.0p-80 },
    226  1.1  christos 	 { 0xa20000.0p-24,  0xea4812.0p-25,  0x1b7be9add7f4d4.0p-80 },
    227  1.1  christos 	 { 0xa18000.0p-24,  0xebdd3d.0p-25,  0x1b3cfb3f7511dd.0p-79 },
    228  1.1  christos 	 { 0xa08000.0p-24,  0xef0b5b.0p-25, -0x1220de1f730190.0p-79 },
    229  1.1  christos 	 { 0xa00000.0p-24,  0xf0a451.0p-25, -0x176364c9ac81cd.0p-80 },
    230  1.1  christos 	 { 0x9f0000.0p-24,  0xf3da16.0p-25,  0x1eed6b9aafac8d.0p-81 },
    231  1.1  christos 	 { 0x9e8000.0p-24,  0xf576e9.0p-25,  0x1d593218675af2.0p-79 },
    232  1.1  christos 	 { 0x9d8000.0p-24,  0xf8b47c.0p-25, -0x13e8eb7da053e0.0p-84 },
    233  1.1  christos 	 { 0x9d0000.0p-24,  0xfa553f.0p-25,  0x1c063259bcade0.0p-79 },
    234  1.1  christos 	 { 0x9c0000.0p-24,  0xfd9ac5.0p-25,  0x1ef491085fa3c1.0p-79 },
    235  1.1  christos 	 { 0x9b8000.0p-24,  0xff3f8c.0p-25,  0x1d607a7c2b8c53.0p-79 },
    236  1.1  christos 	 { 0x9a8000.0p-24,  0x814697.0p-24, -0x12ad3817004f3f.0p-78 },
    237  1.1  christos 	 { 0x9a0000.0p-24,  0x821b06.0p-24, -0x189fc53117f9e5.0p-81 },
    238  1.1  christos 	 { 0x990000.0p-24,  0x83c5f8.0p-24,  0x14cf15a048907b.0p-79 },
    239  1.1  christos 	 { 0x988000.0p-24,  0x849c7d.0p-24,  0x1cbb1d35fb8287.0p-78 },
    240  1.1  christos 	 { 0x978000.0p-24,  0x864ba6.0p-24,  0x1128639b814f9c.0p-78 },
    241  1.1  christos 	 { 0x970000.0p-24,  0x87244c.0p-24,  0x184733853300f0.0p-79 },
    242  1.1  christos 	 { 0x968000.0p-24,  0x87fdaa.0p-24,  0x109d23aef77dd6.0p-80 },
    243  1.1  christos 	 { 0x958000.0p-24,  0x89b293.0p-24, -0x1a81ef367a59de.0p-78 },
    244  1.1  christos 	 { 0x950000.0p-24,  0x8a8e20.0p-24, -0x121ad3dbb2f452.0p-78 },
    245  1.1  christos 	 { 0x948000.0p-24,  0x8b6a6a.0p-24, -0x1cfb981628af72.0p-79 },
    246  1.1  christos 	 { 0x938000.0p-24,  0x8d253a.0p-24, -0x1d21730ea76cfe.0p-79 },
    247  1.1  christos 	 { 0x930000.0p-24,  0x8e03c2.0p-24,  0x135cc00e566f77.0p-78 },
    248  1.1  christos 	 { 0x928000.0p-24,  0x8ee30d.0p-24, -0x10fcb5df257a26.0p-80 },
    249  1.1  christos 	 { 0x918000.0p-24,  0x90a3ee.0p-24, -0x16e171b15433d7.0p-79 },
    250  1.1  christos 	 { 0x910000.0p-24,  0x918587.0p-24, -0x1d050da07f3237.0p-79 },
    251  1.1  christos 	 { 0x908000.0p-24,  0x9267e7.0p-24,  0x1be03669a5268d.0p-79 },
    252  1.1  christos 	 { 0x8f8000.0p-24,  0x942f04.0p-24,  0x10b28e0e26c337.0p-79 },
    253  1.1  christos 	 { 0x8f0000.0p-24,  0x9513c3.0p-24,  0x1a1d820da57cf3.0p-78 },
    254  1.1  christos 	 { 0x8e8000.0p-24,  0x95f950.0p-24, -0x19ef8f13ae3cf1.0p-79 },
    255  1.1  christos 	 { 0x8e0000.0p-24,  0x96dfab.0p-24, -0x109e417a6e507c.0p-78 },
    256  1.1  christos 	 { 0x8d0000.0p-24,  0x98aed2.0p-24,  0x10d01a2c5b0e98.0p-79 },
    257  1.1  christos 	 { 0x8c8000.0p-24,  0x9997a2.0p-24, -0x1d6a50d4b61ea7.0p-78 },
    258  1.1  christos 	 { 0x8c0000.0p-24,  0x9a8145.0p-24,  0x1b3b190b83f952.0p-78 },
    259  1.1  christos 	 { 0x8b8000.0p-24,  0x9b6bbf.0p-24,  0x13a69fad7e7abe.0p-78 },
    260  1.1  christos 	 { 0x8b0000.0p-24,  0x9c5711.0p-24, -0x11cd12316f576b.0p-78 },
    261  1.1  christos 	 { 0x8a8000.0p-24,  0x9d433b.0p-24,  0x1c95c444b807a2.0p-79 },
    262  1.1  christos 	 { 0x898000.0p-24,  0x9f1e22.0p-24, -0x1b9c224ea698c3.0p-79 },
    263  1.1  christos 	 { 0x890000.0p-24,  0xa00ce1.0p-24,  0x125ca93186cf0f.0p-81 },
    264  1.1  christos 	 { 0x888000.0p-24,  0xa0fc80.0p-24, -0x1ee38a7bc228b3.0p-79 },
    265  1.1  christos 	 { 0x880000.0p-24,  0xa1ed00.0p-24, -0x1a0db876613d20.0p-78 },
    266  1.1  christos 	 { 0x878000.0p-24,  0xa2de62.0p-24,  0x193224e8516c01.0p-79 },
    267  1.1  christos 	 { 0x870000.0p-24,  0xa3d0a9.0p-24,  0x1fa28b4d2541ad.0p-79 },
    268  1.1  christos 	 { 0x868000.0p-24,  0xa4c3d6.0p-24,  0x1c1b5760fb4572.0p-78 },
    269  1.1  christos 	 { 0x858000.0p-24,  0xa6acea.0p-24,  0x1fed5d0f65949c.0p-80 },
    270  1.1  christos 	 { 0x850000.0p-24,  0xa7a2d4.0p-24,  0x1ad270c9d74936.0p-80 },
    271  1.1  christos 	 { 0x848000.0p-24,  0xa899ab.0p-24,  0x199ff15ce53266.0p-79 },
    272  1.1  christos 	 { 0x840000.0p-24,  0xa99171.0p-24,  0x1a19e15ccc45d2.0p-79 },
    273  1.1  christos 	 { 0x838000.0p-24,  0xaa8a28.0p-24, -0x121a14ec532b36.0p-80 },
    274  1.1  christos 	 { 0x830000.0p-24,  0xab83d1.0p-24,  0x1aee319980bff3.0p-79 },
    275  1.1  christos 	 { 0x828000.0p-24,  0xac7e6f.0p-24, -0x18ffd9e3900346.0p-80 },
    276  1.1  christos 	 { 0x820000.0p-24,  0xad7a03.0p-24, -0x1e4db102ce29f8.0p-80 },
    277  1.1  christos 	 { 0x818000.0p-24,  0xae768f.0p-24,  0x17c35c55a04a83.0p-81 },
    278  1.1  christos 	 { 0x810000.0p-24,  0xaf7415.0p-24,  0x1448324047019b.0p-78 },
    279  1.1  christos 	 { 0x808000.0p-24,  0xb07298.0p-24, -0x1750ee3915a198.0p-78 },
    280  1.1  christos 	 { 0x800000.0p-24,  0xb17218.0p-24, -0x105c610ca86c39.0p-81 },
    281  1.1  christos };
    282  1.1  christos 
    283  1.1  christos #ifdef USE_UTAB
    284  1.1  christos static const struct {
    285  1.1  christos 	float	H;			/* 1 + i/INTERVALS (exact) */
    286  1.1  christos 	float	E;			/* H(i) * G(i) - 1 (exact) */
    287  1.1  christos } U[TSIZE] = {
    288  1.1  christos 	 { 0x800000.0p-23,  0 },
    289  1.1  christos 	 { 0x810000.0p-23, -0x800000.0p-37 },
    290  1.1  christos 	 { 0x820000.0p-23, -0x800000.0p-35 },
    291  1.1  christos 	 { 0x830000.0p-23, -0x900000.0p-34 },
    292  1.1  christos 	 { 0x840000.0p-23, -0x800000.0p-33 },
    293  1.1  christos 	 { 0x850000.0p-23, -0xc80000.0p-33 },
    294  1.1  christos 	 { 0x860000.0p-23, -0xa00000.0p-36 },
    295  1.1  christos 	 { 0x870000.0p-23,  0x940000.0p-33 },
    296  1.1  christos 	 { 0x880000.0p-23,  0x800000.0p-35 },
    297  1.1  christos 	 { 0x890000.0p-23, -0xc80000.0p-34 },
    298  1.1  christos 	 { 0x8a0000.0p-23,  0xe00000.0p-36 },
    299  1.1  christos 	 { 0x8b0000.0p-23,  0x900000.0p-33 },
    300  1.1  christos 	 { 0x8c0000.0p-23, -0x800000.0p-35 },
    301  1.1  christos 	 { 0x8d0000.0p-23, -0xe00000.0p-33 },
    302  1.1  christos 	 { 0x8e0000.0p-23,  0x880000.0p-33 },
    303  1.1  christos 	 { 0x8f0000.0p-23, -0xa80000.0p-34 },
    304  1.1  christos 	 { 0x900000.0p-23, -0x800000.0p-35 },
    305  1.1  christos 	 { 0x910000.0p-23,  0x800000.0p-37 },
    306  1.1  christos 	 { 0x920000.0p-23,  0x900000.0p-35 },
    307  1.1  christos 	 { 0x930000.0p-23,  0xd00000.0p-35 },
    308  1.1  christos 	 { 0x940000.0p-23,  0xe00000.0p-35 },
    309  1.1  christos 	 { 0x950000.0p-23,  0xc00000.0p-35 },
    310  1.1  christos 	 { 0x960000.0p-23,  0xe00000.0p-36 },
    311  1.1  christos 	 { 0x970000.0p-23, -0x800000.0p-38 },
    312  1.1  christos 	 { 0x980000.0p-23, -0xc00000.0p-35 },
    313  1.1  christos 	 { 0x990000.0p-23, -0xd00000.0p-34 },
    314  1.1  christos 	 { 0x9a0000.0p-23,  0x880000.0p-33 },
    315  1.1  christos 	 { 0x9b0000.0p-23,  0xe80000.0p-35 },
    316  1.1  christos 	 { 0x9c0000.0p-23, -0x800000.0p-35 },
    317  1.1  christos 	 { 0x9d0000.0p-23,  0xb40000.0p-33 },
    318  1.1  christos 	 { 0x9e0000.0p-23,  0x880000.0p-34 },
    319  1.1  christos 	 { 0x9f0000.0p-23, -0xe00000.0p-35 },
    320  1.1  christos 	 { 0xa00000.0p-23,  0x800000.0p-33 },
    321  1.1  christos 	 { 0xa10000.0p-23, -0x900000.0p-36 },
    322  1.1  christos 	 { 0xa20000.0p-23, -0xb00000.0p-33 },
    323  1.1  christos 	 { 0xa30000.0p-23, -0xa00000.0p-36 },
    324  1.1  christos 	 { 0xa40000.0p-23,  0x800000.0p-33 },
    325  1.1  christos 	 { 0xa50000.0p-23, -0xf80000.0p-35 },
    326  1.1  christos 	 { 0xa60000.0p-23,  0x880000.0p-34 },
    327  1.1  christos 	 { 0xa70000.0p-23, -0x900000.0p-33 },
    328  1.1  christos 	 { 0xa80000.0p-23, -0x800000.0p-35 },
    329  1.1  christos 	 { 0xa90000.0p-23,  0x900000.0p-34 },
    330  1.1  christos 	 { 0xaa0000.0p-23,  0xa80000.0p-33 },
    331  1.1  christos 	 { 0xab0000.0p-23, -0xac0000.0p-34 },
    332  1.1  christos 	 { 0xac0000.0p-23, -0x800000.0p-37 },
    333  1.1  christos 	 { 0xad0000.0p-23,  0xf80000.0p-35 },
    334  1.1  christos 	 { 0xae0000.0p-23,  0xf80000.0p-34 },
    335  1.1  christos 	 { 0xaf0000.0p-23, -0xac0000.0p-33 },
    336  1.1  christos 	 { 0xb00000.0p-23, -0x800000.0p-33 },
    337  1.1  christos 	 { 0xb10000.0p-23, -0xb80000.0p-34 },
    338  1.1  christos 	 { 0xb20000.0p-23, -0x800000.0p-34 },
    339  1.1  christos 	 { 0xb30000.0p-23, -0xb00000.0p-35 },
    340  1.1  christos 	 { 0xb40000.0p-23, -0x800000.0p-35 },
    341  1.1  christos 	 { 0xb50000.0p-23, -0xe00000.0p-36 },
    342  1.1  christos 	 { 0xb60000.0p-23, -0x800000.0p-35 },
    343  1.1  christos 	 { 0xb70000.0p-23, -0xb00000.0p-35 },
    344  1.1  christos 	 { 0xb80000.0p-23, -0x800000.0p-34 },
    345  1.1  christos 	 { 0xb90000.0p-23, -0xb80000.0p-34 },
    346  1.1  christos 	 { 0xba0000.0p-23, -0x800000.0p-33 },
    347  1.1  christos 	 { 0xbb0000.0p-23, -0xac0000.0p-33 },
    348  1.1  christos 	 { 0xbc0000.0p-23,  0x980000.0p-33 },
    349  1.1  christos 	 { 0xbd0000.0p-23,  0xbc0000.0p-34 },
    350  1.1  christos 	 { 0xbe0000.0p-23,  0xe00000.0p-36 },
    351  1.1  christos 	 { 0xbf0000.0p-23, -0xb80000.0p-35 },
    352  1.1  christos 	 { 0xc00000.0p-23, -0x800000.0p-33 },
    353  1.1  christos 	 { 0xc10000.0p-23,  0xa80000.0p-33 },
    354  1.1  christos 	 { 0xc20000.0p-23,  0x900000.0p-34 },
    355  1.1  christos 	 { 0xc30000.0p-23, -0x800000.0p-35 },
    356  1.1  christos 	 { 0xc40000.0p-23, -0x900000.0p-33 },
    357  1.1  christos 	 { 0xc50000.0p-23,  0x820000.0p-33 },
    358  1.1  christos 	 { 0xc60000.0p-23,  0x800000.0p-38 },
    359  1.1  christos 	 { 0xc70000.0p-23, -0x820000.0p-33 },
    360  1.1  christos 	 { 0xc80000.0p-23,  0x800000.0p-33 },
    361  1.1  christos 	 { 0xc90000.0p-23, -0xa00000.0p-36 },
    362  1.1  christos 	 { 0xca0000.0p-23, -0xb00000.0p-33 },
    363  1.1  christos 	 { 0xcb0000.0p-23,  0x840000.0p-34 },
    364  1.1  christos 	 { 0xcc0000.0p-23, -0xd00000.0p-34 },
    365  1.1  christos 	 { 0xcd0000.0p-23,  0x800000.0p-33 },
    366  1.1  christos 	 { 0xce0000.0p-23, -0xe00000.0p-35 },
    367  1.1  christos 	 { 0xcf0000.0p-23,  0xa60000.0p-33 },
    368  1.1  christos 	 { 0xd00000.0p-23, -0x800000.0p-35 },
    369  1.1  christos 	 { 0xd10000.0p-23,  0xb40000.0p-33 },
    370  1.1  christos 	 { 0xd20000.0p-23, -0x800000.0p-35 },
    371  1.1  christos 	 { 0xd30000.0p-23,  0xaa0000.0p-33 },
    372  1.1  christos 	 { 0xd40000.0p-23, -0xe00000.0p-35 },
    373  1.1  christos 	 { 0xd50000.0p-23,  0x880000.0p-33 },
    374  1.1  christos 	 { 0xd60000.0p-23, -0xd00000.0p-34 },
    375  1.1  christos 	 { 0xd70000.0p-23,  0x9c0000.0p-34 },
    376  1.1  christos 	 { 0xd80000.0p-23, -0xb00000.0p-33 },
    377  1.1  christos 	 { 0xd90000.0p-23, -0x800000.0p-38 },
    378  1.1  christos 	 { 0xda0000.0p-23,  0xa40000.0p-33 },
    379  1.1  christos 	 { 0xdb0000.0p-23, -0xdc0000.0p-34 },
    380  1.1  christos 	 { 0xdc0000.0p-23,  0xc00000.0p-35 },
    381  1.1  christos 	 { 0xdd0000.0p-23,  0xca0000.0p-33 },
    382  1.1  christos 	 { 0xde0000.0p-23, -0xb80000.0p-34 },
    383  1.1  christos 	 { 0xdf0000.0p-23,  0xd00000.0p-35 },
    384  1.1  christos 	 { 0xe00000.0p-23,  0xc00000.0p-33 },
    385  1.1  christos 	 { 0xe10000.0p-23, -0xf40000.0p-34 },
    386  1.1  christos 	 { 0xe20000.0p-23,  0x800000.0p-37 },
    387  1.1  christos 	 { 0xe30000.0p-23,  0x860000.0p-33 },
    388  1.1  christos 	 { 0xe40000.0p-23, -0xc80000.0p-33 },
    389  1.1  christos 	 { 0xe50000.0p-23, -0xa80000.0p-34 },
    390  1.1  christos 	 { 0xe60000.0p-23,  0xe00000.0p-36 },
    391  1.1  christos 	 { 0xe70000.0p-23,  0x880000.0p-33 },
    392  1.1  christos 	 { 0xe80000.0p-23, -0xe00000.0p-33 },
    393  1.1  christos 	 { 0xe90000.0p-23, -0xfc0000.0p-34 },
    394  1.1  christos 	 { 0xea0000.0p-23, -0x800000.0p-35 },
    395  1.1  christos 	 { 0xeb0000.0p-23,  0xe80000.0p-35 },
    396  1.1  christos 	 { 0xec0000.0p-23,  0x900000.0p-33 },
    397  1.1  christos 	 { 0xed0000.0p-23,  0xe20000.0p-33 },
    398  1.1  christos 	 { 0xee0000.0p-23, -0xac0000.0p-33 },
    399  1.1  christos 	 { 0xef0000.0p-23, -0xc80000.0p-34 },
    400  1.1  christos 	 { 0xf00000.0p-23, -0x800000.0p-35 },
    401  1.1  christos 	 { 0xf10000.0p-23,  0x800000.0p-35 },
    402  1.1  christos 	 { 0xf20000.0p-23,  0xb80000.0p-34 },
    403  1.1  christos 	 { 0xf30000.0p-23,  0x940000.0p-33 },
    404  1.1  christos 	 { 0xf40000.0p-23,  0xc80000.0p-33 },
    405  1.1  christos 	 { 0xf50000.0p-23, -0xf20000.0p-33 },
    406  1.1  christos 	 { 0xf60000.0p-23, -0xc80000.0p-33 },
    407  1.1  christos 	 { 0xf70000.0p-23, -0xa20000.0p-33 },
    408  1.1  christos 	 { 0xf80000.0p-23, -0x800000.0p-33 },
    409  1.1  christos 	 { 0xf90000.0p-23, -0xc40000.0p-34 },
    410  1.1  christos 	 { 0xfa0000.0p-23, -0x900000.0p-34 },
    411  1.1  christos 	 { 0xfb0000.0p-23, -0xc80000.0p-35 },
    412  1.1  christos 	 { 0xfc0000.0p-23, -0x800000.0p-35 },
    413  1.1  christos 	 { 0xfd0000.0p-23, -0x900000.0p-36 },
    414  1.1  christos 	 { 0xfe0000.0p-23, -0x800000.0p-37 },
    415  1.1  christos 	 { 0xff0000.0p-23, -0x800000.0p-39 },
    416  1.1  christos 	 { 0x800000.0p-22,  0 },
    417  1.1  christos };
    418  1.1  christos #endif /* USE_UTAB */
    419  1.1  christos 
    420  1.1  christos #ifdef STRUCT_RETURN
    421  1.1  christos #define	RETURN1(rp, v) do {	\
    422  1.1  christos 	(rp)->hi = (v);		\
    423  1.1  christos 	(rp)->lo_set = 0;	\
    424  1.1  christos 	return;			\
    425  1.1  christos } while (0)
    426  1.1  christos 
    427  1.1  christos #define	RETURN2(rp, h, l) do {	\
    428  1.1  christos 	(rp)->hi = (h);		\
    429  1.1  christos 	(rp)->lo = (l);		\
    430  1.1  christos 	(rp)->lo_set = 1;	\
    431  1.1  christos 	return;			\
    432  1.1  christos } while (0)
    433  1.1  christos 
    434  1.1  christos struct ld {
    435  1.1  christos 	long double hi;
    436  1.1  christos 	long double lo;
    437  1.1  christos 	int	lo_set;
    438  1.1  christos };
    439  1.1  christos #else
    440  1.1  christos #define	RETURN1(rp, v)	RETURNF(v)
    441  1.1  christos #define	RETURN2(rp, h, l)	RETURNI((h) + (l))
    442  1.1  christos #endif
    443  1.1  christos 
    444  1.1  christos #ifdef STRUCT_RETURN
    445  1.1  christos static inline __always_inline void
    446  1.1  christos k_logl(long double x, struct ld *rp)
    447  1.1  christos #else
    448  1.1  christos long double
    449  1.1  christos logl(long double x)
    450  1.1  christos #endif
    451  1.1  christos {
    452  1.1  christos 	long double d, dk, val_hi, val_lo, z;
    453  1.1  christos 	uint64_t ix, lx;
    454  1.1  christos 	int i, k;
    455  1.1  christos 	uint16_t hx;
    456  1.1  christos 
    457  1.1  christos 	EXTRACT_LDBL80_WORDS(hx, lx, x);
    458  1.1  christos 	k = -16383;
    459  1.1  christos #if 0 /* Hard to do efficiently.  Don't do it until we support all modes. */
    460  1.1  christos 	if (x == 1)
    461  1.1  christos 		RETURN1(rp, 0);		/* log(1) = +0 in all rounding modes */
    462  1.1  christos #endif
    463  1.1  christos 	if (hx == 0 || hx >= 0x8000) {	/* zero, negative or subnormal? */
    464  1.1  christos 		if (((hx & 0x7fff) | lx) == 0)
    465  1.1  christos 			RETURN1(rp, -1 / zero);	/* log(+-0) = -Inf */
    466  1.1  christos 		if (hx != 0)
    467  1.1  christos 			/* log(neg or [pseudo-]NaN) = qNaN: */
    468  1.1  christos 			RETURN1(rp, (x - x) / zero);
    469  1.1  christos 		x *= 0x1.0p65;		/* subnormal; scale up x */
    470  1.1  christos 					/* including pseudo-subnormals */
    471  1.1  christos 		EXTRACT_LDBL80_WORDS(hx, lx, x);
    472  1.1  christos 		k = -16383 - 65;
    473  1.1  christos 	} else if (hx >= 0x7fff || (lx & 0x8000000000000000ULL) == 0)
    474  1.1  christos 		RETURN1(rp, x + x);	/* log(Inf or NaN) = Inf or qNaN */
    475  1.1  christos 					/* log(pseudo-Inf) = qNaN */
    476  1.1  christos 					/* log(pseudo-NaN) = qNaN */
    477  1.1  christos 					/* log(unnormal) = qNaN */
    478  1.1  christos #ifndef STRUCT_RETURN
    479  1.1  christos 	ENTERI();
    480  1.1  christos #endif
    481  1.1  christos 	k += hx;
    482  1.1  christos 	ix = lx & 0x7fffffffffffffffULL;
    483  1.1  christos 	dk = k;
    484  1.1  christos 
    485  1.1  christos 	/* Scale x to be in [1, 2). */
    486  1.1  christos 	SET_LDBL_EXPSIGN(x, 0x3fff);
    487  1.1  christos 
    488  1.1  christos 	/* 0 <= i <= INTERVALS: */
    489  1.1  christos #define	L2I	(64 - LOG2_INTERVALS)
    490  1.1  christos 	i = (ix + (1LL << (L2I - 2))) >> (L2I - 1);
    491  1.1  christos 
    492  1.1  christos 	/*
    493  1.1  christos 	 * -0.005280 < d < 0.004838.  In particular, the infinite-
    494  1.1  christos 	 * precision |d| is <= 2**-7.  Rounding of G(i) to 8 bits
    495  1.1  christos 	 * ensures that d is representable without extra precision for
    496  1.1  christos 	 * this bound on |d| (since when this calculation is expressed
    497  1.1  christos 	 * as x*G(i)-1, the multiplication needs as many extra bits as
    498  1.1  christos 	 * G(i) has and the subtraction cancels 8 bits).  But for
    499  1.1  christos 	 * most i (107 cases out of 129), the infinite-precision |d|
    500  1.1  christos 	 * is <= 2**-8.  G(i) is rounded to 9 bits for such i to give
    501  1.1  christos 	 * better accuracy (this works by improving the bound on |d|,
    502  1.1  christos 	 * which in turn allows rounding to 9 bits in more cases).
    503  1.1  christos 	 * This is only important when the original x is near 1 -- it
    504  1.1  christos 	 * lets us avoid using a special method to give the desired
    505  1.1  christos 	 * accuracy for such x.
    506  1.1  christos 	 */
    507  1.1  christos 	if (/*CONSTCOND*/0)
    508  1.1  christos 		/*NOTREACHED*/
    509  1.1  christos 		d = x * G(i) - 1;
    510  1.1  christos 	else {
    511  1.1  christos #ifdef USE_UTAB
    512  1.1  christos 		d = (x - H(i)) * G(i) + E(i);
    513  1.1  christos #else
    514  1.1  christos 		long double x_hi, x_lo;
    515  1.1  christos 		float fx_hi;
    516  1.1  christos 
    517  1.1  christos 		/*
    518  1.1  christos 		 * Split x into x_hi + x_lo to calculate x*G(i)-1 exactly.
    519  1.1  christos 		 * G(i) has at most 9 bits, so the splitting point is not
    520  1.1  christos 		 * critical.
    521  1.1  christos 		 */
    522  1.1  christos 		SET_FLOAT_WORD(fx_hi, (lx >> 40) | 0x3f800000);
    523  1.1  christos 		x_hi = fx_hi;
    524  1.1  christos 		x_lo = x - x_hi;
    525  1.1  christos 		d = x_hi * G(i) - 1 + x_lo * G(i);
    526  1.1  christos #endif
    527  1.1  christos 	}
    528  1.1  christos 
    529  1.1  christos 	/*
    530  1.1  christos 	 * Our algorithm depends on exact cancellation of F_lo(i) and
    531  1.1  christos 	 * F_hi(i) with dk*ln_2_lo and dk*ln2_hi when k is -1 and i is
    532  1.1  christos 	 * at the end of the table.  This and other technical complications
    533  1.1  christos 	 * make it difficult to avoid the double scaling in (dk*ln2) *
    534  1.1  christos 	 * log(base) for base != e without losing more accuracy and/or
    535  1.1  christos 	 * efficiency than is gained.
    536  1.1  christos 	 */
    537  1.1  christos 	z = d * d;
    538  1.1  christos 	val_lo = z * d * z * (z * (d * P8 + P7) + (d * P6 + P5)) +
    539  1.1  christos 	    (F_lo(i) + dk * ln2_lo + z * d * (d * P4 + P3)) + z * P2;
    540  1.1  christos 	val_hi = d;
    541  1.1  christos #ifdef DEBUG
    542  1.1  christos 	if (fetestexcept(FE_UNDERFLOW))
    543  1.1  christos 		breakpoint();
    544  1.1  christos #endif
    545  1.1  christos 
    546  1.1  christos 	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
    547  1.1  christos 	RETURN2(rp, val_hi, val_lo);
    548  1.1  christos }
    549  1.1  christos 
    550  1.1  christos long double
    551  1.1  christos log1pl(long double x)
    552  1.1  christos {
    553  1.1  christos 	long double d, d_hi, d_lo, dk, f_lo, val_hi, val_lo, z;
    554  1.1  christos 	long double f_hi, twopminusk;
    555  1.1  christos 	uint64_t ix, lx;
    556  1.1  christos 	int i, k;
    557  1.1  christos 	int16_t ax, hx;
    558  1.1  christos 
    559  1.1  christos 	EXTRACT_LDBL80_WORDS(hx, lx, x);
    560  1.1  christos 	if (hx < 0x3fff) {		/* x < 1, or x neg NaN */
    561  1.1  christos 		ax = hx & 0x7fff;
    562  1.1  christos 		if (ax >= 0x3fff) {	/* x <= -1, or x neg NaN */
    563  1.1  christos 			if (ax == 0x3fff && lx == 0x8000000000000000ULL)
    564  1.1  christos 				RETURNF(-1 / zero);	/* log1p(-1) = -Inf */
    565  1.1  christos 			/* log1p(x < 1, or x [pseudo-]NaN) = qNaN: */
    566  1.1  christos 			RETURNF((x - x) / (x - x));
    567  1.1  christos 		}
    568  1.1  christos 		if (ax <= 0x3fbe) {	/* |x| < 2**-64 */
    569  1.1  christos 			if ((int)x == 0)
    570  1.1  christos 				RETURNF(x);	/* x with inexact if x != 0 */
    571  1.1  christos 		}
    572  1.1  christos 		f_hi = 1;
    573  1.1  christos 		f_lo = x;
    574  1.1  christos 	} else if (hx >= 0x7fff) {	/* x +Inf or non-neg NaN */
    575  1.1  christos 		RETURNF(x + x);		/* log1p(Inf or NaN) = Inf or qNaN */
    576  1.1  christos 					/* log1p(pseudo-Inf) = qNaN */
    577  1.1  christos 					/* log1p(pseudo-NaN) = qNaN */
    578  1.1  christos 					/* log1p(unnormal) = qNaN */
    579  1.1  christos 	} else if (hx < 0x407f) {	/* 1 <= x < 2**128 */
    580  1.1  christos 		f_hi = x;
    581  1.1  christos 		f_lo = 1;
    582  1.1  christos 	} else {			/* 2**128 <= x < +Inf */
    583  1.1  christos 		f_hi = x;
    584  1.1  christos 		f_lo = 0;		/* avoid underflow of the P5 term */
    585  1.1  christos 	}
    586  1.1  christos 	ENTERI();
    587  1.1  christos 	x = f_hi + f_lo;
    588  1.1  christos 	f_lo = (f_hi - x) + f_lo;
    589  1.1  christos 
    590  1.1  christos 	EXTRACT_LDBL80_WORDS(hx, lx, x);
    591  1.1  christos 	k = -16383;
    592  1.1  christos 
    593  1.1  christos 	k += hx;
    594  1.1  christos 	ix = lx & 0x7fffffffffffffffULL;
    595  1.1  christos 	dk = k;
    596  1.1  christos 
    597  1.1  christos 	SET_LDBL_EXPSIGN(x, 0x3fff);
    598  1.1  christos 	twopminusk = 1;
    599  1.1  christos 	SET_LDBL_EXPSIGN(twopminusk, 0x7ffe - (hx & 0x7fff));
    600  1.1  christos 	f_lo *= twopminusk;
    601  1.1  christos 
    602  1.1  christos 	i = (ix + (1LL << (L2I - 2))) >> (L2I - 1);
    603  1.1  christos 
    604  1.1  christos 	/*
    605  1.1  christos 	 * x*G(i)-1 (with a reduced x) can be represented exactly, as
    606  1.1  christos 	 * above, but now we need to evaluate the polynomial on d =
    607  1.1  christos 	 * (x+f_lo)*G(i)-1 and extra precision is needed for that.
    608  1.1  christos 	 * Since x+x_lo is a hi+lo decomposition and subtracting 1
    609  1.1  christos 	 * doesn't lose too many bits, an inexact calculation for
    610  1.1  christos 	 * f_lo*G(i) is good enough.
    611  1.1  christos 	 */
    612  1.1  christos 	if (/*CONSTCOND*/0)
    613  1.1  christos 		/*NOTREACHED*/
    614  1.1  christos 		d_hi = x * G(i) - 1;
    615  1.1  christos 	else {
    616  1.1  christos #ifdef USE_UTAB
    617  1.1  christos 		d_hi = (x - H(i)) * G(i) + E(i);
    618  1.1  christos #else
    619  1.1  christos 		long double x_hi, x_lo;
    620  1.1  christos 		float fx_hi;
    621  1.1  christos 
    622  1.1  christos 		SET_FLOAT_WORD(fx_hi, (lx >> 40) | 0x3f800000);
    623  1.1  christos 		x_hi = fx_hi;
    624  1.1  christos 		x_lo = x - x_hi;
    625  1.1  christos 		d_hi = x_hi * G(i) - 1 + x_lo * G(i);
    626  1.1  christos #endif
    627  1.1  christos 	}
    628  1.1  christos 	d_lo = f_lo * G(i);
    629  1.1  christos 
    630  1.1  christos 	/*
    631  1.1  christos 	 * This is _2sumF(d_hi, d_lo) inlined.  The condition
    632  1.1  christos 	 * (d_hi == 0 || |d_hi| >= |d_lo|) for using _2sumF() is not
    633  1.1  christos 	 * always satisifed, so it is not clear that this works, but
    634  1.1  christos 	 * it works in practice.  It works even if it gives a wrong
    635  1.1  christos 	 * normalized d_lo, since |d_lo| > |d_hi| implies that i is
    636  1.1  christos 	 * nonzero and d is tiny, so the F(i) term dominates d_lo.
    637  1.1  christos 	 * In float precision:
    638  1.1  christos 	 * (By exhaustive testing, the worst case is d_hi = 0x1.bp-25.
    639  1.1  christos 	 * And if d is only a little tinier than that, we would have
    640  1.1  christos 	 * another underflow problem for the P3 term; this is also ruled
    641  1.1  christos 	 * out by exhaustive testing.)
    642  1.1  christos 	 */
    643  1.1  christos 	d = d_hi + d_lo;
    644  1.1  christos 	d_lo = d_hi - d + d_lo;
    645  1.1  christos 	d_hi = d;
    646  1.1  christos 
    647  1.1  christos 	z = d * d;
    648  1.1  christos 	val_lo = z * d * z * (z * (d * P8 + P7) + (d * P6 + P5)) +
    649  1.1  christos 	    (F_lo(i) + dk * ln2_lo + d_lo + z * d * (d * P4 + P3)) + z * P2;
    650  1.1  christos 	val_hi = d_hi;
    651  1.1  christos #ifdef DEBUG
    652  1.1  christos 	if (fetestexcept(FE_UNDERFLOW))
    653  1.1  christos 		breakpoint();
    654  1.1  christos #endif
    655  1.1  christos 
    656  1.1  christos 	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
    657  1.1  christos 	RETURNI(val_hi + val_lo);
    658  1.1  christos }
    659  1.1  christos 
    660  1.1  christos #ifdef STRUCT_RETURN
    661  1.1  christos 
    662  1.1  christos long double
    663  1.1  christos logl(long double x)
    664  1.1  christos {
    665  1.1  christos 	struct ld r;
    666  1.1  christos 
    667  1.1  christos 	ENTERI();
    668  1.1  christos 	k_logl(x, &r);
    669  1.1  christos 	RETURNSPI(&r);
    670  1.1  christos }
    671  1.1  christos 
    672  1.1  christos /* Use macros since GCC < 8 rejects static const expressions in initializers. */
    673  1.1  christos #define	invln10_hi	4.3429448190317999e-1	/*  0x1bcb7b1526e000.0p-54 */
    674  1.1  christos #define	invln10_lo	7.1842412889749798e-14	/*  0x1438ca9aadd558.0p-96 */
    675  1.1  christos #define	invln2_hi	1.4426950408887933e0	/*  0x171547652b8000.0p-52 */
    676  1.1  christos #define	invln2_lo	1.7010652264631490e-13	/*  0x17f0bbbe87fed0.0p-95 */
    677  1.1  christos /* Let the compiler pre-calculate this sum to avoid FE_INEXACT at run time. */
    678  1.1  christos static const double invln10_lo_plus_hi = invln10_lo + invln10_hi;
    679  1.1  christos static const double invln2_lo_plus_hi = invln2_lo + invln2_hi;
    680  1.1  christos 
    681  1.1  christos long double
    682  1.1  christos log10l(long double x)
    683  1.1  christos {
    684  1.1  christos 	struct ld r;
    685  1.1  christos 	long double hi, lo;
    686  1.1  christos 
    687  1.1  christos 	ENTERI();
    688  1.1  christos 	k_logl(x, &r);
    689  1.1  christos 	if (!r.lo_set)
    690  1.1  christos 		RETURNI(r.hi);
    691  1.1  christos 	_2sumF(r.hi, r.lo);
    692  1.1  christos 	hi = (float)r.hi;
    693  1.1  christos 	lo = r.lo + (r.hi - hi);
    694  1.1  christos 	RETURNI(invln10_hi * hi +
    695  1.1  christos 	    (invln10_lo_plus_hi * lo + invln10_lo * hi));
    696  1.1  christos }
    697  1.1  christos 
    698  1.1  christos long double
    699  1.1  christos log2l(long double x)
    700  1.1  christos {
    701  1.1  christos 	struct ld r;
    702  1.1  christos 	long double hi, lo;
    703  1.1  christos 
    704  1.1  christos 	ENTERI();
    705  1.1  christos 	k_logl(x, &r);
    706  1.1  christos 	if (!r.lo_set)
    707  1.1  christos 		RETURNI(r.hi);
    708  1.1  christos 	_2sumF(r.hi, r.lo);
    709  1.1  christos 	hi = (float)r.hi;
    710  1.1  christos 	lo = r.lo + (r.hi - hi);
    711  1.1  christos 	RETURNI(invln2_hi * hi +
    712  1.1  christos 	    (invln2_lo_plus_hi * lo + invln2_lo * hi));
    713  1.1  christos }
    714  1.1  christos 
    715  1.1  christos #endif /* STRUCT_RETURN */
    716