diff --git a/05/Makefile b/05/Makefile
index d5bd3bd..c41f79f 100644
--- a/05/Makefile
+++ b/05/Makefile
@@ -2,8 +2,8 @@ TCCDIR=tcc-0.9.27
 TCC=$(TCCDIR)/tcc
 TCC0=$(TCC)0
 TCCINST=tcc-bootstrap
-all: out04 a.out $(TCCDIR)/lib/libtcc1.a
-in04: *.b ../04a/out04
+all: out04 a.out tcc
+in04: *.b ../04a/out04 
 	../04a/out04 main.b in04
 out04: in04 ../04/out03
 	../04/out03 in04 out04
@@ -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
diff --git a/05/musl-0.6.0/Makefile b/05/musl-0.6.0/Makefile
index 460ff31..f18fa74 100644
--- a/05/musl-0.6.0/Makefile
+++ b/05/musl-0.6.0/Makefile
@@ -17,7 +17,7 @@ includedir = $(prefix)/include
 libdir = $(prefix)/lib
 
 SRCS = $(sort $(wildcard src/*/*.c))
-OBJS = $(SRCS:.c=.o) src/alloca86_64-bt.o src/alloca86_64.o src/libtcc1.o src/va_list.o src/syscall.o
+OBJS = $(SRCS:.c=.o) src/alloca86_64-bt.o src/alloca86_64.o src/libtcc1.o src/va_list.o src/syscall.o 
 LOBJS = $(OBJS:.o=.lo)
 GENH = include/bits/alltypes.h
 
@@ -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 $@ $<
 
diff --git a/05/musl-0.6.0/include/dlfcn.h b/05/musl-0.6.0/include/dlfcn.h
deleted file mode 100644
index 81b829c..0000000
--- a/05/musl-0.6.0/include/dlfcn.h
+++ /dev/null
@@ -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
diff --git a/05/musl-0.6.0/src/math/__log1p.h b/05/musl-0.6.0/src/math/__log1p.h
new file mode 100644
index 0000000..ec2c77b
--- /dev/null
+++ b/05/musl-0.6.0/src/math/__log1p.h
@@ -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);
+}
diff --git a/05/musl-0.6.0/src/math/__log1pf.h b/05/musl-0.6.0/src/math/__log1pf.h
new file mode 100644
index 0000000..99492c5
--- /dev/null
+++ b/05/musl-0.6.0/src/math/__log1pf.h
@@ -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);
+}
diff --git a/05/musl-0.6.0/src/math/log.c b/05/musl-0.6.0/src/math/log.c
new file mode 100644
index 0000000..71ad73d
--- /dev/null
+++ b/05/musl-0.6.0/src/math/log.c
@@ -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);
+	}
+}
diff --git a/05/musl-0.6.0/src/math/log10.c b/05/musl-0.6.0/src/math/log10.c
new file mode 100644
index 0000000..fa2b1e3
--- /dev/null
+++ b/05/musl-0.6.0/src/math/log10.c
@@ -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;
+}
diff --git a/05/musl-0.6.0/src/math/log10f.c b/05/musl-0.6.0/src/math/log10f.c
new file mode 100644
index 0000000..804359f
--- /dev/null
+++ b/05/musl-0.6.0/src/math/log10f.c
@@ -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;
+}
diff --git a/05/musl-0.6.0/src/math/log10l.c b/05/musl-0.6.0/src/math/log10l.c
new file mode 100644
index 0000000..a378fbb
--- /dev/null
+++ b/05/musl-0.6.0/src/math/log10l.c
@@ -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
diff --git a/05/musl-0.6.0/src/math/log1p.c b/05/musl-0.6.0/src/math/log1p.c
new file mode 100644
index 0000000..8f22e87
--- /dev/null
+++ b/05/musl-0.6.0/src/math/log1p.c
@@ -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);
+// }
diff --git a/05/musl-0.6.0/src/math/log1pf.c b/05/musl-0.6.0/src/math/log1pf.c
new file mode 100644
index 0000000..d746e67
--- /dev/null
+++ b/05/musl-0.6.0/src/math/log1pf.c
@@ -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);
+// }
diff --git a/05/musl-0.6.0/src/math/log1pl.c b/05/musl-0.6.0/src/math/log1pl.c
new file mode 100644
index 0000000..2eda3a6
--- /dev/null
+++ b/05/musl-0.6.0/src/math/log1pl.c
@@ -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
diff --git a/05/musl-0.6.0/src/math/log2.c b/05/musl-0.6.0/src/math/log2.c
new file mode 100644
index 0000000..21290d2
--- /dev/null
+++ b/05/musl-0.6.0/src/math/log2.c
@@ -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;
+}
diff --git a/05/musl-0.6.0/src/math/log2f.c b/05/musl-0.6.0/src/math/log2f.c
new file mode 100644
index 0000000..81c017f
--- /dev/null
+++ b/05/musl-0.6.0/src/math/log2f.c
@@ -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;
+}
diff --git a/05/musl-0.6.0/src/math/log2l.c b/05/musl-0.6.0/src/math/log2l.c
new file mode 100644
index 0000000..f490500
--- /dev/null
+++ b/05/musl-0.6.0/src/math/log2l.c
@@ -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
diff --git a/05/musl-0.6.0/src/math/logb.c b/05/musl-0.6.0/src/math/logb.c
new file mode 100644
index 0000000..2eb81f4
--- /dev/null
+++ b/05/musl-0.6.0/src/math/logb.c
@@ -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;
+}
diff --git a/05/musl-0.6.0/src/math/logbf.c b/05/musl-0.6.0/src/math/logbf.c
new file mode 100644
index 0000000..c04eb29
--- /dev/null
+++ b/05/musl-0.6.0/src/math/logbf.c
@@ -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;
+}
diff --git a/05/musl-0.6.0/src/math/logbl.c b/05/musl-0.6.0/src/math/logbl.c
new file mode 100644
index 0000000..704c853
--- /dev/null
+++ b/05/musl-0.6.0/src/math/logbl.c
@@ -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
diff --git a/05/musl-0.6.0/src/math/logf.c b/05/musl-0.6.0/src/math/logf.c
new file mode 100644
index 0000000..9e49b7a
--- /dev/null
+++ b/05/musl-0.6.0/src/math/logf.c
@@ -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);
+	}
+}
diff --git a/05/musl-0.6.0/src/math/logl.c b/05/musl-0.6.0/src/math/logl.c
new file mode 100644
index 0000000..15d7083
--- /dev/null
+++ b/05/musl-0.6.0/src/math/logl.c
@@ -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
diff --git a/Makefile b/Makefile
index 2df976a..e66cb69 100644
--- a/Makefile
+++ b/Makefile
@@ -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