/*-
* SPDX-License-Identifier: BSD-3-Clause
*
* Copyright (c) 1992, 1993
* The Regents of the University of California. All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
/*
* See bsdsrc/b_log.c for implementation details.
*
* bsdrc/b_log.c converted to long double by Steven G. Kargl.
*/
#define N 128
/*
* Coefficients in the polynomial approximation of log(1+f/F).
* Domain of x is [0,1./256] with 2**(-84.48) precision.
*/
static const union ieee_ext_u
a1u = LD80C(0xaaaaaaaaaaaaaaab, -4, 8.33333333333333333356e-02L),
a2u = LD80C(0xcccccccccccccd29, -7, 1.25000000000000000781e-02L),
a3u = LD80C(0x9249249241ed3764, -9, 2.23214285711721994134e-03L),
a4u = LD80C(0xe38e959e1e7e01cf, -12, 4.34030476540000360640e-04L);
#define A1 (a1u.extu_ld)
#define A2 (a2u.extu_ld)
#define A3 (a3u.extu_ld)
#define A4 (a4u.extu_ld)
/*
* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128.
* Used for generation of extend precision logarithms.
* The constant 35184372088832 is 2^45, so the divide is exact.
* It ensures correct reading of logF_head, even for inaccurate
* decimal-to-binary conversion routines. (Everybody gets the
* right answer for integers less than 2^53.)
* Values for log(F) were generated using error < 10^-57 absolute
* with the bc -l package.
*/
/*
* Argument reduction: 1 <= g < 2; x/2^m = g;
* y = F*(1 + f/F) for |f| <= 2^-8
*/
g = frexpl(x, &m);
g *= 2;
m--;
if (m == DBL_MIN_EXP - 1) {
j = ilogbl(g);
m += j;
g = ldexpl(g, -j);
}
j = N * (g - 1) + 0.5L;
F = (1.L / N) * j + 1;
f = g - F;
g = 1 / (2 * F + f);
u = 2 * f * g;
v = u * u;
q = u * v * (A1 + v * (A2 + v * (A3 + v * A4)));
if (m | j) {
u1 = u + 513;
u1 -= 513;
} else {
u1 = (float)u;
}
u2 = (2 * (f - F * u1) - u1 * f) * g;
u1 += m * (long double)logF_head[N] + logF_head[j];
u2 += logF_tail[j];
u2 += q;
u2 += logF_tail[N] * m;
r.a = (float)(u1 + u2); /* Only difference is here. */
r.b = (u1 - r.a) + u2;
return (r);
}
