dawids/lang-bootstrap
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Patched musl with log functions

Browse source 
I pulled in the implementation of a few log functions from musl 0.8.7
because lua needed log2 to compile. Since lua needed only log2 I
commented out the ones that didn't immediatly compile like log1p.c.
This commit is contained in:
Dawid Sobczak 2023-07-06 11:48:59 +01:00
parent 36a81aa0b9
commit 2ae045cf8a
21 changed files with 1773 additions and 27 deletions
Download patch file Download diff file Expand all files Collapse all files

10
05/Makefile
View file

@ -2,7 +2,7 @@ TCCDIR=tcc-0.9.27
TCC=$(TCCDIR)/tcc
TCC0=$(TCC)0
TCCINST=tcc-bootstrap
all: out04 a.out $(TCCDIR)/lib/libtcc1.a
all: out04 a.out tcc
in04: *.b ../04a/out04
../04a/out04 main.b in04
out04: in04 ../04/out03
@ -28,15 +28,15 @@ $(TCCDIR)/lib/libtcc1.a: $(TCC0) $(TCCDIR)/lib/*.[cS]
$(TCC0) -c $(TCCDIR)/lib/libtcc1.c -o $(TCCDIR)/lib/libtcc1.o
$(TCC0) -ar $(TCCDIR)/lib/libtcc1.a $(TCCDIR)/lib/*.o
musl: tcc-files
$(MAKE) -C musl-0.6.0
$(MAKE) -j8 -C musl-0.6.0
$(MAKE) -C musl-0.6.0 install
tcc-files: $(TCCDIR)/lib/libtcc1.a $(TCCDIR)/include/*.h
mkdir -p $(TCCINST)/include
cp -r $(TCCDIR)/include/*.h $(TCCINST)/include/
cp -r $(TCCDIR)/lib/libtcc1.a $(TCCINST)/
$(TCC)1: $(TCC0) $(TCCINST)/libtcc1.a
cd $(TCCDIR) && ./tcc0 -nostdinc -nostdlib -B ../tcc-boostrap -L../musl-bootstrap/lib -lc -I ../musl-bootstrap/include tcc.c -o tcc1
tcc: $(TCC)1
$(TCC): $(TCC0) musl
cd $(TCCDIR) && ./tcc0 -nostdinc -nostdlib -B ../tcc-boostrap -I ../musl-bootstrap/include tcc.c ../musl-bootstrap/lib/*.[oa] -o tcc
tcc: $(TCC)
$(TCC)2: $(TCC)1
cd $(TCCDIR) && ./tcc1 tcc.c -o tcc2

3
05/musl-0.6.0/Makefile
View file

@ -70,6 +70,9 @@ include/bits/alltypes.h: include/bits/alltypes.h.sh
%.o: $(ARCH)/%.s
$(CC) $(CFLAGS) $(INC) -c -o $@ $<
%.o: $(ARCH)/%.s
$(CC) $(CFLAGS) $(INC) -c -o $@ $<
%.o: %.c $(GENH)
$(CC) $(CFLAGS) $(INC) -c -o $@ $<

19
05/musl-0.6.0/include/dlfcn.h
View file

@ -1,19 +0,0 @@
#ifndef _DLFCN_H
#define _DLFCN_H
#define RTLD_LAZY 0x10000
#define RTLD_NOW 0x20000
#define RTLD_GLOBAL 0x40000
#define RTLD_LOCAL 0x80000
#if 1
#define RTLD_NEXT ((void *) -1l)
#define RTLD_DEFAULT ((void *) 0)
#endif
int dlclose(void *);
char *dlerror(void);
void *dlopen(const char *, int);
void *dlsym(void *, const char *);
#endif

94
05/musl-0.6.0/src/math/__log1p.h Normal file
View file

@ -0,0 +1,94 @@
/* origin: FreeBSD /usr/src/lib/msun/src/k_log.h */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* __log1p(f):
* Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)].
*
* The following describes the overall strategy for computing
* logarithms in base e. The argument reduction and adding the final
* term of the polynomial are done by the caller for increased accuracy
* when different bases are used.
*
* Method :
* 1. Argument Reduction: find k and f such that
* x = 2^k * (1+f),
* where sqrt(2)/2 < 1+f < sqrt(2) .
*
* 2. Approximation of log(1+f).
* Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
* = 2s + s*R
* We use a special Reme algorithm on [0,0.1716] to generate
* a polynomial of degree 14 to approximate R The maximum error
* of this polynomial approximation is bounded by 2**-58.45. In
* other words,
* 2 4 6 8 10 12 14
* R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
* (the values of Lg1 to Lg7 are listed in the program)
* and
* | 2 14 | -58.45
* | Lg1*s +...+Lg7*s - R(z) | <= 2
* | |
* Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
* In order to guarantee error in log below 1ulp, we compute log
* by
* log(1+f) = f - s*(f - R) (if f is not too large)
* log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
*
* 3. Finally, log(x) = k*ln2 + log(1+f).
* = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
* Here ln2 is split into two floating point number:
* ln2_hi + ln2_lo,
* where n*ln2_hi is always exact for |n| < 2000.
*
* Special cases:
* log(x) is NaN with signal if x < 0 (including -INF) ;
* log(+INF) is +INF; log(0) is -INF with signal;
* log(NaN) is that NaN with no signal.
*
* Accuracy:
* according to an error analysis, the error is always less than
* 1 ulp (unit in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
static const double
Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
/*
* We always inline __log1p(), since doing so produces a
* substantial performance improvement (~40% on amd64).
*/
static inline double __log1p(double f)
{
double hfsq,s,z,R,w,t1,t2;
s = f/(2.0+f);
z = s*s;
w = z*z;
t1= w*(Lg2+w*(Lg4+w*Lg6));
t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
R = t2+t1;
hfsq = 0.5*f*f;
return s*(hfsq+R);
}

35
05/musl-0.6.0/src/math/__log1pf.h Normal file
View file

@ -0,0 +1,35 @@
/* origin: FreeBSD /usr/src/lib/msun/src/k_logf.h */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* See comments in __log1p.h.
*/
/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
static const float
Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
static inline float __log1pf(float f)
{
float hfsq,s,z,R,w,t1,t2;
s = f/(2.0f + f);
z = s*s;
w = z*z;
t1 = w*(Lg2+w*Lg4);
t2 = z*(Lg1+w*Lg3);
R = t2+t1;
hfsq = 0.5f * f * f;
return s*(hfsq+R);
}

142
05/musl-0.6.0/src/math/log.c Normal file
View file

@ -0,0 +1,142 @@
/* origin: FreeBSD /usr/src/lib/msun/src/e_log.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* log(x)
* Return the logrithm of x
*
* Method :
* 1. Argument Reduction: find k and f such that
* x = 2^k * (1+f),
* where sqrt(2)/2 < 1+f < sqrt(2) .
*
* 2. Approximation of log(1+f).
* Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
* = 2s + s*R
* We use a special Remez algorithm on [0,0.1716] to generate
* a polynomial of degree 14 to approximate R The maximum error
* of this polynomial approximation is bounded by 2**-58.45. In
* other words,
* 2 4 6 8 10 12 14
* R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
* (the values of Lg1 to Lg7 are listed in the program)
* and
* | 2 14 | -58.45
* | Lg1*s +...+Lg7*s - R(z) | <= 2
* | |
* Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
* In order to guarantee error in log below 1ulp, we compute log
* by
* log(1+f) = f - s*(f - R) (if f is not too large)
* log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
*
* 3. Finally, log(x) = k*ln2 + log(1+f).
* = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
* Here ln2 is split into two floating point number:
* ln2_hi + ln2_lo,
* where n*ln2_hi is always exact for |n| < 2000.
*
* Special cases:
* log(x) is NaN with signal if x < 0 (including -INF) ;
* log(+INF) is +INF; log(0) is -INF with signal;
* log(NaN) is that NaN with no signal.
*
* Accuracy:
* according to an error analysis, the error is always less than
* 1 ulp (unit in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
#include "math_private.h"
#include "math.h"
#include <stdint.h>
static const double
ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
double log(double x)
{
double hfsq,f,s,z,R,w,t1,t2,dk;
int32_t k,hx,i,j;
uint32_t lx;
EXTRACT_WORDS(hx, lx, x);
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
return -two54/zero; /* log(+-0)=-inf */
if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 54;
x *= two54;
GET_HIGH_WORD(hx,x);
}
if (hx >= 0x7ff00000)
return x+x;
k += (hx>>20) - 1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
SET_HIGH_WORD(x, hx|(i^0x3ff00000)); /* normalize x or x/2 */
k += i>>20;
f = x - 1.0;
if ((0x000fffff&(2+hx)) < 3) { /* -2**-20 <= f < 2**-20 */
if (f == zero) {
if (k == 0) {
return zero;
}
dk = (double)k;
return dk*ln2_hi + dk*ln2_lo;
}
R = f*f*(0.5-0.33333333333333333*f);
if (k == 0)
return f - R;
dk = (double)k;
return dk*ln2_hi - ((R-dk*ln2_lo)-f);
}
s = f/(2.0+f);
dk = (double)k;
z = s*s;
i = hx - 0x6147a;
w = z*z;
j = 0x6b851 - hx;
t1 = w*(Lg2+w*(Lg4+w*Lg6));
t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
i |= j;
R = t2 + t1;
if (i > 0) {
hfsq = 0.5*f*f;
if (k == 0)
return f - (hfsq-s*(hfsq+R));
return dk*ln2_hi - ((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
} else {
if (k == 0)
return f - s*(f-R);
return dk*ln2_hi - ((s*(f-R)-dk*ln2_lo)-f);
}
}

86
05/musl-0.6.0/src/math/log10.c Normal file
View file

@ -0,0 +1,86 @@
/* origin: FreeBSD /usr/src/lib/msun/src/e_log10.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* Return the base 10 logarithm of x. See e_log.c and k_log.h for most
* comments.
*
* log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2)
* in not-quite-routine extra precision.
*/
#include "math_private.h"
#include "math.h"
#include "__log1p.h"
#include <stdint.h>
static const double
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
ivln10hi = 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */
ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
static const double zero = 0.0;
double log10(double x)
{
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
int32_t i,k,hx;
uint32_t lx;
EXTRACT_WORDS(hx, lx, x);
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
return -two54/zero; /* log(+-0)=-inf */
if (hx<0)
return (x-x)/zero; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 54;
x *= two54;
GET_HIGH_WORD(hx, x);
}
if (hx >= 0x7ff00000)
return x+x;
if (hx == 0x3ff00000 && lx == 0)
return zero; /* log(1) = +0 */
k += (hx>>20) - 1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
SET_HIGH_WORD(x, hx|(i^0x3ff00000)); /* normalize x or x/2 */
k += i>>20;
y = (double)k;
f = x - 1.0;
hfsq = 0.5*f*f;
r = __log1p(f);
/* See log2.c for details. */
hi = f - hfsq;
SET_LOW_WORD(hi, 0);
lo = (f - hi) - hfsq + r;
val_hi = hi*ivln10hi;
y2 = y*log10_2hi;
val_lo = y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi;
/*
* Extra precision in for adding y*log10_2hi is not strictly needed
* since there is no very large cancellation near x = sqrt(2) or
* x = 1/sqrt(2), but we do it anyway since it costs little on CPUs
* with some parallelism and it reduces the error for many args.
*/
w = y2 + val_hi;
val_lo += (y2 - w) + val_hi;
val_hi = w;
return val_lo + val_hi;
}

72
05/musl-0.6.0/src/math/log10f.c Normal file
View file

@ -0,0 +1,72 @@
/* origin: FreeBSD /usr/src/lib/msun/src/e_log10f.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* See comments in log10.c.
*/
#include "math.h"
#include "__log1pf.h"
#include <stdint.h>
static const float
two25 = 3.3554432000e+07, /* 0x4c000000 */
ivln10hi = 4.3432617188e-01, /* 0x3ede6000 */
ivln10lo = -3.1689971365e-05, /* 0xb804ead9 */
log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */
log10_2lo = 7.9034151668e-07; /* 0x355427db */
static const float zero = 0.0;
float log10f(float x)
{
float f,hfsq,hi,lo,r,y;
int32_t i,k,hx;
GET_FLOAT_WORD(hx, x);
k = 0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff) == 0)
return -two25/zero; /* log(+-0)=-inf */
if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 25;
x *= two25;
GET_FLOAT_WORD(hx, x);
}
if (hx >= 0x7f800000)
return x+x;
if (hx == 0x3f800000)
return zero; /* log(1) = +0 */
k += (hx>>23) - 127;
hx &= 0x007fffff;
i = (hx+(0x4afb0d))&0x800000;
SET_FLOAT_WORD(x, hx|(i^0x3f800000)); /* normalize x or x/2 */
k += i>>23;
y = (float)k;
f = x - 1.0f;
hfsq = 0.5f * f * f;
r = __log1pf(f);
// FIXME
// /* See log2f.c and log2.c for details. */
// if (sizeof(float_t) > sizeof(float))
// return (r - hfsq + f) * ((float_t)ivln10lo + ivln10hi) +
// y * ((float_t)log10_2lo + log10_2hi);
hi = f - hfsq;
GET_FLOAT_WORD(hx, hi);
SET_FLOAT_WORD(hi, hx&0xfffff000);
lo = (f - hi) - hfsq + r;
return y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi +
hi*ivln10hi + y*log10_2hi;
}

186
05/musl-0.6.0/src/math/log10l.c Normal file
View file

@ -0,0 +1,186 @@
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_log10l.c */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/*
* Common logarithm, long double precision
*
*
* SYNOPSIS:
*
* long double x, y, log10l();
*
* y = log10l( x );
*
*
* DESCRIPTION:
*
* Returns the base 10 logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 30000 9.0e-20 2.6e-20
* IEEE exp(+-10000) 30000 6.0e-20 2.3e-20
*
* In the tests over the interval exp(+-10000), the logarithms
* of the random arguments were uniformly distributed over
* [-10000, +10000].
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns MINLOG
* log domain: x < 0; returns MINLOG
*/
#include "math.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double log10l(long double x)
{
return log10(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.2e-22
*/
static long double P[] = {
4.9962495940332550844739E-1L,
1.0767376367209449010438E1L,
7.7671073698359539859595E1L,
2.5620629828144409632571E2L,
4.2401812743503691187826E2L,
3.4258224542413922935104E2L,
1.0747524399916215149070E2L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0,*/
2.3479774160285863271658E1L,
1.9444210022760132894510E2L,
7.7952888181207260646090E2L,
1.6911722418503949084863E3L,
2.0307734695595183428202E3L,
1.2695660352705325274404E3L,
3.2242573199748645407652E2L,
};
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.16e-22
*/
static long double R[4] = {
1.9757429581415468984296E-3L,
-7.1990767473014147232598E-1L,
1.0777257190312272158094E1L,
-3.5717684488096787370998E1L,
};
static long double S[4] = {
/* 1.00000000000000000000E0L,*/
-2.6201045551331104417768E1L,
1.9361891836232102174846E2L,
-4.2861221385716144629696E2L,
};
/* log10(2) */
#define L102A 0.3125L
#define L102B -1.1470004336018804786261e-2L
/* log10(e) */
#define L10EA 0.5L
#define L10EB -6.5705518096748172348871e-2L
#define SQRTH 0.70710678118654752440L
long double log10l(long double x)
{
long double y;
volatile long double z;
int e;
if (isnan(x))
return x;
if(x <= 0.0L) {
if(x == 0.0L)
return -1.0L / (x - x);
return (x - x) / (x - x);
}
if (x == INFINITY)
return INFINITY;
/* separate mantissa from exponent */
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
x = frexpl(x, &e);
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
} else { /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x*x;
y = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3));
goto done;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH) {
e -= 1;
x = ldexpl(x, 1) - 1.0L; /* 2x - 1 */
} else {
x = x - 1.0L;
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
y = y - ldexpl(z, -1); /* -0.5x^2 + ... */
done:
/* Multiply log of fraction by log10(e)
* and base 2 exponent by log10(2).
*
* ***CAUTION***
*
* This sequence of operations is critical and it may
* be horribly defeated by some compiler optimizers.
*/
z = y * (L10EB);
z += x * (L10EB);
z += e * (L102B);
z += y * (L10EA);
z += x * (L10EA);
z += e * (L102A);
return z;
}
#endif

172
05/musl-0.6.0/src/math/log1p.c Normal file
View file

@ -0,0 +1,172 @@
/* origin: FreeBSD /usr/src/lib/msun/src/s_log1p.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* double log1p(double x)
*
* Method :
* 1. Argument Reduction: find k and f such that
* 1+x = 2^k * (1+f),
* where sqrt(2)/2 < 1+f < sqrt(2) .
*
* Note. If k=0, then f=x is exact. However, if k!=0, then f
* may not be representable exactly. In that case, a correction
* term is need. Let u=1+x rounded. Let c = (1+x)-u, then
* log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
* and add back the correction term c/u.
* (Note: when x > 2**53, one can simply return log(x))
*
* 2. Approximation of log1p(f).
* Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
* = 2s + s*R
* We use a special Reme algorithm on [0,0.1716] to generate
* a polynomial of degree 14 to approximate R The maximum error
* of this polynomial approximation is bounded by 2**-58.45. In
* other words,
* 2 4 6 8 10 12 14
* R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s
* (the values of Lp1 to Lp7 are listed in the program)
* and
* | 2 14 | -58.45
* | Lp1*s +...+Lp7*s - R(z) | <= 2
* | |
* Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
* In order to guarantee error in log below 1ulp, we compute log
* by
* log1p(f) = f - (hfsq - s*(hfsq+R)).
*
* 3. Finally, log1p(x) = k*ln2 + log1p(f).
* = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
* Here ln2 is split into two floating point number:
* ln2_hi + ln2_lo,
* where n*ln2_hi is always exact for |n| < 2000.
*
* Special cases:
* log1p(x) is NaN with signal if x < -1 (including -INF) ;
* log1p(+INF) is +INF; log1p(-1) is -INF with signal;
* log1p(NaN) is that NaN with no signal.
*
* Accuracy:
* according to an error analysis, the error is always less than
* 1 ulp (unit in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*
* Note: Assuming log() return accurate answer, the following
* algorithm can be used to compute log1p(x) to within a few ULP:
*
* u = 1+x;
* if(u==1.0) return x ; else
* return log(u)*(x/(u-1.0));
*
* See HP-15C Advanced Functions Handbook, p.193.
*/
#include "math.h"
#include <stdint.h>
static const double
ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
Lp1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
Lp2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
Lp3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
Lp4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
Lp5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
// double log1p(double x)
// {
// double hfsq,f,c,s,z,R,u;
// int32_t k,hx,hu,ax;
// GET_HIGH_WORD(hx, x);
// ax = hx & 0x7fffffff;
// k = 1;
// if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */
// if (ax >= 0x3ff00000) { /* x <= -1.0 */
// if (x == -1.0)
// return -two54/zero; /* log1p(-1)=+inf */
// return (x-x)/(x-x); /* log1p(x<-1)=NaN */
// }
// if (ax < 0x3e200000) { /* |x| < 2**-29 */
// /* raise inexact */
// if (two54 + x > zero && ax < 0x3c900000) /* |x| < 2**-54 */
// return x;
// return x - x*x*0.5;
// }
// if (hx > 0 || hx <= (int32_t)0xbfd2bec4) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
// k = 0;
// f = x;
// hu = 1;
// }
// }
// if (hx >= 0x7ff00000)
// return x+x;
// if (k != 0) {
// if (hx < 0x43400000) {
// STRICT_ASSIGN(double, u, 1.0 + x);
// GET_HIGH_WORD(hu, u);
// k = (hu>>20) - 1023;
// c = k > 0 ? 1.0-(u-x) : x-(u-1.0); /* correction term */
// c /= u;
// } else {
// u = x;
// GET_HIGH_WORD(hu,u);
// k = (hu>>20) - 1023;
// c = 0;
// }
// hu &= 0x000fffff;
// /*
// * The approximation to sqrt(2) used in thresholds is not
// * critical. However, the ones used above must give less
// * strict bounds than the one here so that the k==0 case is
// * never reached from here, since here we have committed to
// * using the correction term but don't use it if k==0.
// */
// if (hu < 0x6a09e) { /* u ~< sqrt(2) */
// SET_HIGH_WORD(u, hu|0x3ff00000); /* normalize u */
// } else {
// k += 1;
// SET_HIGH_WORD(u, hu|0x3fe00000); /* normalize u/2 */
// hu = (0x00100000-hu)>>2;
// }
// f = u - 1.0;
// }
// hfsq = 0.5*f*f;
// if (hu == 0) { /* |f| < 2**-20 */
// if (f == zero) {
// if(k == 0)
// return zero;
// c += k*ln2_lo;
// return k*ln2_hi + c;
// }
// R = hfsq*(1.0 - 0.66666666666666666*f);
// if (k == 0)
// return f - R;
// return k*ln2_hi - ((R-(k*ln2_lo+c))-f);
// }
// s = f/(2.0+f);
// z = s*s;
// R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
// if (k == 0)
// return f - (hfsq-s*(hfsq+R));
// return k*ln2_hi - ((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
// }

112
05/musl-0.6.0/src/math/log1pf.c Normal file
View file

@ -0,0 +1,112 @@
/* origin: FreeBSD /usr/src/lib/msun/src/s_log1pf.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "math.h"
#include <stdint.h>
static const float
ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
two25 = 3.355443200e+07, /* 0x4c000000 */
Lp1 = 6.6666668653e-01, /* 3F2AAAAB */
Lp2 = 4.0000000596e-01, /* 3ECCCCCD */
Lp3 = 2.8571429849e-01, /* 3E924925 */
Lp4 = 2.2222198546e-01, /* 3E638E29 */
Lp5 = 1.8183572590e-01, /* 3E3A3325 */
Lp6 = 1.5313838422e-01, /* 3E1CD04F */
Lp7 = 1.4798198640e-01; /* 3E178897 */
static const float zero = 0.0;
// float log1pf(float x)
// {
// float hfsq,f,c,s,z,R,u;
// int32_t k,hx,hu,ax;
// GET_FLOAT_WORD(hx, x);
// ax = hx & 0x7fffffff;
// k = 1;
// if (hx < 0x3ed413d0) { /* 1+x < sqrt(2)+ */
// if (ax >= 0x3f800000) { /* x <= -1.0 */
// if (x == -1.0f)
// return -two25/zero; /* log1p(-1)=+inf */
// return (x-x)/(x-x); /* log1p(x<-1)=NaN */
// }
// if (ax < 0x38000000) { /* |x| < 2**-15 */
// /* raise inexact */
// if (two25 + x > zero && ax < 0x33800000) /* |x| < 2**-24 */
// return x;
// return x - x*x*0.5f;
// }
// if (hx > 0 || hx <= (int32_t)0xbe95f619) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
// k = 0;
// f = x;
// hu = 1;
// }
// }
// if (hx >= 0x7f800000)
// return x+x;
// if (k != 0) {
// if (hx < 0x5a000000) {
// STRICT_ASSIGN(float, u, 1.0f + x);
// GET_FLOAT_WORD(hu, u);
// k = (hu>>23) - 127;
// /* correction term */
// c = k > 0 ? 1.0f-(u-x) : x-(u-1.0f);
// c /= u;
// } else {
// u = x;
// GET_FLOAT_WORD(hu,u);
// k = (hu>>23) - 127;
// c = 0;
// }
// hu &= 0x007fffff;
// /*
// * The approximation to sqrt(2) used in thresholds is not
// * critical. However, the ones used above must give less
// * strict bounds than the one here so that the k==0 case is
// * never reached from here, since here we have committed to
// * using the correction term but don't use it if k==0.
// */
// if (hu < 0x3504f4) { /* u < sqrt(2) */
// SET_FLOAT_WORD(u, hu|0x3f800000); /* normalize u */
// } else {
// k += 1;
// SET_FLOAT_WORD(u, hu|0x3f000000); /* normalize u/2 */
// hu = (0x00800000-hu)>>2;
// }
// f = u - 1.0f;
// }
// hfsq = 0.5f * f * f;
// if (hu == 0) { /* |f| < 2**-20 */
// if (f == zero) {
// if (k == 0)
// return zero;
// c += k*ln2_lo;
// return k*ln2_hi+c;
// }
// R = hfsq*(1.0f - 0.66666666666666666f * f);
// if (k == 0)
// return f - R;
// return k*ln2_hi - ((R-(k*ln2_lo+c))-f);
// }
// s = f/(2.0f + f);
// z = s*s;
// R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
// if (k == 0)
// return f - (hfsq-s*(hfsq+R));
// return k*ln2_hi - ((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
// }

176
05/musl-0.6.0/src/math/log1pl.c Normal file
View file

@ -0,0 +1,176 @@
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_log1pl.c */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/*
* Relative error logarithm
* Natural logarithm of 1+x, long double precision
*
*
* SYNOPSIS:
*
* long double x, y, log1pl();
*
* y = log1pl( x );
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of 1+x.
*
* The argument 1+x is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z^3 P(z)/Q(z).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1.0, 9.0 100000 8.2e-20 2.5e-20
*
* ERROR MESSAGES:
*
* log singularity: x-1 = 0; returns -INFINITY
* log domain: x-1 < 0; returns NAN
*/
#include "math.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double log1pl(long double x)
{
return log1p(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 2.32e-20
*/
static long double P[] = {
4.5270000862445199635215E-5L,
4.9854102823193375972212E-1L,
6.5787325942061044846969E0L,
2.9911919328553073277375E1L,
6.0949667980987787057556E1L,
5.7112963590585538103336E1L,
2.0039553499201281259648E1L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0,*/
1.5062909083469192043167E1L,
8.3047565967967209469434E1L,
2.2176239823732856465394E2L,
3.0909872225312059774938E2L,
2.1642788614495947685003E2L,
6.0118660497603843919306E1L,
};
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.16e-22
*/
static long double R[4] = {
1.9757429581415468984296E-3L,
-7.1990767473014147232598E-1L,
1.0777257190312272158094E1L,
-3.5717684488096787370998E1L,
};
static long double S[4] = {
/* 1.00000000000000000000E0L,*/
-2.6201045551331104417768E1L,
1.9361891836232102174846E2L,
-4.2861221385716144629696E2L,
};
static const long double C1 = 6.9314575195312500000000E-1L;
static const long double C2 = 1.4286068203094172321215E-6L;
#define SQRTH 0.70710678118654752440L
long double log1pl(long double xm1)
{
long double x, y, z;
int e;
if (isnan(xm1))
return xm1;
if (xm1 == INFINITY)
return xm1;
if (xm1 == 0.0)
return xm1;
x = xm1 + 1.0L;
/* Test for domain errors. */
if (x <= 0.0L) {
if (x == 0.0L)
return -INFINITY;
return NAN;
}
/* Separate mantissa from exponent.
Use frexp so that denormal numbers will be handled properly. */
x = frexpl(x, &e);
/* logarithm using log(x) = z + z^3 P(z)/Q(z),
where z = 2(x-1)/x+1) */
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
} else { /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x*x;
z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3));
z = z + e * C2;
z = z + x;
z = z + e * C1;
return z;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH) {
e -= 1;
if (e != 0)
x = 2.0 * x - 1.0L;
else
x = xm1;
} else {
if (e != 0)
x = x - 1.0L;
else
x = xm1;
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6));
y = y + e * C2;
z = y - 0.5 * z;
z = z + x;
z = z + e * C1;
return z;
}
#endif

109
05/musl-0.6.0/src/math/log2.c Normal file
View file

@ -0,0 +1,109 @@
/* origin: FreeBSD /usr/src/lib/msun/src/e_log2.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* Return the base 2 logarithm of x. See log.c and __log1p.h for most
* comments.
*
* This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel,
* then does the combining and scaling steps
* log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k
* in not-quite-routine extra precision.
*/
#include "math_private.h"
#include "math.h"
#include "__log1p.h"
#include <stdint.h>
static const double
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
static const double zero = 0.0;
double log2(double x)
{
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
int32_t i,k,hx;
uint32_t lx;
EXTRACT_WORDS(hx, lx, x);
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
return -two54/zero; /* log(+-0)=-inf */
if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 54;
x *= two54;
GET_HIGH_WORD(hx, x);
}
if (hx >= 0x7ff00000)
return x+x;
if (hx == 0x3ff00000 && lx == 0)
return zero; /* log(1) = +0 */
k += (hx>>20) - 1023;
hx &= 0x000fffff;
i = (hx+0x95f64) & 0x100000;
SET_HIGH_WORD(x, hx|(i^0x3ff00000)); /* normalize x or x/2 */
k += i>>20;
y = (double)k;
f = x - 1.0;
hfsq = 0.5*f*f;
r = __log1p(f);
/*
* f-hfsq must (for args near 1) be evaluated in extra precision
* to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
* This is fairly efficient since f-hfsq only depends on f, so can
* be evaluated in parallel with R. Not combining hfsq with R also
* keeps R small (though not as small as a true `lo' term would be),
* so that extra precision is not needed for terms involving R.
*
* Compiler bugs involving extra precision used to break Dekker's
* theorem for spitting f-hfsq as hi+lo, unless double_t was used
* or the multi-precision calculations were avoided when double_t
* has extra precision. These problems are now automatically
* avoided as a side effect of the optimization of combining the
* Dekker splitting step with the clear-low-bits step.
*
* y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
* precision to avoid a very large cancellation when x is very near
* these values. Unlike the above cancellations, this problem is
* specific to base 2. It is strange that adding +-1 is so much
* harder than adding +-ln2 or +-log10_2.
*
* This uses Dekker's theorem to normalize y+val_hi, so the
* compiler bugs are back in some configurations, sigh. And I
* don't want to used double_t to avoid them, since that gives a
* pessimization and the support for avoiding the pessimization
* is not yet available.
*
* The multi-precision calculations for the multiplications are
* routine.
*/
hi = f - hfsq;
SET_LOW_WORD(hi, 0);
lo = (f - hi) - hfsq + r;
val_hi = hi*ivln2hi;
val_lo = (lo+hi)*ivln2lo + lo*ivln2hi;
/* spadd(val_hi, val_lo, y), except for not using double_t: */
w = y + val_hi;
val_lo += (y - w) + val_hi;
val_hi = w;
return val_lo + val_hi;
}

82
05/musl-0.6.0/src/math/log2f.c Normal file
View file

@ -0,0 +1,82 @@
/* origin: FreeBSD /usr/src/lib/msun/src/e_log2f.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* See comments in log2.c.
*/
#include "math.h"
#include "__log1pf.h"
#include <stdint.h>
static const float
two25 = 3.3554432000e+07, /* 0x4c000000 */
ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */
ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */
static const float zero = 0.0;
float log2f(float x)
{
float f,hfsq,hi,lo,r,y;
int32_t i,k,hx;
GET_FLOAT_WORD(hx, x);
k = 0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff) == 0)
return -two25/zero; /* log(+-0)=-inf */
if (hx < 0)
return (x-x)/zero; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 25;
x *= two25;
GET_FLOAT_WORD(hx, x);
}
if (hx >= 0x7f800000)
return x+x;
if (hx == 0x3f800000)
return zero; /* log(1) = +0 */
k += (hx>>23) - 127;
hx &= 0x007fffff;
i = (hx+(0x4afb0d))&0x800000;
SET_FLOAT_WORD(x, hx|(i^0x3f800000)); /* normalize x or x/2 */
k += i>>23;
y = (float)k;
f = x - 1.0f;
hfsq = 0.5f * f * f;
r = __log1pf(f);
/*
* We no longer need to avoid falling into the multi-precision
* calculations due to compiler bugs breaking Dekker's theorem.
* Keep avoiding this as an optimization. See log2.c for more
* details (some details are here only because the optimization
* is not yet available in double precision).
*
* Another compiler bug turned up. With gcc on i386,
* (ivln2lo + ivln2hi) would be evaluated in float precision
* despite runtime evaluations using double precision. So we
* must cast one of its terms to float_t. This makes the whole
* expression have type float_t, so return is forced to waste
* time clobbering its extra precision.
*/
// FIXME
// if (sizeof(float_t) > sizeof(float))
// return (r - hfsq + f) * ((float_t)ivln2lo + ivln2hi) + y;
hi = f - hfsq;
GET_FLOAT_WORD(hx,hi);
SET_FLOAT_WORD(hi,hx&0xfffff000);
lo = (f - hi) - hfsq + r;
return (lo+hi)*ivln2lo + lo*ivln2hi + hi*ivln2hi + y;
}

182
05/musl-0.6.0/src/math/log2l.c Normal file
View file

@ -0,0 +1,182 @@
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_log2l.c */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/*
* Base 2 logarithm, long double precision
*
*
* SYNOPSIS:
*
* long double x, y, log2l();
*
* y = log2l( x );
*
*
* DESCRIPTION:
*
* Returns the base 2 logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the (natural)
* logarithm of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 30000 9.8e-20 2.7e-20
* IEEE exp(+-10000) 70000 5.4e-20 2.3e-20
*
* In the tests over the interval exp(+-10000), the logarithms
* of the random arguments were uniformly distributed over
* [-10000, +10000].
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns -INFINITY
* log domain: x < 0; returns NAN
*/
#include "math.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double log2l(long double x)
{
return log2(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.2e-22
*/
static long double P[] = {
4.9962495940332550844739E-1L,
1.0767376367209449010438E1L,
7.7671073698359539859595E1L,
2.5620629828144409632571E2L,
4.2401812743503691187826E2L,
3.4258224542413922935104E2L,
1.0747524399916215149070E2L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0,*/
2.3479774160285863271658E1L,
1.9444210022760132894510E2L,
7.7952888181207260646090E2L,
1.6911722418503949084863E3L,
2.0307734695595183428202E3L,
1.2695660352705325274404E3L,
3.2242573199748645407652E2L,
};
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.16e-22
*/
static long double R[4] = {
1.9757429581415468984296E-3L,
-7.1990767473014147232598E-1L,
1.0777257190312272158094E1L,
-3.5717684488096787370998E1L,
};
static long double S[4] = {
/* 1.00000000000000000000E0L,*/
-2.6201045551331104417768E1L,
1.9361891836232102174846E2L,
-4.2861221385716144629696E2L,
};
/* log2(e) - 1 */
#define LOG2EA 4.4269504088896340735992e-1L
#define SQRTH 0.70710678118654752440L
long double log2l(long double x)
{
volatile long double z;
long double y;
int e;
if (isnan(x))
return x;
if (x == INFINITY)
return x;
if (x <= 0.0L) {
if (x == 0.0L)
return -INFINITY;
return NAN;
}
/* separate mantissa from exponent */
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
x = frexpl(x, &e);
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
} else { /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x*x;
y = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3));
goto done;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH) {
e -= 1;
x = ldexpl(x, 1) - 1.0L; /* 2x - 1 */
} else {
x = x - 1.0L;
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
y = y - ldexpl(z, -1); /* -0.5x^2 + ... */
done:
/* Multiply log of fraction by log2(e)
* and base 2 exponent by 1
*
* ***CAUTION***
*
* This sequence of operations is critical and it may
* be horribly defeated by some compiler optimizers.
*/
z = y * LOG2EA;
z += x * LOG2EA;
z += y;
z += x;
z += e;
return z;
}
#endif

20
05/musl-0.6.0/src/math/logb.c Normal file
View file

@ -0,0 +1,20 @@
#include <limits.h>
#include "math.h"
/*
special cases:
logb(+-0) = -inf
logb(+-inf) = +inf
logb(nan) = nan
these are calculated at runtime to raise fp exceptions
*/
double logb(double x) {
int i = ilogb(x);
if (i == FP_ILOGB0)
return -1.0/fabs(x);
if (i == FP_ILOGBNAN || i == INT_MAX)
return x * x;
return i;
}

12
05/musl-0.6.0/src/math/logbf.c Normal file
View file

@ -0,0 +1,12 @@
#include <limits.h>
#include "math.h"
float logbf(float x) {
int i = ilogbf(x);
if (i == FP_ILOGB0)
return -1.0f/fabsf(x);
if (i == FP_ILOGBNAN || i == INT_MAX)
return x * x;
return i;
}

19
05/musl-0.6.0/src/math/logbl.c Normal file
View file

@ -0,0 +1,19 @@
#include <limits.h>
#include "math.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double logbl(long double x)
{
return logb(x);
}
#else
long double logbl(long double x)
{
int i = ilogbl(x);
if (i == FP_ILOGB0)
return -1.0/fabsl(x);
if (i == FP_ILOGBNAN || i == INT_MAX)
return x * x;
return i;
}
#endif

90
05/musl-0.6.0/src/math/logf.c Normal file
View file

@ -0,0 +1,90 @@
/* origin: FreeBSD /usr/src/lib/msun/src/e_logf.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "math.h"
#include <stdint.h>
static const float
ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
two25 = 3.355443200e+07, /* 0x4c000000 */
/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
static const float zero = 0.0;
float logf(float x)
{
float hfsq,f,s,z,R,w,t1,t2,dk;
int32_t k,ix,i,j;
GET_FLOAT_WORD(ix, x);
k = 0;
if (ix < 0x00800000) { /* x < 2**-126 */
if ((ix & 0x7fffffff) == 0)
return -two25/zero; /* log(+-0)=-inf */
if (ix < 0)
return (x-x)/zero; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 25;
x *= two25;
GET_FLOAT_WORD(ix, x);
}
if (ix >= 0x7f800000)
return x+x;
k += (ix>>23) - 127;
ix &= 0x007fffff;
i = (ix + (0x95f64<<3)) & 0x800000;
SET_FLOAT_WORD(x, ix|(i^0x3f800000)); /* normalize x or x/2 */
k += i>>23;
f = x - 1.0f;
if ((0x007fffff & (0x8000 + ix)) < 0xc000) { /* -2**-9 <= f < 2**-9 */
if (f == zero) {
if (k == 0)
return zero;
dk = (float)k;
return dk*ln2_hi + dk*ln2_lo;
}
R = f*f*(0.5f - 0.33333333333333333f*f);
if (k == 0)
return f-R;
dk = (float)k;
return dk*ln2_hi - ((R-dk*ln2_lo)-f);
}
s = f/(2.0f + f);
dk = (float)k;
z = s*s;
i = ix-(0x6147a<<3);
w = z*z;
j = (0x6b851<<3)-ix;
t1= w*(Lg2+w*Lg4);
t2= z*(Lg1+w*Lg3);
i |= j;
R = t2 + t1;
if (i > 0) {
hfsq = 0.5f * f * f;
if (k == 0)
return f - (hfsq-s*(hfsq+R));
return dk*ln2_hi - ((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
} else {
if (k == 0)
return f - s*(f-R);
return dk*ln2_hi - ((s*(f-R)-dk*ln2_lo)-f);
}
}

174
05/musl-0.6.0/src/math/logl.c Normal file
View file

@ -0,0 +1,174 @@
/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_logl.c */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
/*
* Natural logarithm, long double precision
*
*
* SYNOPSIS:
*
* long double x, y, logl();
*
* y = logl( x );
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 150000 8.71e-20 2.75e-20
* IEEE exp(+-10000) 100000 5.39e-20 2.34e-20
*
* In the tests over the interval exp(+-10000), the logarithms
* of the random arguments were uniformly distributed over
* [-10000, +10000].
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns -INFINITY
* log domain: x < 0; returns NAN
*/
#include "math.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double logl(long double x)
{
return log(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 2.32e-20
*/
static long double P[] = {
4.5270000862445199635215E-5L,
4.9854102823193375972212E-1L,
6.5787325942061044846969E0L,
2.9911919328553073277375E1L,
6.0949667980987787057556E1L,
5.7112963590585538103336E1L,
2.0039553499201281259648E1L,
};
static long double Q[] = {
/* 1.0000000000000000000000E0,*/
1.5062909083469192043167E1L,
8.3047565967967209469434E1L,
2.2176239823732856465394E2L,
3.0909872225312059774938E2L,
2.1642788614495947685003E2L,
6.0118660497603843919306E1L,
};
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
* Theoretical peak relative error = 6.16e-22
*/
static long double R[4] = {
1.9757429581415468984296E-3L,
-7.1990767473014147232598E-1L,
1.0777257190312272158094E1L,
-3.5717684488096787370998E1L,
};
static long double S[4] = {
/* 1.00000000000000000000E0L,*/
-2.6201045551331104417768E1L,
1.9361891836232102174846E2L,
-4.2861221385716144629696E2L,
};
static const long double C1 = 6.9314575195312500000000E-1L;
static const long double C2 = 1.4286068203094172321215E-6L;
#define SQRTH 0.70710678118654752440L
long double logl(long double x)
{
long double y, z;
int e;
if (isnan(x))
return x;
if (x == INFINITY)
return x;
if (x <= 0.0L) {
if (x == 0.0L)
return -INFINITY;
return NAN;
}
/* separate mantissa from exponent */
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
x = frexpl(x, &e);
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
z = x - 0.5L;
y = 0.5L * z + 0.5L;
} else { /* 2 (x-1)/(x+1) */
z = x - 0.5L;
z -= 0.5L;
y = 0.5L * x + 0.5L;
}
x = z / y;
z = x*x;
z = x * (z * __polevll(z, R, 3) / __p1evll(z, S, 3));
z = z + e * C2;
z = z + x;
z = z + e * C1;
return z;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH) {
e -= 1;
x = ldexpl(x, 1) - 1.0L; /* 2x - 1 */
} else {
x = x - 1.0L;
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6));
y = y + e * C2;
z = y - ldexpl(z, -1); /* y - 0.5 * z */
/* Note, the sum of above terms does not exceed x/4,
* so it contributes at most about 1/4 lsb to the error.
*/
z = z + x;
z = z + e * C1; /* This sum has an error of 1/2 lsb. */
return z;
}
#endif

1
Makefile
View file

@ -7,7 +7,6 @@ all: markdown README.html
$(MAKE) -C 04a
# don't compile all of 05 because it takes a while
$(MAKE) -C 05
$(MAKE) -C 05 install
$(MAKE) -C 05 tcc
clean:
$(MAKE) -C 00 clean