ability to build musl from source

bundles musl 1.1.21

See #514

libc/musl/src/aio/aio.c

@@ -0,0 +1,400 @@
+#include <aio.h>
+#include <pthread.h>
+#include <semaphore.h>
+#include <limits.h>
+#include <errno.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <sys/auxv.h>
+#include "syscall.h"
+#include "atomic.h"
+#include "pthread_impl.h"
+
+/* The following is a threads-based implementation of AIO with minimal
+ * dependence on implementation details. Most synchronization is
+ * performed with pthread primitives, but atomics and futex operations
+ * are used for notification in a couple places where the pthread
+ * primitives would be inefficient or impractical.
+ *
+ * For each fd with outstanding aio operations, an aio_queue structure
+ * is maintained. These are reference-counted and destroyed by the last
+ * aio worker thread to exit. Accessing any member of the aio_queue
+ * structure requires a lock on the aio_queue. Adding and removing aio
+ * queues themselves requires a write lock on the global map object,
+ * a 4-level table mapping file descriptor numbers to aio queues. A
+ * read lock on the map is used to obtain locks on existing queues by
+ * excluding destruction of the queue by a different thread while it is
+ * being locked.
+ *
+ * Each aio queue has a list of active threads/operations. Presently there
+ * is a one to one relationship between threads and operations. The only
+ * members of the aio_thread structure which are accessed by other threads
+ * are the linked list pointers, op (which is immutable), running (which
+ * is updated atomically), and err (which is synchronized via running),
+ * so no locking is necessary. Most of the other other members are used
+ * for sharing data between the main flow of execution and cancellation
+ * cleanup handler.
+ *
+ * Taking any aio locks requires having all signals blocked. This is
+ * necessary because aio_cancel is needed by close, and close is required
+ * to be async-signal safe. All aio worker threads run with all signals
+ * blocked permanently.
+ */
+
+struct aio_thread {
+	pthread_t td;
+	struct aiocb *cb;
+	struct aio_thread *next, *prev;
+	struct aio_queue *q;
+	volatile int running;
+	int err, op;
+	ssize_t ret;
+};
+
+struct aio_queue {
+	int fd, seekable, append, ref, init;
+	pthread_mutex_t lock;
+	pthread_cond_t cond;
+	struct aio_thread *head;
+};
+
+struct aio_args {
+	struct aiocb *cb;
+	struct aio_queue *q;
+	int op;
+	sem_t sem;
+};
+
+static pthread_rwlock_t maplock = PTHREAD_RWLOCK_INITIALIZER;
+static struct aio_queue *****map;
+static volatile int aio_fd_cnt;
+volatile int __aio_fut;
+
+static struct aio_queue *__aio_get_queue(int fd, int need)
+{
+	if (fd < 0) {
+		errno = EBADF;
+		return 0;
+	}
+	int a=fd>>24;
+	unsigned char b=fd>>16, c=fd>>8, d=fd;
+	struct aio_queue *q = 0;
+	pthread_rwlock_rdlock(&maplock);
+	if ((!map || !map[a] || !map[a][b] || !map[a][b][c] || !(q=map[a][b][c][d])) && need) {
+		pthread_rwlock_unlock(&maplock);
+		if (fcntl(fd, F_GETFD) < 0) return 0;
+		pthread_rwlock_wrlock(&maplock);
+		if (!map) map = calloc(sizeof *map, (-1U/2+1)>>24);
+		if (!map) goto out;
+		if (!map[a]) map[a] = calloc(sizeof **map, 256);
+		if (!map[a]) goto out;
+		if (!map[a][b]) map[a][b] = calloc(sizeof ***map, 256);
+		if (!map[a][b]) goto out;
+		if (!map[a][b][c]) map[a][b][c] = calloc(sizeof ****map, 256);
+		if (!map[a][b][c]) goto out;
+		if (!(q = map[a][b][c][d])) {
+			map[a][b][c][d] = q = calloc(sizeof *****map, 1);
+			if (q) {
+				q->fd = fd;
+				pthread_mutex_init(&q->lock, 0);
+				pthread_cond_init(&q->cond, 0);
+				a_inc(&aio_fd_cnt);
+			}
+		}
+	}
+	if (q) pthread_mutex_lock(&q->lock);
+out:
+	pthread_rwlock_unlock(&maplock);
+	return q;
+}
+
+static void __aio_unref_queue(struct aio_queue *q)
+{
+	if (q->ref > 1) {
+		q->ref--;
+		pthread_mutex_unlock(&q->lock);
+		return;
+	}
+
+	/* This is potentially the last reference, but a new reference
+	 * may arrive since we cannot free the queue object without first
+	 * taking the maplock, which requires releasing the queue lock. */
+	pthread_mutex_unlock(&q->lock);
+	pthread_rwlock_wrlock(&maplock);
+	pthread_mutex_lock(&q->lock);
+	if (q->ref == 1) {
+		int fd=q->fd;
+		int a=fd>>24;
+		unsigned char b=fd>>16, c=fd>>8, d=fd;
+		map[a][b][c][d] = 0;
+		a_dec(&aio_fd_cnt);
+		pthread_rwlock_unlock(&maplock);
+		pthread_mutex_unlock(&q->lock);
+		free(q);
+	} else {
+		q->ref--;
+		pthread_rwlock_unlock(&maplock);
+		pthread_mutex_unlock(&q->lock);
+	}
+}
+
+static void cleanup(void *ctx)
+{
+	struct aio_thread *at = ctx;
+	struct aio_queue *q = at->q;
+	struct aiocb *cb = at->cb;
+	struct sigevent sev = cb->aio_sigevent;
+
+	/* There are four potential types of waiters we could need to wake:
+	 *   1. Callers of aio_cancel/close.
+	 *   2. Callers of aio_suspend with a single aiocb.
+	 *   3. Callers of aio_suspend with a list.
+	 *   4. AIO worker threads waiting for sequenced operations.
+	 * Types 1-3 are notified via atomics/futexes, mainly for AS-safety
+	 * considerations. Type 4 is notified later via a cond var. */
+
+	cb->__ret = at->ret;
+	if (a_swap(&at->running, 0) < 0)
+		__wake(&at->running, -1, 1);
+	if (a_swap(&cb->__err, at->err) != EINPROGRESS)
+		__wake(&cb->__err, -1, 1);
+	if (a_swap(&__aio_fut, 0))
+		__wake(&__aio_fut, -1, 1);
+
+	pthread_mutex_lock(&q->lock);
+
+	if (at->next) at->next->prev = at->prev;
+	if (at->prev) at->prev->next = at->next;
+	else q->head = at->next;
+
+	/* Signal aio worker threads waiting for sequenced operations. */
+	pthread_cond_broadcast(&q->cond);
+
+	__aio_unref_queue(q);
+
+	if (sev.sigev_notify == SIGEV_SIGNAL) {
+		siginfo_t si = {
+			.si_signo = sev.sigev_signo,
+			.si_value = sev.sigev_value,
+			.si_code = SI_ASYNCIO,
+			.si_pid = getpid(),
+			.si_uid = getuid()
+		};
+		__syscall(SYS_rt_sigqueueinfo, si.si_pid, si.si_signo, &si);
+	}
+	if (sev.sigev_notify == SIGEV_THREAD) {
+		a_store(&__pthread_self()->cancel, 0);
+		sev.sigev_notify_function(sev.sigev_value);
+	}
+}
+
+static void *io_thread_func(void *ctx)
+{
+	struct aio_thread at, *p;
+
+	struct aio_args *args = ctx;
+	struct aiocb *cb = args->cb;
+	int fd = cb->aio_fildes;
+	int op = args->op;
+	void *buf = (void *)cb->aio_buf;
+	size_t len = cb->aio_nbytes;
+	off_t off = cb->aio_offset;
+
+	struct aio_queue *q = args->q;
+	ssize_t ret;
+
+	pthread_mutex_lock(&q->lock);
+	sem_post(&args->sem);
+
+	at.op = op;
+	at.running = 1;
+	at.ret = -1;
+	at.err = ECANCELED;
+	at.q = q;
+	at.td = __pthread_self();
+	at.cb = cb;
+	at.prev = 0;
+	if ((at.next = q->head)) at.next->prev = &at;
+	q->head = &at;
+
+	if (!q->init) {
+		int seekable = lseek(fd, 0, SEEK_CUR) >= 0;
+		q->seekable = seekable;
+		q->append = !seekable || (fcntl(fd, F_GETFL) & O_APPEND);
+		q->init = 1;
+	}
+
+	pthread_cleanup_push(cleanup, &at);
+
+	/* Wait for sequenced operations. */
+	if (op!=LIO_READ && (op!=LIO_WRITE || q->append)) {
+		for (;;) {
+			for (p=at.next; p && p->op!=LIO_WRITE; p=p->next);
+			if (!p) break;
+			pthread_cond_wait(&q->cond, &q->lock);
+		}
+	}
+
+	pthread_mutex_unlock(&q->lock);
+
+	switch (op) {
+	case LIO_WRITE:
+		ret = q->append ? write(fd, buf, len) : pwrite(fd, buf, len, off);
+		break;
+	case LIO_READ:
+		ret = !q->seekable ? read(fd, buf, len) : pread(fd, buf, len, off);
+		break;
+	case O_SYNC:
+		ret = fsync(fd);
+		break;
+	case O_DSYNC:
+		ret = fdatasync(fd);
+		break;
+	}
+	at.ret = ret;
+	at.err = ret<0 ? errno : 0;
+	
+	pthread_cleanup_pop(1);
+
+	return 0;
+}
+
+static size_t io_thread_stack_size = MINSIGSTKSZ+2048;
+static pthread_once_t init_stack_size_once;
+
+static void init_stack_size()
+{
+	unsigned long val = __getauxval(AT_MINSIGSTKSZ);
+	if (val > MINSIGSTKSZ) io_thread_stack_size = val + 512;
+}
+
+static int submit(struct aiocb *cb, int op)
+{
+	int ret = 0;
+	pthread_attr_t a;
+	sigset_t allmask, origmask;
+	pthread_t td;
+	struct aio_queue *q = __aio_get_queue(cb->aio_fildes, 1);
+	struct aio_args args = { .cb = cb, .op = op, .q = q };
+	sem_init(&args.sem, 0, 0);
+
+	if (!q) {
+		if (errno != EBADF) errno = EAGAIN;
+		cb->__ret = -1;
+		cb->__err = errno;
+		return -1;
+	}
+	q->ref++;
+	pthread_mutex_unlock(&q->lock);
+
+	if (cb->aio_sigevent.sigev_notify == SIGEV_THREAD) {
+		if (cb->aio_sigevent.sigev_notify_attributes)
+			a = *cb->aio_sigevent.sigev_notify_attributes;
+		else
+			pthread_attr_init(&a);
+	} else {
+		pthread_once(&init_stack_size_once, init_stack_size);
+		pthread_attr_init(&a);
+		pthread_attr_setstacksize(&a, io_thread_stack_size);
+		pthread_attr_setguardsize(&a, 0);
+	}
+	pthread_attr_setdetachstate(&a, PTHREAD_CREATE_DETACHED);
+	sigfillset(&allmask);
+	pthread_sigmask(SIG_BLOCK, &allmask, &origmask);
+	cb->__err = EINPROGRESS;
+	if (pthread_create(&td, &a, io_thread_func, &args)) {
+		pthread_mutex_lock(&q->lock);
+		__aio_unref_queue(q);
+		cb->__err = errno = EAGAIN;
+		cb->__ret = ret = -1;
+	}
+	pthread_sigmask(SIG_SETMASK, &origmask, 0);
+
+	if (!ret) {
+		while (sem_wait(&args.sem));
+	}
+
+	return ret;
+}
+
+int aio_read(struct aiocb *cb)
+{
+	return submit(cb, LIO_READ);
+}
+
+int aio_write(struct aiocb *cb)
+{
+	return submit(cb, LIO_WRITE);
+}
+
+int aio_fsync(int op, struct aiocb *cb)
+{
+	if (op != O_SYNC && op != O_DSYNC) {
+		errno = EINVAL;
+		return -1;
+	}
+	return submit(cb, op);
+}
+
+ssize_t aio_return(struct aiocb *cb)
+{
+	return cb->__ret;
+}
+
+int aio_error(const struct aiocb *cb)
+{
+	a_barrier();
+	return cb->__err & 0x7fffffff;
+}
+
+int aio_cancel(int fd, struct aiocb *cb)
+{
+	sigset_t allmask, origmask;
+	int ret = AIO_ALLDONE;
+	struct aio_thread *p;
+	struct aio_queue *q;
+
+	/* Unspecified behavior case. Report an error. */
+	if (cb && fd != cb->aio_fildes) {
+		errno = EINVAL;
+		return -1;
+	}
+
+	sigfillset(&allmask);
+	pthread_sigmask(SIG_BLOCK, &allmask, &origmask);
+
+	errno = ENOENT;
+	if (!(q = __aio_get_queue(fd, 0))) {
+		if (errno == EBADF) ret = -1;
+		goto done;
+	}
+
+	for (p = q->head; p; p = p->next) {
+		if (cb && cb != p->cb) continue;
+		/* Transition target from running to running-with-waiters */
+		if (a_cas(&p->running, 1, -1)) {
+			pthread_cancel(p->td);
+			__wait(&p->running, 0, -1, 1);
+			if (p->err == ECANCELED) ret = AIO_CANCELED;
+		}
+	}
+
+	pthread_mutex_unlock(&q->lock);
+done:
+	pthread_sigmask(SIG_SETMASK, &origmask, 0);
+	return ret;
+}
+
+int __aio_close(int fd)
+{
+	a_barrier();
+	if (aio_fd_cnt) aio_cancel(fd, 0);
+	return fd;
+}
+
+weak_alias(aio_cancel, aio_cancel64);
+weak_alias(aio_error, aio_error64);
+weak_alias(aio_fsync, aio_fsync64);
+weak_alias(aio_read, aio_read64);
+weak_alias(aio_write, aio_write64);
+weak_alias(aio_return, aio_return64);

libc/musl/src/aio/aio_suspend.c

@@ -0,0 +1,76 @@
+#include <aio.h>
+#include <errno.h>
+#include <time.h>
+#include "atomic.h"
+#include "pthread_impl.h"
+
+int aio_suspend(const struct aiocb *const cbs[], int cnt, const struct timespec *ts)
+{
+	int i, tid = 0, ret, expect = 0;
+	struct timespec at;
+	volatile int dummy_fut, *pfut;
+	int nzcnt = 0;
+	const struct aiocb *cb = 0;
+
+	pthread_testcancel();
+
+	if (cnt<0) {
+		errno = EINVAL;
+		return -1;
+	}
+
+	for (i=0; i<cnt; i++) if (cbs[i]) {
+		if (aio_error(cbs[i]) != EINPROGRESS) return 0;
+		nzcnt++;
+		cb = cbs[i];
+	}
+
+	if (ts) {
+		clock_gettime(CLOCK_MONOTONIC, &at);
+		at.tv_sec += ts->tv_sec;
+		if ((at.tv_nsec += ts->tv_nsec) >= 1000000000) {
+			at.tv_nsec -= 1000000000;
+			at.tv_sec++;
+		}
+	}
+
+	for (;;) {
+		for (i=0; i<cnt; i++)
+			if (cbs[i] && aio_error(cbs[i]) != EINPROGRESS)
+				return 0;
+
+		switch (nzcnt) {
+		case 0:
+			pfut = &dummy_fut;
+			break;
+		case 1:
+			pfut = (void *)&cb->__err;
+			expect = EINPROGRESS | 0x80000000;
+			a_cas(pfut, EINPROGRESS, expect);
+			break;
+		default:
+			pfut = &__aio_fut;
+			if (!tid) tid = __pthread_self()->tid;
+			expect = a_cas(pfut, 0, tid);
+			if (!expect) expect = tid;
+			/* Need to recheck the predicate before waiting. */
+			for (i=0; i<cnt; i++)
+				if (cbs[i] && aio_error(cbs[i]) != EINPROGRESS)
+					return 0;
+			break;
+		}
+
+		ret = __timedwait_cp(pfut, expect, CLOCK_MONOTONIC, ts?&at:0, 1);
+
+		switch (ret) {
+		case ETIMEDOUT:
+			ret = EAGAIN;
+		case ECANCELED:
+		case EINTR:
+			errno = ret;
+			return -1;
+		}
+	}
+}
+
+weak_alias(aio_suspend, aio_suspend64);

libc/musl/src/aio/lio_listio.c

@@ -0,0 +1,143 @@
+#include <aio.h>
+#include <errno.h>
+#include <unistd.h>
+#include <string.h>
+#include "pthread_impl.h"
+
+struct lio_state {
+	struct sigevent *sev;
+	int cnt;
+	struct aiocb *cbs[];
+};
+
+static int lio_wait(struct lio_state *st)
+{
+	int i, err, got_err = 0;
+	int cnt = st->cnt;
+	struct aiocb **cbs = st->cbs;
+
+	for (;;) {
+		for (i=0; i<cnt; i++) {
+			if (!cbs[i]) continue;
+			err = aio_error(cbs[i]);
+			if (err==EINPROGRESS)
+				break;
+			if (err) got_err=1;
+			cbs[i] = 0;
+		}
+		if (i==cnt) {
+			if (got_err) {
+				errno = EIO;
+				return -1;
+			}
+			return 0;
+		}
+		if (aio_suspend((void *)cbs, cnt, 0))
+			return -1;
+	}
+}
+
+static void notify_signal(struct sigevent *sev)
+{
+	siginfo_t si = {
+		.si_signo = sev->sigev_signo,
+		.si_value = sev->sigev_value,
+		.si_code = SI_ASYNCIO,
+		.si_pid = getpid(),
+		.si_uid = getuid()
+	};
+	__syscall(SYS_rt_sigqueueinfo, si.si_pid, si.si_signo, &si);
+}
+
+static void *wait_thread(void *p)
+{
+	struct lio_state *st = p;
+	struct sigevent *sev = st->sev;
+	lio_wait(st);
+	free(st);
+	switch (sev->sigev_notify) {
+	case SIGEV_SIGNAL:
+		notify_signal(sev);
+		break;
+	case SIGEV_THREAD:
+		sev->sigev_notify_function(sev->sigev_value);
+		break;
+	}
+	return 0;
+}
+
+int lio_listio(int mode, struct aiocb *restrict const *restrict cbs, int cnt, struct sigevent *restrict sev)
+{
+	int i, ret;
+	struct lio_state *st=0;
+
+	if (cnt < 0) {
+		errno = EINVAL;
+		return -1;
+	}
+
+	if (mode == LIO_WAIT || (sev && sev->sigev_notify != SIGEV_NONE)) {
+		if (!(st = malloc(sizeof *st + cnt*sizeof *cbs))) {
+			errno = EAGAIN;
+			return -1;
+		}
+		st->cnt = cnt;
+		st->sev = sev;
+		memcpy(st->cbs, (void*) cbs, cnt*sizeof *cbs);
+	}
+
+	for (i=0; i<cnt; i++) {
+		if (!cbs[i]) continue;
+		switch (cbs[i]->aio_lio_opcode) {
+		case LIO_READ:
+			ret = aio_read(cbs[i]);
+			break;
+		case LIO_WRITE:
+			ret = aio_write(cbs[i]);
+			break;
+		default:
+			continue;
+		}
+		if (ret) {
+			free(st);
+			errno = EAGAIN;
+			return -1;
+		}
+	}
+
+	if (mode == LIO_WAIT) {
+		ret = lio_wait(st);
+		free(st);
+		return ret;
+	}
+
+	if (st) {
+		pthread_attr_t a;
+		sigset_t set;
+		pthread_t td;
+
+		if (sev->sigev_notify == SIGEV_THREAD) {
+			if (sev->sigev_notify_attributes)
+				a = *sev->sigev_notify_attributes;
+			else
+				pthread_attr_init(&a);
+		} else {
+			pthread_attr_init(&a);
+			pthread_attr_setstacksize(&a, PAGE_SIZE);
+			pthread_attr_setguardsize(&a, 0);
+		}
+		pthread_attr_setdetachstate(&a, PTHREAD_CREATE_DETACHED);
+		sigfillset(&set);
+		pthread_sigmask(SIG_BLOCK, &set, &set);
+		if (pthread_create(&td, &a, wait_thread, st)) {
+			free(st);
+			errno = EAGAIN;
+			return -1;
+		}
+		pthread_sigmask(SIG_SETMASK, &set, 0);
+	}
+
+	return 0;
+}
+
+weak_alias(lio_listio, lio_listio64);

libc/musl/src/complex/__cexp.c

@@ -0,0 +1,87 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/k_exp.c */
+/*-
+ * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include "libm.h"
+
+static const uint32_t k = 1799; /* constant for reduction */
+static const double kln2 = 1246.97177782734161156; /* k * ln2 */
+
+/*
+ * Compute exp(x), scaled to avoid spurious overflow.  An exponent is
+ * returned separately in 'expt'.
+ *
+ * Input:  ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91
+ * Output: 2**1023 <= y < 2**1024
+ */
+static double __frexp_exp(double x, int *expt)
+{
+	double exp_x;
+	uint32_t hx;
+
+	/*
+	 * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to
+	 * minimize |exp(kln2) - 2**k|.  We also scale the exponent of
+	 * exp_x to MAX_EXP so that the result can be multiplied by
+	 * a tiny number without losing accuracy due to denormalization.
+	 */
+	exp_x = exp(x - kln2);
+	GET_HIGH_WORD(hx, exp_x);
+	*expt = (hx >> 20) - (0x3ff + 1023) + k;
+	SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20));
+	return exp_x;
+}
+
+/*
+ * __ldexp_cexp(x, expt) compute exp(x) * 2**expt.
+ * It is intended for large arguments (real part >= ln(DBL_MAX))
+ * where care is needed to avoid overflow.
+ *
+ * The present implementation is narrowly tailored for our hyperbolic and
+ * exponential functions.  We assume expt is small (0 or -1), and the caller
+ * has filtered out very large x, for which overflow would be inevitable.
+ */
+double complex __ldexp_cexp(double complex z, int expt)
+{
+	double x, y, exp_x, scale1, scale2;
+	int ex_expt, half_expt;
+
+	x = creal(z);
+	y = cimag(z);
+	exp_x = __frexp_exp(x, &ex_expt);
+	expt += ex_expt;
+
+	/*
+	 * Arrange so that scale1 * scale2 == 2**expt.  We use this to
+	 * compensate for scalbn being horrendously slow.
+	 */
+	half_expt = expt / 2;
+	INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0);
+	half_expt = expt - half_expt;
+	INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
+
+	return CMPLX(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2);
+}

libc/musl/src/complex/__cexpf.c

@@ -0,0 +1,68 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/k_expf.c */
+/*-
+ * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include "libm.h"
+
+static const uint32_t k = 235; /* constant for reduction */
+static const float kln2 = 162.88958740F; /* k * ln2 */
+
+/*
+ * See __cexp.c for details.
+ *
+ * Input:  ln(FLT_MAX) <= x < ln(2 * FLT_MAX / FLT_MIN_DENORM) ~= 192.7
+ * Output: 2**127 <= y < 2**128
+ */
+static float __frexp_expf(float x, int *expt)
+{
+	float exp_x;
+	uint32_t hx;
+
+	exp_x = expf(x - kln2);
+	GET_FLOAT_WORD(hx, exp_x);
+	*expt = (hx >> 23) - (0x7f + 127) + k;
+	SET_FLOAT_WORD(exp_x, (hx & 0x7fffff) | ((0x7f + 127) << 23));
+	return exp_x;
+}
+
+float complex __ldexp_cexpf(float complex z, int expt)
+{
+	float x, y, exp_x, scale1, scale2;
+	int ex_expt, half_expt;
+
+	x = crealf(z);
+	y = cimagf(z);
+	exp_x = __frexp_expf(x, &ex_expt);
+	expt += ex_expt;
+
+	half_expt = expt / 2;
+	SET_FLOAT_WORD(scale1, (0x7f + half_expt) << 23);
+	half_expt = expt - half_expt;
+	SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23);
+
+	return CMPLXF(cosf(y) * exp_x * scale1 * scale2,
+	  sinf(y) * exp_x * scale1 * scale2);
+}

libc/musl/src/complex/cabs.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+double cabs(double complex z)
+{
+	return hypot(creal(z), cimag(z));
+}

libc/musl/src/complex/cabsf.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+float cabsf(float complex z)
+{
+	return hypotf(crealf(z), cimagf(z));
+}

libc/musl/src/complex/cabsl.c

@@ -0,0 +1,13 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double cabsl(long double complex z)
+{
+	return cabs(z);
+}
+#else
+long double cabsl(long double complex z)
+{
+	return hypotl(creall(z), cimagl(z));
+}
+#endif

libc/musl/src/complex/cacos.c

@@ -0,0 +1,11 @@
+#include "libm.h"
+
+// FIXME: Hull et al. "Implementing the complex arcsine and arccosine functions using exception handling" 1997
+
+/* acos(z) = pi/2 - asin(z) */
+
+double complex cacos(double complex z)
+{
+	z = casin(z);
+	return CMPLX(M_PI_2 - creal(z), -cimag(z));
+}

libc/musl/src/complex/cacosf.c

@@ -0,0 +1,9 @@
+#include "libm.h"
+
+// FIXME
+
+float complex cacosf(float complex z)
+{
+	z = casinf(z);
+	return CMPLXF((float)M_PI_2 - crealf(z), -cimagf(z));
+}

libc/musl/src/complex/cacosh.c

@@ -0,0 +1,9 @@
+#include "libm.h"
+
+/* acosh(z) = i acos(z) */
+
+double complex cacosh(double complex z)
+{
+	z = cacos(z);
+	return CMPLX(-cimag(z), creal(z));
+}

libc/musl/src/complex/cacoshf.c

@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float complex cacoshf(float complex z)
+{
+	z = cacosf(z);
+	return CMPLXF(-cimagf(z), crealf(z));
+}

libc/musl/src/complex/cacoshl.c

@@ -0,0 +1,14 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex cacoshl(long double complex z)
+{
+	return cacosh(z);
+}
+#else
+long double complex cacoshl(long double complex z)
+{
+	z = cacosl(z);
+	return CMPLXL(-cimagl(z), creall(z));
+}
+#endif

libc/musl/src/complex/cacosl.c

@@ -0,0 +1,16 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex cacosl(long double complex z)
+{
+	return cacos(z);
+}
+#else
+// FIXME
+#define PI_2 1.57079632679489661923132169163975144L
+long double complex cacosl(long double complex z)
+{
+	z = casinl(z);
+	return CMPLXL(PI_2 - creall(z), -cimagl(z));
+}
+#endif

libc/musl/src/complex/carg.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+double carg(double complex z)
+{
+	return atan2(cimag(z), creal(z));
+}

libc/musl/src/complex/cargf.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+float cargf(float complex z)
+{
+	return atan2f(cimagf(z), crealf(z));
+}

libc/musl/src/complex/cargl.c

@@ -0,0 +1,13 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double cargl(long double complex z)
+{
+	return carg(z);
+}
+#else
+long double cargl(long double complex z)
+{
+	return atan2l(cimagl(z), creall(z));
+}
+#endif

libc/musl/src/complex/casin.c

@@ -0,0 +1,17 @@
+#include "libm.h"
+
+// FIXME
+
+/* asin(z) = -i log(i z + sqrt(1 - z*z)) */
+
+double complex casin(double complex z)
+{
+	double complex w;
+	double x, y;
+
+	x = creal(z);
+	y = cimag(z);
+	w = CMPLX(1.0 - (x - y)*(x + y), -2.0*x*y);
+	double complex r = clog(CMPLX(-y, x) + csqrt(w));
+	return CMPLX(cimag(r), -creal(r));
+}

libc/musl/src/complex/casinf.c

@@ -0,0 +1,15 @@
+#include "libm.h"
+
+// FIXME
+
+float complex casinf(float complex z)
+{
+	float complex w;
+	float x, y;
+
+	x = crealf(z);
+	y = cimagf(z);
+	w = CMPLXF(1.0 - (x - y)*(x + y), -2.0*x*y);
+	float complex r = clogf(CMPLXF(-y, x) + csqrtf(w));
+	return CMPLXF(cimagf(r), -crealf(r));
+}

libc/musl/src/complex/casinh.c

@@ -0,0 +1,9 @@
+#include "libm.h"
+
+/* asinh(z) = -i asin(i z) */
+
+double complex casinh(double complex z)
+{
+	z = casin(CMPLX(-cimag(z), creal(z)));
+	return CMPLX(cimag(z), -creal(z));
+}

libc/musl/src/complex/casinhf.c

@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float complex casinhf(float complex z)
+{
+	z = casinf(CMPLXF(-cimagf(z), crealf(z)));
+	return CMPLXF(cimagf(z), -crealf(z));
+}

libc/musl/src/complex/casinhl.c

@@ -0,0 +1,14 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex casinhl(long double complex z)
+{
+	return casinh(z);
+}
+#else
+long double complex casinhl(long double complex z)
+{
+	z = casinl(CMPLXL(-cimagl(z), creall(z)));
+	return CMPLXL(cimagl(z), -creall(z));
+}
+#endif

libc/musl/src/complex/casinl.c

@@ -0,0 +1,21 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex casinl(long double complex z)
+{
+	return casin(z);
+}
+#else
+// FIXME
+long double complex casinl(long double complex z)
+{
+	long double complex w;
+	long double x, y;
+
+	x = creall(z);
+	y = cimagl(z);
+	w = CMPLXL(1.0 - (x - y)*(x + y), -2.0*x*y);
+	long double complex r = clogl(CMPLXL(-y, x) + csqrtl(w));
+	return CMPLXL(cimagl(r), -creall(r));
+}
+#endif

libc/musl/src/complex/catan.c

@@ -0,0 +1,107 @@
+/* origin: OpenBSD /usr/src/lib/libm/src/s_catan.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.
+ */
+/*
+ *      Complex circular arc tangent
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double complex catan();
+ * double complex z, w;
+ *
+ * w = catan (z);
+ *
+ *
+ * DESCRIPTION:
+ *
+ * If
+ *     z = x + iy,
+ *
+ * then
+ *          1       (    2x     )
+ * Re w  =  - arctan(-----------)  +  k PI
+ *          2       (     2    2)
+ *                  (1 - x  - y )
+ *
+ *               ( 2         2)
+ *          1    (x  +  (y+1) )
+ * Im w  =  - log(------------)
+ *          4    ( 2         2)
+ *               (x  +  (y-1) )
+ *
+ * Where k is an arbitrary integer.
+ *
+ * catan(z) = -i catanh(iz).
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10      5900       1.3e-16     7.8e-18
+ *    IEEE      -10,+10     30000       2.3e-15     8.5e-17
+ * The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,
+ * had peak relative error 1.5e-16, rms relative error
+ * 2.9e-17.  See also clog().
+ */
+
+#include "libm.h"
+
+#define MAXNUM 1.0e308
+
+static const double DP1 = 3.14159265160560607910E0;
+static const double DP2 = 1.98418714791870343106E-9;
+static const double DP3 = 1.14423774522196636802E-17;
+
+static double _redupi(double x)
+{
+	double t;
+	long i;
+
+	t = x/M_PI;
+	if (t >= 0.0)
+		t += 0.5;
+	else
+		t -= 0.5;
+
+	i = t;  /* the multiple */
+	t = i;
+	t = ((x - t * DP1) - t * DP2) - t * DP3;
+	return t;
+}
+
+double complex catan(double complex z)
+{
+	double complex w;
+	double a, t, x, x2, y;
+
+	x = creal(z);
+	y = cimag(z);
+
+	x2 = x * x;
+	a = 1.0 - x2 - (y * y);
+
+	t = 0.5 * atan2(2.0 * x, a);
+	w = _redupi(t);
+
+	t = y - 1.0;
+	a = x2 + (t * t);
+
+	t = y + 1.0;
+	a = (x2 + t * t)/a;
+	w = CMPLX(w, 0.25 * log(a));
+	return w;
+}

libc/musl/src/complex/catanf.c

@@ -0,0 +1,115 @@
+/* origin: OpenBSD /usr/src/lib/libm/src/s_catanf.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.
+ */
+/*
+ *      Complex circular arc tangent
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float complex catanf();
+ * float complex z, w;
+ *
+ * w = catanf( z );
+ *
+ *
+ * DESCRIPTION:
+ *
+ * If
+ *     z = x + iy,
+ *
+ * then
+ *          1       (    2x     )
+ * Re w  =  - arctan(-----------)  +  k PI
+ *          2       (     2    2)
+ *                  (1 - x  - y )
+ *
+ *               ( 2         2)
+ *          1    (x  +  (y+1) )
+ * Im w  =  - log(------------)
+ *          4    ( 2         2)
+ *               (x  +  (y-1) )
+ *
+ * Where k is an arbitrary integer.
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10,+10     30000        2.3e-6      5.2e-8
+ */
+
+#include "libm.h"
+
+#define MAXNUMF 1.0e38F
+
+static const double DP1 = 3.140625;
+static const double DP2 = 9.67502593994140625E-4;
+static const double DP3 = 1.509957990978376432E-7;
+
+static float _redupif(float xx)
+{
+	float x, t;
+	long i;
+
+	x = xx;
+	t = x/(float)M_PI;
+	if (t >= 0.0f)
+		t += 0.5f;
+	else
+		t -= 0.5f;
+
+	i = t;  /* the multiple */
+	t = i;
+	t = ((x - t * DP1) - t * DP2) - t * DP3;
+	return t;
+}
+
+float complex catanf(float complex z)
+{
+	float complex w;
+	float a, t, x, x2, y;
+
+	x = crealf(z);
+	y = cimagf(z);
+
+	if ((x == 0.0f) && (y > 1.0f))
+		goto ovrf;
+
+	x2 = x * x;
+	a = 1.0f - x2 - (y * y);
+	if (a == 0.0f)
+		goto ovrf;
+
+	t = 0.5f * atan2f(2.0f * x, a);
+	w = _redupif(t);
+
+	t = y - 1.0f;
+	a = x2 + (t * t);
+	if (a == 0.0f)
+		goto ovrf;
+
+	t = y + 1.0f;
+	a = (x2 + (t * t))/a;
+	w = w + (0.25f * logf (a)) * I;
+	return w;
+
+ovrf:
+	// FIXME
+	w = MAXNUMF + MAXNUMF * I;
+	return w;
+}

libc/musl/src/complex/catanh.c

@@ -0,0 +1,9 @@
+#include "libm.h"
+
+/* atanh = -i atan(i z) */
+
+double complex catanh(double complex z)
+{
+	z = catan(CMPLX(-cimag(z), creal(z)));
+	return CMPLX(cimag(z), -creal(z));
+}

libc/musl/src/complex/catanhf.c

@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float complex catanhf(float complex z)
+{
+	z = catanf(CMPLXF(-cimagf(z), crealf(z)));
+	return CMPLXF(cimagf(z), -crealf(z));
+}

libc/musl/src/complex/catanhl.c

@@ -0,0 +1,14 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex catanhl(long double complex z)
+{
+	return catanh(z);
+}
+#else
+long double complex catanhl(long double complex z)
+{
+	z = catanl(CMPLXL(-cimagl(z), creall(z)));
+	return CMPLXL(cimagl(z), -creall(z));
+}
+#endif

libc/musl/src/complex/catanl.c

@@ -0,0 +1,126 @@
+/* origin: OpenBSD /usr/src/lib/libm/src/s_catanl.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.
+ */
+/*
+ *      Complex circular arc tangent
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double complex catanl();
+ * long double complex z, w;
+ *
+ * w = catanl( z );
+ *
+ *
+ * DESCRIPTION:
+ *
+ * If
+ *     z = x + iy,
+ *
+ * then
+ *          1       (    2x     )
+ * Re w  =  - arctan(-----------)  +  k PI
+ *          2       (     2    2)
+ *                  (1 - x  - y )
+ *
+ *               ( 2         2)
+ *          1    (x  +  (y+1) )
+ * Im w  =  - log(------------)
+ *          4    ( 2         2)
+ *               (x  +  (y-1) )
+ *
+ * Where k is an arbitrary integer.
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       -10,+10      5900       1.3e-16     7.8e-18
+ *    IEEE      -10,+10     30000       2.3e-15     8.5e-17
+ * The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,
+ * had peak relative error 1.5e-16, rms relative error
+ * 2.9e-17.  See also clog().
+ */
+
+#include <complex.h>
+#include <float.h>
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex catanl(long double complex z)
+{
+	return catan(z);
+}
+#else
+static const long double PIL = 3.141592653589793238462643383279502884197169L;
+static const long double DP1 = 3.14159265358979323829596852490908531763125L;
+static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L;
+static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L;
+
+static long double redupil(long double x)
+{
+	long double t;
+	long i;
+
+	t = x / PIL;
+	if (t >= 0.0L)
+		t += 0.5L;
+	else
+		t -= 0.5L;
+
+	i = t;  /* the multiple */
+	t = i;
+	t = ((x - t * DP1) - t * DP2) - t * DP3;
+	return t;
+}
+
+long double complex catanl(long double complex z)
+{
+	long double complex w;
+	long double a, t, x, x2, y;
+
+	x = creall(z);
+	y = cimagl(z);
+
+	if ((x == 0.0L) && (y > 1.0L))
+		goto ovrf;
+
+	x2 = x * x;
+	a = 1.0L - x2 - (y * y);
+	if (a == 0.0L)
+		goto ovrf;
+
+	t = atan2l(2.0L * x, a) * 0.5L;
+	w = redupil(t);
+
+	t = y - 1.0L;
+	a = x2 + (t * t);
+	if (a == 0.0L)
+		goto ovrf;
+
+	t = y + 1.0L;
+	a = (x2 + (t * t)) / a;
+	w = w + (0.25L * logl(a)) * I;
+	return w;
+
+ovrf:
+	// FIXME
+	w = LDBL_MAX + LDBL_MAX * I;
+	return w;
+}
+#endif

libc/musl/src/complex/ccos.c

@@ -0,0 +1,8 @@
+#include "libm.h"
+
+/* cos(z) = cosh(i z) */
+
+double complex ccos(double complex z)
+{
+	return ccosh(CMPLX(-cimag(z), creal(z)));
+}

libc/musl/src/complex/ccosf.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+float complex ccosf(float complex z)
+{
+	return ccoshf(CMPLXF(-cimagf(z), crealf(z)));
+}

libc/musl/src/complex/ccosh.c

@@ -0,0 +1,140 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_ccosh.c */
+/*-
+ * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+/*
+ * Hyperbolic cosine of a complex argument z = x + i y.
+ *
+ * cosh(z) = cosh(x+iy)
+ *         = cosh(x) cos(y) + i sinh(x) sin(y).
+ *
+ * Exceptional values are noted in the comments within the source code.
+ * These values and the return value were taken from n1124.pdf.
+ */
+
+#include "libm.h"
+
+static const double huge = 0x1p1023;
+
+double complex ccosh(double complex z)
+{
+	double x, y, h;
+	int32_t hx, hy, ix, iy, lx, ly;
+
+	x = creal(z);
+	y = cimag(z);
+
+	EXTRACT_WORDS(hx, lx, x);
+	EXTRACT_WORDS(hy, ly, y);
+
+	ix = 0x7fffffff & hx;
+	iy = 0x7fffffff & hy;
+
+	/* Handle the nearly-non-exceptional cases where x and y are finite. */
+	if (ix < 0x7ff00000 && iy < 0x7ff00000) {
+		if ((iy | ly) == 0)
+			return CMPLX(cosh(x), x * y);
+		if (ix < 0x40360000)    /* small x: normal case */
+			return CMPLX(cosh(x) * cos(y), sinh(x) * sin(y));
+
+		/* |x| >= 22, so cosh(x) ~= exp(|x|) */
+		if (ix < 0x40862e42) {
+			/* x < 710: exp(|x|) won't overflow */
+			h = exp(fabs(x)) * 0.5;
+			return CMPLX(h * cos(y), copysign(h, x) * sin(y));
+		} else if (ix < 0x4096bbaa) {
+			/* x < 1455: scale to avoid overflow */
+			z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
+			return CMPLX(creal(z), cimag(z) * copysign(1, x));
+		} else {
+			/* x >= 1455: the result always overflows */
+			h = huge * x;
+			return CMPLX(h * h * cos(y), h * sin(y));
+		}
+	}
+
+	/*
+	 * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0.
+	 * The sign of 0 in the result is unspecified.  Choice = normally
+	 * the same as dNaN.  Raise the invalid floating-point exception.
+	 *
+	 * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0.
+	 * The sign of 0 in the result is unspecified.  Choice = normally
+	 * the same as d(NaN).
+	 */
+	if ((ix | lx) == 0 && iy >= 0x7ff00000)
+		return CMPLX(y - y, copysign(0, x * (y - y)));
+
+	/*
+	 * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
+	 *
+	 * cosh(NaN +- I 0)   = d(NaN) + I sign(d(NaN, +-0))0.
+	 * The sign of 0 in the result is unspecified.
+	 */
+	if ((iy | ly) == 0 && ix >= 0x7ff00000) {
+		if (((hx & 0xfffff) | lx) == 0)
+			return CMPLX(x * x, copysign(0, x) * y);
+		return CMPLX(x * x, copysign(0, (x + x) * y));
+	}
+
+	/*
+	 * cosh(x +- I Inf) = dNaN + I dNaN.
+	 * Raise the invalid floating-point exception for finite nonzero x.
+	 *
+	 * cosh(x + I NaN) = d(NaN) + I d(NaN).
+	 * Optionally raises the invalid floating-point exception for finite
+	 * nonzero x.  Choice = don't raise (except for signaling NaNs).
+	 */
+	if (ix < 0x7ff00000 && iy >= 0x7ff00000)
+		return CMPLX(y - y, x * (y - y));
+
+	/*
+	 * cosh(+-Inf + I NaN)  = +Inf + I d(NaN).
+	 *
+	 * cosh(+-Inf +- I Inf) = +Inf + I dNaN.
+	 * The sign of Inf in the result is unspecified.  Choice = always +.
+	 * Raise the invalid floating-point exception.
+	 *
+	 * cosh(+-Inf + I y)   = +Inf cos(y) +- I Inf sin(y)
+	 */
+	if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
+		if (iy >= 0x7ff00000)
+			return CMPLX(x * x, x * (y - y));
+		return CMPLX((x * x) * cos(y), x * sin(y));
+	}
+
+	/*
+	 * cosh(NaN + I NaN)  = d(NaN) + I d(NaN).
+	 *
+	 * cosh(NaN +- I Inf) = d(NaN) + I d(NaN).
+	 * Optionally raises the invalid floating-point exception.
+	 * Choice = raise.
+	 *
+	 * cosh(NaN + I y)    = d(NaN) + I d(NaN).
+	 * Optionally raises the invalid floating-point exception for finite
+	 * nonzero y.  Choice = don't raise (except for signaling NaNs).
+	 */
+	return CMPLX((x * x) * (y - y), (x + x) * (y - y));
+}

libc/musl/src/complex/ccoshf.c

@@ -0,0 +1,90 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_ccoshf.c */
+/*-
+ * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+/*
+ * Hyperbolic cosine of a complex argument.  See s_ccosh.c for details.
+ */
+
+#include "libm.h"
+
+static const float huge = 0x1p127;
+
+float complex ccoshf(float complex z)
+{
+	float x, y, h;
+	int32_t hx, hy, ix, iy;
+
+	x = crealf(z);
+	y = cimagf(z);
+
+	GET_FLOAT_WORD(hx, x);
+	GET_FLOAT_WORD(hy, y);
+
+	ix = 0x7fffffff & hx;
+	iy = 0x7fffffff & hy;
+
+	if (ix < 0x7f800000 && iy < 0x7f800000) {
+		if (iy == 0)
+			return CMPLXF(coshf(x), x * y);
+		if (ix < 0x41100000)    /* small x: normal case */
+			return CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y));
+
+		/* |x| >= 9, so cosh(x) ~= exp(|x|) */
+		if (ix < 0x42b17218) {
+			/* x < 88.7: expf(|x|) won't overflow */
+			h = expf(fabsf(x)) * 0.5f;
+			return CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y));
+		} else if (ix < 0x4340b1e7) {
+			/* x < 192.7: scale to avoid overflow */
+			z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
+			return CMPLXF(crealf(z), cimagf(z) * copysignf(1, x));
+		} else {
+			/* x >= 192.7: the result always overflows */
+			h = huge * x;
+			return CMPLXF(h * h * cosf(y), h * sinf(y));
+		}
+	}
+
+	if (ix == 0 && iy >= 0x7f800000)
+		return CMPLXF(y - y, copysignf(0, x * (y - y)));
+
+	if (iy == 0 && ix >= 0x7f800000) {
+		if ((hx & 0x7fffff) == 0)
+			return CMPLXF(x * x, copysignf(0, x) * y);
+		return CMPLXF(x * x, copysignf(0, (x + x) * y));
+	}
+
+	if (ix < 0x7f800000 && iy >= 0x7f800000)
+		return CMPLXF(y - y, x * (y - y));
+
+	if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
+		if (iy >= 0x7f800000)
+			return CMPLXF(x * x, x * (y - y));
+		return CMPLXF((x * x) * cosf(y), x * sinf(y));
+	}
+
+	return CMPLXF((x * x) * (y - y), (x + x) * (y - y));
+}

libc/musl/src/complex/ccoshl.c

@@ -0,0 +1,7 @@
+#include "libm.h"
+
+//FIXME
+long double complex ccoshl(long double complex z)
+{
+	return ccosh(z);
+}

libc/musl/src/complex/ccosl.c

@@ -0,0 +1,13 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex ccosl(long double complex z)
+{
+	return ccos(z);
+}
+#else
+long double complex ccosl(long double complex z)
+{
+	return ccoshl(CMPLXL(-cimagl(z), creall(z)));
+}
+#endif

libc/musl/src/complex/cexp.c

@@ -0,0 +1,83 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_cexp.c */
+/*-
+ * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include "libm.h"
+
+static const uint32_t
+exp_ovfl  = 0x40862e42,  /* high bits of MAX_EXP * ln2 ~= 710 */
+cexp_ovfl = 0x4096b8e4;  /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
+
+double complex cexp(double complex z)
+{
+	double x, y, exp_x;
+	uint32_t hx, hy, lx, ly;
+
+	x = creal(z);
+	y = cimag(z);
+
+	EXTRACT_WORDS(hy, ly, y);
+	hy &= 0x7fffffff;
+
+	/* cexp(x + I 0) = exp(x) + I 0 */
+	if ((hy | ly) == 0)
+		return CMPLX(exp(x), y);
+	EXTRACT_WORDS(hx, lx, x);
+	/* cexp(0 + I y) = cos(y) + I sin(y) */
+	if (((hx & 0x7fffffff) | lx) == 0)
+		return CMPLX(cos(y), sin(y));
+
+	if (hy >= 0x7ff00000) {
+		if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
+			/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
+			return CMPLX(y - y, y - y);
+		} else if (hx & 0x80000000) {
+			/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
+			return CMPLX(0.0, 0.0);
+		} else {
+			/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
+			return CMPLX(x, y - y);
+		}
+	}
+
+	if (hx >= exp_ovfl && hx <= cexp_ovfl) {
+		/*
+		 * x is between 709.7 and 1454.3, so we must scale to avoid
+		 * overflow in exp(x).
+		 */
+		return __ldexp_cexp(z, 0);
+	} else {
+		/*
+		 * Cases covered here:
+		 *  -  x < exp_ovfl and exp(x) won't overflow (common case)
+		 *  -  x > cexp_ovfl, so exp(x) * s overflows for all s > 0
+		 *  -  x = +-Inf (generated by exp())
+		 *  -  x = NaN (spurious inexact exception from y)
+		 */
+		exp_x = exp(x);
+		return CMPLX(exp_x * cos(y), exp_x * sin(y));
+	}
+}

libc/musl/src/complex/cexpf.c

@@ -0,0 +1,83 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_cexpf.c */
+/*-
+ * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include "libm.h"
+
+static const uint32_t
+exp_ovfl  = 0x42b17218,  /* MAX_EXP * ln2 ~= 88.722839355 */
+cexp_ovfl = 0x43400074;  /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
+
+float complex cexpf(float complex z)
+{
+	float x, y, exp_x;
+	uint32_t hx, hy;
+
+	x = crealf(z);
+	y = cimagf(z);
+
+	GET_FLOAT_WORD(hy, y);
+	hy &= 0x7fffffff;
+
+	/* cexp(x + I 0) = exp(x) + I 0 */
+	if (hy == 0)
+		return CMPLXF(expf(x), y);
+	GET_FLOAT_WORD(hx, x);
+	/* cexp(0 + I y) = cos(y) + I sin(y) */
+	if ((hx & 0x7fffffff) == 0)
+		return CMPLXF(cosf(y), sinf(y));
+
+	if (hy >= 0x7f800000) {
+		if ((hx & 0x7fffffff) != 0x7f800000) {
+			/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
+			return CMPLXF(y - y, y - y);
+		} else if (hx & 0x80000000) {
+			/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
+			return CMPLXF(0.0, 0.0);
+		} else {
+			/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
+			return CMPLXF(x, y - y);
+		}
+	}
+
+	if (hx >= exp_ovfl && hx <= cexp_ovfl) {
+		/*
+		 * x is between 88.7 and 192, so we must scale to avoid
+		 * overflow in expf(x).
+		 */
+		return __ldexp_cexpf(z, 0);
+	} else {
+		/*
+		 * Cases covered here:
+		 *  -  x < exp_ovfl and exp(x) won't overflow (common case)
+		 *  -  x > cexp_ovfl, so exp(x) * s overflows for all s > 0
+		 *  -  x = +-Inf (generated by exp())
+		 *  -  x = NaN (spurious inexact exception from y)
+		 */
+		exp_x = expf(x);
+		return CMPLXF(exp_x * cosf(y), exp_x * sinf(y));
+	}
+}

libc/musl/src/complex/cexpl.c

@@ -0,0 +1,7 @@
+#include "libm.h"
+
+//FIXME
+long double complex cexpl(long double complex z)
+{
+	return cexp(z);
+}

libc/musl/src/complex/cimag.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+double (cimag)(double complex z)
+{
+	return cimag(z);
+}

libc/musl/src/complex/cimagf.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+float (cimagf)(float complex z)
+{
+	return cimagf(z);
+}

libc/musl/src/complex/cimagl.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+long double (cimagl)(long double complex z)
+{
+	return cimagl(z);
+}

libc/musl/src/complex/clog.c

@@ -0,0 +1,14 @@
+#include "libm.h"
+
+// FIXME
+
+/* log(z) = log(|z|) + i arg(z) */
+
+double complex clog(double complex z)
+{
+	double r, phi;
+
+	r = cabs(z);
+	phi = carg(z);
+	return CMPLX(log(r), phi);
+}

libc/musl/src/complex/clogf.c

@@ -0,0 +1,12 @@
+#include "libm.h"
+
+// FIXME
+
+float complex clogf(float complex z)
+{
+	float r, phi;
+
+	r = cabsf(z);
+	phi = cargf(z);
+	return CMPLXF(logf(r), phi);
+}

libc/musl/src/complex/clogl.c

@@ -0,0 +1,18 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex clogl(long double complex z)
+{
+	return clog(z);
+}
+#else
+// FIXME
+long double complex clogl(long double complex z)
+{
+	long double r, phi;
+
+	r = cabsl(z);
+	phi = cargl(z);
+	return CMPLXL(logl(r), phi);
+}
+#endif

libc/musl/src/complex/conj.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+double complex conj(double complex z)
+{
+	return CMPLX(creal(z), -cimag(z));
+}

libc/musl/src/complex/conjf.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+float complex conjf(float complex z)
+{
+	return CMPLXF(crealf(z), -cimagf(z));
+}

libc/musl/src/complex/conjl.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+long double complex conjl(long double complex z)
+{
+	return CMPLXL(creall(z), -cimagl(z));
+}

libc/musl/src/complex/cpow.c

@@ -0,0 +1,8 @@
+#include "libm.h"
+
+/* pow(z, c) = exp(c log(z)), See C99 G.6.4.1 */
+
+double complex cpow(double complex z, double complex c)
+{
+	return cexp(c * clog(z));
+}

libc/musl/src/complex/cpowf.c

@@ -0,0 +1,6 @@
+#include "libm.h"
+
+float complex cpowf(float complex z, float complex c)
+{
+	return cexpf(c * clogf(z));
+}

libc/musl/src/complex/cpowl.c

@@ -0,0 +1,13 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex cpowl(long double complex z, long double complex c)
+{
+	return cpow(z, c);
+}
+#else
+long double complex cpowl(long double complex z, long double complex c)
+{
+	return cexpl(c * clogl(z));
+}
+#endif

libc/musl/src/complex/cproj.c

@@ -0,0 +1,8 @@
+#include "libm.h"
+
+double complex cproj(double complex z)
+{
+	if (isinf(creal(z)) || isinf(cimag(z)))
+		return CMPLX(INFINITY, copysign(0.0, creal(z)));
+	return z;
+}

libc/musl/src/complex/cprojf.c

@@ -0,0 +1,8 @@
+#include "libm.h"
+
+float complex cprojf(float complex z)
+{
+	if (isinf(crealf(z)) || isinf(cimagf(z)))
+		return CMPLXF(INFINITY, copysignf(0.0, crealf(z)));
+	return z;
+}

libc/musl/src/complex/cprojl.c

@@ -0,0 +1,15 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex cprojl(long double complex z)
+{
+	return cproj(z);
+}
+#else
+long double complex cprojl(long double complex z)
+{
+	if (isinf(creall(z)) || isinf(cimagl(z)))
+		return CMPLXL(INFINITY, copysignl(0.0, creall(z)));
+	return z;
+}
+#endif

libc/musl/src/complex/creal.c

@@ -0,0 +1,6 @@
+#include <complex.h>
+
+double (creal)(double complex z)
+{
+	return creal(z);
+}

libc/musl/src/complex/crealf.c

@@ -0,0 +1,6 @@
+#include <complex.h>
+
+float (crealf)(float complex z)
+{
+	return crealf(z);
+}

libc/musl/src/complex/creall.c

@@ -0,0 +1,6 @@
+#include <complex.h>
+
+long double (creall)(long double complex z)
+{
+	return creall(z);
+}

libc/musl/src/complex/csin.c

@@ -0,0 +1,9 @@
+#include "libm.h"
+
+/* sin(z) = -i sinh(i z) */
+
+double complex csin(double complex z)
+{
+	z = csinh(CMPLX(-cimag(z), creal(z)));
+	return CMPLX(cimag(z), -creal(z));
+}

libc/musl/src/complex/csinf.c

@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float complex csinf(float complex z)
+{
+	z = csinhf(CMPLXF(-cimagf(z), crealf(z)));
+	return CMPLXF(cimagf(z), -crealf(z));
+}

libc/musl/src/complex/csinh.c

@@ -0,0 +1,141 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_csinh.c */
+/*-
+ * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+/*
+ * Hyperbolic sine of a complex argument z = x + i y.
+ *
+ * sinh(z) = sinh(x+iy)
+ *         = sinh(x) cos(y) + i cosh(x) sin(y).
+ *
+ * Exceptional values are noted in the comments within the source code.
+ * These values and the return value were taken from n1124.pdf.
+ */
+
+#include "libm.h"
+
+static const double huge = 0x1p1023;
+
+double complex csinh(double complex z)
+{
+	double x, y, h;
+	int32_t hx, hy, ix, iy, lx, ly;
+
+	x = creal(z);
+	y = cimag(z);
+
+	EXTRACT_WORDS(hx, lx, x);
+	EXTRACT_WORDS(hy, ly, y);
+
+	ix = 0x7fffffff & hx;
+	iy = 0x7fffffff & hy;
+
+	/* Handle the nearly-non-exceptional cases where x and y are finite. */
+	if (ix < 0x7ff00000 && iy < 0x7ff00000) {
+		if ((iy | ly) == 0)
+			return CMPLX(sinh(x), y);
+		if (ix < 0x40360000)    /* small x: normal case */
+			return CMPLX(sinh(x) * cos(y), cosh(x) * sin(y));
+
+		/* |x| >= 22, so cosh(x) ~= exp(|x|) */
+		if (ix < 0x40862e42) {
+			/* x < 710: exp(|x|) won't overflow */
+			h = exp(fabs(x)) * 0.5;
+			return CMPLX(copysign(h, x) * cos(y), h * sin(y));
+		} else if (ix < 0x4096bbaa) {
+			/* x < 1455: scale to avoid overflow */
+			z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
+			return CMPLX(creal(z) * copysign(1, x), cimag(z));
+		} else {
+			/* x >= 1455: the result always overflows */
+			h = huge * x;
+			return CMPLX(h * cos(y), h * h * sin(y));
+		}
+	}
+
+	/*
+	 * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN.
+	 * The sign of 0 in the result is unspecified.  Choice = normally
+	 * the same as dNaN.  Raise the invalid floating-point exception.
+	 *
+	 * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN).
+	 * The sign of 0 in the result is unspecified.  Choice = normally
+	 * the same as d(NaN).
+	 */
+	if ((ix | lx) == 0 && iy >= 0x7ff00000)
+		return CMPLX(copysign(0, x * (y - y)), y - y);
+
+	/*
+	 * sinh(+-Inf +- I 0) = +-Inf + I +-0.
+	 *
+	 * sinh(NaN +- I 0)   = d(NaN) + I +-0.
+	 */
+	if ((iy | ly) == 0 && ix >= 0x7ff00000) {
+		if (((hx & 0xfffff) | lx) == 0)
+			return CMPLX(x, y);
+		return CMPLX(x, copysign(0, y));
+	}
+
+	/*
+	 * sinh(x +- I Inf) = dNaN + I dNaN.
+	 * Raise the invalid floating-point exception for finite nonzero x.
+	 *
+	 * sinh(x + I NaN) = d(NaN) + I d(NaN).
+	 * Optionally raises the invalid floating-point exception for finite
+	 * nonzero x.  Choice = don't raise (except for signaling NaNs).
+	 */
+	if (ix < 0x7ff00000 && iy >= 0x7ff00000)
+		return CMPLX(y - y, x * (y - y));
+
+	/*
+	 * sinh(+-Inf + I NaN)  = +-Inf + I d(NaN).
+	 * The sign of Inf in the result is unspecified.  Choice = normally
+	 * the same as d(NaN).
+	 *
+	 * sinh(+-Inf +- I Inf) = +Inf + I dNaN.
+	 * The sign of Inf in the result is unspecified.  Choice = always +.
+	 * Raise the invalid floating-point exception.
+	 *
+	 * sinh(+-Inf + I y)   = +-Inf cos(y) + I Inf sin(y)
+	 */
+	if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
+		if (iy >= 0x7ff00000)
+			return CMPLX(x * x, x * (y - y));
+		return CMPLX(x * cos(y), INFINITY * sin(y));
+	}
+
+	/*
+	 * sinh(NaN + I NaN)  = d(NaN) + I d(NaN).
+	 *
+	 * sinh(NaN +- I Inf) = d(NaN) + I d(NaN).
+	 * Optionally raises the invalid floating-point exception.
+	 * Choice = raise.
+	 *
+	 * sinh(NaN + I y)    = d(NaN) + I d(NaN).
+	 * Optionally raises the invalid floating-point exception for finite
+	 * nonzero y.  Choice = don't raise (except for signaling NaNs).
+	 */
+	return CMPLX((x * x) * (y - y), (x + x) * (y - y));
+}

libc/musl/src/complex/csinhf.c

@@ -0,0 +1,90 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_csinhf.c */
+/*-
+ * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+/*
+ * Hyperbolic sine of a complex argument z.  See s_csinh.c for details.
+ */
+
+#include "libm.h"
+
+static const float huge = 0x1p127;
+
+float complex csinhf(float complex z)
+{
+	float x, y, h;
+	int32_t hx, hy, ix, iy;
+
+	x = crealf(z);
+	y = cimagf(z);
+
+	GET_FLOAT_WORD(hx, x);
+	GET_FLOAT_WORD(hy, y);
+
+	ix = 0x7fffffff & hx;
+	iy = 0x7fffffff & hy;
+
+	if (ix < 0x7f800000 && iy < 0x7f800000) {
+		if (iy == 0)
+			return CMPLXF(sinhf(x), y);
+		if (ix < 0x41100000)    /* small x: normal case */
+			return CMPLXF(sinhf(x) * cosf(y), coshf(x) * sinf(y));
+
+		/* |x| >= 9, so cosh(x) ~= exp(|x|) */
+		if (ix < 0x42b17218) {
+			/* x < 88.7: expf(|x|) won't overflow */
+			h = expf(fabsf(x)) * 0.5f;
+			return CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y));
+		} else if (ix < 0x4340b1e7) {
+			/* x < 192.7: scale to avoid overflow */
+			z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
+			return CMPLXF(crealf(z) * copysignf(1, x), cimagf(z));
+		} else {
+			/* x >= 192.7: the result always overflows */
+			h = huge * x;
+			return CMPLXF(h * cosf(y), h * h * sinf(y));
+		}
+	}
+
+	if (ix == 0 && iy >= 0x7f800000)
+		return CMPLXF(copysignf(0, x * (y - y)), y - y);
+
+	if (iy == 0 && ix >= 0x7f800000) {
+		if ((hx & 0x7fffff) == 0)
+			return CMPLXF(x, y);
+		return CMPLXF(x, copysignf(0, y));
+	}
+
+	if (ix < 0x7f800000 && iy >= 0x7f800000)
+		return CMPLXF(y - y, x * (y - y));
+
+	if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
+		if (iy >= 0x7f800000)
+			return CMPLXF(x * x, x * (y - y));
+		return CMPLXF(x * cosf(y), INFINITY * sinf(y));
+	}
+
+	return CMPLXF((x * x) * (y - y), (x + x) * (y - y));
+}

libc/musl/src/complex/csinhl.c

@@ -0,0 +1,7 @@
+#include "libm.h"
+
+//FIXME
+long double complex csinhl(long double complex z)
+{
+	return csinh(z);
+}

libc/musl/src/complex/csinl.c

@@ -0,0 +1,14 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex csinl(long double complex z)
+{
+	return csin(z);
+}
+#else
+long double complex csinl(long double complex z)
+{
+	z = csinhl(CMPLXL(-cimagl(z), creall(z)));
+	return CMPLXL(cimagl(z), -creall(z));
+}
+#endif

libc/musl/src/complex/csqrt.c

@@ -0,0 +1,100 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_csqrt.c */
+/*-
+ * Copyright (c) 2007 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include "libm.h"
+
+/*
+ * gcc doesn't implement complex multiplication or division correctly,
+ * so we need to handle infinities specially. We turn on this pragma to
+ * notify conforming c99 compilers that the fast-but-incorrect code that
+ * gcc generates is acceptable, since the special cases have already been
+ * handled.
+ */
+#pragma STDC CX_LIMITED_RANGE ON
+
+/* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */
+#define THRESH  0x1.a827999fcef32p+1022
+
+double complex csqrt(double complex z)
+{
+	double complex result;
+	double a, b;
+	double t;
+	int scale;
+
+	a = creal(z);
+	b = cimag(z);
+
+	/* Handle special cases. */
+	if (z == 0)
+		return CMPLX(0, b);
+	if (isinf(b))
+		return CMPLX(INFINITY, b);
+	if (isnan(a)) {
+		t = (b - b) / (b - b);  /* raise invalid if b is not a NaN */
+		return CMPLX(a, t);   /* return NaN + NaN i */
+	}
+	if (isinf(a)) {
+		/*
+		 * csqrt(inf + NaN i)  = inf +  NaN i
+		 * csqrt(inf + y i)    = inf +  0 i
+		 * csqrt(-inf + NaN i) = NaN +- inf i
+		 * csqrt(-inf + y i)   = 0   +  inf i
+		 */
+		if (signbit(a))
+			return CMPLX(fabs(b - b), copysign(a, b));
+		else
+			return CMPLX(a, copysign(b - b, b));
+	}
+	/*
+	 * The remaining special case (b is NaN) is handled just fine by
+	 * the normal code path below.
+	 */
+
+	/* Scale to avoid overflow. */
+	if (fabs(a) >= THRESH || fabs(b) >= THRESH) {
+		a *= 0.25;
+		b *= 0.25;
+		scale = 1;
+	} else {
+		scale = 0;
+	}
+
+	/* Algorithm 312, CACM vol 10, Oct 1967. */
+	if (a >= 0) {
+		t = sqrt((a + hypot(a, b)) * 0.5);
+		result = CMPLX(t, b / (2 * t));
+	} else {
+		t = sqrt((-a + hypot(a, b)) * 0.5);
+		result = CMPLX(fabs(b) / (2 * t), copysign(t, b));
+	}
+
+	/* Rescale. */
+	if (scale)
+		result *= 2;
+	return result;
+}

libc/musl/src/complex/csqrtf.c

@@ -0,0 +1,82 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_csqrtf.c */
+/*-
+ * Copyright (c) 2007 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include "libm.h"
+
+/*
+ * gcc doesn't implement complex multiplication or division correctly,
+ * so we need to handle infinities specially. We turn on this pragma to
+ * notify conforming c99 compilers that the fast-but-incorrect code that
+ * gcc generates is acceptable, since the special cases have already been
+ * handled.
+ */
+#pragma STDC CX_LIMITED_RANGE ON
+
+float complex csqrtf(float complex z)
+{
+	float a = crealf(z), b = cimagf(z);
+	double t;
+
+	/* Handle special cases. */
+	if (z == 0)
+		return CMPLXF(0, b);
+	if (isinf(b))
+		return CMPLXF(INFINITY, b);
+	if (isnan(a)) {
+		t = (b - b) / (b - b);  /* raise invalid if b is not a NaN */
+		return CMPLXF(a, t);  /* return NaN + NaN i */
+	}
+	if (isinf(a)) {
+		/*
+		 * csqrtf(inf + NaN i)  = inf +  NaN i
+		 * csqrtf(inf + y i)    = inf +  0 i
+		 * csqrtf(-inf + NaN i) = NaN +- inf i
+		 * csqrtf(-inf + y i)   = 0   +  inf i
+		 */
+		if (signbit(a))
+			return CMPLXF(fabsf(b - b), copysignf(a, b));
+		else
+			return CMPLXF(a, copysignf(b - b, b));
+	}
+	/*
+	 * The remaining special case (b is NaN) is handled just fine by
+	 * the normal code path below.
+	 */
+
+	/*
+	 * We compute t in double precision to avoid overflow and to
+	 * provide correct rounding in nearly all cases.
+	 * This is Algorithm 312, CACM vol 10, Oct 1967.
+	 */
+	if (a >= 0) {
+		t = sqrt((a + hypot(a, b)) * 0.5);
+		return CMPLXF(t, b / (2.0 * t));
+	} else {
+		t = sqrt((-a + hypot(a, b)) * 0.5);
+		return CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b));
+	}
+}

libc/musl/src/complex/csqrtl.c

@@ -0,0 +1,7 @@
+#include "libm.h"
+
+//FIXME
+long double complex csqrtl(long double complex z)
+{
+	return csqrt(z);
+}

libc/musl/src/complex/ctan.c

@@ -0,0 +1,9 @@
+#include "libm.h"
+
+/* tan(z) = -i tanh(i z) */
+
+double complex ctan(double complex z)
+{
+	z = ctanh(CMPLX(-cimag(z), creal(z)));
+	return CMPLX(cimag(z), -creal(z));
+}

libc/musl/src/complex/ctanf.c

@@ -0,0 +1,7 @@
+#include "libm.h"
+
+float complex ctanf(float complex z)
+{
+	z = ctanhf(CMPLXF(-cimagf(z), crealf(z)));
+	return CMPLXF(cimagf(z), -crealf(z));
+}

libc/musl/src/complex/ctanh.c

@@ -0,0 +1,129 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_ctanh.c */
+/*-
+ * Copyright (c) 2011 David Schultz
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+/*
+ * Hyperbolic tangent of a complex argument z = x + i y.
+ *
+ * The algorithm is from:
+ *
+ *   W. Kahan.  Branch Cuts for Complex Elementary Functions or Much
+ *   Ado About Nothing's Sign Bit.  In The State of the Art in
+ *   Numerical Analysis, pp. 165 ff.  Iserles and Powell, eds., 1987.
+ *
+ * Method:
+ *
+ *   Let t    = tan(x)
+ *       beta = 1/cos^2(y)
+ *       s    = sinh(x)
+ *       rho  = cosh(x)
+ *
+ *   We have:
+ *
+ *   tanh(z) = sinh(z) / cosh(z)
+ *
+ *             sinh(x) cos(y) + i cosh(x) sin(y)
+ *           = ---------------------------------
+ *             cosh(x) cos(y) + i sinh(x) sin(y)
+ *
+ *             cosh(x) sinh(x) / cos^2(y) + i tan(y)
+ *           = -------------------------------------
+ *                    1 + sinh^2(x) / cos^2(y)
+ *
+ *             beta rho s + i t
+ *           = ----------------
+ *               1 + beta s^2
+ *
+ * Modifications:
+ *
+ *   I omitted the original algorithm's handling of overflow in tan(x) after
+ *   verifying with nearpi.c that this can't happen in IEEE single or double
+ *   precision.  I also handle large x differently.
+ */
+
+#include "libm.h"
+
+double complex ctanh(double complex z)
+{
+	double x, y;
+	double t, beta, s, rho, denom;
+	uint32_t hx, ix, lx;
+
+	x = creal(z);
+	y = cimag(z);
+
+	EXTRACT_WORDS(hx, lx, x);
+	ix = hx & 0x7fffffff;
+
+	/*
+	 * ctanh(NaN + i 0) = NaN + i 0
+	 *
+	 * ctanh(NaN + i y) = NaN + i NaN               for y != 0
+	 *
+	 * The imaginary part has the sign of x*sin(2*y), but there's no
+	 * special effort to get this right.
+	 *
+	 * ctanh(+-Inf +- i Inf) = +-1 +- 0
+	 *
+	 * ctanh(+-Inf + i y) = +-1 + 0 sin(2y)         for y finite
+	 *
+	 * The imaginary part of the sign is unspecified.  This special
+	 * case is only needed to avoid a spurious invalid exception when
+	 * y is infinite.
+	 */
+	if (ix >= 0x7ff00000) {
+		if ((ix & 0xfffff) | lx)        /* x is NaN */
+			return CMPLX(x, (y == 0 ? y : x * y));
+		SET_HIGH_WORD(x, hx - 0x40000000);      /* x = copysign(1, x) */
+		return CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)));
+	}
+
+	/*
+	 * ctanh(+-0 + i NAN) = +-0 + i NaN
+	 * ctanh(+-0 +- i Inf) = +-0 + i NaN
+	 * ctanh(x + i NAN) = NaN + i NaN
+	 * ctanh(x +- i Inf) = NaN + i NaN
+	 */
+	if (!isfinite(y))
+		return CMPLX(x ? y - y : x, y - y);
+
+	/*
+	 * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the
+	 * approximation sinh^2(huge) ~= exp(2*huge) / 4.
+	 * We use a modified formula to avoid spurious overflow.
+	 */
+	if (ix >= 0x40360000) { /* x >= 22 */
+		double exp_mx = exp(-fabs(x));
+		return CMPLX(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx);
+	}
+
+	/* Kahan's algorithm */
+	t = tan(y);
+	beta = 1.0 + t * t;     /* = 1 / cos^2(y) */
+	s = sinh(x);
+	rho = sqrt(1 + s * s);  /* = cosh(x) */
+	denom = 1 + beta * s * s;
+	return CMPLX((beta * rho * s) / denom, t / denom);
+}

libc/musl/src/complex/ctanhf.c

@@ -0,0 +1,66 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_ctanhf.c */
+/*-
+ * Copyright (c) 2011 David Schultz
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+/*
+ * Hyperbolic tangent of a complex argument z.  See s_ctanh.c for details.
+ */
+
+#include "libm.h"
+
+float complex ctanhf(float complex z)
+{
+	float x, y;
+	float t, beta, s, rho, denom;
+	uint32_t hx, ix;
+
+	x = crealf(z);
+	y = cimagf(z);
+
+	GET_FLOAT_WORD(hx, x);
+	ix = hx & 0x7fffffff;
+
+	if (ix >= 0x7f800000) {
+		if (ix & 0x7fffff)
+			return CMPLXF(x, (y == 0 ? y : x * y));
+		SET_FLOAT_WORD(x, hx - 0x40000000);
+		return CMPLXF(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)));
+	}
+
+	if (!isfinite(y))
+		return CMPLXF(ix ? y - y : x, y - y);
+
+	if (ix >= 0x41300000) { /* x >= 11 */
+		float exp_mx = expf(-fabsf(x));
+		return CMPLXF(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx);
+	}
+
+	t = tanf(y);
+	beta = 1.0 + t * t;
+	s = sinhf(x);
+	rho = sqrtf(1 + s * s);
+	denom = 1 + beta * s * s;
+	return CMPLXF((beta * rho * s) / denom, t / denom);
+}

libc/musl/src/complex/ctanhl.c

@@ -0,0 +1,7 @@
+#include "libm.h"
+
+//FIXME
+long double complex ctanhl(long double complex z)
+{
+	return ctanh(z);
+}

libc/musl/src/complex/ctanl.c

@@ -0,0 +1,14 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double complex ctanl(long double complex z)
+{
+	return ctan(z);
+}
+#else
+long double complex ctanl(long double complex z)
+{
+	z = ctanhl(CMPLXL(-cimagl(z), creall(z)));
+	return CMPLXL(cimagl(z), -creall(z));
+}
+#endif

libc/musl/src/conf/confstr.c

@@ -0,0 +1,17 @@
+#include <unistd.h>
+#include <stdio.h>
+#include <errno.h>
+
+size_t confstr(int name, char *buf, size_t len)
+{
+	const char *s = "";
+	if (!name) {
+		s = "/bin:/usr/bin";
+	} else if ((name&~4U)!=1 && name-_CS_POSIX_V6_ILP32_OFF32_CFLAGS>33U) {
+		errno = EINVAL;
+		return 0;
+	}
+	// snprintf is overkill but avoid wasting code size to implement
+	// this completely useless function and its truncation semantics
+	return snprintf(buf, len, "%s", s) + 1;
+}

libc/musl/src/conf/fpathconf.c

@@ -0,0 +1,35 @@
+#include <unistd.h>
+#include <limits.h>
+#include <errno.h>
+
+long fpathconf(int fd, int name)
+{
+	static const short values[] = {
+		[_PC_LINK_MAX] = _POSIX_LINK_MAX,
+		[_PC_MAX_CANON] = _POSIX_MAX_CANON,
+		[_PC_MAX_INPUT] = _POSIX_MAX_INPUT,
+		[_PC_NAME_MAX] = NAME_MAX,
+		[_PC_PATH_MAX] = PATH_MAX,
+		[_PC_PIPE_BUF] = PIPE_BUF,
+		[_PC_CHOWN_RESTRICTED] = 1,
+		[_PC_NO_TRUNC] = 1,
+		[_PC_VDISABLE] = 0,
+		[_PC_SYNC_IO] = 1,
+		[_PC_ASYNC_IO] = -1,
+		[_PC_PRIO_IO] = -1,
+		[_PC_SOCK_MAXBUF] = -1,
+		[_PC_FILESIZEBITS] = FILESIZEBITS,
+		[_PC_REC_INCR_XFER_SIZE] = 4096,
+		[_PC_REC_MAX_XFER_SIZE] = 4096,
+		[_PC_REC_MIN_XFER_SIZE] = 4096,
+		[_PC_REC_XFER_ALIGN] = 4096,
+		[_PC_ALLOC_SIZE_MIN] = 4096,
+		[_PC_SYMLINK_MAX] = -1,
+		[_PC_2_SYMLINKS] = 1
+	};
+	if (name >= sizeof(values)/sizeof(values[0])) {
+		errno = EINVAL;
+		return -1;
+	}
+	return values[name];
+}

libc/musl/src/conf/legacy.c

@@ -0,0 +1,22 @@
+#include <sys/sysinfo.h>
+#include <unistd.h>
+
+int get_nprocs_conf()
+{
+	return sysconf(_SC_NPROCESSORS_CONF);
+}
+
+int get_nprocs()
+{
+	return sysconf(_SC_NPROCESSORS_ONLN);
+}
+
+long get_phys_pages()
+{
+	return sysconf(_SC_PHYS_PAGES);	
+}
+
+long get_avphys_pages()
+{
+	return sysconf(_SC_AVPHYS_PAGES);	
+}

libc/musl/src/conf/pathconf.c

@@ -0,0 +1,6 @@
+#include <unistd.h>
+
+long pathconf(const char *path, int name)
+{
+	return fpathconf(-1, name);
+}

libc/musl/src/conf/sysconf.c

@@ -0,0 +1,219 @@
+#include <unistd.h>
+#include <limits.h>
+#include <errno.h>
+#include <sys/resource.h>
+#include <signal.h>
+#include <sys/sysinfo.h>
+#include "syscall.h"
+#include "libc.h"
+
+#define JT(x) (-256|(x))
+#define VER JT(1)
+#define JT_ARG_MAX JT(2)
+#define JT_MQ_PRIO_MAX JT(3)
+#define JT_PAGE_SIZE JT(4)
+#define JT_SEM_VALUE_MAX JT(5)
+#define JT_NPROCESSORS_CONF JT(6)
+#define JT_NPROCESSORS_ONLN JT(7)
+#define JT_PHYS_PAGES JT(8)
+#define JT_AVPHYS_PAGES JT(9)
+#define JT_ZERO JT(10)
+#define JT_DELAYTIMER_MAX JT(11)
+
+#define RLIM(x) (-32768|(RLIMIT_ ## x))
+
+long sysconf(int name)
+{
+	static const short values[] = {
+		[_SC_ARG_MAX] = JT_ARG_MAX,
+		[_SC_CHILD_MAX] = RLIM(NPROC),
+		[_SC_CLK_TCK] = 100,
+		[_SC_NGROUPS_MAX] = 32,
+		[_SC_OPEN_MAX] = RLIM(NOFILE),
+		[_SC_STREAM_MAX] = -1,
+		[_SC_TZNAME_MAX] = TZNAME_MAX,
+		[_SC_JOB_CONTROL] = 1,
+		[_SC_SAVED_IDS] = 1,
+		[_SC_REALTIME_SIGNALS] = VER,
+		[_SC_PRIORITY_SCHEDULING] = -1,
+		[_SC_TIMERS] = VER,
+		[_SC_ASYNCHRONOUS_IO] = VER,
+		[_SC_PRIORITIZED_IO] = -1,
+		[_SC_SYNCHRONIZED_IO] = -1,
+		[_SC_FSYNC] = VER,
+		[_SC_MAPPED_FILES] = VER,
+		[_SC_MEMLOCK] = VER,
+		[_SC_MEMLOCK_RANGE] = VER,
+		[_SC_MEMORY_PROTECTION] = VER,
+		[_SC_MESSAGE_PASSING] = VER,
+		[_SC_SEMAPHORES] = VER,
+		[_SC_SHARED_MEMORY_OBJECTS] = VER,
+		[_SC_AIO_LISTIO_MAX] = -1,
+		[_SC_AIO_MAX] = -1,
+		[_SC_AIO_PRIO_DELTA_MAX] = JT_ZERO, /* ?? */
+		[_SC_DELAYTIMER_MAX] = JT_DELAYTIMER_MAX,
+		[_SC_MQ_OPEN_MAX] = -1,
+		[_SC_MQ_PRIO_MAX] = JT_MQ_PRIO_MAX,
+		[_SC_VERSION] = VER,
+		[_SC_PAGE_SIZE] = JT_PAGE_SIZE,
+		[_SC_RTSIG_MAX] = _NSIG - 1 - 31 - 3,
+		[_SC_SEM_NSEMS_MAX] = SEM_NSEMS_MAX,
+		[_SC_SEM_VALUE_MAX] = JT_SEM_VALUE_MAX,
+		[_SC_SIGQUEUE_MAX] = -1,
+		[_SC_TIMER_MAX] = -1,
+		[_SC_BC_BASE_MAX] = _POSIX2_BC_BASE_MAX,
+		[_SC_BC_DIM_MAX] = _POSIX2_BC_DIM_MAX,
+		[_SC_BC_SCALE_MAX] = _POSIX2_BC_SCALE_MAX,
+		[_SC_BC_STRING_MAX] = _POSIX2_BC_STRING_MAX,
+		[_SC_COLL_WEIGHTS_MAX] = COLL_WEIGHTS_MAX,
+		[_SC_EXPR_NEST_MAX] = -1,
+		[_SC_LINE_MAX] = -1,
+		[_SC_RE_DUP_MAX] = RE_DUP_MAX,
+		[_SC_2_VERSION] = VER,
+		[_SC_2_C_BIND] = VER,
+		[_SC_2_C_DEV] = -1,
+		[_SC_2_FORT_DEV] = -1,
+		[_SC_2_FORT_RUN] = -1,
+		[_SC_2_SW_DEV] = -1,
+		[_SC_2_LOCALEDEF] = -1,
+		[_SC_IOV_MAX] = IOV_MAX,
+		[_SC_THREADS] = VER,
+		[_SC_THREAD_SAFE_FUNCTIONS] = VER,
+		[_SC_GETGR_R_SIZE_MAX] = -1,
+		[_SC_GETPW_R_SIZE_MAX] = -1,
+		[_SC_LOGIN_NAME_MAX] = 256,
+		[_SC_TTY_NAME_MAX] = TTY_NAME_MAX,
+		[_SC_THREAD_DESTRUCTOR_ITERATIONS] = PTHREAD_DESTRUCTOR_ITERATIONS,
+		[_SC_THREAD_KEYS_MAX] = PTHREAD_KEYS_MAX,
+		[_SC_THREAD_STACK_MIN] = PTHREAD_STACK_MIN,
+		[_SC_THREAD_THREADS_MAX] = -1,
+		[_SC_THREAD_ATTR_STACKADDR] = VER,
+		[_SC_THREAD_ATTR_STACKSIZE] = VER,
+		[_SC_THREAD_PRIORITY_SCHEDULING] = VER,
+		[_SC_THREAD_PRIO_INHERIT] = -1,
+		[_SC_THREAD_PRIO_PROTECT] = -1,
+		[_SC_THREAD_PROCESS_SHARED] = VER,
+		[_SC_NPROCESSORS_CONF] = JT_NPROCESSORS_CONF,
+		[_SC_NPROCESSORS_ONLN] = JT_NPROCESSORS_ONLN,
+		[_SC_PHYS_PAGES] = JT_PHYS_PAGES,
+		[_SC_AVPHYS_PAGES] = JT_AVPHYS_PAGES,
+		[_SC_ATEXIT_MAX] = -1,
+		[_SC_PASS_MAX] = -1,
+		[_SC_XOPEN_VERSION] = _XOPEN_VERSION,
+		[_SC_XOPEN_XCU_VERSION] = _XOPEN_VERSION,
+		[_SC_XOPEN_UNIX] = 1,
+		[_SC_XOPEN_CRYPT] = -1,
+		[_SC_XOPEN_ENH_I18N] = 1,
+		[_SC_XOPEN_SHM] = 1,
+		[_SC_2_CHAR_TERM] = -1,
+		[_SC_2_UPE] = -1,
+		[_SC_XOPEN_XPG2] = -1,
+		[_SC_XOPEN_XPG3] = -1,
+		[_SC_XOPEN_XPG4] = -1,
+		[_SC_NZERO] = NZERO,
+		[_SC_XBS5_ILP32_OFF32] = -1,
+		[_SC_XBS5_ILP32_OFFBIG] = sizeof(long)==4 ? 1 : -1,
+		[_SC_XBS5_LP64_OFF64] = sizeof(long)==8 ? 1 : -1,
+		[_SC_XBS5_LPBIG_OFFBIG] = -1,
+		[_SC_XOPEN_LEGACY] = -1,
+		[_SC_XOPEN_REALTIME] = -1,
+		[_SC_XOPEN_REALTIME_THREADS] = -1,
+		[_SC_ADVISORY_INFO] = VER,
+		[_SC_BARRIERS] = VER,
+		[_SC_CLOCK_SELECTION] = VER,
+		[_SC_CPUTIME] = VER,
+		[_SC_THREAD_CPUTIME] = VER,
+		[_SC_MONOTONIC_CLOCK] = VER,
+		[_SC_READER_WRITER_LOCKS] = VER,
+		[_SC_SPIN_LOCKS] = VER,
+		[_SC_REGEXP] = 1,
+		[_SC_SHELL] = 1,
+		[_SC_SPAWN] = VER,
+		[_SC_SPORADIC_SERVER] = -1,
+		[_SC_THREAD_SPORADIC_SERVER] = -1,
+		[_SC_TIMEOUTS] = VER,
+		[_SC_TYPED_MEMORY_OBJECTS] = -1,
+		[_SC_2_PBS] = -1,
+		[_SC_2_PBS_ACCOUNTING] = -1,
+		[_SC_2_PBS_LOCATE] = -1,
+		[_SC_2_PBS_MESSAGE] = -1,
+		[_SC_2_PBS_TRACK] = -1,
+		[_SC_SYMLOOP_MAX] = SYMLOOP_MAX,
+		[_SC_STREAMS] = JT_ZERO,
+		[_SC_2_PBS_CHECKPOINT] = -1,
+		[_SC_V6_ILP32_OFF32] = -1,
+		[_SC_V6_ILP32_OFFBIG] = sizeof(long)==4 ? 1 : -1,
+		[_SC_V6_LP64_OFF64] = sizeof(long)==8 ? 1 : -1,
+		[_SC_V6_LPBIG_OFFBIG] = -1,
+		[_SC_HOST_NAME_MAX] = HOST_NAME_MAX,
+		[_SC_TRACE] = -1,
+		[_SC_TRACE_EVENT_FILTER] = -1,
+		[_SC_TRACE_INHERIT] = -1,
+		[_SC_TRACE_LOG] = -1,
+
+		[_SC_IPV6] = VER,
+		[_SC_RAW_SOCKETS] = VER,
+		[_SC_V7_ILP32_OFF32] = -1,
+		[_SC_V7_ILP32_OFFBIG] = sizeof(long)==4 ? 1 : -1,
+		[_SC_V7_LP64_OFF64] = sizeof(long)==8 ? 1 : -1,
+		[_SC_V7_LPBIG_OFFBIG] = -1,
+		[_SC_SS_REPL_MAX] = -1,
+		[_SC_TRACE_EVENT_NAME_MAX] = -1,
+		[_SC_TRACE_NAME_MAX] = -1,
+		[_SC_TRACE_SYS_MAX] = -1,
+		[_SC_TRACE_USER_EVENT_MAX] = -1,
+		[_SC_XOPEN_STREAMS] = JT_ZERO,
+		[_SC_THREAD_ROBUST_PRIO_INHERIT] = -1,
+		[_SC_THREAD_ROBUST_PRIO_PROTECT] = -1,
+	};
+
+	if (name >= sizeof(values)/sizeof(values[0]) || !values[name]) {
+		errno = EINVAL;
+		return -1;
+	} else if (values[name] >= -1) {
+		return values[name];
+	} else if (values[name] < -256) {
+		struct rlimit lim;
+		getrlimit(values[name]&16383, &lim);
+		if (lim.rlim_cur == RLIM_INFINITY)
+			return -1;
+		return lim.rlim_cur > LONG_MAX ? LONG_MAX : lim.rlim_cur;
+	}
+
+	switch ((unsigned char)values[name]) {
+	case VER & 255:
+		return _POSIX_VERSION;
+	case JT_ARG_MAX & 255:
+		return ARG_MAX;
+	case JT_MQ_PRIO_MAX & 255:
+		return MQ_PRIO_MAX;
+	case JT_PAGE_SIZE & 255:
+		return PAGE_SIZE;
+	case JT_SEM_VALUE_MAX & 255:
+		return SEM_VALUE_MAX;
+	case JT_DELAYTIMER_MAX & 255:
+		return DELAYTIMER_MAX;
+	case JT_NPROCESSORS_CONF & 255:
+	case JT_NPROCESSORS_ONLN & 255: ;
+		unsigned char set[128] = {1};
+		int i, cnt;
+		__syscall(SYS_sched_getaffinity, 0, sizeof set, set);
+		for (i=cnt=0; i<sizeof set; i++)
+			for (; set[i]; set[i]&=set[i]-1, cnt++);
+		return cnt;
+	case JT_PHYS_PAGES & 255:
+	case JT_AVPHYS_PAGES & 255: ;
+		unsigned long long mem;
+		struct sysinfo si;
+		__lsysinfo(&si);
+		if (!si.mem_unit) si.mem_unit = 1;
+		if (name==_SC_PHYS_PAGES) mem = si.totalram;
+		else mem = si.freeram + si.bufferram;
+		mem *= si.mem_unit;
+		mem /= PAGE_SIZE;
+		return (mem > LONG_MAX) ? LONG_MAX : mem;
+	case JT_ZERO & 255:
+		return 0;
+	}
+	return values[name];
+}

libc/musl/src/crypt/crypt.c

@@ -0,0 +1,14 @@
+#include <unistd.h>
+#include <crypt.h>
+
+char *crypt(const char *key, const char *salt)
+{
+	/* This buffer is sufficiently large for all
+	 * currently-supported hash types. It needs to be updated if
+	 * longer hashes are added. The cast to struct crypt_data * is
+	 * purely to meet the public API requirements of the crypt_r
+	 * function; the implementation of crypt_r uses the object
+	 * purely as a char buffer. */
+	static char buf[128];
+	return __crypt_r(key, salt, (struct crypt_data *)buf);
+}

libc/musl/src/crypt/crypt_blowfish.c

@@ -0,0 +1,798 @@
+/* Modified by Rich Felker in for inclusion in musl libc, based on
+ * Solar Designer's second size-optimized version sent to the musl
+ * mailing list. */
+
+/*
+ * The crypt_blowfish homepage is:
+ *
+ *	http://www.openwall.com/crypt/
+ *
+ * This code comes from John the Ripper password cracker, with reentrant
+ * and crypt(3) interfaces added, but optimizations specific to password
+ * cracking removed.
+ *
+ * Written by Solar Designer <solar at openwall.com> in 1998-2012.
+ * No copyright is claimed, and the software is hereby placed in the public
+ * domain.  In case this attempt to disclaim copyright and place the software
+ * in the public domain is deemed null and void, then the software is
+ * Copyright (c) 1998-2012 Solar Designer and it is hereby released to the
+ * general public under the following terms:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted.
+ *
+ * There's ABSOLUTELY NO WARRANTY, express or implied.
+ *
+ * It is my intent that you should be able to use this on your system,
+ * as part of a software package, or anywhere else to improve security,
+ * ensure compatibility, or for any other purpose.  I would appreciate
+ * it if you give credit where it is due and keep your modifications in
+ * the public domain as well, but I don't require that in order to let
+ * you place this code and any modifications you make under a license
+ * of your choice.
+ *
+ * This implementation is mostly compatible with OpenBSD's bcrypt.c (prefix
+ * "$2a$") by Niels Provos <provos at citi.umich.edu>, and uses some of his
+ * ideas.  The password hashing algorithm was designed by David Mazieres
+ * <dm at lcs.mit.edu>.  For more information on the level of compatibility,
+ * please refer to the comments in BF_set_key() below and to the included
+ * crypt(3) man page.
+ *
+ * There's a paper on the algorithm that explains its design decisions:
+ *
+ *	http://www.usenix.org/events/usenix99/provos.html
+ *
+ * Some of the tricks in BF_ROUND might be inspired by Eric Young's
+ * Blowfish library (I can't be sure if I would think of something if I
+ * hadn't seen his code).
+ */
+
+#include <string.h>
+#include <stdint.h>
+
+typedef uint32_t BF_word;
+typedef int32_t BF_word_signed;
+
+/* Number of Blowfish rounds, this is also hardcoded into a few places */
+#define BF_N				16
+
+typedef BF_word BF_key[BF_N + 2];
+
+typedef union {
+	struct {
+		BF_key P;
+		BF_word S[4][0x100];
+	} s;
+	BF_word PS[BF_N + 2 + 4 * 0x100];
+} BF_ctx;
+
+/*
+ * Magic IV for 64 Blowfish encryptions that we do at the end.
+ * The string is "OrpheanBeholderScryDoubt" on big-endian.
+ */
+static const BF_word BF_magic_w[6] = {
+	0x4F727068, 0x65616E42, 0x65686F6C,
+	0x64657253, 0x63727944, 0x6F756274
+};
+
+/*
+ * P-box and S-box tables initialized with digits of Pi.
+ */
+static const BF_ctx BF_init_state = {{
+	{
+		0x243f6a88, 0x85a308d3, 0x13198a2e, 0x03707344,
+		0xa4093822, 0x299f31d0, 0x082efa98, 0xec4e6c89,
+		0x452821e6, 0x38d01377, 0xbe5466cf, 0x34e90c6c,
+		0xc0ac29b7, 0xc97c50dd, 0x3f84d5b5, 0xb5470917,
+		0x9216d5d9, 0x8979fb1b
+	}, {
+		{
+			0xd1310ba6, 0x98dfb5ac, 0x2ffd72db, 0xd01adfb7,
+			0xb8e1afed, 0x6a267e96, 0xba7c9045, 0xf12c7f99,
+			0x24a19947, 0xb3916cf7, 0x0801f2e2, 0x858efc16,
+			0x636920d8, 0x71574e69, 0xa458fea3, 0xf4933d7e,
+			0x0d95748f, 0x728eb658, 0x718bcd58, 0x82154aee,
+			0x7b54a41d, 0xc25a59b5, 0x9c30d539, 0x2af26013,
+			0xc5d1b023, 0x286085f0, 0xca417918, 0xb8db38ef,
+			0x8e79dcb0, 0x603a180e, 0x6c9e0e8b, 0xb01e8a3e,
+			0xd71577c1, 0xbd314b27, 0x78af2fda, 0x55605c60,
+			0xe65525f3, 0xaa55ab94, 0x57489862, 0x63e81440,
+			0x55ca396a, 0x2aab10b6, 0xb4cc5c34, 0x1141e8ce,
+			0xa15486af, 0x7c72e993, 0xb3ee1411, 0x636fbc2a,
+			0x2ba9c55d, 0x741831f6, 0xce5c3e16, 0x9b87931e,
+			0xafd6ba33, 0x6c24cf5c, 0x7a325381, 0x28958677,
+			0x3b8f4898, 0x6b4bb9af, 0xc4bfe81b, 0x66282193,
+			0x61d809cc, 0xfb21a991, 0x487cac60, 0x5dec8032,
+			0xef845d5d, 0xe98575b1, 0xdc262302, 0xeb651b88,
+			0x23893e81, 0xd396acc5, 0x0f6d6ff3, 0x83f44239,
+			0x2e0b4482, 0xa4842004, 0x69c8f04a, 0x9e1f9b5e,
+			0x21c66842, 0xf6e96c9a, 0x670c9c61, 0xabd388f0,
+			0x6a51a0d2, 0xd8542f68, 0x960fa728, 0xab5133a3,
+			0x6eef0b6c, 0x137a3be4, 0xba3bf050, 0x7efb2a98,
+			0xa1f1651d, 0x39af0176, 0x66ca593e, 0x82430e88,
+			0x8cee8619, 0x456f9fb4, 0x7d84a5c3, 0x3b8b5ebe,
+			0xe06f75d8, 0x85c12073, 0x401a449f, 0x56c16aa6,
+			0x4ed3aa62, 0x363f7706, 0x1bfedf72, 0x429b023d,
+			0x37d0d724, 0xd00a1248, 0xdb0fead3, 0x49f1c09b,
+			0x075372c9, 0x80991b7b, 0x25d479d8, 0xf6e8def7,
+			0xe3fe501a, 0xb6794c3b, 0x976ce0bd, 0x04c006ba,
+			0xc1a94fb6, 0x409f60c4, 0x5e5c9ec2, 0x196a2463,
+			0x68fb6faf, 0x3e6c53b5, 0x1339b2eb, 0x3b52ec6f,
+			0x6dfc511f, 0x9b30952c, 0xcc814544, 0xaf5ebd09,
+			0xbee3d004, 0xde334afd, 0x660f2807, 0x192e4bb3,
+			0xc0cba857, 0x45c8740f, 0xd20b5f39, 0xb9d3fbdb,
+			0x5579c0bd, 0x1a60320a, 0xd6a100c6, 0x402c7279,
+			0x679f25fe, 0xfb1fa3cc, 0x8ea5e9f8, 0xdb3222f8,
+			0x3c7516df, 0xfd616b15, 0x2f501ec8, 0xad0552ab,
+			0x323db5fa, 0xfd238760, 0x53317b48, 0x3e00df82,
+			0x9e5c57bb, 0xca6f8ca0, 0x1a87562e, 0xdf1769db,
+			0xd542a8f6, 0x287effc3, 0xac6732c6, 0x8c4f5573,
+			0x695b27b0, 0xbbca58c8, 0xe1ffa35d, 0xb8f011a0,
+			0x10fa3d98, 0xfd2183b8, 0x4afcb56c, 0x2dd1d35b,
+			0x9a53e479, 0xb6f84565, 0xd28e49bc, 0x4bfb9790,
+			0xe1ddf2da, 0xa4cb7e33, 0x62fb1341, 0xcee4c6e8,
+			0xef20cada, 0x36774c01, 0xd07e9efe, 0x2bf11fb4,
+			0x95dbda4d, 0xae909198, 0xeaad8e71, 0x6b93d5a0,
+			0xd08ed1d0, 0xafc725e0, 0x8e3c5b2f, 0x8e7594b7,
+			0x8ff6e2fb, 0xf2122b64, 0x8888b812, 0x900df01c,
+			0x4fad5ea0, 0x688fc31c, 0xd1cff191, 0xb3a8c1ad,
+			0x2f2f2218, 0xbe0e1777, 0xea752dfe, 0x8b021fa1,
+			0xe5a0cc0f, 0xb56f74e8, 0x18acf3d6, 0xce89e299,
+			0xb4a84fe0, 0xfd13e0b7, 0x7cc43b81, 0xd2ada8d9,
+			0x165fa266, 0x80957705, 0x93cc7314, 0x211a1477,
+			0xe6ad2065, 0x77b5fa86, 0xc75442f5, 0xfb9d35cf,
+			0xebcdaf0c, 0x7b3e89a0, 0xd6411bd3, 0xae1e7e49,
+			0x00250e2d, 0x2071b35e, 0x226800bb, 0x57b8e0af,
+			0x2464369b, 0xf009b91e, 0x5563911d, 0x59dfa6aa,
+			0x78c14389, 0xd95a537f, 0x207d5ba2, 0x02e5b9c5,
+			0x83260376, 0x6295cfa9, 0x11c81968, 0x4e734a41,
+			0xb3472dca, 0x7b14a94a, 0x1b510052, 0x9a532915,
+			0xd60f573f, 0xbc9bc6e4, 0x2b60a476, 0x81e67400,
+			0x08ba6fb5, 0x571be91f, 0xf296ec6b, 0x2a0dd915,
+			0xb6636521, 0xe7b9f9b6, 0xff34052e, 0xc5855664,
+			0x53b02d5d, 0xa99f8fa1, 0x08ba4799, 0x6e85076a
+		}, {
+			0x4b7a70e9, 0xb5b32944, 0xdb75092e, 0xc4192623,
+			0xad6ea6b0, 0x49a7df7d, 0x9cee60b8, 0x8fedb266,
+			0xecaa8c71, 0x699a17ff, 0x5664526c, 0xc2b19ee1,
+			0x193602a5, 0x75094c29, 0xa0591340, 0xe4183a3e,
+			0x3f54989a, 0x5b429d65, 0x6b8fe4d6, 0x99f73fd6,
+			0xa1d29c07, 0xefe830f5, 0x4d2d38e6, 0xf0255dc1,
+			0x4cdd2086, 0x8470eb26, 0x6382e9c6, 0x021ecc5e,
+			0x09686b3f, 0x3ebaefc9, 0x3c971814, 0x6b6a70a1,
+			0x687f3584, 0x52a0e286, 0xb79c5305, 0xaa500737,
+			0x3e07841c, 0x7fdeae5c, 0x8e7d44ec, 0x5716f2b8,
+			0xb03ada37, 0xf0500c0d, 0xf01c1f04, 0x0200b3ff,
+			0xae0cf51a, 0x3cb574b2, 0x25837a58, 0xdc0921bd,
+			0xd19113f9, 0x7ca92ff6, 0x94324773, 0x22f54701,
+			0x3ae5e581, 0x37c2dadc, 0xc8b57634, 0x9af3dda7,
+			0xa9446146, 0x0fd0030e, 0xecc8c73e, 0xa4751e41,
+			0xe238cd99, 0x3bea0e2f, 0x3280bba1, 0x183eb331,
+			0x4e548b38, 0x4f6db908, 0x6f420d03, 0xf60a04bf,
+			0x2cb81290, 0x24977c79, 0x5679b072, 0xbcaf89af,
+			0xde9a771f, 0xd9930810, 0xb38bae12, 0xdccf3f2e,
+			0x5512721f, 0x2e6b7124, 0x501adde6, 0x9f84cd87,
+			0x7a584718, 0x7408da17, 0xbc9f9abc, 0xe94b7d8c,
+			0xec7aec3a, 0xdb851dfa, 0x63094366, 0xc464c3d2,
+			0xef1c1847, 0x3215d908, 0xdd433b37, 0x24c2ba16,
+			0x12a14d43, 0x2a65c451, 0x50940002, 0x133ae4dd,
+			0x71dff89e, 0x10314e55, 0x81ac77d6, 0x5f11199b,
+			0x043556f1, 0xd7a3c76b, 0x3c11183b, 0x5924a509,
+			0xf28fe6ed, 0x97f1fbfa, 0x9ebabf2c, 0x1e153c6e,
+			0x86e34570, 0xeae96fb1, 0x860e5e0a, 0x5a3e2ab3,
+			0x771fe71c, 0x4e3d06fa, 0x2965dcb9, 0x99e71d0f,
+			0x803e89d6, 0x5266c825, 0x2e4cc978, 0x9c10b36a,
+			0xc6150eba, 0x94e2ea78, 0xa5fc3c53, 0x1e0a2df4,
+			0xf2f74ea7, 0x361d2b3d, 0x1939260f, 0x19c27960,
+			0x5223a708, 0xf71312b6, 0xebadfe6e, 0xeac31f66,
+			0xe3bc4595, 0xa67bc883, 0xb17f37d1, 0x018cff28,
+			0xc332ddef, 0xbe6c5aa5, 0x65582185, 0x68ab9802,
+			0xeecea50f, 0xdb2f953b, 0x2aef7dad, 0x5b6e2f84,
+			0x1521b628, 0x29076170, 0xecdd4775, 0x619f1510,
+			0x13cca830, 0xeb61bd96, 0x0334fe1e, 0xaa0363cf,
+			0xb5735c90, 0x4c70a239, 0xd59e9e0b, 0xcbaade14,
+			0xeecc86bc, 0x60622ca7, 0x9cab5cab, 0xb2f3846e,
+			0x648b1eaf, 0x19bdf0ca, 0xa02369b9, 0x655abb50,
+			0x40685a32, 0x3c2ab4b3, 0x319ee9d5, 0xc021b8f7,
+			0x9b540b19, 0x875fa099, 0x95f7997e, 0x623d7da8,
+			0xf837889a, 0x97e32d77, 0x11ed935f, 0x16681281,
+			0x0e358829, 0xc7e61fd6, 0x96dedfa1, 0x7858ba99,
+			0x57f584a5, 0x1b227263, 0x9b83c3ff, 0x1ac24696,
+			0xcdb30aeb, 0x532e3054, 0x8fd948e4, 0x6dbc3128,
+			0x58ebf2ef, 0x34c6ffea, 0xfe28ed61, 0xee7c3c73,
+			0x5d4a14d9, 0xe864b7e3, 0x42105d14, 0x203e13e0,
+			0x45eee2b6, 0xa3aaabea, 0xdb6c4f15, 0xfacb4fd0,
+			0xc742f442, 0xef6abbb5, 0x654f3b1d, 0x41cd2105,
+			0xd81e799e, 0x86854dc7, 0xe44b476a, 0x3d816250,
+			0xcf62a1f2, 0x5b8d2646, 0xfc8883a0, 0xc1c7b6a3,
+			0x7f1524c3, 0x69cb7492, 0x47848a0b, 0x5692b285,
+			0x095bbf00, 0xad19489d, 0x1462b174, 0x23820e00,
+			0x58428d2a, 0x0c55f5ea, 0x1dadf43e, 0x233f7061,
+			0x3372f092, 0x8d937e41, 0xd65fecf1, 0x6c223bdb,
+			0x7cde3759, 0xcbee7460, 0x4085f2a7, 0xce77326e,
+			0xa6078084, 0x19f8509e, 0xe8efd855, 0x61d99735,
+			0xa969a7aa, 0xc50c06c2, 0x5a04abfc, 0x800bcadc,
+			0x9e447a2e, 0xc3453484, 0xfdd56705, 0x0e1e9ec9,
+			0xdb73dbd3, 0x105588cd, 0x675fda79, 0xe3674340,
+			0xc5c43465, 0x713e38d8, 0x3d28f89e, 0xf16dff20,
+			0x153e21e7, 0x8fb03d4a, 0xe6e39f2b, 0xdb83adf7
+		}, {
+			0xe93d5a68, 0x948140f7, 0xf64c261c, 0x94692934,
+			0x411520f7, 0x7602d4f7, 0xbcf46b2e, 0xd4a20068,
+			0xd4082471, 0x3320f46a, 0x43b7d4b7, 0x500061af,
+			0x1e39f62e, 0x97244546, 0x14214f74, 0xbf8b8840,
+			0x4d95fc1d, 0x96b591af, 0x70f4ddd3, 0x66a02f45,
+			0xbfbc09ec, 0x03bd9785, 0x7fac6dd0, 0x31cb8504,
+			0x96eb27b3, 0x55fd3941, 0xda2547e6, 0xabca0a9a,
+			0x28507825, 0x530429f4, 0x0a2c86da, 0xe9b66dfb,
+			0x68dc1462, 0xd7486900, 0x680ec0a4, 0x27a18dee,
+			0x4f3ffea2, 0xe887ad8c, 0xb58ce006, 0x7af4d6b6,
+			0xaace1e7c, 0xd3375fec, 0xce78a399, 0x406b2a42,
+			0x20fe9e35, 0xd9f385b9, 0xee39d7ab, 0x3b124e8b,
+			0x1dc9faf7, 0x4b6d1856, 0x26a36631, 0xeae397b2,
+			0x3a6efa74, 0xdd5b4332, 0x6841e7f7, 0xca7820fb,
+			0xfb0af54e, 0xd8feb397, 0x454056ac, 0xba489527,
+			0x55533a3a, 0x20838d87, 0xfe6ba9b7, 0xd096954b,
+			0x55a867bc, 0xa1159a58, 0xcca92963, 0x99e1db33,
+			0xa62a4a56, 0x3f3125f9, 0x5ef47e1c, 0x9029317c,
+			0xfdf8e802, 0x04272f70, 0x80bb155c, 0x05282ce3,
+			0x95c11548, 0xe4c66d22, 0x48c1133f, 0xc70f86dc,
+			0x07f9c9ee, 0x41041f0f, 0x404779a4, 0x5d886e17,
+			0x325f51eb, 0xd59bc0d1, 0xf2bcc18f, 0x41113564,
+			0x257b7834, 0x602a9c60, 0xdff8e8a3, 0x1f636c1b,
+			0x0e12b4c2, 0x02e1329e, 0xaf664fd1, 0xcad18115,
+			0x6b2395e0, 0x333e92e1, 0x3b240b62, 0xeebeb922,
+			0x85b2a20e, 0xe6ba0d99, 0xde720c8c, 0x2da2f728,
+			0xd0127845, 0x95b794fd, 0x647d0862, 0xe7ccf5f0,
+			0x5449a36f, 0x877d48fa, 0xc39dfd27, 0xf33e8d1e,
+			0x0a476341, 0x992eff74, 0x3a6f6eab, 0xf4f8fd37,
+			0xa812dc60, 0xa1ebddf8, 0x991be14c, 0xdb6e6b0d,
+			0xc67b5510, 0x6d672c37, 0x2765d43b, 0xdcd0e804,
+			0xf1290dc7, 0xcc00ffa3, 0xb5390f92, 0x690fed0b,
+			0x667b9ffb, 0xcedb7d9c, 0xa091cf0b, 0xd9155ea3,
+			0xbb132f88, 0x515bad24, 0x7b9479bf, 0x763bd6eb,
+			0x37392eb3, 0xcc115979, 0x8026e297, 0xf42e312d,
+			0x6842ada7, 0xc66a2b3b, 0x12754ccc, 0x782ef11c,
+			0x6a124237, 0xb79251e7, 0x06a1bbe6, 0x4bfb6350,
+			0x1a6b1018, 0x11caedfa, 0x3d25bdd8, 0xe2e1c3c9,
+			0x44421659, 0x0a121386, 0xd90cec6e, 0xd5abea2a,
+			0x64af674e, 0xda86a85f, 0xbebfe988, 0x64e4c3fe,
+			0x9dbc8057, 0xf0f7c086, 0x60787bf8, 0x6003604d,
+			0xd1fd8346, 0xf6381fb0, 0x7745ae04, 0xd736fccc,
+			0x83426b33, 0xf01eab71, 0xb0804187, 0x3c005e5f,
+			0x77a057be, 0xbde8ae24, 0x55464299, 0xbf582e61,
+			0x4e58f48f, 0xf2ddfda2, 0xf474ef38, 0x8789bdc2,
+			0x5366f9c3, 0xc8b38e74, 0xb475f255, 0x46fcd9b9,
+			0x7aeb2661, 0x8b1ddf84, 0x846a0e79, 0x915f95e2,
+			0x466e598e, 0x20b45770, 0x8cd55591, 0xc902de4c,
+			0xb90bace1, 0xbb8205d0, 0x11a86248, 0x7574a99e,
+			0xb77f19b6, 0xe0a9dc09, 0x662d09a1, 0xc4324633,
+			0xe85a1f02, 0x09f0be8c, 0x4a99a025, 0x1d6efe10,
+			0x1ab93d1d, 0x0ba5a4df, 0xa186f20f, 0x2868f169,
+			0xdcb7da83, 0x573906fe, 0xa1e2ce9b, 0x4fcd7f52,
+			0x50115e01, 0xa70683fa, 0xa002b5c4, 0x0de6d027,
+			0x9af88c27, 0x773f8641, 0xc3604c06, 0x61a806b5,
+			0xf0177a28, 0xc0f586e0, 0x006058aa, 0x30dc7d62,
+			0x11e69ed7, 0x2338ea63, 0x53c2dd94, 0xc2c21634,
+			0xbbcbee56, 0x90bcb6de, 0xebfc7da1, 0xce591d76,
+			0x6f05e409, 0x4b7c0188, 0x39720a3d, 0x7c927c24,
+			0x86e3725f, 0x724d9db9, 0x1ac15bb4, 0xd39eb8fc,
+			0xed545578, 0x08fca5b5, 0xd83d7cd3, 0x4dad0fc4,
+			0x1e50ef5e, 0xb161e6f8, 0xa28514d9, 0x6c51133c,
+			0x6fd5c7e7, 0x56e14ec4, 0x362abfce, 0xddc6c837,
+			0xd79a3234, 0x92638212, 0x670efa8e, 0x406000e0
+		}, {
+			0x3a39ce37, 0xd3faf5cf, 0xabc27737, 0x5ac52d1b,
+			0x5cb0679e, 0x4fa33742, 0xd3822740, 0x99bc9bbe,
+			0xd5118e9d, 0xbf0f7315, 0xd62d1c7e, 0xc700c47b,
+			0xb78c1b6b, 0x21a19045, 0xb26eb1be, 0x6a366eb4,
+			0x5748ab2f, 0xbc946e79, 0xc6a376d2, 0x6549c2c8,
+			0x530ff8ee, 0x468dde7d, 0xd5730a1d, 0x4cd04dc6,
+			0x2939bbdb, 0xa9ba4650, 0xac9526e8, 0xbe5ee304,
+			0xa1fad5f0, 0x6a2d519a, 0x63ef8ce2, 0x9a86ee22,
+			0xc089c2b8, 0x43242ef6, 0xa51e03aa, 0x9cf2d0a4,
+			0x83c061ba, 0x9be96a4d, 0x8fe51550, 0xba645bd6,
+			0x2826a2f9, 0xa73a3ae1, 0x4ba99586, 0xef5562e9,
+			0xc72fefd3, 0xf752f7da, 0x3f046f69, 0x77fa0a59,
+			0x80e4a915, 0x87b08601, 0x9b09e6ad, 0x3b3ee593,
+			0xe990fd5a, 0x9e34d797, 0x2cf0b7d9, 0x022b8b51,
+			0x96d5ac3a, 0x017da67d, 0xd1cf3ed6, 0x7c7d2d28,
+			0x1f9f25cf, 0xadf2b89b, 0x5ad6b472, 0x5a88f54c,
+			0xe029ac71, 0xe019a5e6, 0x47b0acfd, 0xed93fa9b,
+			0xe8d3c48d, 0x283b57cc, 0xf8d56629, 0x79132e28,
+			0x785f0191, 0xed756055, 0xf7960e44, 0xe3d35e8c,
+			0x15056dd4, 0x88f46dba, 0x03a16125, 0x0564f0bd,
+			0xc3eb9e15, 0x3c9057a2, 0x97271aec, 0xa93a072a,
+			0x1b3f6d9b, 0x1e6321f5, 0xf59c66fb, 0x26dcf319,
+			0x7533d928, 0xb155fdf5, 0x03563482, 0x8aba3cbb,
+			0x28517711, 0xc20ad9f8, 0xabcc5167, 0xccad925f,
+			0x4de81751, 0x3830dc8e, 0x379d5862, 0x9320f991,
+			0xea7a90c2, 0xfb3e7bce, 0x5121ce64, 0x774fbe32,
+			0xa8b6e37e, 0xc3293d46, 0x48de5369, 0x6413e680,
+			0xa2ae0810, 0xdd6db224, 0x69852dfd, 0x09072166,
+			0xb39a460a, 0x6445c0dd, 0x586cdecf, 0x1c20c8ae,
+			0x5bbef7dd, 0x1b588d40, 0xccd2017f, 0x6bb4e3bb,
+			0xdda26a7e, 0x3a59ff45, 0x3e350a44, 0xbcb4cdd5,
+			0x72eacea8, 0xfa6484bb, 0x8d6612ae, 0xbf3c6f47,
+			0xd29be463, 0x542f5d9e, 0xaec2771b, 0xf64e6370,
+			0x740e0d8d, 0xe75b1357, 0xf8721671, 0xaf537d5d,
+			0x4040cb08, 0x4eb4e2cc, 0x34d2466a, 0x0115af84,
+			0xe1b00428, 0x95983a1d, 0x06b89fb4, 0xce6ea048,
+			0x6f3f3b82, 0x3520ab82, 0x011a1d4b, 0x277227f8,
+			0x611560b1, 0xe7933fdc, 0xbb3a792b, 0x344525bd,
+			0xa08839e1, 0x51ce794b, 0x2f32c9b7, 0xa01fbac9,
+			0xe01cc87e, 0xbcc7d1f6, 0xcf0111c3, 0xa1e8aac7,
+			0x1a908749, 0xd44fbd9a, 0xd0dadecb, 0xd50ada38,
+			0x0339c32a, 0xc6913667, 0x8df9317c, 0xe0b12b4f,
+			0xf79e59b7, 0x43f5bb3a, 0xf2d519ff, 0x27d9459c,
+			0xbf97222c, 0x15e6fc2a, 0x0f91fc71, 0x9b941525,
+			0xfae59361, 0xceb69ceb, 0xc2a86459, 0x12baa8d1,
+			0xb6c1075e, 0xe3056a0c, 0x10d25065, 0xcb03a442,
+			0xe0ec6e0e, 0x1698db3b, 0x4c98a0be, 0x3278e964,
+			0x9f1f9532, 0xe0d392df, 0xd3a0342b, 0x8971f21e,
+			0x1b0a7441, 0x4ba3348c, 0xc5be7120, 0xc37632d8,
+			0xdf359f8d, 0x9b992f2e, 0xe60b6f47, 0x0fe3f11d,
+			0xe54cda54, 0x1edad891, 0xce6279cf, 0xcd3e7e6f,
+			0x1618b166, 0xfd2c1d05, 0x848fd2c5, 0xf6fb2299,
+			0xf523f357, 0xa6327623, 0x93a83531, 0x56cccd02,
+			0xacf08162, 0x5a75ebb5, 0x6e163697, 0x88d273cc,
+			0xde966292, 0x81b949d0, 0x4c50901b, 0x71c65614,
+			0xe6c6c7bd, 0x327a140a, 0x45e1d006, 0xc3f27b9a,
+			0xc9aa53fd, 0x62a80f00, 0xbb25bfe2, 0x35bdd2f6,
+			0x71126905, 0xb2040222, 0xb6cbcf7c, 0xcd769c2b,
+			0x53113ec0, 0x1640e3d3, 0x38abbd60, 0x2547adf0,
+			0xba38209c, 0xf746ce76, 0x77afa1c5, 0x20756060,
+			0x85cbfe4e, 0x8ae88dd8, 0x7aaaf9b0, 0x4cf9aa7e,
+			0x1948c25c, 0x02fb8a8c, 0x01c36ae4, 0xd6ebe1f9,
+			0x90d4f869, 0xa65cdea0, 0x3f09252d, 0xc208e69f,
+			0xb74e6132, 0xce77e25b, 0x578fdfe3, 0x3ac372e6
+		}
+	}
+}};
+
+static const unsigned char BF_itoa64[64 + 1] =
+	"./ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789";
+
+static const unsigned char BF_atoi64[0x60] = {
+	64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 0, 1,
+	54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 64, 64, 64, 64, 64,
+	64, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
+	17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 64, 64, 64, 64, 64,
+	64, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
+	43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 64, 64, 64, 64, 64
+};
+
+#define BF_safe_atoi64(dst, src) \
+{ \
+	tmp = (unsigned char)(src); \
+	if ((unsigned int)(tmp -= 0x20) >= 0x60) return -1; \
+	tmp = BF_atoi64[tmp]; \
+	if (tmp > 63) return -1; \
+	(dst) = tmp; \
+}
+
+static int BF_decode(BF_word *dst, const char *src, int size)
+{
+	unsigned char *dptr = (unsigned char *)dst;
+	unsigned char *end = dptr + size;
+	const unsigned char *sptr = (const unsigned char *)src;
+	unsigned int tmp, c1, c2, c3, c4;
+
+	do {
+		BF_safe_atoi64(c1, *sptr++);
+		BF_safe_atoi64(c2, *sptr++);
+		*dptr++ = (c1 << 2) | ((c2 & 0x30) >> 4);
+		if (dptr >= end) break;
+
+		BF_safe_atoi64(c3, *sptr++);
+		*dptr++ = ((c2 & 0x0F) << 4) | ((c3 & 0x3C) >> 2);
+		if (dptr >= end) break;
+
+		BF_safe_atoi64(c4, *sptr++);
+		*dptr++ = ((c3 & 0x03) << 6) | c4;
+	} while (dptr < end);
+
+	return 0;
+}
+
+static void BF_encode(char *dst, const BF_word *src, int size)
+{
+	const unsigned char *sptr = (const unsigned char *)src;
+	const unsigned char *end = sptr + size;
+	unsigned char *dptr = (unsigned char *)dst;
+	unsigned int c1, c2;
+
+	do {
+		c1 = *sptr++;
+		*dptr++ = BF_itoa64[c1 >> 2];
+		c1 = (c1 & 0x03) << 4;
+		if (sptr >= end) {
+			*dptr++ = BF_itoa64[c1];
+			break;
+		}
+
+		c2 = *sptr++;
+		c1 |= c2 >> 4;
+		*dptr++ = BF_itoa64[c1];
+		c1 = (c2 & 0x0f) << 2;
+		if (sptr >= end) {
+			*dptr++ = BF_itoa64[c1];
+			break;
+		}
+
+		c2 = *sptr++;
+		c1 |= c2 >> 6;
+		*dptr++ = BF_itoa64[c1];
+		*dptr++ = BF_itoa64[c2 & 0x3f];
+	} while (sptr < end);
+}
+
+static void BF_swap(BF_word *x, int count)
+{
+	if ((union { int i; char c; }){1}.c)
+	do {
+		BF_word tmp = *x;
+		tmp = (tmp << 16) | (tmp >> 16);
+		*x++ = ((tmp & 0x00FF00FF) << 8) | ((tmp >> 8) & 0x00FF00FF);
+	} while (--count);
+}
+
+#define BF_ROUND(L, R, N) \
+	tmp1 = L & 0xFF; \
+	tmp2 = L >> 8; \
+	tmp2 &= 0xFF; \
+	tmp3 = L >> 16; \
+	tmp3 &= 0xFF; \
+	tmp4 = L >> 24; \
+	tmp1 = ctx->s.S[3][tmp1]; \
+	tmp2 = ctx->s.S[2][tmp2]; \
+	tmp3 = ctx->s.S[1][tmp3]; \
+	tmp3 += ctx->s.S[0][tmp4]; \
+	tmp3 ^= tmp2; \
+	R ^= ctx->s.P[N + 1]; \
+	tmp3 += tmp1; \
+	R ^= tmp3;
+
+static BF_word BF_encrypt(BF_ctx *ctx,
+    BF_word L, BF_word R,
+    BF_word *start, BF_word *end)
+{
+	BF_word tmp1, tmp2, tmp3, tmp4;
+	BF_word *ptr = start;
+
+	do {
+		L ^= ctx->s.P[0];
+#if 0
+		BF_ROUND(L, R, 0);
+		BF_ROUND(R, L, 1);
+		BF_ROUND(L, R, 2);
+		BF_ROUND(R, L, 3);
+		BF_ROUND(L, R, 4);
+		BF_ROUND(R, L, 5);
+		BF_ROUND(L, R, 6);
+		BF_ROUND(R, L, 7);
+		BF_ROUND(L, R, 8);
+		BF_ROUND(R, L, 9);
+		BF_ROUND(L, R, 10);
+		BF_ROUND(R, L, 11);
+		BF_ROUND(L, R, 12);
+		BF_ROUND(R, L, 13);
+		BF_ROUND(L, R, 14);
+		BF_ROUND(R, L, 15);
+#else
+		for (int i=0; i<16; i+=2) {
+			BF_ROUND(L, R, i);
+			BF_ROUND(R, L, i+1);
+		}
+#endif
+		tmp4 = R;
+		R = L;
+		L = tmp4 ^ ctx->s.P[BF_N + 1];
+		*ptr++ = L;
+		*ptr++ = R;
+	} while (ptr < end);
+
+	return L;
+}
+
+static void BF_set_key(const char *key, BF_key expanded, BF_key initial,
+    unsigned char flags)
+{
+	const char *ptr = key;
+	unsigned int bug, i, j;
+	BF_word safety, sign, diff, tmp[2];
+
+/*
+ * There was a sign extension bug in older revisions of this function.  While
+ * we would have liked to simply fix the bug and move on, we have to provide
+ * a backwards compatibility feature (essentially the bug) for some systems and
+ * a safety measure for some others.  The latter is needed because for certain
+ * multiple inputs to the buggy algorithm there exist easily found inputs to
+ * the correct algorithm that produce the same hash.  Thus, we optionally
+ * deviate from the correct algorithm just enough to avoid such collisions.
+ * While the bug itself affected the majority of passwords containing
+ * characters with the 8th bit set (although only a percentage of those in a
+ * collision-producing way), the anti-collision safety measure affects
+ * only a subset of passwords containing the '\xff' character (not even all of
+ * those passwords, just some of them).  This character is not found in valid
+ * UTF-8 sequences and is rarely used in popular 8-bit character encodings.
+ * Thus, the safety measure is unlikely to cause much annoyance, and is a
+ * reasonable tradeoff to use when authenticating against existing hashes that
+ * are not reliably known to have been computed with the correct algorithm.
+ *
+ * We use an approach that tries to minimize side-channel leaks of password
+ * information - that is, we mostly use fixed-cost bitwise operations instead
+ * of branches or table lookups.  (One conditional branch based on password
+ * length remains.  It is not part of the bug aftermath, though, and is
+ * difficult and possibly unreasonable to avoid given the use of C strings by
+ * the caller, which results in similar timing leaks anyway.)
+ *
+ * For actual implementation, we set an array index in the variable "bug"
+ * (0 means no bug, 1 means sign extension bug emulation) and a flag in the
+ * variable "safety" (bit 16 is set when the safety measure is requested).
+ * Valid combinations of settings are:
+ *
+ * Prefix "$2a$": bug = 0, safety = 0x10000
+ * Prefix "$2x$": bug = 1, safety = 0
+ * Prefix "$2y$": bug = 0, safety = 0
+ */
+	bug = flags & 1;
+	safety = ((BF_word)flags & 2) << 15;
+
+	sign = diff = 0;
+
+	for (i = 0; i < BF_N + 2; i++) {
+		tmp[0] = tmp[1] = 0;
+		for (j = 0; j < 4; j++) {
+			tmp[0] <<= 8;
+			tmp[0] |= (unsigned char)*ptr; /* correct */
+			tmp[1] <<= 8;
+			tmp[1] |= (signed char)*ptr; /* bug */
+/*
+ * Sign extension in the first char has no effect - nothing to overwrite yet,
+ * and those extra 24 bits will be fully shifted out of the 32-bit word.  For
+ * chars 2, 3, 4 in each four-char block, we set bit 7 of "sign" if sign
+ * extension in tmp[1] occurs.  Once this flag is set, it remains set.
+ */
+			if (j)
+				sign |= tmp[1] & 0x80;
+			if (!*ptr)
+				ptr = key;
+			else
+				ptr++;
+		}
+		diff |= tmp[0] ^ tmp[1]; /* Non-zero on any differences */
+
+		expanded[i] = tmp[bug];
+		initial[i] = BF_init_state.s.P[i] ^ tmp[bug];
+	}
+
+/*
+ * At this point, "diff" is zero iff the correct and buggy algorithms produced
+ * exactly the same result.  If so and if "sign" is non-zero, which indicates
+ * that there was a non-benign sign extension, this means that we have a
+ * collision between the correctly computed hash for this password and a set of
+ * passwords that could be supplied to the buggy algorithm.  Our safety measure
+ * is meant to protect from such many-buggy to one-correct collisions, by
+ * deviating from the correct algorithm in such cases.  Let's check for this.
+ */
+	diff |= diff >> 16; /* still zero iff exact match */
+	diff &= 0xffff; /* ditto */
+	diff += 0xffff; /* bit 16 set iff "diff" was non-zero (on non-match) */
+	sign <<= 9; /* move the non-benign sign extension flag to bit 16 */
+	sign &= ~diff & safety; /* action needed? */
+
+/*
+ * If we have determined that we need to deviate from the correct algorithm,
+ * flip bit 16 in initial expanded key.  (The choice of 16 is arbitrary, but
+ * let's stick to it now.  It came out of the approach we used above, and it's
+ * not any worse than any other choice we could make.)
+ *
+ * It is crucial that we don't do the same to the expanded key used in the main
+ * Eksblowfish loop.  By doing it to only one of these two, we deviate from a
+ * state that could be directly specified by a password to the buggy algorithm
+ * (and to the fully correct one as well, but that's a side-effect).
+ */
+	initial[0] ^= sign;
+}
+
+static char *BF_crypt(const char *key, const char *setting,
+	char *output, BF_word min)
+{
+	static const unsigned char flags_by_subtype[26] =
+		{2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 0};
+	struct {
+		BF_ctx ctx;
+		BF_key expanded_key;
+		union {
+			BF_word salt[4];
+			BF_word output[6];
+		} binary;
+	} data;
+	BF_word count;
+	int i;
+
+	if (setting[0] != '$' ||
+	    setting[1] != '2' ||
+	    setting[2] - 'a' > 25U ||
+	    !flags_by_subtype[setting[2] - 'a'] ||
+	    setting[3] != '$' ||
+	    setting[4] - '0' > 1U ||
+	    setting[5] - '0' > 9U ||
+	    setting[6] != '$') {
+		return NULL;
+	}
+
+	count = (BF_word)1 << ((setting[4] - '0') * 10 + (setting[5] - '0'));
+	if (count < min || BF_decode(data.binary.salt, &setting[7], 16)) {
+		return NULL;
+	}
+	BF_swap(data.binary.salt, 4);
+
+	BF_set_key(key, data.expanded_key, data.ctx.s.P,
+	    flags_by_subtype[setting[2] - 'a']);
+
+	memcpy(data.ctx.s.S, BF_init_state.s.S, sizeof(data.ctx.s.S));
+
+	{
+		BF_word L = 0, R = 0;
+		BF_word *ptr = &data.ctx.PS[0];
+		do {
+			L = BF_encrypt(&data.ctx,
+			    L ^ data.binary.salt[0], R ^ data.binary.salt[1],
+			    ptr, ptr);
+			R = *(ptr + 1);
+			ptr += 2;
+
+			if (ptr >= &data.ctx.PS[BF_N + 2 + 4 * 0x100])
+				break;
+
+			L = BF_encrypt(&data.ctx,
+			    L ^ data.binary.salt[2], R ^ data.binary.salt[3],
+			    ptr, ptr);
+			R = *(ptr + 1);
+			ptr += 2;
+		} while (1);
+	}
+
+	do {
+		int done;
+
+		for (i = 0; i < BF_N + 2; i += 2) {
+			data.ctx.s.P[i] ^= data.expanded_key[i];
+			data.ctx.s.P[i + 1] ^= data.expanded_key[i + 1];
+		}
+
+		done = 0;
+		do {
+			BF_encrypt(&data.ctx, 0, 0,
+			    &data.ctx.PS[0],
+			    &data.ctx.PS[BF_N + 2 + 4 * 0x100]);
+
+			if (done)
+				break;
+			done = 1;
+
+			{
+				BF_word tmp1, tmp2, tmp3, tmp4;
+
+				tmp1 = data.binary.salt[0];
+				tmp2 = data.binary.salt[1];
+				tmp3 = data.binary.salt[2];
+				tmp4 = data.binary.salt[3];
+				for (i = 0; i < BF_N; i += 4) {
+					data.ctx.s.P[i] ^= tmp1;
+					data.ctx.s.P[i + 1] ^= tmp2;
+					data.ctx.s.P[i + 2] ^= tmp3;
+					data.ctx.s.P[i + 3] ^= tmp4;
+				}
+				data.ctx.s.P[16] ^= tmp1;
+				data.ctx.s.P[17] ^= tmp2;
+			}
+		} while (1);
+	} while (--count);
+
+	for (i = 0; i < 6; i += 2) {
+		BF_word L, LR[2];
+
+		L = BF_magic_w[i];
+		LR[1] = BF_magic_w[i + 1];
+
+		count = 64;
+		do {
+			L = BF_encrypt(&data.ctx, L, LR[1],
+			    &LR[0], &LR[0]);
+		} while (--count);
+
+		data.binary.output[i] = L;
+		data.binary.output[i + 1] = LR[1];
+	}
+
+	memcpy(output, setting, 7 + 22 - 1);
+	output[7 + 22 - 1] = BF_itoa64[
+		BF_atoi64[setting[7 + 22 - 1] - 0x20] & 0x30];
+
+/* This has to be bug-compatible with the original implementation, so
+ * only encode 23 of the 24 bytes. :-) */
+	BF_swap(data.binary.output, 6);
+	BF_encode(&output[7 + 22], data.binary.output, 23);
+	output[7 + 22 + 31] = '\0';
+
+	return output;
+}
+
+/*
+ * Please preserve the runtime self-test.  It serves two purposes at once:
+ *
+ * 1. We really can't afford the risk of producing incompatible hashes e.g.
+ * when there's something like gcc bug 26587 again, whereas an application or
+ * library integrating this code might not also integrate our external tests or
+ * it might not run them after every build.  Even if it does, the miscompile
+ * might only occur on the production build, but not on a testing build (such
+ * as because of different optimization settings).  It is painful to recover
+ * from incorrectly-computed hashes - merely fixing whatever broke is not
+ * enough.  Thus, a proactive measure like this self-test is needed.
+ *
+ * 2. We don't want to leave sensitive data from our actual password hash
+ * computation on the stack or in registers.  Previous revisions of the code
+ * would do explicit cleanups, but simply running the self-test after hash
+ * computation is more reliable.
+ *
+ * The performance cost of this quick self-test is around 0.6% at the "$2a$08"
+ * setting.
+ */
+char *__crypt_blowfish(const char *key, const char *setting, char *output)
+{
+	const char *test_key = "8b \xd0\xc1\xd2\xcf\xcc\xd8";
+	const char *test_setting = "$2a$00$abcdefghijklmnopqrstuu";
+	static const char test_hash[2][34] =
+		{"VUrPmXD6q/nVSSp7pNDhCR9071IfIRe\0\x55", /* $2x$ */
+		"i1D709vfamulimlGcq0qq3UvuUasvEa\0\x55"}; /* $2a$, $2y$ */
+	char *retval;
+	const char *p;
+	int ok;
+	struct {
+		char s[7 + 22 + 1];
+		char o[7 + 22 + 31 + 1 + 1 + 1];
+	} buf;
+
+/* Hash the supplied password */
+	retval = BF_crypt(key, setting, output, 16);
+
+/*
+ * Do a quick self-test.  It is important that we make both calls to BF_crypt()
+ * from the same scope such that they likely use the same stack locations,
+ * which makes the second call overwrite the first call's sensitive data on the
+ * stack and makes it more likely that any alignment related issues would be
+ * detected by the self-test.
+ */
+	memcpy(buf.s, test_setting, sizeof(buf.s));
+	if (retval)
+		buf.s[2] = setting[2];
+	memset(buf.o, 0x55, sizeof(buf.o));
+	buf.o[sizeof(buf.o) - 1] = 0;
+	p = BF_crypt(test_key, buf.s, buf.o, 1);
+
+	ok = (p == buf.o &&
+	    !memcmp(p, buf.s, 7 + 22) &&
+	    !memcmp(p + (7 + 22),
+	    test_hash[buf.s[2] & 1],
+	    31 + 1 + 1 + 1));
+
+	{
+		const char *k = "\xff\xa3" "34" "\xff\xff\xff\xa3" "345";
+		BF_key ae, ai, ye, yi;
+		BF_set_key(k, ae, ai, 2); /* $2a$ */
+		BF_set_key(k, ye, yi, 4); /* $2y$ */
+		ai[0] ^= 0x10000; /* undo the safety (for comparison) */
+		ok = ok && ai[0] == 0xdb9c59bc && ye[17] == 0x33343500 &&
+		    !memcmp(ae, ye, sizeof(ae)) &&
+		    !memcmp(ai, yi, sizeof(ai));
+	}
+
+	if (ok && retval)
+		return retval;
+
+	return "*";
+}

libc/musl/src/crypt/crypt_des.h

@@ -0,0 +1,14 @@
+#ifndef CRYPT_DES_H
+#define CRYPT_DES_H
+
+#include <stdint.h>
+
+struct expanded_key {
+	uint32_t l[16], r[16];
+};
+
+hidden void __des_setkey(const unsigned char *, struct expanded_key *);
+hidden void __do_des(uint32_t, uint32_t, uint32_t *, uint32_t *,
+                     uint32_t, uint32_t, const struct expanded_key *);
+
+#endif

libc/musl/src/crypt/crypt_md5.c

@@ -0,0 +1,285 @@
+/*
+ * md5 crypt implementation
+ *
+ * original md5 crypt design is from Poul-Henning Kamp
+ * this implementation was created based on the code in freebsd
+ * at least 32bit int is assumed, key is limited and $1$ prefix is mandatory,
+ * on error "*" is returned
+ */
+#include <string.h>
+#include <stdint.h>
+
+/* public domain md5 implementation based on rfc1321 and libtomcrypt */
+
+struct md5 {
+	uint64_t len;    /* processed message length */
+	uint32_t h[4];   /* hash state */
+	uint8_t buf[64]; /* message block buffer */
+};
+
+static uint32_t rol(uint32_t n, int k) { return (n << k) | (n >> (32-k)); }
+#define F(x,y,z) (z ^ (x & (y ^ z)))
+#define G(x,y,z) (y ^ (z & (y ^ x)))
+#define H(x,y,z) (x ^ y ^ z)
+#define I(x,y,z) (y ^ (x | ~z))
+#define FF(a,b,c,d,w,s,t) a += F(b,c,d) + w + t; a = rol(a,s) + b
+#define GG(a,b,c,d,w,s,t) a += G(b,c,d) + w + t; a = rol(a,s) + b
+#define HH(a,b,c,d,w,s,t) a += H(b,c,d) + w + t; a = rol(a,s) + b
+#define II(a,b,c,d,w,s,t) a += I(b,c,d) + w + t; a = rol(a,s) + b
+
+static const uint32_t tab[64] = {
+0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, 0xa8304613, 0xfd469501,
+0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e, 0x49b40821,
+0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x02441453, 0xd8a1e681, 0xe7d3fbc8,
+0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8, 0x676f02d9, 0x8d2a4c8a,
+0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60, 0xbebfbc70,
+0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x04881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8, 0xc4ac5665,
+0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, 0xffeff47d, 0x85845dd1,
+0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1, 0xf7537e82, 0xbd3af235, 0x2ad7d2bb, 0xeb86d391
+};
+
+static void processblock(struct md5 *s, const uint8_t *buf)
+{
+	uint32_t i, W[16], a, b, c, d;
+
+	for (i = 0; i < 16; i++) {
+		W[i] = buf[4*i];
+		W[i] |= (uint32_t)buf[4*i+1]<<8;
+		W[i] |= (uint32_t)buf[4*i+2]<<16;
+		W[i] |= (uint32_t)buf[4*i+3]<<24;
+	}
+
+	a = s->h[0];
+	b = s->h[1];
+	c = s->h[2];
+	d = s->h[3];
+
+	i = 0;
+	while (i < 16) {
+		FF(a,b,c,d, W[i],  7, tab[i]); i++;
+		FF(d,a,b,c, W[i], 12, tab[i]); i++;
+		FF(c,d,a,b, W[i], 17, tab[i]); i++;
+		FF(b,c,d,a, W[i], 22, tab[i]); i++;
+	}
+	while (i < 32) {
+		GG(a,b,c,d, W[(5*i+1)%16],  5, tab[i]); i++;
+		GG(d,a,b,c, W[(5*i+1)%16],  9, tab[i]); i++;
+		GG(c,d,a,b, W[(5*i+1)%16], 14, tab[i]); i++;
+		GG(b,c,d,a, W[(5*i+1)%16], 20, tab[i]); i++;
+	}
+	while (i < 48) {
+		HH(a,b,c,d, W[(3*i+5)%16],  4, tab[i]); i++;
+		HH(d,a,b,c, W[(3*i+5)%16], 11, tab[i]); i++;
+		HH(c,d,a,b, W[(3*i+5)%16], 16, tab[i]); i++;
+		HH(b,c,d,a, W[(3*i+5)%16], 23, tab[i]); i++;
+	}
+	while (i < 64) {
+		II(a,b,c,d, W[7*i%16],  6, tab[i]); i++;
+		II(d,a,b,c, W[7*i%16], 10, tab[i]); i++;
+		II(c,d,a,b, W[7*i%16], 15, tab[i]); i++;
+		II(b,c,d,a, W[7*i%16], 21, tab[i]); i++;
+	}
+
+	s->h[0] += a;
+	s->h[1] += b;
+	s->h[2] += c;
+	s->h[3] += d;
+}
+
+static void pad(struct md5 *s)
+{
+	unsigned r = s->len % 64;
+
+	s->buf[r++] = 0x80;
+	if (r > 56) {
+		memset(s->buf + r, 0, 64 - r);
+		r = 0;
+		processblock(s, s->buf);
+	}
+	memset(s->buf + r, 0, 56 - r);
+	s->len *= 8;
+	s->buf[56] = s->len;
+	s->buf[57] = s->len >> 8;
+	s->buf[58] = s->len >> 16;
+	s->buf[59] = s->len >> 24;
+	s->buf[60] = s->len >> 32;
+	s->buf[61] = s->len >> 40;
+	s->buf[62] = s->len >> 48;
+	s->buf[63] = s->len >> 56;
+	processblock(s, s->buf);
+}
+
+static void md5_init(struct md5 *s)
+{
+	s->len = 0;
+	s->h[0] = 0x67452301;
+	s->h[1] = 0xefcdab89;
+	s->h[2] = 0x98badcfe;
+	s->h[3] = 0x10325476;
+}
+
+static void md5_sum(struct md5 *s, uint8_t *md)
+{
+	int i;
+
+	pad(s);
+	for (i = 0; i < 4; i++) {
+		md[4*i] = s->h[i];
+		md[4*i+1] = s->h[i] >> 8;
+		md[4*i+2] = s->h[i] >> 16;
+		md[4*i+3] = s->h[i] >> 24;
+	}
+}
+
+static void md5_update(struct md5 *s, const void *m, unsigned long len)
+{
+	const uint8_t *p = m;
+	unsigned r = s->len % 64;
+
+	s->len += len;
+	if (r) {
+		if (len < 64 - r) {
+			memcpy(s->buf + r, p, len);
+			return;
+		}
+		memcpy(s->buf + r, p, 64 - r);
+		len -= 64 - r;
+		p += 64 - r;
+		processblock(s, s->buf);
+	}
+	for (; len >= 64; len -= 64, p += 64)
+		processblock(s, p);
+	memcpy(s->buf, p, len);
+}
+
+/*-
+ * Copyright (c) 2003 Poul-Henning Kamp
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+/* key limit is not part of the original design, added for DoS protection */
+#define KEY_MAX 30000
+#define SALT_MAX 8
+
+static const unsigned char b64[] =
+"./0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
+
+static char *to64(char *s, unsigned int u, int n)
+{
+	while (--n >= 0) {
+		*s++ = b64[u % 64];
+		u /= 64;
+	}
+	return s;
+}
+
+static char *md5crypt(const char *key, const char *setting, char *output)
+{
+	struct md5 ctx;
+	unsigned char md[16];
+	unsigned int i, klen, slen;
+	const char *salt;
+	char *p;
+
+	/* reject large keys */
+	klen = strnlen(key, KEY_MAX+1);
+	if (klen > KEY_MAX)
+		return 0;
+
+	/* setting: $1$salt$ (closing $ is optional) */
+	if (strncmp(setting, "$1$", 3) != 0)
+		return 0;
+	salt = setting + 3;
+	for (i = 0; i < SALT_MAX && salt[i] && salt[i] != '$'; i++);
+	slen = i;
+
+	/* md5(key salt key) */
+	md5_init(&ctx);
+	md5_update(&ctx, key, klen);
+	md5_update(&ctx, salt, slen);
+	md5_update(&ctx, key, klen);
+	md5_sum(&ctx, md);
+
+	/* md5(key $1$ salt repeated-md weird-key[0]-0) */
+	md5_init(&ctx);
+	md5_update(&ctx, key, klen);
+	md5_update(&ctx, setting, 3 + slen);
+	for (i = klen; i > sizeof md; i -= sizeof md)
+		md5_update(&ctx, md, sizeof md);
+	md5_update(&ctx, md, i);
+	md[0] = 0;
+	for (i = klen; i; i >>= 1)
+		if (i & 1)
+			md5_update(&ctx, md, 1);
+		else
+			md5_update(&ctx, key, 1);
+	md5_sum(&ctx, md);
+
+	/* md = f(md, key, salt) iteration */
+	for (i = 0; i < 1000; i++) {
+		md5_init(&ctx);
+		if (i % 2)
+			md5_update(&ctx, key, klen);
+		else
+			md5_update(&ctx, md, sizeof md);
+		if (i % 3)
+			md5_update(&ctx, salt, slen);
+		if (i % 7)
+			md5_update(&ctx, key, klen);
+		if (i % 2)
+			md5_update(&ctx, md, sizeof md);
+		else
+			md5_update(&ctx, key, klen);
+		md5_sum(&ctx, md);
+	}
+
+	/* output is $1$salt$hash */
+	memcpy(output, setting, 3 + slen);
+	p = output + 3 + slen;
+	*p++ = '$';
+	static const unsigned char perm[][3] = {
+		0,6,12,1,7,13,2,8,14,3,9,15,4,10,5 };
+	for (i=0; i<5; i++) p = to64(p,
+		(md[perm[i][0]]<<16)|(md[perm[i][1]]<<8)|md[perm[i][2]], 4);
+	p = to64(p, md[11], 2);
+	*p = 0;
+
+	return output;
+}
+
+char *__crypt_md5(const char *key, const char *setting, char *output)
+{
+	static const char testkey[] = "Xy01@#\x01\x02\x80\x7f\xff\r\n\x81\t !";
+	static const char testsetting[] = "$1$abcd0123$";
+	static const char testhash[] = "$1$abcd0123$9Qcg8DyviekV3tDGMZynJ1";
+	char testbuf[64];
+	char *p, *q;
+
+	p = md5crypt(key, setting, output);
+	/* self test and stack cleanup */
+	q = md5crypt(testkey, testsetting, testbuf);
+	if (!p || q != testbuf || memcmp(testbuf, testhash, sizeof testhash))
+		return "*";
+	return p;
+}

libc/musl/src/crypt/crypt_r.c

@@ -0,0 +1,23 @@
+#include <crypt.h>
+
+char *__crypt_r(const char *key, const char *salt, struct crypt_data *data)
+{
+	/* Per the crypt_r API, the caller has provided a pointer to
+	 * struct crypt_data; however, this implementation does not
+	 * use the structure to store any internal state, and treats
+	 * it purely as a char buffer for storing the result. */
+	char *output = (char *)data;
+	if (salt[0] == '$' && salt[1] && salt[2]) {
+		if (salt[1] == '1' && salt[2] == '$')
+			return __crypt_md5(key, salt, output);
+		if (salt[1] == '2' && salt[3] == '$')
+			return __crypt_blowfish(key, salt, output);
+		if (salt[1] == '5' && salt[2] == '$')
+			return __crypt_sha256(key, salt, output);
+		if (salt[1] == '6' && salt[2] == '$')
+			return __crypt_sha512(key, salt, output);
+	}
+	return __crypt_des(key, salt, output);
+}
+
+weak_alias(__crypt_r, crypt_r);

libc/musl/src/crypt/crypt_sha256.c

@@ -0,0 +1,322 @@
+/*
+ * public domain sha256 crypt implementation
+ *
+ * original sha crypt design: http://people.redhat.com/drepper/SHA-crypt.txt
+ * in this implementation at least 32bit int is assumed,
+ * key length is limited, the $5$ prefix is mandatory, '\n' and ':' is rejected
+ * in the salt and rounds= setting must contain a valid iteration count,
+ * on error "*" is returned.
+ */
+#include <ctype.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>
+
+/* public domain sha256 implementation based on fips180-3 */
+
+struct sha256 {
+	uint64_t len;    /* processed message length */
+	uint32_t h[8];   /* hash state */
+	uint8_t buf[64]; /* message block buffer */
+};
+
+static uint32_t ror(uint32_t n, int k) { return (n >> k) | (n << (32-k)); }
+#define Ch(x,y,z)  (z ^ (x & (y ^ z)))
+#define Maj(x,y,z) ((x & y) | (z & (x | y)))
+#define S0(x)      (ror(x,2) ^ ror(x,13) ^ ror(x,22))
+#define S1(x)      (ror(x,6) ^ ror(x,11) ^ ror(x,25))
+#define R0(x)      (ror(x,7) ^ ror(x,18) ^ (x>>3))
+#define R1(x)      (ror(x,17) ^ ror(x,19) ^ (x>>10))
+
+static const uint32_t K[64] = {
+0x428a2f98, 0x71374491, 0xb5c0fbcf, 0xe9b5dba5, 0x3956c25b, 0x59f111f1, 0x923f82a4, 0xab1c5ed5,
+0xd807aa98, 0x12835b01, 0x243185be, 0x550c7dc3, 0x72be5d74, 0x80deb1fe, 0x9bdc06a7, 0xc19bf174,
+0xe49b69c1, 0xefbe4786, 0x0fc19dc6, 0x240ca1cc, 0x2de92c6f, 0x4a7484aa, 0x5cb0a9dc, 0x76f988da,
+0x983e5152, 0xa831c66d, 0xb00327c8, 0xbf597fc7, 0xc6e00bf3, 0xd5a79147, 0x06ca6351, 0x14292967,
+0x27b70a85, 0x2e1b2138, 0x4d2c6dfc, 0x53380d13, 0x650a7354, 0x766a0abb, 0x81c2c92e, 0x92722c85,
+0xa2bfe8a1, 0xa81a664b, 0xc24b8b70, 0xc76c51a3, 0xd192e819, 0xd6990624, 0xf40e3585, 0x106aa070,
+0x19a4c116, 0x1e376c08, 0x2748774c, 0x34b0bcb5, 0x391c0cb3, 0x4ed8aa4a, 0x5b9cca4f, 0x682e6ff3,
+0x748f82ee, 0x78a5636f, 0x84c87814, 0x8cc70208, 0x90befffa, 0xa4506ceb, 0xbef9a3f7, 0xc67178f2
+};
+
+static void processblock(struct sha256 *s, const uint8_t *buf)
+{
+	uint32_t W[64], t1, t2, a, b, c, d, e, f, g, h;
+	int i;
+
+	for (i = 0; i < 16; i++) {
+		W[i] = (uint32_t)buf[4*i]<<24;
+		W[i] |= (uint32_t)buf[4*i+1]<<16;
+		W[i] |= (uint32_t)buf[4*i+2]<<8;
+		W[i] |= buf[4*i+3];
+	}
+	for (; i < 64; i++)
+		W[i] = R1(W[i-2]) + W[i-7] + R0(W[i-15]) + W[i-16];
+	a = s->h[0];
+	b = s->h[1];
+	c = s->h[2];
+	d = s->h[3];
+	e = s->h[4];
+	f = s->h[5];
+	g = s->h[6];
+	h = s->h[7];
+	for (i = 0; i < 64; i++) {
+		t1 = h + S1(e) + Ch(e,f,g) + K[i] + W[i];
+		t2 = S0(a) + Maj(a,b,c);
+		h = g;
+		g = f;
+		f = e;
+		e = d + t1;
+		d = c;
+		c = b;
+		b = a;
+		a = t1 + t2;
+	}
+	s->h[0] += a;
+	s->h[1] += b;
+	s->h[2] += c;
+	s->h[3] += d;
+	s->h[4] += e;
+	s->h[5] += f;
+	s->h[6] += g;
+	s->h[7] += h;
+}
+
+static void pad(struct sha256 *s)
+{
+	unsigned r = s->len % 64;
+
+	s->buf[r++] = 0x80;
+	if (r > 56) {
+		memset(s->buf + r, 0, 64 - r);
+		r = 0;
+		processblock(s, s->buf);
+	}
+	memset(s->buf + r, 0, 56 - r);
+	s->len *= 8;
+	s->buf[56] = s->len >> 56;
+	s->buf[57] = s->len >> 48;
+	s->buf[58] = s->len >> 40;
+	s->buf[59] = s->len >> 32;
+	s->buf[60] = s->len >> 24;
+	s->buf[61] = s->len >> 16;
+	s->buf[62] = s->len >> 8;
+	s->buf[63] = s->len;
+	processblock(s, s->buf);
+}
+
+static void sha256_init(struct sha256 *s)
+{
+	s->len = 0;
+	s->h[0] = 0x6a09e667;
+	s->h[1] = 0xbb67ae85;
+	s->h[2] = 0x3c6ef372;
+	s->h[3] = 0xa54ff53a;
+	s->h[4] = 0x510e527f;
+	s->h[5] = 0x9b05688c;
+	s->h[6] = 0x1f83d9ab;
+	s->h[7] = 0x5be0cd19;
+}
+
+static void sha256_sum(struct sha256 *s, uint8_t *md)
+{
+	int i;
+
+	pad(s);
+	for (i = 0; i < 8; i++) {
+		md[4*i] = s->h[i] >> 24;
+		md[4*i+1] = s->h[i] >> 16;
+		md[4*i+2] = s->h[i] >> 8;
+		md[4*i+3] = s->h[i];
+	}
+}
+
+static void sha256_update(struct sha256 *s, const void *m, unsigned long len)
+{
+	const uint8_t *p = m;
+	unsigned r = s->len % 64;
+
+	s->len += len;
+	if (r) {
+		if (len < 64 - r) {
+			memcpy(s->buf + r, p, len);
+			return;
+		}
+		memcpy(s->buf + r, p, 64 - r);
+		len -= 64 - r;
+		p += 64 - r;
+		processblock(s, s->buf);
+	}
+	for (; len >= 64; len -= 64, p += 64)
+		processblock(s, p);
+	memcpy(s->buf, p, len);
+}
+
+static const unsigned char b64[] =
+"./0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
+
+static char *to64(char *s, unsigned int u, int n)
+{
+	while (--n >= 0) {
+		*s++ = b64[u % 64];
+		u /= 64;
+	}
+	return s;
+}
+
+/* key limit is not part of the original design, added for DoS protection.
+ * rounds limit has been lowered (versus the reference/spec), also for DoS
+ * protection. runtime is O(klen^2 + klen*rounds) */
+#define KEY_MAX 256
+#define SALT_MAX 16
+#define ROUNDS_DEFAULT 5000
+#define ROUNDS_MIN 1000
+#define ROUNDS_MAX 9999999
+
+/* hash n bytes of the repeated md message digest */
+static void hashmd(struct sha256 *s, unsigned int n, const void *md)
+{
+	unsigned int i;
+
+	for (i = n; i > 32; i -= 32)
+		sha256_update(s, md, 32);
+	sha256_update(s, md, i);
+}
+
+static char *sha256crypt(const char *key, const char *setting, char *output)
+{
+	struct sha256 ctx;
+	unsigned char md[32], kmd[32], smd[32];
+	unsigned int i, r, klen, slen;
+	char rounds[20] = "";
+	const char *salt;
+	char *p;
+
+	/* reject large keys */
+	klen = strnlen(key, KEY_MAX+1);
+	if (klen > KEY_MAX)
+		return 0;
+
+	/* setting: $5$rounds=n$salt$ (rounds=n$ and closing $ are optional) */
+	if (strncmp(setting, "$5$", 3) != 0)
+		return 0;
+	salt = setting + 3;
+
+	r = ROUNDS_DEFAULT;
+	if (strncmp(salt, "rounds=", sizeof "rounds=" - 1) == 0) {
+		unsigned long u;
+		char *end;
+
+		/*
+		 * this is a deviation from the reference:
+		 * bad rounds setting is rejected if it is
+		 * - empty
+		 * - unterminated (missing '$')
+		 * - begins with anything but a decimal digit
+		 * the reference implementation treats these bad
+		 * rounds as part of the salt or parse them with
+		 * strtoul semantics which may cause problems
+		 * including non-portable hashes that depend on
+		 * the host's value of ULONG_MAX.
+		 */
+		salt += sizeof "rounds=" - 1;
+		if (!isdigit(*salt))
+			return 0;
+		u = strtoul(salt, &end, 10);
+		if (*end != '$')
+			return 0;
+		salt = end+1;
+		if (u < ROUNDS_MIN)
+			r = ROUNDS_MIN;
+		else if (u > ROUNDS_MAX)
+			return 0;
+		else
+			r = u;
+		/* needed when rounds is zero prefixed or out of bounds */
+		sprintf(rounds, "rounds=%u$", r);
+	}
+
+	for (i = 0; i < SALT_MAX && salt[i] && salt[i] != '$'; i++)
+		/* reject characters that interfere with /etc/shadow parsing */
+		if (salt[i] == '\n' || salt[i] == ':')
+			return 0;
+	slen = i;
+
+	/* B = sha(key salt key) */
+	sha256_init(&ctx);
+	sha256_update(&ctx, key, klen);
+	sha256_update(&ctx, salt, slen);
+	sha256_update(&ctx, key, klen);
+	sha256_sum(&ctx, md);
+
+	/* A = sha(key salt repeat-B alternate-B-key) */
+	sha256_init(&ctx);
+	sha256_update(&ctx, key, klen);
+	sha256_update(&ctx, salt, slen);
+	hashmd(&ctx, klen, md);
+	for (i = klen; i > 0; i >>= 1)
+		if (i & 1)
+			sha256_update(&ctx, md, sizeof md);
+		else
+			sha256_update(&ctx, key, klen);
+	sha256_sum(&ctx, md);
+
+	/* DP = sha(repeat-key), this step takes O(klen^2) time */
+	sha256_init(&ctx);
+	for (i = 0; i < klen; i++)
+		sha256_update(&ctx, key, klen);
+	sha256_sum(&ctx, kmd);
+
+	/* DS = sha(repeat-salt) */
+	sha256_init(&ctx);
+	for (i = 0; i < 16 + md[0]; i++)
+		sha256_update(&ctx, salt, slen);
+	sha256_sum(&ctx, smd);
+
+	/* iterate A = f(A,DP,DS), this step takes O(rounds*klen) time */
+	for (i = 0; i < r; i++) {
+		sha256_init(&ctx);
+		if (i % 2)
+			hashmd(&ctx, klen, kmd);
+		else
+			sha256_update(&ctx, md, sizeof md);
+		if (i % 3)
+			sha256_update(&ctx, smd, slen);
+		if (i % 7)
+			hashmd(&ctx, klen, kmd);
+		if (i % 2)
+			sha256_update(&ctx, md, sizeof md);
+		else
+			hashmd(&ctx, klen, kmd);
+		sha256_sum(&ctx, md);
+	}
+
+	/* output is $5$rounds=n$salt$hash */
+	p = output;
+	p += sprintf(p, "$5$%s%.*s$", rounds, slen, salt);
+	static const unsigned char perm[][3] = {
+		0,10,20,21,1,11,12,22,2,3,13,23,24,4,14,
+		15,25,5,6,16,26,27,7,17,18,28,8,9,19,29 };
+	for (i=0; i<10; i++) p = to64(p,
+		(md[perm[i][0]]<<16)|(md[perm[i][1]]<<8)|md[perm[i][2]], 4);
+	p = to64(p, (md[31]<<8)|md[30], 3);
+	*p = 0;
+	return output;
+}
+
+char *__crypt_sha256(const char *key, const char *setting, char *output)
+{
+	static const char testkey[] = "Xy01@#\x01\x02\x80\x7f\xff\r\n\x81\t !";
+	static const char testsetting[] = "$5$rounds=1234$abc0123456789$";
+	static const char testhash[] = "$5$rounds=1234$abc0123456789$3VfDjPt05VHFn47C/ojFZ6KRPYrOjj1lLbH.dkF3bZ6";
+	char testbuf[128];
+	char *p, *q;
+
+	p = sha256crypt(key, setting, output);
+	/* self test and stack cleanup */
+	q = sha256crypt(testkey, testsetting, testbuf);
+	if (!p || q != testbuf || memcmp(testbuf, testhash, sizeof testhash))
+		return "*";
+	return p;
+}

libc/musl/src/crypt/crypt_sha512.c

@@ -0,0 +1,371 @@
+/*
+ * public domain sha512 crypt implementation
+ *
+ * original sha crypt design: http://people.redhat.com/drepper/SHA-crypt.txt
+ * in this implementation at least 32bit int is assumed,
+ * key length is limited, the $6$ prefix is mandatory, '\n' and ':' is rejected
+ * in the salt and rounds= setting must contain a valid iteration count,
+ * on error "*" is returned.
+ */
+#include <ctype.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>
+
+/* public domain sha512 implementation based on fips180-3 */
+/* >=2^64 bits messages are not supported (about 2000 peta bytes) */
+
+struct sha512 {
+	uint64_t len;     /* processed message length */
+	uint64_t h[8];    /* hash state */
+	uint8_t buf[128]; /* message block buffer */
+};
+
+static uint64_t ror(uint64_t n, int k) { return (n >> k) | (n << (64-k)); }
+#define Ch(x,y,z)  (z ^ (x & (y ^ z)))
+#define Maj(x,y,z) ((x & y) | (z & (x | y)))
+#define S0(x)      (ror(x,28) ^ ror(x,34) ^ ror(x,39))
+#define S1(x)      (ror(x,14) ^ ror(x,18) ^ ror(x,41))
+#define R0(x)      (ror(x,1) ^ ror(x,8) ^ (x>>7))
+#define R1(x)      (ror(x,19) ^ ror(x,61) ^ (x>>6))
+
+static const uint64_t K[80] = {
+0x428a2f98d728ae22ULL, 0x7137449123ef65cdULL, 0xb5c0fbcfec4d3b2fULL, 0xe9b5dba58189dbbcULL,
+0x3956c25bf348b538ULL, 0x59f111f1b605d019ULL, 0x923f82a4af194f9bULL, 0xab1c5ed5da6d8118ULL,
+0xd807aa98a3030242ULL, 0x12835b0145706fbeULL, 0x243185be4ee4b28cULL, 0x550c7dc3d5ffb4e2ULL,
+0x72be5d74f27b896fULL, 0x80deb1fe3b1696b1ULL, 0x9bdc06a725c71235ULL, 0xc19bf174cf692694ULL,
+0xe49b69c19ef14ad2ULL, 0xefbe4786384f25e3ULL, 0x0fc19dc68b8cd5b5ULL, 0x240ca1cc77ac9c65ULL,
+0x2de92c6f592b0275ULL, 0x4a7484aa6ea6e483ULL, 0x5cb0a9dcbd41fbd4ULL, 0x76f988da831153b5ULL,
+0x983e5152ee66dfabULL, 0xa831c66d2db43210ULL, 0xb00327c898fb213fULL, 0xbf597fc7beef0ee4ULL,
+0xc6e00bf33da88fc2ULL, 0xd5a79147930aa725ULL, 0x06ca6351e003826fULL, 0x142929670a0e6e70ULL,
+0x27b70a8546d22ffcULL, 0x2e1b21385c26c926ULL, 0x4d2c6dfc5ac42aedULL, 0x53380d139d95b3dfULL,
+0x650a73548baf63deULL, 0x766a0abb3c77b2a8ULL, 0x81c2c92e47edaee6ULL, 0x92722c851482353bULL,
+0xa2bfe8a14cf10364ULL, 0xa81a664bbc423001ULL, 0xc24b8b70d0f89791ULL, 0xc76c51a30654be30ULL,
+0xd192e819d6ef5218ULL, 0xd69906245565a910ULL, 0xf40e35855771202aULL, 0x106aa07032bbd1b8ULL,
+0x19a4c116b8d2d0c8ULL, 0x1e376c085141ab53ULL, 0x2748774cdf8eeb99ULL, 0x34b0bcb5e19b48a8ULL,
+0x391c0cb3c5c95a63ULL, 0x4ed8aa4ae3418acbULL, 0x5b9cca4f7763e373ULL, 0x682e6ff3d6b2b8a3ULL,
+0x748f82ee5defb2fcULL, 0x78a5636f43172f60ULL, 0x84c87814a1f0ab72ULL, 0x8cc702081a6439ecULL,
+0x90befffa23631e28ULL, 0xa4506cebde82bde9ULL, 0xbef9a3f7b2c67915ULL, 0xc67178f2e372532bULL,
+0xca273eceea26619cULL, 0xd186b8c721c0c207ULL, 0xeada7dd6cde0eb1eULL, 0xf57d4f7fee6ed178ULL,
+0x06f067aa72176fbaULL, 0x0a637dc5a2c898a6ULL, 0x113f9804bef90daeULL, 0x1b710b35131c471bULL,
+0x28db77f523047d84ULL, 0x32caab7b40c72493ULL, 0x3c9ebe0a15c9bebcULL, 0x431d67c49c100d4cULL,
+0x4cc5d4becb3e42b6ULL, 0x597f299cfc657e2aULL, 0x5fcb6fab3ad6faecULL, 0x6c44198c4a475817ULL
+};
+
+static void processblock(struct sha512 *s, const uint8_t *buf)
+{
+	uint64_t W[80], t1, t2, a, b, c, d, e, f, g, h;
+	int i;
+
+	for (i = 0; i < 16; i++) {
+		W[i] = (uint64_t)buf[8*i]<<56;
+		W[i] |= (uint64_t)buf[8*i+1]<<48;
+		W[i] |= (uint64_t)buf[8*i+2]<<40;
+		W[i] |= (uint64_t)buf[8*i+3]<<32;
+		W[i] |= (uint64_t)buf[8*i+4]<<24;
+		W[i] |= (uint64_t)buf[8*i+5]<<16;
+		W[i] |= (uint64_t)buf[8*i+6]<<8;
+		W[i] |= buf[8*i+7];
+	}
+	for (; i < 80; i++)
+		W[i] = R1(W[i-2]) + W[i-7] + R0(W[i-15]) + W[i-16];
+	a = s->h[0];
+	b = s->h[1];
+	c = s->h[2];
+	d = s->h[3];
+	e = s->h[4];
+	f = s->h[5];
+	g = s->h[6];
+	h = s->h[7];
+	for (i = 0; i < 80; i++) {
+		t1 = h + S1(e) + Ch(e,f,g) + K[i] + W[i];
+		t2 = S0(a) + Maj(a,b,c);
+		h = g;
+		g = f;
+		f = e;
+		e = d + t1;
+		d = c;
+		c = b;
+		b = a;
+		a = t1 + t2;
+	}
+	s->h[0] += a;
+	s->h[1] += b;
+	s->h[2] += c;
+	s->h[3] += d;
+	s->h[4] += e;
+	s->h[5] += f;
+	s->h[6] += g;
+	s->h[7] += h;
+}
+
+static void pad(struct sha512 *s)
+{
+	unsigned r = s->len % 128;
+
+	s->buf[r++] = 0x80;
+	if (r > 112) {
+		memset(s->buf + r, 0, 128 - r);
+		r = 0;
+		processblock(s, s->buf);
+	}
+	memset(s->buf + r, 0, 120 - r);
+	s->len *= 8;
+	s->buf[120] = s->len >> 56;
+	s->buf[121] = s->len >> 48;
+	s->buf[122] = s->len >> 40;
+	s->buf[123] = s->len >> 32;
+	s->buf[124] = s->len >> 24;
+	s->buf[125] = s->len >> 16;
+	s->buf[126] = s->len >> 8;
+	s->buf[127] = s->len;
+	processblock(s, s->buf);
+}
+
+static void sha512_init(struct sha512 *s)
+{
+	s->len = 0;
+	s->h[0] = 0x6a09e667f3bcc908ULL;
+	s->h[1] = 0xbb67ae8584caa73bULL;
+	s->h[2] = 0x3c6ef372fe94f82bULL;
+	s->h[3] = 0xa54ff53a5f1d36f1ULL;
+	s->h[4] = 0x510e527fade682d1ULL;
+	s->h[5] = 0x9b05688c2b3e6c1fULL;
+	s->h[6] = 0x1f83d9abfb41bd6bULL;
+	s->h[7] = 0x5be0cd19137e2179ULL;
+}
+
+static void sha512_sum(struct sha512 *s, uint8_t *md)
+{
+	int i;
+
+	pad(s);
+	for (i = 0; i < 8; i++) {
+		md[8*i] = s->h[i] >> 56;
+		md[8*i+1] = s->h[i] >> 48;
+		md[8*i+2] = s->h[i] >> 40;
+		md[8*i+3] = s->h[i] >> 32;
+		md[8*i+4] = s->h[i] >> 24;
+		md[8*i+5] = s->h[i] >> 16;
+		md[8*i+6] = s->h[i] >> 8;
+		md[8*i+7] = s->h[i];
+	}
+}
+
+static void sha512_update(struct sha512 *s, const void *m, unsigned long len)
+{
+	const uint8_t *p = m;
+	unsigned r = s->len % 128;
+
+	s->len += len;
+	if (r) {
+		if (len < 128 - r) {
+			memcpy(s->buf + r, p, len);
+			return;
+		}
+		memcpy(s->buf + r, p, 128 - r);
+		len -= 128 - r;
+		p += 128 - r;
+		processblock(s, s->buf);
+	}
+	for (; len >= 128; len -= 128, p += 128)
+		processblock(s, p);
+	memcpy(s->buf, p, len);
+}
+
+static const unsigned char b64[] =
+"./0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
+
+static char *to64(char *s, unsigned int u, int n)
+{
+	while (--n >= 0) {
+		*s++ = b64[u % 64];
+		u /= 64;
+	}
+	return s;
+}
+
+/* key limit is not part of the original design, added for DoS protection.
+ * rounds limit has been lowered (versus the reference/spec), also for DoS
+ * protection. runtime is O(klen^2 + klen*rounds) */
+#define KEY_MAX 256
+#define SALT_MAX 16
+#define ROUNDS_DEFAULT 5000
+#define ROUNDS_MIN 1000
+#define ROUNDS_MAX 9999999
+
+/* hash n bytes of the repeated md message digest */
+static void hashmd(struct sha512 *s, unsigned int n, const void *md)
+{
+	unsigned int i;
+
+	for (i = n; i > 64; i -= 64)
+		sha512_update(s, md, 64);
+	sha512_update(s, md, i);
+}
+
+static char *sha512crypt(const char *key, const char *setting, char *output)
+{
+	struct sha512 ctx;
+	unsigned char md[64], kmd[64], smd[64];
+	unsigned int i, r, klen, slen;
+	char rounds[20] = "";
+	const char *salt;
+	char *p;
+
+	/* reject large keys */
+	for (i = 0; i <= KEY_MAX && key[i]; i++);
+	if (i > KEY_MAX)
+		return 0;
+	klen = i;
+
+	/* setting: $6$rounds=n$salt$ (rounds=n$ and closing $ are optional) */
+	if (strncmp(setting, "$6$", 3) != 0)
+		return 0;
+	salt = setting + 3;
+
+	r = ROUNDS_DEFAULT;
+	if (strncmp(salt, "rounds=", sizeof "rounds=" - 1) == 0) {
+		unsigned long u;
+		char *end;
+
+		/*
+		 * this is a deviation from the reference:
+		 * bad rounds setting is rejected if it is
+		 * - empty
+		 * - unterminated (missing '$')
+		 * - begins with anything but a decimal digit
+		 * the reference implementation treats these bad
+		 * rounds as part of the salt or parse them with
+		 * strtoul semantics which may cause problems
+		 * including non-portable hashes that depend on
+		 * the host's value of ULONG_MAX.
+		 */
+		salt += sizeof "rounds=" - 1;
+		if (!isdigit(*salt))
+			return 0;
+		u = strtoul(salt, &end, 10);
+		if (*end != '$')
+			return 0;
+		salt = end+1;
+		if (u < ROUNDS_MIN)
+			r = ROUNDS_MIN;
+		else if (u > ROUNDS_MAX)
+			return 0;
+		else
+			r = u;
+		/* needed when rounds is zero prefixed or out of bounds */
+		sprintf(rounds, "rounds=%u$", r);
+	}
+
+	for (i = 0; i < SALT_MAX && salt[i] && salt[i] != '$'; i++)
+		/* reject characters that interfere with /etc/shadow parsing */
+		if (salt[i] == '\n' || salt[i] == ':')
+			return 0;
+	slen = i;
+
+	/* B = sha(key salt key) */
+	sha512_init(&ctx);
+	sha512_update(&ctx, key, klen);
+	sha512_update(&ctx, salt, slen);
+	sha512_update(&ctx, key, klen);
+	sha512_sum(&ctx, md);
+
+	/* A = sha(key salt repeat-B alternate-B-key) */
+	sha512_init(&ctx);
+	sha512_update(&ctx, key, klen);
+	sha512_update(&ctx, salt, slen);
+	hashmd(&ctx, klen, md);
+	for (i = klen; i > 0; i >>= 1)
+		if (i & 1)
+			sha512_update(&ctx, md, sizeof md);
+		else
+			sha512_update(&ctx, key, klen);
+	sha512_sum(&ctx, md);
+
+	/* DP = sha(repeat-key), this step takes O(klen^2) time */
+	sha512_init(&ctx);
+	for (i = 0; i < klen; i++)
+		sha512_update(&ctx, key, klen);
+	sha512_sum(&ctx, kmd);
+
+	/* DS = sha(repeat-salt) */
+	sha512_init(&ctx);
+	for (i = 0; i < 16 + md[0]; i++)
+		sha512_update(&ctx, salt, slen);
+	sha512_sum(&ctx, smd);
+
+	/* iterate A = f(A,DP,DS), this step takes O(rounds*klen) time */
+	for (i = 0; i < r; i++) {
+		sha512_init(&ctx);
+		if (i % 2)
+			hashmd(&ctx, klen, kmd);
+		else
+			sha512_update(&ctx, md, sizeof md);
+		if (i % 3)
+			sha512_update(&ctx, smd, slen);
+		if (i % 7)
+			hashmd(&ctx, klen, kmd);
+		if (i % 2)
+			sha512_update(&ctx, md, sizeof md);
+		else
+			hashmd(&ctx, klen, kmd);
+		sha512_sum(&ctx, md);
+	}
+
+	/* output is $6$rounds=n$salt$hash */
+	p = output;
+	p += sprintf(p, "$6$%s%.*s$", rounds, slen, salt);
+#if 1
+	static const unsigned char perm[][3] = {
+		0,21,42,22,43,1,44,2,23,3,24,45,25,46,4,
+		47,5,26,6,27,48,28,49,7,50,8,29,9,30,51,
+		31,52,10,53,11,32,12,33,54,34,55,13,56,14,35,
+		15,36,57,37,58,16,59,17,38,18,39,60,40,61,19,
+		62,20,41 };
+	for (i=0; i<21; i++) p = to64(p,
+		(md[perm[i][0]]<<16)|(md[perm[i][1]]<<8)|md[perm[i][2]], 4);
+#else
+	p = to64(p, (md[0]<<16)|(md[21]<<8)|md[42], 4);
+	p = to64(p, (md[22]<<16)|(md[43]<<8)|md[1], 4);
+	p = to64(p, (md[44]<<16)|(md[2]<<8)|md[23], 4);
+	p = to64(p, (md[3]<<16)|(md[24]<<8)|md[45], 4);
+	p = to64(p, (md[25]<<16)|(md[46]<<8)|md[4], 4);
+	p = to64(p, (md[47]<<16)|(md[5]<<8)|md[26], 4);
+	p = to64(p, (md[6]<<16)|(md[27]<<8)|md[48], 4);
+	p = to64(p, (md[28]<<16)|(md[49]<<8)|md[7], 4);
+	p = to64(p, (md[50]<<16)|(md[8]<<8)|md[29], 4);
+	p = to64(p, (md[9]<<16)|(md[30]<<8)|md[51], 4);
+	p = to64(p, (md[31]<<16)|(md[52]<<8)|md[10], 4);
+	p = to64(p, (md[53]<<16)|(md[11]<<8)|md[32], 4);
+	p = to64(p, (md[12]<<16)|(md[33]<<8)|md[54], 4);
+	p = to64(p, (md[34]<<16)|(md[55]<<8)|md[13], 4);
+	p = to64(p, (md[56]<<16)|(md[14]<<8)|md[35], 4);
+	p = to64(p, (md[15]<<16)|(md[36]<<8)|md[57], 4);
+	p = to64(p, (md[37]<<16)|(md[58]<<8)|md[16], 4);
+	p = to64(p, (md[59]<<16)|(md[17]<<8)|md[38], 4);
+	p = to64(p, (md[18]<<16)|(md[39]<<8)|md[60], 4);
+	p = to64(p, (md[40]<<16)|(md[61]<<8)|md[19], 4);
+	p = to64(p, (md[62]<<16)|(md[20]<<8)|md[41], 4);
+#endif
+	p = to64(p, md[63], 2);
+	*p = 0;
+	return output;
+}
+
+char *__crypt_sha512(const char *key, const char *setting, char *output)
+{
+	static const char testkey[] = "Xy01@#\x01\x02\x80\x7f\xff\r\n\x81\t !";
+	static const char testsetting[] = "$6$rounds=1234$abc0123456789$";
+	static const char testhash[] = "$6$rounds=1234$abc0123456789$BCpt8zLrc/RcyuXmCDOE1ALqMXB2MH6n1g891HhFj8.w7LxGv.FTkqq6Vxc/km3Y0jE0j24jY5PIv/oOu6reg1";
+	char testbuf[128];
+	char *p, *q;
+
+	p = sha512crypt(key, setting, output);
+	/* self test and stack cleanup */
+	q = sha512crypt(testkey, testsetting, testbuf);
+	if (!p || q != testbuf || memcmp(testbuf, testhash, sizeof testhash))
+		return "*";
+	return p;
+}

libc/musl/src/crypt/encrypt.c

@@ -0,0 +1,52 @@
+#include <stdint.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+#include "crypt_des.h"
+
+static struct expanded_key __encrypt_key;
+
+void setkey(const char *key)
+{
+	unsigned char bkey[8];
+	int i, j;
+
+	for (i = 0; i < 8; i++) {
+		bkey[i] = 0;
+		for (j = 7; j >= 0; j--, key++)
+			bkey[i] |= (uint32_t)(*key & 1) << j;
+	}
+
+	__des_setkey(bkey, &__encrypt_key);
+}
+
+void encrypt(char *block, int edflag)
+{
+	struct expanded_key decrypt_key, *key;
+	uint32_t b[2];
+	int i, j;
+	char *p;
+
+	p = block;
+	for (i = 0; i < 2; i++) {
+		b[i] = 0;
+		for (j = 31; j >= 0; j--, p++)
+			b[i] |= (uint32_t)(*p & 1) << j;
+	}
+
+	key = &__encrypt_key;
+	if (edflag) {
+		key = &decrypt_key;
+		for (i = 0; i < 16; i++) {
+			decrypt_key.l[i] = __encrypt_key.l[15-i];
+			decrypt_key.r[i] = __encrypt_key.r[15-i];
+		}
+	}
+
+	__do_des(b[0], b[1], b, b + 1, 1, 0, key);
+
+	p = block;
+	for (i = 0; i < 2; i++)
+		for (j = 31; j >= 0; j--)
+			*p++ = b[i]>>j & 1;
+}

libc/musl/src/ctype/__ctype_b_loc.c

@@ -0,0 +1,41 @@
+#include <endian.h>
+
+#if __BYTE_ORDER == __BIG_ENDIAN
+#define X(x) x
+#else
+#define X(x) (((x)/256 | (x)*256) % 65536)
+#endif
+
+static const unsigned short table[] = {
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),
+X(0x200),X(0x320),X(0x220),X(0x220),X(0x220),X(0x220),X(0x200),X(0x200),
+X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),
+X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),X(0x200),
+X(0x160),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),
+X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),
+X(0x8d8),X(0x8d8),X(0x8d8),X(0x8d8),X(0x8d8),X(0x8d8),X(0x8d8),X(0x8d8),
+X(0x8d8),X(0x8d8),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),
+X(0x4c0),X(0x8d5),X(0x8d5),X(0x8d5),X(0x8d5),X(0x8d5),X(0x8d5),X(0x8c5),
+X(0x8c5),X(0x8c5),X(0x8c5),X(0x8c5),X(0x8c5),X(0x8c5),X(0x8c5),X(0x8c5),
+X(0x8c5),X(0x8c5),X(0x8c5),X(0x8c5),X(0x8c5),X(0x8c5),X(0x8c5),X(0x8c5),
+X(0x8c5),X(0x8c5),X(0x8c5),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),
+X(0x4c0),X(0x8d6),X(0x8d6),X(0x8d6),X(0x8d6),X(0x8d6),X(0x8d6),X(0x8c6),
+X(0x8c6),X(0x8c6),X(0x8c6),X(0x8c6),X(0x8c6),X(0x8c6),X(0x8c6),X(0x8c6),
+X(0x8c6),X(0x8c6),X(0x8c6),X(0x8c6),X(0x8c6),X(0x8c6),X(0x8c6),X(0x8c6),
+X(0x8c6),X(0x8c6),X(0x8c6),X(0x4c0),X(0x4c0),X(0x4c0),X(0x4c0),X(0x200),
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+};
+
+static const unsigned short *const ptable = table+128;
+
+const unsigned short **__ctype_b_loc(void)
+{
+	return (void *)&ptable;
+}

libc/musl/src/ctype/__ctype_get_mb_cur_max.c

@@ -0,0 +1,7 @@
+#include <stdlib.h>
+#include "locale_impl.h"
+
+size_t __ctype_get_mb_cur_max()
+{
+	return MB_CUR_MAX;
+}

libc/musl/src/ctype/__ctype_tolower_loc.c

@@ -0,0 +1,30 @@
+#include <stdint.h>
+
+static const int32_t table[] = {
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,
+16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,
+32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,
+48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,
+64,
+'a','b','c','d','e','f','g','h','i','j','k','l','m',
+'n','o','p','q','r','s','t','u','v','w','x','y','z',
+91,92,93,94,95,96,
+'a','b','c','d','e','f','g','h','i','j','k','l','m',
+'n','o','p','q','r','s','t','u','v','w','x','y','z',
+123,124,125,126,127,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+};
+
+static const int32_t *const ptable = table+128;
+
+const int32_t **__ctype_tolower_loc(void)
+{
+	return (void *)&ptable;
+}

libc/musl/src/ctype/__ctype_toupper_loc.c

@@ -0,0 +1,30 @@
+#include <stdint.h>
+
+static const int32_t table[] = {
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,
+16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,
+32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,
+48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,
+64,
+'A','B','C','D','E','F','G','H','I','J','K','L','M',
+'N','O','P','Q','R','S','T','U','V','W','X','Y','Z',
+91,92,93,94,95,96,
+'A','B','C','D','E','F','G','H','I','J','K','L','M',
+'N','O','P','Q','R','S','T','U','V','W','X','Y','Z',
+123,124,125,126,127,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+};
+
+static const int32_t *const ptable = table+128;
+
+const int32_t **__ctype_toupper_loc(void)
+{
+	return (void *)&ptable;
+}

libc/musl/src/ctype/alpha.h

@@ -0,0 +1,163 @@
+18,17,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,17,34,35,36,17,37,38,39,40,
+41,42,43,44,17,45,46,47,16,16,48,16,16,16,16,16,16,16,49,50,51,16,52,53,16,16,
+17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,54,
+17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,
+17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,
+17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,
+17,17,17,55,17,17,17,17,56,17,57,58,59,60,61,62,17,17,17,17,17,17,17,17,17,17,
+17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,
+17,17,17,17,17,17,17,63,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,
+16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,17,64,65,17,66,67,
+68,69,70,71,72,73,74,17,75,76,77,78,79,80,16,16,16,81,82,83,84,85,86,87,88,89,
+16,90,16,91,92,16,16,17,17,17,93,94,95,16,16,16,16,16,16,16,16,16,16,17,17,17,
+17,96,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,17,17,97,16,16,16,16,16,16,
+16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,
+16,17,17,98,99,16,16,16,100,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,
+17,17,17,17,17,17,17,101,17,17,102,16,16,16,16,16,16,16,16,16,16,16,16,16,16,
+16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,17,103,
+104,16,16,16,16,16,16,16,16,16,105,16,16,16,16,16,16,16,16,16,16,16,16,16,16,
+16,16,16,16,16,16,16,16,16,106,107,108,109,16,16,16,16,16,16,16,16,110,16,16,
+16,16,16,16,16,111,112,16,16,16,16,113,16,16,114,16,16,16,16,16,16,16,16,16,
+16,16,16,16,
+16,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,0,0,0,0,0,0,0,0,254,255,255,7,254,
+255,255,7,0,0,0,0,0,4,32,4,255,255,127,255,255,255,127,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
+255,195,255,3,0,31,80,0,0,0,0,0,0,0,0,0,0,32,0,0,0,0,0,223,188,64,215,255,255,
+251,255,255,255,255,255,255,255,255,255,191,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,3,252,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,254,255,255,255,127,2,254,255,255,
+255,255,0,0,0,0,0,255,191,182,0,255,255,255,7,7,0,0,0,255,7,255,255,255,255,
+255,255,255,254,255,195,255,255,255,255,255,255,255,255,255,255,255,255,239,
+31,254,225,255,
+159,0,0,255,255,255,255,255,255,0,224,255,255,255,255,255,255,255,255,255,255,
+255,255,3,0,255,255,255,255,255,7,48,4,255,255,255,252,255,31,0,0,255,255,255,
+1,255,7,0,0,0,0,0,0,255,255,223,63,0,0,240,255,248,3,255,255,255,255,255,255,
+255,255,255,239,255,223,225,255,207,255,254,255,239,159,249,255,255,253,197,
+227,159,89,128,176,207,255,3,16,238,135,249,255,255,253,109,195,135,25,2,94,
+192,255,63,0,238,191,251,255,255,253,237,227,191,27,1,0,207,255,0,30,238,159,
+249,255,255,253,237,227,159,25,192,176,207,255,2,0,236,199,61,214,24,199,255,
+195,199,29,129,0,192,255,0,0,239,223,253,255,255,253,255,227,223,29,96,7,207,
+255,0,0,239,223,253,255,255,253,239,227,223,29,96,64,207,255,6,0,239,223,253,
+255,255,255,255,231,223,93,240,128,207,255,0,252,236,255,127,252,255,255,251,
+47,127,128,95,255,192,255,12,0,254,255,255,255,255,127,255,7,63,32,255,3,0,0,
+0,0,150,37,240,254,174,236,255,59,95,32,255,243,0,0,0,
+0,1,0,0,0,255,3,0,0,255,254,255,255,255,31,254,255,3,255,255,254,255,255,255,
+31,0,0,0,0,0,0,0,0,255,255,255,255,255,255,127,249,255,3,255,255,231,193,255,
+255,127,64,255,51,255,255,255,255,191,32,255,255,255,255,255,247,255,255,255,
+255,255,255,255,255,255,61,127,61,255,255,255,255,255,61,255,255,255,255,61,
+127,61,255,127,255,255,255,255,255,255,255,61,255,255,255,255,255,255,255,255,
+135,0,0,0,0,255,255,0,0,255,255,255,255,255,255,255,255,255,255,63,63,254,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,159,255,255,254,255,255,7,255,255,255,255,255,255,255,255,
+255,199,255,1,255,223,15,0,255,255,15,0,255,255,15,0,255,223,13,0,255,255,255,
+255,255,255,207,255,255,1,128,16,255,3,0,0,0,0,255,3,255,255,255,255,255,255,
+255,255,255,255,255,0,255,255,255,255,255,7,255,255,255,255,255,255,255,255,
+63,
+0,255,255,255,127,255,15,255,1,192,255,255,255,255,63,31,0,255,255,255,255,
+255,15,255,255,255,3,255,3,0,0,0,0,255,255,255,15,255,255,255,255,255,255,255,
+127,254,255,31,0,255,3,255,3,128,0,0,0,0,0,0,0,0,0,0,0,255,255,255,255,255,
+255,239,255,239,15,255,3,0,0,0,0,255,255,255,255,255,243,255,255,255,255,255,
+255,191,255,3,0,255,255,255,255,255,255,63,0,255,227,255,255,255,255,255,63,
+255,1,0,0,0,0,0,0,0,0,0,0,0,222,111,0,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,0,0,0,0,128,255,31,0,
+255,255,63,63,255,255,255,255,63,63,255,170,255,255,255,63,255,255,255,255,
+255,255,223,95,220,31,207,15,255,31,220,31,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,128,
+0,0,255,31,0,0,0,0,0,0,0,0,0,0,0,0,132,252,47,62,80,189,255,243,224,67,0,0,
+255,255,255,255,255,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,192,255,255,255,255,255,255,3,0,
+0,255,255,255,255,255,127,255,255,255,255,255,127,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,31,120,12,0,255,255,255,255,191,32,255,
+255,255,255,255,255,255,128,0,0,255,255,127,0,127,127,127,127,127,127,127,127,
+255,255,255,255,0,0,0,0,0,128,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,224,0,0,0,254,3,62,31,254,255,255,255,255,255,255,255,255,255,127,224,254,
+255,255,255,255,255,255,255,255,255,255,247,224,255,255,255,255,127,254,255,
+255,255,255,255,255,255,255,255,255,127,0,0,255,255,255,7,0,0,0,0,0,0,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,63,0,0,0,0,0,0,0,0,0,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,7,0,
+0,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,31,0,0,
+0,0,0,0,0,0,255,255,255,255,255,63,255,31,255,255,255,15,0,0,255,255,255,255,
+255,127,240,143,255,255,255,255,255,255,255,255,255,255,255,255,255,255,0,0,0,
+0,128,255,252,255,255,255,255,255,255,255,255,255,255,255,255,249,255,255,255,
+127,255,0,0,0,0,0,0,0,128,255,187,247,255,255,255,0,0,0,255,255,255,255,255,
+255,15,0,255,255,255,255,255,255,255,255,47,0,255,3,0,0,252,40,255,255,255,
+255,255,7,255,255,255,255,7,0,255,255,255,31,255,255,255,255,255,255,247,255,
+0,128,255,3,223,255,255,127,255,255,255,255,255,255,127,0,255,63,255,3,255,
+255,127,196,255,255,255,255,255,255,255,127,5,0,0,56,255,255,60,0,126,126,126,
+0,127,127,255,255,255,255,255,247,63,0,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,7,255,3,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,15,0,255,255,127,248,255,255,255,255,255,
+15,255,255,255,255,255,255,255,255,255,255,255,255,255,63,255,255,255,255,255,
+255,255,255,255,255,255,255,255,3,0,0,0,0,127,0,248,224,255,253,127,95,219,
+255,255,255,255,255,255,255,255,255,255,255,255,255,3,0,0,0,248,255,255,255,
+255,255,255,255,255,255,255,255,255,63,0,0,255,255,255,255,255,255,255,255,
+252,255,255,255,255,255,255,0,0,0,0,0,255,15,0,0,0,0,0,0,0,0,0,0,0,0,0,0,223,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,31,0,0,255,3,
+254,255,255,7,254,255,255,7,192,255,255,255,255,255,255,255,255,255,255,127,
+252,252,252,28,0,0,0,0,255,239,255,255,127,255,255,183,255,63,255,63,0,0,0,0,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,7,0,0,0,0,0,0,0,0,
+255,255,255,255,255,255,31,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,255,255,255,31,255,255,255,255,255,255,1,0,0,0,0,
+0,255,255,255,255,0,224,255,255,255,7,255,255,255,255,255,7,255,255,255,63,
+255,255,255,255,15,255,62,0,0,0,0,0,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,63,255,3,255,255,255,255,15,255,255,255,
+255,15,255,255,255,255,255,0,255,255,255,255,255,255,15,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,255,255,255,255,255,255,127,0,255,255,63,0,255,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,63,253,255,255,255,255,191,145,255,255,63,0,255,255,
+127,0,255,255,255,127,0,0,0,0,0,0,0,0,255,255,55,0,255,255,63,0,255,255,255,3,
+0,0,0,0,0,0,0,0,255,255,255,255,255,255,255,192,0,0,0,0,0,0,0,0,111,240,239,
+254,255,255,15,0,0,0,0,0,255,255,255,31,255,255,255,31,0,0,0,0,255,254,255,
+255,31,0,0,0,255,255,255,255,255,255,63,0,255,255,63,0,255,255,7,0,255,255,3,
+0,0,0,0,0,0,0,0,0,0,0,0,
+0,255,255,255,255,255,255,255,255,255,1,0,0,0,0,0,0,255,255,255,255,255,255,7,
+0,255,255,255,255,255,255,7,0,255,255,255,255,255,255,255,255,63,0,0,0,192,
+255,0,0,252,255,255,255,255,255,255,1,0,0,255,255,255,1,255,3,255,255,255,255,
+255,255,199,255,0,0,255,255,255,255,71,0,255,255,255,255,255,255,255,255,30,0,
+255,23,0,0,0,0,255,255,251,255,255,255,159,64,0,0,0,0,0,0,0,0,127,189,255,191,
+255,1,255,255,255,255,255,255,255,1,255,3,239,159,249,255,255,253,237,227,159,
+25,129,224,15,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,255,255,255,255,255,255,
+255,255,187,7,255,3,0,0,0,0,255,255,255,255,255,255,255,255,179,0,255,3,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,255,255,255,255,255,255,63,127,0,0,0,63,0,0,
+0,0,255,255,255,255,255,255,255,127,17,0,255,3,0,0,0,0,255,255,255,255,255,
+255,63,0,255,3,0,0,0,0,0,
+0,255,255,255,227,255,7,255,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,255,255,255,255,255,255,255,255,255,3,
+0,128,255,255,255,255,255,255,231,127,0,0,255,255,255,255,255,255,207,255,255,
+0,0,0,0,0,255,255,255,255,255,255,255,1,255,253,255,255,255,255,127,127,1,0,
+255,3,0,0,252,255,255,255,252,255,255,254,127,0,0,0,0,0,0,0,0,0,127,251,255,
+255,255,255,127,180,203,0,255,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,3,0,0,
+0,0,0,0,0,0,0,0,0,0,255,255,255,255,255,255,255,255,255,255,255,255,255,127,0,
+0,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,15,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,255,255,255,255,255,127,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+255,255,255,255,255,255,255,255,127,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,255,255,255,255,255,255,255,1,255,255,255,127,255,3,0,0,0,0,0,0,0,0,0,0,0,
+0,255,255,255,63,0,0,255,255,255,255,255,255,127,0,15,0,255,3,248,255,255,224,
+255,255,0,0,0,0,0,0,0,0,0,0,0,0,0,0,255,255,255,255,255,255,255,255,31,0,255,
+255,255,255,255,127,0,0,248,255,0,0,0,0,0,0,0,0,3,0,0,0,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,31,0,0,255,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,7,0,
+255,255,255,127,0,0,0,0,0,0,0,0,0,0,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,15,255,255,255,255,255,
+255,255,255,255,255,255,255,255,7,255,31,255,1,255,67,0,0,0,0,0,0,0,0,0,0,0,0,
+255,255,255,255,255,255,255,255,255,255,223,255,255,255,255,255,255,255,255,
+223,100,222,255,235,239,255,255,255,255,255,255,255,191,231,223,223,255,255,
+255,123,95,252,253,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,63,255,255,255,253,255,255,247,255,255,255,
+247,255,255,223,255,255,255,223,255,255,127,255,255,255,127,255,255,255,253,
+255,255,255,253,255,255,247,207,255,255,255,255,255,255,127,255,255,249,219,7,
+0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,255,255,255,255,255,255,
+255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,31,0,
+0,0,0,0,0,
+0,255,255,255,255,255,255,255,255,143,0,255,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0,0,239,255,255,255,150,254,247,10,132,234,150,170,150,247,247,94,255,251,
+255,15,238,251,255,15,0,0,0,0,0,0,0,0,0,0,0,0,0,0,255,255,255,3,255,255,255,3,
+255,255,255,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

libc/musl/src/ctype/isalnum.c

@@ -0,0 +1,13 @@
+#include <ctype.h>
+
+int isalnum(int c)
+{
+	return isalpha(c) || isdigit(c);
+}
+
+int __isalnum_l(int c, locale_t l)
+{
+	return isalnum(c);
+}
+
+weak_alias(__isalnum_l, isalnum_l);

libc/musl/src/ctype/isalpha.c

@@ -0,0 +1,14 @@
+#include <ctype.h>
+#undef isalpha
+
+int isalpha(int c)
+{
+	return ((unsigned)c|32)-'a' < 26;
+}
+
+int __isalpha_l(int c, locale_t l)
+{
+	return isalpha(c);
+}
+
+weak_alias(__isalpha_l, isalpha_l);

libc/musl/src/ctype/isascii.c

@@ -0,0 +1,7 @@
+#include <ctype.h>
+#undef isascii
+
+int isascii(int c)
+{
+	return !(c&~0x7f);
+}

libc/musl/src/ctype/isblank.c

@@ -0,0 +1,13 @@
+#include <ctype.h>
+
+int isblank(int c)
+{
+	return (c == ' ' || c == '\t');
+}
+
+int __isblank_l(int c, locale_t l)
+{
+	return isblank(c);
+}
+
+weak_alias(__isblank_l, isblank_l);

libc/musl/src/ctype/iscntrl.c

@@ -0,0 +1,13 @@
+#include <ctype.h>
+
+int iscntrl(int c)
+{
+	return (unsigned)c < 0x20 || c == 0x7f;
+}
+
+int __iscntrl_l(int c, locale_t l)
+{
+	return iscntrl(c);
+}
+
+weak_alias(__iscntrl_l, iscntrl_l);

libc/musl/src/ctype/isdigit.c

@@ -0,0 +1,14 @@
+#include <ctype.h>
+#undef isdigit
+
+int isdigit(int c)
+{
+	return (unsigned)c-'0' < 10;
+}
+
+int __isdigit_l(int c, locale_t l)
+{
+	return isdigit(c);
+}
+
+weak_alias(__isdigit_l, isdigit_l);

libc/musl/src/ctype/isgraph.c

@@ -0,0 +1,14 @@
+#include <ctype.h>
+#undef isgraph
+
+int isgraph(int c)
+{
+	return (unsigned)c-0x21 < 0x5e;
+}
+
+int __isgraph_l(int c, locale_t l)
+{
+	return isgraph(c);
+}
+
+weak_alias(__isgraph_l, isgraph_l);

libc/musl/src/ctype/islower.c

@@ -0,0 +1,14 @@
+#include <ctype.h>
+#undef islower
+
+int islower(int c)
+{
+	return (unsigned)c-'a' < 26;
+}
+
+int __islower_l(int c, locale_t l)
+{
+	return islower(c);
+}
+
+weak_alias(__islower_l, islower_l);

libc/musl/src/ctype/isprint.c

@@ -0,0 +1,14 @@
+#include <ctype.h>
+#undef isprint
+
+int isprint(int c)
+{
+	return (unsigned)c-0x20 < 0x5f;
+}
+
+int __isprint_l(int c, locale_t l)
+{
+	return isprint(c);
+}
+
+weak_alias(__isprint_l, isprint_l);