Get more Python tests passing (#141)

examples/rusage.c

@@ -106,6 +106,7 @@ int main(int argc, char *argv[]) {
     sigaction(SIGQUIT, &dflt, 0);
     sigprocmask(SIG_SETMASK, &savemask, 0);
     execvp(argv[1], argv + 1);
+    fprintf(stderr, "exec failed %d\n", errno);
     _Exit(127);
   }
   while (wait4(pid, &wstatus, 0, &rusage) == -1) {

libc/calls/calls.h

@@ -117,6 +117,7 @@ int fsync(int);
 int ftruncate(int, int64_t);
 int getdomainname(char *, size_t);
 int gethostname(char *, size_t);
+int getpgid(int);
 int getpid(void);
 int getppid(void);
 int getpriority(int, unsigned);
@@ -151,8 +152,8 @@ int posix_fadvise(int, uint64_t, uint64_t, int);
 int posix_madvise(void *, uint64_t, int);
 int prctl();
 int raise(int);
-int reboot(int);
 int readlink(const char *, char *, size_t);
+int reboot(int);
 int remove(const char *);
 int rename(const char *, const char *);
 int renameat(int, const char *, int, const char *);

libc/calls/commandv.c

@@ -26,8 +26,9 @@ static bool AccessCommand(char path[hasatleast PATH_MAX], const char *name,
                           size_t namelen, size_t pathlen) {
   if (pathlen + 1 + namelen + 1 + 4 + 1 > PATH_MAX) return -1;
   if (pathlen && (path[pathlen - 1] != '/' && path[pathlen - 1] != '\\')) {
-    path[pathlen] =
+    path[pathlen] = !IsWindows()                  ? '/'
-        !IsWindows() ? '/' : memchr(path, '\\', pathlen) ? '\\' : '/';
+                    : memchr(path, '\\', pathlen) ? '\\'
+                                                  : '/';
     pathlen++;
   }
   memcpy(path + pathlen, name, namelen + 1);
@@ -93,7 +94,9 @@ char *commandv(const char *name, char pathbuf[hasatleast PATH_MAX]) {
        (AccessCommand(pathbuf, name, namelen,
                       stpcpy(pathbuf, kNtSystemDirectory) - pathbuf) ||
         AccessCommand(pathbuf, name, namelen,
-                      stpcpy(pathbuf, kNtWindowsDirectory) - pathbuf))) ||
+                      stpcpy(pathbuf, kNtWindowsDirectory) - pathbuf) ||
+        AccessCommand(pathbuf, name, namelen,
+                      stpcpy(pathbuf, ".") - pathbuf))) ||
       SearchPath(pathbuf, name, namelen)) {
     errno = olderr;
     return pathbuf;

libc/calls/getconsolectrlevent.h

@@ -0,0 +1,23 @@
+#ifndef COSMOPOLITAN_LIBC_CALLS_GETCONSOLECTRLEVENT_H_
+#define COSMOPOLITAN_LIBC_CALLS_GETCONSOLECTRLEVENT_H_
+#include "libc/nt/enum/ctrlevent.h"
+#include "libc/sysv/consts/sig.h"
+#if !(__ASSEMBLER__ + __LINKER__ + 0)
+COSMOPOLITAN_C_START_
+
+static inline int GetConsoleCtrlEvent(int sig) {
+  switch (sig) {
+    case SIGINT:
+      return kNtCtrlCEvent;
+    case SIGHUP:
+      return kNtCtrlCloseEvent;
+    case SIGQUIT:
+      return kNtCtrlBreakEvent;
+    default:
+      return -1;
+  }
+}
+
+COSMOPOLITAN_C_END_
+#endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */
+#endif /* COSMOPOLITAN_LIBC_CALLS_GETCONSOLECTRLEVENT_H_ */

libc/calls/getpgid.c

@@ -0,0 +1,32 @@
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+╞══════════════════════════════════════════════════════════════════════════════╡
+│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
+│                                                                              │
+│ Permission to use, copy, modify, and/or 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.                                                │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/calls/calls.h"
+#include "libc/calls/internal.h"
+#include "libc/dce.h"
+
+/**
+ * Returns process group id.
+ */
+int getpgid(int pid) {
+  if (!IsWindows()) {
+    return sys_getpgid(pid);
+  } else {
+    return getpid();
+  }
+}

libc/calls/getrusage-nt.c

@@ -38,8 +38,8 @@ textwindows int sys_getrusage_nt(int who, struct rusage *usage) {
   if ((who == RUSAGE_SELF ? GetProcessTimes : GetThreadTimes)(
           (who == RUSAGE_SELF ? GetCurrentProcess : GetCurrentThread)(),
           &CreationFileTime, &ExitFileTime, &KernelFileTime, &UserFileTime)) {
-    FileTimeToTimeVal(&usage->ru_utime, UserFileTime);
+    usage->ru_utime = FileTimeToTimeVal(UserFileTime);
-    FileTimeToTimeVal(&usage->ru_stime, KernelFileTime);
+    usage->ru_stime = FileTimeToTimeVal(KernelFileTime);
     return 0;
   } else {
     return __winerr();

libc/calls/gettimeofday-nt.c

@@ -28,7 +28,7 @@
 int sys_gettimeofday_nt(struct timeval *tv, struct timezone *tz) {
   struct NtFileTime ft;
   GetSystemTimeAsFileTime(&ft);
-  FileTimeToTimeVal(tv, ft);
+  *tv = FileTimeToTimeVal(ft);
   if (tz) memset(tz, 0, sizeof(*tz));
   return 0;
 }

libc/calls/internal.h

@@ -142,6 +142,7 @@ i32 sys_ftruncate(i32, i64, i64) hidden;
 i32 sys_futimes(i32, const struct timeval *) hidden;
 i32 sys_futimesat(i32, const char *, const struct timeval *) hidden;
 i32 sys_getitimer(i32, struct itimerval *) hidden;
+i32 sys_getpgid(i32) hidden;
 i32 sys_getppid(void) hidden;
 i32 sys_getpriority(i32, u32) hidden;
 i32 sys_getrlimit(i32, struct rlimit *) hidden;

libc/calls/kill-nt.c

@@ -16,41 +16,37 @@
 │ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR             │
 │ PERFORMANCE OF THIS SOFTWARE.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/calls/getconsolectrlevent.h"
 #include "libc/calls/internal.h"
 #include "libc/dce.h"
 #include "libc/macros.internal.h"
 #include "libc/nt/console.h"
 #include "libc/nt/enum/ctrlevent.h"
+#include "libc/nt/enum/processaccess.h"
 #include "libc/nt/process.h"
+#include "libc/nt/runtime.h"
 #include "libc/sysv/errfuns.h"
 
 textwindows int sys_kill_nt(int pid, int sig) {
-  int target;
+  bool ok;
-  uint32_t event;
+  int event;
-  if (!pid) return raise(sig);
+  int64_t handle;
-  if ((pid > 0 && __isfdkind(pid, kFdProcess)) ||
+  if (pid) {
-      (pid < 0 && __isfdkind(-pid, kFdProcess))) {
+    pid = ABS(pid);
-    target = GetProcessId(g_fds.p[ABS(pid)].handle);
+    if ((event = GetConsoleCtrlEvent(sig)) != -1) {
-  } else {
+      ok = !!GenerateConsoleCtrlEvent(
-    target = pid;
+          event, __isfdkind(pid, kFdProcess) ? GetProcessId(g_fds.p[pid].handle)
-  }
+                                             : pid);
-  if (target == GetCurrentProcessId()) {
+    } else if (__isfdkind(pid, kFdProcess)) {
-    return raise(sig);
+      ok = !!TerminateProcess(g_fds.p[pid].handle, 128 + sig);
-  } else {
+    } else if ((handle = OpenProcess(kNtProcessAllAccess, false, pid))) {
-    switch (sig) {
+      ok = !!TerminateProcess(handle, 128 + sig);
-      case SIGINT:
+      CloseHandle(handle);
-        event = kNtCtrlCEvent;
-      case SIGHUP:
-        event = kNtCtrlCloseEvent;
-      case SIGQUIT:
-        event = kNtCtrlBreakEvent;
-      default:
-        return einval();
-    }
-    if (GenerateConsoleCtrlEvent(event, target)) {
-      return 0;
     } else {
-      return __winerr();
+      ok = false;
     }
+    return ok ? 0 : __winerr();
+  } else {
+    return raise(sig);
   }
 }

libc/calls/makedev.h

@@ -0,0 +1,56 @@
+#ifndef COSMOPOLITAN_LIBC_CALLS_MAKEDEV_H_
+#define COSMOPOLITAN_LIBC_CALLS_MAKEDEV_H_
+#include "libc/dce.h"
+#if !(__ASSEMBLER__ + __LINKER__ + 0)
+
+static inline uint64_t major(uint64_t x) {
+  if (IsXnu()) {
+    return (x >> 24) & 0xff;
+  } else if (IsNetbsd()) {
+    return (x & 0x000fff00) >> 8;
+  } else if (IsOpenbsd()) {
+    return (x >> 8) & 0xff;
+  } else if (IsFreebsd()) {
+    return ((x >> 32) & 0xffffff00) | ((x >> 8) & 0x000000ff);
+  } else {
+    return ((x >> 32) & 0xfffff000) | ((x >> 8) & 0x00000fff);
+  }
+}
+
+static inline uint64_t minor(uint64_t x) {
+  if (IsXnu()) {
+    return x & 0x00ffffff;
+  } else if (IsNetbsd()) {
+    return (x & 0x000000ff) | (x & 0xfff00000) >> 12;
+  } else if (IsOpenbsd()) {
+    return (x & 0x000000ff) | (x & 0x0ffff000) >> 8;
+  } else if (IsFreebsd()) {
+    return ((x >> 24) & 0x0000ff00) | (x & 0xffff00ff);
+  } else {
+    return ((x >> 12) & 0xffffff00) | (x & 0x000000ff);
+  }
+}
+
+static inline uint64_t makedev(uint64_t x, uint64_t y) {
+  if (IsXnu()) {
+    return x << 24 | y;
+  } else if (IsNetbsd()) {
+    return ((x << 8) & 0x000fff00) | ((y << 12) & 0xfff00000u) |
+           (y & 0x000000ff);
+  } else if (IsOpenbsd()) {
+    return (x & 0xff) << 8 | (y & 0xff) | (y & 0xffff00) << 8;
+  } else if (IsFreebsd()) {
+    return (x & 0xffffff00) << 32 | (x & 0x000000ff) << 8 |
+           (y & 0x0000ff00) << 24 | (y & 0xffff00ff);
+  } else {
+    return (x & 0xfffff000) << 32 | (x & 0x00000fff) << 8 |
+           (y & 0xffffff00) << 12 | (y & 0x000000ff);
+  }
+}
+
+#define major(x)      major(x)
+#define minor(x)      minor(x)
+#define makedev(x, y) makedev(x, y)
+
+#endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */
+#endif /* COSMOPOLITAN_LIBC_CALLS_MAKEDEV_H_ */

libc/calls/raise.c

@@ -17,25 +17,12 @@
 │ PERFORMANCE OF THIS SOFTWARE.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/calls/calls.h"
+#include "libc/calls/getconsolectrlevent.h"
 #include "libc/calls/internal.h"
 #include "libc/nt/console.h"
-#include "libc/nt/enum/ctrlevent.h"
 #include "libc/nt/runtime.h"
 #include "libc/sysv/consts/sig.h"
 
-static textwindows uint32_t GetCtrlEvent(int sig) {
-  switch (sig) {
-    case SIGINT:
-      return kNtCtrlCEvent;
-    case SIGHUP:
-      return kNtCtrlCloseEvent;
-    case SIGQUIT:
-      return kNtCtrlBreakEvent;
-    default:
-      ExitProcess(128 + sig);
-  }
-}
-
 /**
  * Sends signal to this process.
  *
@@ -44,6 +31,7 @@ static textwindows uint32_t GetCtrlEvent(int sig) {
  * @asyncsignalsafe
  */
 int raise(int sig) {
+  int event;
   if (sig == SIGTRAP) {
     DebugBreak();
     return 0;
@@ -55,11 +43,13 @@ int raise(int sig) {
   }
   if (!IsWindows()) {
     return sys_kill(getpid(), sig, 1);
-  } else {
+  } else if ((event = GetConsoleCtrlEvent(sig))) {
-    if (GenerateConsoleCtrlEvent(GetCtrlEvent(sig), 0)) {
+    if (GenerateConsoleCtrlEvent(event, 0)) {
       return 0;
     } else {
       return __winerr();
     }
+  } else {
+    ExitProcess(128 + sig);
   }
 }

libc/calls/symlinkat-nt.c

@@ -18,13 +18,25 @@
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/calls/calls.h"
 #include "libc/calls/internal.h"
+#include "libc/errno.h"
+#include "libc/nt/enum/accessmask.h"
+#include "libc/nt/enum/tokeninformationclass.h"
+#include "libc/nt/errors.h"
 #include "libc/nt/files.h"
+#include "libc/nt/privilege.h"
+#include "libc/nt/runtime.h"
+#include "libc/nt/struct/luid.h"
+#include "libc/nt/struct/tokenprivileges.h"
+#include "libc/sysv/errfuns.h"
+
+static bool g_can_symlink;
 
 textwindows int sys_symlinkat_nt(const char *target, int newdirfd,
-                             const char *linkpath) {
+                                 const char *linkpath) {
   uint32_t flags;
   char16_t target16[PATH_MAX];
   char16_t linkpath16[PATH_MAX];
+  if (!g_can_symlink) return eacces();
   flags = isdirectory(target) ? kNtSymbolicLinkFlagDirectory : 0;
   if (__mkntpathat(newdirfd, linkpath, 0, linkpath16) == -1) return -1;
   if (__mkntpath(target, target16) == -1) return -1;
@@ -34,3 +46,24 @@ textwindows int sys_symlinkat_nt(const char *target, int newdirfd,
     return __winerr();
   }
 }
+
+static textstartup bool EnableSymlink(void) {
+  int64_t tok;
+  struct NtLuid id;
+  struct NtTokenPrivileges tp;
+  if (!OpenProcessToken(GetCurrentProcess(), kNtTokenAllAccess, &tok)) return 0;
+  if (!LookupPrivilegeValue(0, u"SeCreateSymbolicLinkPrivilege", &id)) return 0;
+  tp.PrivilegeCount = 1;
+  tp.Privileges[0].Luid = id;
+  tp.Privileges[0].Attributes = kNtSePrivilegeEnabled;
+  if (!AdjustTokenPrivileges(tok, 0, &tp, sizeof(tp), 0, 0)) return 0;
+  return GetLastError() != kNtErrorNotAllAssigned;
+}
+
+static textstartup void g_can_symlink_init() {
+  g_can_symlink = EnableSymlink();
+}
+
+const void *const g_can_symlink_ctor[] initarray = {
+    g_can_symlink_init,
+};

libc/calls/wait4-nt.c

@@ -74,10 +74,17 @@ textwindows int sys_wait4_nt(int pid, int *opt_out_wstatus, int options,
     }
     if (opt_out_rusage) {
       memset(opt_out_rusage, 0, sizeof(*opt_out_rusage));
-      GetProcessTimes(GetCurrentProcess(), &createfiletime, &exitfiletime,
+      if (GetProcessTimes(g_fds.p[pids[i]].handle, &createfiletime,
-                      &kernelfiletime, &userfiletime);
+                          &exitfiletime, &kernelfiletime, &userfiletime)) {
-      FileTimeToTimeVal(&opt_out_rusage->ru_utime, userfiletime);
+        opt_out_rusage->ru_utime.tv_sec =
-      FileTimeToTimeVal(&opt_out_rusage->ru_stime, kernelfiletime);
+            ReadFileTime(userfiletime) / HECTONANOSECONDS;
+        opt_out_rusage->ru_utime.tv_usec =
+            ReadFileTime(userfiletime) % HECTONANOSECONDS;
+        opt_out_rusage->ru_stime.tv_sec =
+            ReadFileTime(kernelfiletime) / HECTONANOSECONDS;
+        opt_out_rusage->ru_stime.tv_usec =
+            ReadFileTime(kernelfiletime) % HECTONANOSECONDS;
+      }
     }
     CloseHandle(g_fds.p[pids[i]].handle);
     g_fds.p[pids[i]].kind = kFdEmpty;

libc/fmt/conv.h

@@ -39,20 +39,27 @@ long sizetol(const char *, long) paramsnonnull() libcesque;
 │ cosmopolitan § conversion » time                                         ─╬─│┼
 ╚────────────────────────────────────────────────────────────────────────────│*/
 
-struct timespec WindowsTimeToTime(uint64_t);
+int64_t DosDateTimeToUnix(unsigned, unsigned) nothrow;
-int64_t DosDateTimeToUnix(unsigned, unsigned);
+struct timeval WindowsTimeToTimeVal(int64_t) nothrow;
-struct timespec FileTimeToTimeSpec(struct NtFileTime);
+struct timespec WindowsTimeToTimeSpec(int64_t) nothrow;
-struct NtFileTime TimeSpecToFileTime(struct timespec);
+int64_t TimeSpecToWindowsTime(struct timespec) nothrow;
-struct NtFileTime TimeToFileTime(int64_t) nothrow pureconst;
+int64_t TimeValToWindowsTime(struct timeval) nothrow;
-int64_t filetimetotime(struct NtFileTime) nothrow pureconst;
+struct timeval WindowsDurationToTimeVal(int64_t) nothrow;
-void FileTimeToTimeVal(struct timeval *, struct NtFileTime) nothrow;
+struct timespec WindowsDurationToTimeSpec(int64_t) nothrow;
-struct NtFileTime TimeValToFileTime(const struct timeval *) nosideeffect;
-long convertmicros(const struct timeval *, long) paramsnonnull() nosideeffect;
 
-/* forceinline struct timespec WindowsTimeToTime(uint64_t x) { */
+static inline struct NtFileTime MakeFileTime(int64_t x) {
-/*   return (struct timespec){x / HECTONANOSECONDS - MODERNITYSECONDS, */
+  return (struct NtFileTime){x, x >> 32};
-/*                            x % HECTONANOSECONDS * 100}; */
+}
-/* } */
+
+static inline int64_t ReadFileTime(struct NtFileTime t) {
+  uint64_t x = t.dwHighDateTime;
+  return x << 32 | t.dwLowDateTime;
+}
+
+#define FileTimeToTimeSpec(x) WindowsTimeToTimeSpec(ReadFileTime(x))
+#define FileTimeToTimeVal(x)  WindowsTimeToTimeVal(ReadFileTime(x))
+#define TimeSpecToFileTime(x) MakeFileTime(TimeSpecToWindowsTime(x))
+#define TimeValToFileTime(x)  MakeFileTime(TimeValToWindowsTime(x))
 
 /*───────────────────────────────────────────────────────────────────────────│─╗
 │ cosmopolitan § conversion » manipulation                                 ─╬─│┼

libc/fmt/fmt.mk

@@ -58,14 +58,8 @@ $(LIBC_FMT_A_OBJS):				\
 		OVERRIDE_CFLAGS +=		\
 			-fno-jump-tables
 
-o/$(MODE)/libc/fmt/windowstimetotime.o		\
 o/$(MODE)/libc/fmt/dosdatetimetounix.o		\
-o/$(MODE)/libc/fmt/itoa64radix10.greg.o		\
+o/$(MODE)/libc/fmt/itoa64radix10.greg.o:	\
-o/$(MODE)/libc/fmt/timetofiletime.o		\
-o/$(MODE)/libc/fmt/filetimetotime.o		\
-o/$(MODE)/libc/fmt/timespectofiletime.o		\
-o/$(MODE)/libc/fmt/filetimetotimespec.o		\
-o/$(MODE)/libc/fmt/filetimetotimeval.o:		\
 		OVERRIDE_CFLAGS +=		\
 			-O3
 

libc/fmt/timespectofiletime.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
-│ Copyright 2020 Justine Alexandra Roberts Tunney                              │
-│                                                                              │
-│ Permission to use, copy, modify, and/or 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.                                                │
-╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/calls/calls.h"
-#include "libc/fmt/conv.h"
-#include "libc/nexgen32e/nexgen32e.h"
-#include "libc/nt/struct/filetime.h"
-
-/**
- * Converts UNIX nanosecond timestamp to Windows COBOL timestamp.
- */
-struct NtFileTime TimeSpecToFileTime(struct timespec ts) {
-  uint64_t x;
-  x = MODERNITYSECONDS;
-  x += ts.tv_sec * HECTONANOSECONDS;
-  x += div100int64(ts.tv_nsec);
-  return (struct NtFileTime){x, x >> 32};
-}

libc/isystem/sys/sysmacros.h

@@ -0,0 +1,4 @@
+#ifndef COSMOPOLITAN_LIBC_ISYSTEM_SYS_SYSMACROS_H_
+#define COSMOPOLITAN_LIBC_ISYSTEM_SYS_SYSMACROS_H_
+#include "libc/calls/makedev.h"
+#endif /* COSMOPOLITAN_LIBC_ISYSTEM_SYS_SYSMACROS_H_ */

libc/isystem/sys/types.h

@@ -1,5 +1,6 @@
 #ifndef LIBC_ISYSTEM_SYS_TYPES_H_
 #define LIBC_ISYSTEM_SYS_TYPES_H_
+#include "libc/calls/makedev.h"
 #include "libc/calls/weirdtypes.h"
 
 #if defined(_GNU_SOURCE) || defined(_BSD_SOURCE)
@@ -13,8 +14,8 @@ typedef unsigned u_int, uint;
 typedef unsigned long u_long, ulong;
 typedef long long quad_t;
 typedef unsigned long long u_quad_t;
-#include <endian.h>
+#include "libc/isystem/endian.h"
-#include <sys/select.h>
+#include "libc/isystem/sys/select.h"
 #endif
 
 #endif

libc/log/oldtermios.c

@@ -21,16 +21,15 @@
 #include "libc/log/internal.h"
 #include "libc/str/str.h"
 
+#define RESET_COLOR   "\e[0m"
 #define SHOW_CURSOR   "\e[?25h"
 #define DISABLE_MOUSE "\e[?1000;1002;1015;1006l"
-#define ANSI_RESTORE  SHOW_CURSOR DISABLE_MOUSE
+#define ANSI_RESTORE  RESET_COLOR SHOW_CURSOR DISABLE_MOUSE
 
 struct termios g_oldtermios;
 
 static textstartup void g_oldtermios_init() {
-  if (isatty(1)) {
+  tcgetattr(1, &g_oldtermios);
-    tcgetattr(1, &g_oldtermios);
-  }
 }
 
 const void *const g_oldtermios_ctor[] initarray = {
@@ -38,7 +37,7 @@ const void *const g_oldtermios_ctor[] initarray = {
 };
 
 void __restore_tty(void) {
-  if (isatty(1)) {
+  if (g_oldtermios.c_lflag && isatty(1)) {
     write(1, ANSI_RESTORE, strlen(ANSI_RESTORE));
     tcsetattr(1, TCSAFLUSH, &g_oldtermios);
   }

libc/nexgen32e/memrchr32.S

@@ -41,7 +41,7 @@ wmemrchr:
 	vpcmpeqd %ymm1,%ymm0,%ymm1
 	vpmovmskb %ymm1,%eax
 	lzcnt	%eax,%eax
-	shr	%eax
+	shr	$2,%eax
 	mov	%eax,%ecx
 	sub	%rcx,%rdx
 	cmp	$8,%eax

libc/runtime/clktck.c

@@ -18,6 +18,7 @@
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/calls/calls.h"
 #include "libc/dce.h"
+#include "libc/fmt/conv.h"
 #include "libc/runtime/clktck.h"
 #include "libc/sysv/consts/auxv.h"
 
@@ -36,7 +37,9 @@ static noinline int __clk_tck_init(void) {
   int cmd[2];
   size_t len;
   struct clockinfo_netbsd clock;
-  if (IsXnu() || IsOpenbsd()) {
+  if (IsWindows()) {
+    x = HECTONANOSECONDS;
+  } else if (IsXnu() || IsOpenbsd()) {
     x = 100;
   } else if (IsFreebsd()) {
     x = 128;

libc/runtime/getcpucount.c

@@ -33,6 +33,51 @@
 #define HW_NCPUONLINE_NETBSD  16
 #define ALL_PROCESSOR_GROUPS  0xffff
 
+static unsigned GetCpuCountLinux(void) {
+  uint64_t s[16];
+  unsigned i, c, n;
+  if ((n = sched_getaffinity(0, sizeof(s), s)) > 0) {
+    assert(!(n & 7));
+    for (n >>= 3, c = i = 0; i < n; ++i) {
+      c += popcnt(s[i]);
+    }
+    return c;
+  } else {
+    return 0;
+  }
+}
+
+static unsigned GetCpuCountBsd(void) {
+  size_t n;
+  int c, cmd[2];
+  n = sizeof(c);
+  cmd[0] = CTL_HW;
+  if (IsOpenbsd()) {
+    cmd[1] = HW_NCPUONLINE_OPENBSD;
+  } else if (IsNetbsd()) {
+    cmd[1] = HW_NCPUONLINE_NETBSD;
+  } else {
+    cmd[1] = HW_NCPU;
+  }
+  if (!sysctl(cmd, 2, &c, &n, 0, 0)) {
+    return c;
+  } else {
+    return 0;
+  }
+}
+
+static textwindows unsigned GetCpuCountWindows(void) {
+  struct NtSystemInfo si;
+  uint32_t (*f)(uint16_t);
+  if ((f = GetProcAddress(GetModuleHandle("KERNEL32"),
+                          "GetMaximumProcessorCount"))) {
+    return f(ALL_PROCESSOR_GROUPS);
+  } else {
+    GetSystemInfo(&si);
+    return si.dwNumberOfProcessors;
+  }
+}
+
 /**
  * Returns number of CPUs in system.
  *
@@ -42,43 +87,13 @@
  * @return cpu count or 0 if it couldn't be determined
  */
 unsigned GetCpuCount(void) {
-  size_t len;
+  if (!IsWindows()) {
-  int i, c, n, cmd[2];
+    if (!IsBsd()) {
-  uint64_t cpuset[16];
+      return GetCpuCountLinux();
-  uint32_t (*f)(uint16_t);
-  struct NtSystemInfo sysinfo;
-  if (IsWindows()) {
-    if ((f = GetProcAddress(GetModuleHandle("KERNEL32"),
-                            "GetMaximumProcessorCount"))) {
-      return f(ALL_PROCESSOR_GROUPS);
     } else {
-      GetSystemInfo(&sysinfo);
+      return GetCpuCountBsd();
-      return sysinfo.dwNumberOfProcessors;
-    }
-  } else if (IsBsd()) {
-    len = sizeof(c);
-    cmd[0] = CTL_HW;
-    if (IsOpenbsd()) {
-      cmd[1] = HW_NCPUONLINE_OPENBSD;
-    } else if (IsNetbsd()) {
-      cmd[1] = HW_NCPUONLINE_NETBSD;
-    } else {
-      cmd[1] = HW_NCPU;
-    }
-    if (!sysctl(cmd, 2, &c, &len, 0, 0)) {
-      return c;
-    } else {
-      return 0;
     }
   } else {
-    if ((n = sched_getaffinity(0, sizeof(cpuset), cpuset)) > 0) {
+    return GetCpuCountWindows();
-      assert(!(n & 7));
-      for (n >>= 3, c = i = 0; i < ARRAYLEN(cpuset); ++i) {
-        c += popcnt(cpuset[i]);
-      }
-      return c;
-    } else {
-      return 0;
-    }
   }
 }

libc/str/getzipcfiletimestamps.c

@@ -46,15 +46,15 @@ void GetZipCfileTimestamps(const uint8_t *cf, struct timespec *mtim,
         READ16LE(ZIP_EXTRA_CONTENT(p) + 4) == 1 &&
         READ16LE(ZIP_EXTRA_CONTENT(p) + 6) >= 8) {
       if (mtim) {
-        *mtim = WindowsTimeToTime(READ64LE(ZIP_EXTRA_CONTENT(p) + 8));
+        *mtim = WindowsTimeToTimeSpec(READ64LE(ZIP_EXTRA_CONTENT(p) + 8));
       }
       if (atim && ZIP_EXTRA_CONTENTSIZE(p) >= 4 + 4 + 8 * 2 &&
           READ16LE(ZIP_EXTRA_CONTENT(p) + 6) >= 16) {
-        *atim = WindowsTimeToTime(READ64LE(ZIP_EXTRA_CONTENT(p) + 8 * 2));
+        *atim = WindowsTimeToTimeSpec(READ64LE(ZIP_EXTRA_CONTENT(p) + 8 * 2));
       }
       if (ctim && ZIP_EXTRA_CONTENTSIZE(p) >= 4 + 4 + 8 * 3 &&
           READ16LE(ZIP_EXTRA_CONTENT(p) + 6) >= 24) {
-        *ctim = WindowsTimeToTime(READ64LE(ZIP_EXTRA_CONTENT(p) + 8 * 3));
+        *ctim = WindowsTimeToTimeSpec(READ64LE(ZIP_EXTRA_CONTENT(p) + 8 * 3));
       }
       return;
     }

libc/str/str.mk

@@ -77,7 +77,12 @@ o/$(MODE)/libc/str/iswlower.o:					\
 		OVERRIDE_CFLAGS +=				\
 			-fno-jump-tables
 
-o/$(MODE)/libc/str/windowstimetotime.o:				\
+o/$(MODE)/libc/str/windowsdurationtotimeval.o			\
+o/$(MODE)/libc/str/windowsdurationtotimespec.o			\
+o/$(MODE)/libc/str/timevaltowindowstime.o			\
+o/$(MODE)/libc/str/timespectowindowstime.o			\
+o/$(MODE)/libc/str/windowstimetotimeval.o			\
+o/$(MODE)/libc/str/windowstimetotimespec.o:			\
 		OVERRIDE_CFLAGS +=				\
 			-O3
 

libc/str/strcasecmp.c

@@ -18,20 +18,44 @@
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/str/str.h"
 
+static inline noasan uint64_t UncheckedAlignedRead64(const char *p) {
+  return (uint64_t)(255 & p[7]) << 070 | (uint64_t)(255 & p[6]) << 060 |
+         (uint64_t)(255 & p[5]) << 050 | (uint64_t)(255 & p[4]) << 040 |
+         (uint64_t)(255 & p[3]) << 030 | (uint64_t)(255 & p[2]) << 020 |
+         (uint64_t)(255 & p[1]) << 010 | (uint64_t)(255 & p[0]) << 000;
+}
+
 /**
  * Compares NUL-terminated strings case-insensitively.
  *
- * @param a is first non-null NUL-terminated string pointer
+ * @param a is first non-null nul-terminated string pointer
- * @param b is second non-null NUL-terminated string pointer
+ * @param b is second non-null nul-terminated string pointer
- * @return is <0, 0, or >0 based on uint8_t comparison
+ * @return is <0, 0, or >0 based on tolower(uint8_t) comparison
  * @asyncsignalsafe
  */
 int strcasecmp(const char *a, const char *b) {
   int x, y;
   size_t i = 0;
+  uint64_t v, w, d;
   if (a == b) return 0;
-  while ((x = kToLower[a[i] & 0xff]) == (y = kToLower[b[i] & 0xff]) && b[i]) {
+  if (((uintptr_t)a & 7) == ((uintptr_t)b & 7)) {
-    ++i;
+    for (; (uintptr_t)(a + i) & 7; ++i) {
+    CheckEm:
+      if ((x = kToLower[a[i] & 255]) != (y = kToLower[b[i] & 255]) || !y) {
+        return x - y;
+      }
+    }
+    for (;; i += 8) {
+      v = UncheckedAlignedRead64(a + i);
+      w = UncheckedAlignedRead64(b + i);
+      w = (v ^ w) | (~v & (v - 0x0101010101010101) & 0x8080808080808080);
+      if (w) {
+        i += (unsigned)__builtin_ctzll(w) >> 3;
+        goto CheckEm;
+      }
+    }
+  } else {
+    while ((x = kToLower[a[i] & 255]) == (y = kToLower[b[i] & 255]) && y) ++i;
+    return x - y;
   }
-  return x - y;
 }

libc/str/strcat.c

@@ -27,5 +27,6 @@
  * @asyncsignalsafe
  */
 char *strcat(char *d, const char *s) {
-  return strcpy(d + strlen(d), s);
+  strcpy(d + strlen(d), s);
+  return d;
 }

libc/str/strcat16.c

@@ -27,5 +27,6 @@
  * @asyncsignalsafe
  */
 char16_t *strcat16(char16_t *d, const char16_t *s) {
-  return strcpy16(d + strlen16(d), s);
+  strcpy16(d + strlen16(d), s);
+  return d;
 }

libc/str/strcmp.c

@@ -38,7 +38,7 @@ int strcmp(const char *a, const char *b) {
   uint64_t v, w, d;
   if (a == b) return 0;
   if (((uintptr_t)a & 7) == ((uintptr_t)b & 7)) {
-    for (; (uintptr_t)a & 7; ++i) {
+    for (; (uintptr_t)(a + i) & 7; ++i) {
       if (a[i] != b[i] || !b[i]) {
         return (a[i] & 0xff) - (b[i] & 0xff);
       }

libc/str/strncmp.c

@@ -29,6 +29,6 @@
 int strncmp(const char *a, const char *b, size_t n) {
   size_t i = 0;
   if (!n-- || a == b) return 0;
-  while (a[i] == b[i] && b[i] && i < n) ++i;
+  while (i < n && a[i] == b[i] && b[i]) ++i;
   return (a[i] & 0xff) - (b[i] & 0xff);
 }

libc/str/timespectowindowstime.c

@@ -18,7 +18,6 @@
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/fmt/conv.h"
 
-struct timespec WindowsTimeToTime(uint64_t x) {
+int64_t TimeSpecToWindowsTime(struct timespec t) {
-  return (struct timespec){x / HECTONANOSECONDS - MODERNITYSECONDS,
+  return t.tv_nsec / 100 + (t.tv_sec + MODERNITYSECONDS) * HECTONANOSECONDS;
-                           x % HECTONANOSECONDS * 100};
 }

libc/str/timevaltowindowstime.c

@@ -1,7 +1,7 @@
 /*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
 │vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
 ╞══════════════════════════════════════════════════════════════════════════════╡
-│ Copyright 2020 Justine Alexandra Roberts Tunney                              │
+│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
 │ Permission to use, copy, modify, and/or distribute this software for         │
 │ any purpose with or without fee is hereby granted, provided that the         │
@@ -16,10 +16,8 @@
 │ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR             │
 │ PERFORMANCE OF THIS SOFTWARE.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/calls/struct/timeval.h"
 #include "libc/fmt/conv.h"
-#include "libc/time/time.h"
 
-long convertmicros(const struct timeval *tv, long tick) {
+int64_t TimeValToWindowsTime(struct timeval t) {
-  return tv->tv_sec * tick + tv->tv_usec / (1000000 / tick);
+  return t.tv_usec * 10 + (t.tv_sec + MODERNITYSECONDS) * HECTONANOSECONDS;
 }

libc/str/wcscat.c

@@ -27,5 +27,6 @@
  * @asyncsignalsafe
  */
 wchar_t *wcscat(wchar_t *d, const wchar_t *s) {
-  return wcscpy(d + wcslen(d), s);
+  wcscpy(d + wcslen(d), s);
+  return d;
 }

libc/str/wcscpy.c

@@ -29,5 +29,6 @@
  * @asyncsignalsafe
  */
 wchar_t *wcscpy(wchar_t *d, const wchar_t *s) {
-  return memcpy(d, s, (wcslen(s) + 1) * sizeof(wchar_t));
+  memcpy(d, s, (wcslen(s) + 1) * sizeof(wchar_t));
+  return d;
 }

libc/str/windowsdurationtotimespec.c

@@ -1,7 +1,7 @@
 /*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
 │vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
 ╞══════════════════════════════════════════════════════════════════════════════╡
-│ Copyright 2020 Justine Alexandra Roberts Tunney                              │
+│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
 │ Permission to use, copy, modify, and/or distribute this software for         │
 │ any purpose with or without fee is hereby granted, provided that the         │
@@ -17,9 +17,7 @@
 │ PERFORMANCE OF THIS SOFTWARE.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/fmt/conv.h"
-#include "libc/nt/struct/filetime.h"
 
-struct NtFileTime TimeToFileTime(int64_t t) {
+struct timespec WindowsDurationToTimeSpec(int64_t x) {
-  uint64_t t2 = (t + MODERNITYSECONDS) * HECTONANOSECONDS;
+  return (struct timespec){x / HECTONANOSECONDS, x % HECTONANOSECONDS * 100};
-  return (struct NtFileTime){(uint32_t)t2, (uint32_t)(t2 >> 32)};
 }

libc/str/windowsdurationtotimeval.c

@@ -1,7 +1,7 @@
 /*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
 │vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
 ╞══════════════════════════════════════════════════════════════════════════════╡
-│ Copyright 2020 Justine Alexandra Roberts Tunney                              │
+│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
 │ Permission to use, copy, modify, and/or distribute this software for         │
 │ any purpose with or without fee is hereby granted, provided that the         │
@@ -17,9 +17,7 @@
 │ PERFORMANCE OF THIS SOFTWARE.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/fmt/conv.h"
-#include "libc/nt/struct/filetime.h"
 
-textwindows int64_t filetimetotime(struct NtFileTime ft) {
+struct timeval WindowsDurationToTimeVal(int64_t x) {
-  uint64_t t = (uint64_t)ft.dwHighDateTime << 32 | ft.dwLowDateTime;
+  return (struct timeval){x / HECTONANOSECONDS, x % HECTONANOSECONDS / 10};
-  return (t - MODERNITYSECONDS * HECTONANOSECONDS) / HECTONANOSECONDS;
 }

libc/str/windowstimetotimespec.c

@@ -0,0 +1,26 @@
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+╞══════════════════════════════════════════════════════════════════════════════╡
+│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
+│                                                                              │
+│ Permission to use, copy, modify, and/or 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.                                                │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/fmt/conv.h"
+
+/**
+ * Converts Windows COBOL timestamp to UNIX epoch in nanoseconds.
+ */
+struct timespec WindowsTimeToTimeSpec(int64_t x) {
+  return WindowsDurationToTimeSpec(x - MODERNITYSECONDS * HECTONANOSECONDS);
+}

libc/str/windowstimetotimeval.c

@@ -0,0 +1,29 @@
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+╞══════════════════════════════════════════════════════════════════════════════╡
+│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
+│                                                                              │
+│ Permission to use, copy, modify, and/or 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.                                                │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/fmt/conv.h"
+
+/**
+ * Converts Windows COBOL timestamp to UNIX epoch in microseconds.
+ */
+struct timeval WindowsTimeToTimeVal(int64_t x) {
+  /* return WindowsDurationToTimeVal(x - MODERNITYSECONDS * HECTONANOSECONDS);
+   */
+  return (struct timeval){x / HECTONANOSECONDS - MODERNITYSECONDS,
+                          x % HECTONANOSECONDS / 10};
+}

libc/sysv/calls/getpgid.s

@@ -1,2 +0,0 @@
-.include "o/libc/sysv/macros.internal.inc"
-.scall getpgid,0x0cf0cf0cf2097079,globl

libc/sysv/calls/sys_getpgid.s

@@ -0,0 +1,2 @@
+.include "o/libc/sysv/macros.internal.inc"
+.scall sys_getpgid,0x0cf0cf0cf2097079,globl,hidden

libc/sysv/consts.sh

@@ -352,23 +352,23 @@ syscon	stat	S_IRWXO					0000007			0000007			0000007			0000007			0000007			000000
 #
 #	group	name					GNU/Systemd		XNU's Not UNIX!		FreeBSD			OpenBSD			NetBSD			The New Technology	Commentary
 syscon	fcntl2	F_DUPFD					0			0			0			0			0			0			# consensus
-
 syscon	fcntl2	F_GETFD					1			1			1			1			1			1			# unix consensus & faked nt
 syscon	fcntl2	F_SETFD					2			2			2			2			2			2			# unix consensus & faked nt
+syscon	fcntl2	F_GETFL					3			3			3			3			3			3			# unix consensus & faked nt
+syscon	fcntl2	F_SETFL					4			4			4			4			4			4			# unix consensus & faked nt
+syscon	fcntl2	F_SETOWN				8			6			6			6			6			0			# bsd consensus
+syscon	fcntl2	F_GETOWN				9			5			5			5			5			0			# bsd consensus
+syscon	fcntl2	F_FULLFSYNC				0			51			0			0			0			0			#
+syscon	fcntl2	F_NOCACHE				0			48			0			0			0			0			#
 syscon	fcntl3	FD_CLOEXEC				1			1			1			1			1			1			# unix consensus & faked nt
 syscon	fcntl	F_DUPFD_CLOEXEC				0x0406			67			17			10			12			0x0406			# faked nt
 
-syscon	fcntl2	F_GETFL					3			3			3			3			3			3			# unix consensus & faked nt
-syscon	fcntl2	F_SETFL					4			4			4			4			4			4			# unix consensus & faked nt
 #       fcntl3	O_NONBLOCK
 #       fcntl3	O_APPEND
 #       fcntl3	O_ASYNC
 #       fcntl3	O_DIRECT
 #       fcntl3	O_NOATIME
 
-syscon	fcntl2	F_SETOWN				8			6			6			6			6			0			# bsd consensus
-syscon	fcntl2	F_GETOWN				9			5			5			5			5			0			# bsd consensus
-
 #	fcntl() POSIX Advisory Locks
 #
 #	group	name					GNU/Systemd		XNU's Not UNIX!		FreeBSD			OpenBSD			NetBSD			The New Technology	Commentary

libc/sysv/consts/F_FULLFSYNC.S

@@ -0,0 +1,2 @@
+#include "libc/sysv/consts/syscon.internal.h"
+.syscon fcntl2,F_FULLFSYNC,0,51,0,0,0,0

libc/sysv/consts/F_NOCACHE.S

@@ -0,0 +1,2 @@
+#include "libc/sysv/consts/syscon.internal.h"
+.syscon fcntl2,F_NOCACHE,0,48,0,0,0,0

libc/sysv/consts/f.h

@@ -6,6 +6,7 @@ COSMOPOLITAN_C_START_
 
 extern const long F_DUPFD;
 extern const long F_DUPFD_CLOEXEC;
+extern const long F_FULLFSYNC;
 extern const long F_GETFD;
 extern const long F_GETFL;
 extern const long F_GETLEASE;
@@ -16,6 +17,7 @@ extern const long F_GETOWN_EX;
 extern const long F_GETPIPE_SZ;
 extern const long F_GETSIG;
 extern const long F_LOCK;
+extern const long F_NOCACHE;
 extern const long F_NOTIFY;
 extern const long F_OFD_GETLK;
 extern const long F_OFD_SETLK;
@@ -41,39 +43,41 @@ extern const long F_WRLCK;
 COSMOPOLITAN_C_END_
 #endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */
 
+#define F_DUPFD LITERALLY(0)
 #define F_GETFD LITERALLY(1)
 #define F_SETFD LITERALLY(2)
 #define F_GETFL LITERALLY(3)
 #define F_SETFL LITERALLY(4)
 
-#define F_DUPFD SYMBOLIC(F_DUPFD)
 #define F_DUPFD_CLOEXEC SYMBOLIC(F_DUPFD_CLOEXEC)
-#define F_GETLEASE SYMBOLIC(F_GETLEASE)
+#define F_FULLFSYNC     SYMBOLIC(F_FULLFSYNC)
-#define F_GETLK SYMBOLIC(F_GETLK)
+#define F_GETLEASE      SYMBOLIC(F_GETLEASE)
-#define F_GETLK64 SYMBOLIC(F_GETLK64)
+#define F_GETLK         SYMBOLIC(F_GETLK)
-#define F_GETOWN SYMBOLIC(F_GETOWN)
+#define F_GETLK64       SYMBOLIC(F_GETLK64)
-#define F_GETOWN_EX SYMBOLIC(F_GETOWN_EX)
+#define F_GETOWN        SYMBOLIC(F_GETOWN)
-#define F_GETPIPE_SZ SYMBOLIC(F_GETPIPE_SZ)
+#define F_GETOWN_EX     SYMBOLIC(F_GETOWN_EX)
-#define F_GETSIG SYMBOLIC(F_GETSIG)
+#define F_GETPIPE_SZ    SYMBOLIC(F_GETPIPE_SZ)
-#define F_LOCK SYMBOLIC(F_LOCK)
+#define F_GETSIG        SYMBOLIC(F_GETSIG)
-#define F_NOTIFY SYMBOLIC(F_NOTIFY)
+#define F_LOCK          SYMBOLIC(F_LOCK)
-#define F_OFD_GETLK SYMBOLIC(F_OFD_GETLK)
+#define F_NOCACHE       SYMBOLIC(F_NOCACHE)
-#define F_OFD_SETLK SYMBOLIC(F_OFD_SETLK)
+#define F_NOTIFY        SYMBOLIC(F_NOTIFY)
-#define F_OFD_SETLKW SYMBOLIC(F_OFD_SETLKW)
+#define F_OFD_GETLK     SYMBOLIC(F_OFD_GETLK)
-#define F_RDLCK SYMBOLIC(F_RDLCK)
+#define F_OFD_SETLK     SYMBOLIC(F_OFD_SETLK)
-#define F_SETLEASE SYMBOLIC(F_SETLEASE)
+#define F_OFD_SETLKW    SYMBOLIC(F_OFD_SETLKW)
-#define F_SETLK SYMBOLIC(F_SETLK)
+#define F_RDLCK         SYMBOLIC(F_RDLCK)
-#define F_SETLK64 SYMBOLIC(F_SETLK64)
+#define F_SETLEASE      SYMBOLIC(F_SETLEASE)
-#define F_SETLKW SYMBOLIC(F_SETLKW)
+#define F_SETLK         SYMBOLIC(F_SETLK)
-#define F_SETLKW64 SYMBOLIC(F_SETLKW64)
+#define F_SETLK64       SYMBOLIC(F_SETLK64)
-#define F_SETOWN SYMBOLIC(F_SETOWN)
+#define F_SETLKW        SYMBOLIC(F_SETLKW)
-#define F_SETOWN_EX SYMBOLIC(F_SETOWN_EX)
+#define F_SETLKW64      SYMBOLIC(F_SETLKW64)
-#define F_SETPIPE_SZ SYMBOLIC(F_SETPIPE_SZ)
+#define F_SETOWN        SYMBOLIC(F_SETOWN)
-#define F_SETSIG SYMBOLIC(F_SETSIG)
+#define F_SETOWN_EX     SYMBOLIC(F_SETOWN_EX)
-#define F_TEST SYMBOLIC(F_TEST)
+#define F_SETPIPE_SZ    SYMBOLIC(F_SETPIPE_SZ)
-#define F_TLOCK SYMBOLIC(F_TLOCK)
+#define F_SETSIG        SYMBOLIC(F_SETSIG)
-#define F_ULOCK SYMBOLIC(F_ULOCK)
+#define F_TEST          SYMBOLIC(F_TEST)
-#define F_UNLCK SYMBOLIC(F_UNLCK)
+#define F_TLOCK         SYMBOLIC(F_TLOCK)
-#define F_WRLCK SYMBOLIC(F_WRLCK)
+#define F_ULOCK         SYMBOLIC(F_ULOCK)
+#define F_UNLCK         SYMBOLIC(F_UNLCK)
+#define F_WRLCK         SYMBOLIC(F_WRLCK)
 
 #endif /* COSMOPOLITAN_LIBC_SYSV_CONSTS_F_H_ */

libc/sysv/syscalls.sh

@@ -148,7 +148,7 @@ scall	sys_getppid		0xfff027027202706e	globl hidden # see sys_getpid()→edx for
 scall	getpgrp			0x051051051205106f	globl
 scall	sys_setsid		0x0930930932093070	globl hidden
 scall	sys_getsid		0x11e0ff136213607c	globl hidden
-scall	getpgid			0x0cf0cf0cf2097079	globl
+scall	sys_getpgid		0x0cf0cf0cf2097079	globl hidden
 scall	setpgid			0x052052052205206d	globl
 scall	geteuid			0xfff019019201906b	globl
 scall	getegid			0xfff02b02b202b06c	globl

libc/testlib/ezbenchreport.c

@@ -27,7 +27,6 @@ void __testlib_ezbenchreport(const char *form, uint64_t c1, uint64_t c2) {
   uint64_t ns1, ns2;
   ns1 = rintl(ConvertTicksToNanos(c1));
   ns2 = rintl(ConvertTicksToNanos(c2));
-  (fprintf)(stderr,
+  (fprintf)(stderr, VEIL("r", "%-26s l: %,9lu𝑐 %,9lu𝑛𝑠   m: %,9lu𝑐 %,9lu𝑛𝑠\n"),
-            VEIL("r", "%-30s l: %,10lu𝑐 %,10lu𝑛𝑠   m: %,10lu𝑐 %,10lu𝑛𝑠\n"),
             form, c1, ns1, c2, ns2);
 }

libc/time/times.c

@@ -25,36 +25,37 @@
 #include "libc/fmt/conv.h"
 #include "libc/nt/accounting.h"
 #include "libc/nt/runtime.h"
+#include "libc/runtime/clktck.h"
 #include "libc/runtime/sysconf.h"
 #include "libc/sysv/consts/rusage.h"
 #include "libc/time/time.h"
 
+static noinline long ConvertMicros(struct timeval tv) {
+  return tv.tv_sec * CLK_TCK + tv.tv_usec / (1000000 / CLK_TCK);
+}
+
 static noinline long times2(struct tms *out_times, struct rusage *ru) {
-  long tick;
   struct timeval tv;
-  tick = sysconf(_SC_CLK_TCK);
+  struct NtFileTime CreationTime, ExitTime, KernelTime, UserTime;
   if (!IsWindows()) {
     if (getrusage(RUSAGE_SELF, ru) == -1) return -1;
-    out_times->tms_utime = convertmicros(&ru->ru_utime, tick);
+    out_times->tms_utime = ConvertMicros(ru->ru_utime);
-    out_times->tms_stime = convertmicros(&ru->ru_stime, tick);
+    out_times->tms_stime = ConvertMicros(ru->ru_stime);
     if (getrusage(RUSAGE_CHILDREN, ru) == -1) return -1;
-    out_times->tms_cutime = convertmicros(&ru->ru_utime, tick);
+    out_times->tms_cutime = ConvertMicros(ru->ru_utime);
-    out_times->tms_cstime = convertmicros(&ru->ru_stime, tick);
+    out_times->tms_cstime = ConvertMicros(ru->ru_stime);
   } else {
-    struct NtFileTime CreationTime, ExitTime, KernelTime, UserTime;
     if (!GetProcessTimes(GetCurrentProcess(), &CreationTime, &ExitTime,
                          &KernelTime, &UserTime)) {
       return __winerr();
     }
-    FileTimeToTimeVal(&tv, UserTime);
+    out_times->tms_utime = ReadFileTime(UserTime);
-    out_times->tms_utime = convertmicros(&tv, tick);
+    out_times->tms_stime = ReadFileTime(KernelTime);
-    FileTimeToTimeVal(&tv, KernelTime);
-    out_times->tms_stime = convertmicros(&tv, tick);
     out_times->tms_cutime = 0;
     out_times->tms_cstime = 0;
   }
   if (gettimeofday(&tv, NULL) == -1) return -1;
-  return convertmicros(&tv, tick);
+  return ConvertMicros(tv);
 }
 
 /**

libc/tinymath/acosh.c

@@ -1,26 +1,55 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
 
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
+
 /**
  * Returns inverse hyperbolic cosine of 𝑥.
+ * @define acosh(x) = log(x + sqrt(x*x-1))
  */
-double acosh(double x) {
+double acosh(double x)
-  return log(x + sqrt(x * x - 1));
+{
+	union {double f; uint64_t i;} u = {.f = x};
+	unsigned e = u.i >> 52 & 0x7ff;
+
+	/* x < 1 domain error is handled in the called functions */
+
+	if (e < 0x3ff + 1)
+		/* |x| < 2, up to 2ulp error in [1,1.125] */
+		return log1p(x-1 + sqrt((x-1)*(x-1)+2*(x-1)));
+	if (e < 0x3ff + 26)
+		/* |x| < 0x1p26 */
+		return log(2*x - 1/(x+sqrt(x*x-1)));
+	/* |x| >= 0x1p26 or nan */
+	return log(x) + 0.693147180559945309417232121458176568;
 }

libc/tinymath/acoshf.c

@@ -1,26 +1,54 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
 
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
+
 /**
  * Returns inverse hyperbolic cosine of 𝑥.
+ * @define acosh(x) = log(x + sqrt(x*x-1))
  */
-float acoshf(float x) {
+float acoshf(float x)
-  return logf(x + sqrtf(x * x - 1));
+{
+	union {float f; uint32_t i;} u = {x};
+	uint32_t a = u.i & 0x7fffffff;
+
+	if (a < 0x3f800000+(1<<23))
+		/* |x| < 2, invalid if x < 1 */
+		/* up to 2ulp error in [1,1.125] */
+		return log1pf(x-1 + sqrtf((x-1)*(x-1)+2*(x-1)));
+	if (u.i < 0x3f800000+(12<<23))
+		/* 2 <= x < 0x1p12 */
+		return logf(2*x - 1/(x+sqrtf(x*x-1)));
+	/* x >= 0x1p12 or x <= -2 or nan */
+	return logf(x) + 0.693147180559945309417232121458176568f;
 }

libc/tinymath/acoshl.c

@@ -1,26 +1,52 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/ldshape.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns inverse hyperbolic cosine of 𝑥.
  */
-long double acoshl(long double x) {
+long double acoshl(long double x)
-  return logl(x + sqrtl(x * x - 1));
+{
+	union ldshape u = {x};
+	int e = u.i.se & 0x7fff;
+
+	if (e < 0x3fff + 1)
+		/* |x| < 2, invalid if x < 1 or nan */
+		return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1)));
+	if (e < 0x3fff + 32)
+		/* |x| < 0x1p32 */
+		return logl(2*x - 1/(x+sqrtl(x*x-1)));
+	return logl(x) + 0.693147180559945309417232121458176568L;
 }

libc/tinymath/asinh.c

@@ -1,26 +1,65 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns inverse hyperbolic sine of 𝑥.
+ * @define asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5)
  */
-double asinh(double x) {
+double asinh(double x)
-  return copysign(log(fabs(x) + sqrt(x * x + 1)), x);
+{
+	union {double f; uint64_t i;} u = {.f = x};
+	unsigned e = u.i >> 52 & 0x7ff;
+	unsigned s = u.i >> 63;
+
+	/* |x| */
+	u.i &= (uint64_t)-1/2;
+	x = u.f;
+
+	if (e >= 0x3ff + 26) {
+		/* |x| >= 0x1p26 or inf or nan */
+		x = log(x) + 0.693147180559945309417232121458176568;
+	} else if (e >= 0x3ff + 1) {
+		/* |x| >= 2 */
+		x = log(2*x + 1/(sqrt(x*x+1)+x));
+	} else if (e >= 0x3ff - 26) {
+		/* |x| >= 0x1p-26, up to 1.6ulp error in [0.125,0.5] */
+		x = log1p(x + x*x/(sqrt(x*x+1)+1));
+	} else {
+		/* |x| < 0x1p-26, raise inexact if x != 0 */
+		feval(x + 0x1p120f);
+	}
+	return s ? -x : x;
 }

libc/tinymath/asinhf.c

@@ -1,26 +1,65 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns inverse hyperbolic sine of 𝑥.
+ * @define asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5)
  */
-float asinhf(float x) {
+float asinhf(float x)
-  return copysignf(logf(fabsf(x) + sqrtf(x * x + 1)), x);
+{
+	union {float f; uint32_t i;} u = {.f = x};
+	uint32_t i = u.i & 0x7fffffff;
+	unsigned s = u.i >> 31;
+
+	/* |x| */
+	u.i = i;
+	x = u.f;
+
+	if (i >= 0x3f800000 + (12<<23)) {
+		/* |x| >= 0x1p12 or inf or nan */
+		x = logf(x) + 0.693147180559945309417232121458176568f;
+	} else if (i >= 0x3f800000 + (1<<23)) {
+		/* |x| >= 2 */
+		x = logf(2*x + 1/(sqrtf(x*x+1)+x));
+	} else if (i >= 0x3f800000 - (12<<23)) {
+		/* |x| >= 0x1p-12, up to 1.6ulp error in [0.125,0.5] */
+		x = log1pf(x + x*x/(sqrtf(x*x+1)+1));
+	} else {
+		/* |x| < 0x1p-12, raise inexact if x!=0 */
+		feval(x + 0x1p120f);
+	}
+	return s ? -x : x;
 }

libc/tinymath/asinhl.c

@@ -1,26 +1,66 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
+#include "libc/tinymath/ldshape.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns inverse hyperbolic sine of 𝑥.
+ * @define asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5)
  */
-long double asinhl(long double x) {
+long double asinhl(long double x)
-  return copysignl(logl(fabsl(x) + sqrtl(x * x + 1)), x);
+{
+	union ldshape u = {x};
+	unsigned e = u.i.se & 0x7fff;
+	unsigned s = u.i.se >> 15;
+
+	/* |x| */
+	u.i.se = e;
+	x = u.f;
+
+	if (e >= 0x3fff + 32) {
+		/* |x| >= 0x1p32 or inf or nan */
+		x = logl(x) + 0.693147180559945309417232121458176568L;
+	} else if (e >= 0x3fff + 1) {
+		/* |x| >= 2 */
+		x = logl(2*x + 1/(sqrtl(x*x+1)+x));
+	} else if (e >= 0x3fff - 32) {
+		/* |x| >= 0x1p-32 */
+		x = log1pl(x + x*x/(sqrtl(x*x+1)+1));
+	} else {
+		/* |x| < 0x1p-32, raise inexact if x!=0 */
+		fevall(x + 0x1p120f);
+	}
+	return s ? -x : x;
 }

libc/tinymath/atanh.c

@@ -1,26 +1,66 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns inverse hyperbolic tangent of 𝑥.
+ * @define x ? log1p(2 * x / (1 - x)) / 2 : x
  */
-double atanh(double x) {
+double atanh(double x)
-  return x ? log((1 + x) / (1 - x)) / 2 : x;
+{
+	union {double f; uint64_t i;} u = {.f = x};
+	unsigned e = u.i >> 52 & 0x7ff;
+	unsigned s = u.i >> 63;
+	double_t y;
+
+	/* |x| */
+	u.i &= (uint64_t)-1/2;
+	y = u.f;
+
+	if (e < 0x3ff - 1) {
+		if (e < 0x3ff - 32) {
+			/* handle underflow */
+			if (e == 0)
+				fevalf(y);
+		} else {
+			/* |x| < 0.5, up to 1.7ulp error */
+			y = 0.5*log1p(2*y + 2*y*y/(1-y));
+		}
+	} else {
+		/* avoid overflow */
+		y = 0.5*log1p(2*(y/(1-y)));
+	}
+	return s ? -y : y;
 }

libc/tinymath/atanhf.c

@@ -1,26 +1,65 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns inverse hyperbolic tangent of 𝑥.
+ * @define x ? log1p(2 * x / (1 - x)) / 2 : x
  */
-float atanhf(float x) {
+float atanhf(float x)
-  return x ? logf((1 + x) / (1 - x)) / 2 : x;
+{
+	union {float f; uint32_t i;} u = {.f = x};
+	unsigned s = u.i >> 31;
+	float_t y;
+
+	/* |x| */
+	u.i &= 0x7fffffff;
+	y = u.f;
+
+	if (u.i < 0x3f800000 - (1<<23)) {
+		if (u.i < 0x3f800000 - (32<<23)) {
+			/* handle underflow */
+			if (u.i < (1<<23))
+				fevalf(y*y);
+		} else {
+			/* |x| < 0.5, up to 1.7ulp error */
+			y = 0.5f*log1pf(2*y + 2*y*y/(1-y));
+		}
+	} else {
+		/* avoid overflow */
+		y = 0.5f*log1pf(2*(y/(1-y)));
+	}
+	return s ? -y : y;
 }

libc/tinymath/atanhl.c

@@ -1,26 +1,66 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
+#include "libc/tinymath/ldshape.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns inverse hyperbolic tangent of 𝑥.
+ * @define x ? log1p(2 * x / (1 - x)) / 2 : x
  */
-long double atanhl(long double x) {
+long double atanhl(long double x)
-  return x ? logl((1 + x) / (1 - x)) / 2 : x;
+{
+	union ldshape u = {x};
+	unsigned e = u.i.se & 0x7fff;
+	unsigned s = u.i.se >> 15;
+
+	/* |x| */
+	u.i.se = e;
+	x = u.f;
+
+	if (e < 0x3ff - 1) {
+		if (e < 0x3ff - LDBL_MANT_DIG/2) {
+			/* handle underflow */
+			if (e == 0)
+				fevalf(x);
+		} else {
+			/* |x| < 0.5, up to 1.7ulp error */
+			x = 0.5*log1pl(2*x + 2*x*x/(1-x));
+		}
+	} else {
+		/* avoid overflow */
+		x = 0.5*log1pl(2*(x/(1-x)));
+	}
+	return s ? -x : x;
 }

libc/tinymath/cos.S

@@ -1,28 +0,0 @@
-/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8     -*-│
-│vi: set et ft=asm ts=8 tw=8 fenc=utf-8                                     :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
-│ Copyright 2020 Justine Alexandra Roberts Tunney                              │
-│                                                                              │
-│ Permission to use, copy, modify, and/or 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.                                                │
-╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/macros.internal.h"
-
-//	Returns cosine of 𝑥.
-//
-//	@param	𝑥 is double scalar in low half of %xmm0
-//	@return	double scalar in low half of %xmm0
-//	@domain	-(3π/8) < 𝑥 < 3π/8 for best accuracy
-cos:	ezlea	cosl,ax
-	jmp	_d2ld2
-	.endfn	cos,globl

libc/tinymath/cos.c

@@ -0,0 +1,118 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│  Musl Libc                                                                   │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
+#include "libc/tinymath/kernel.internal.h"
+
+asm(".ident\t\"\\n\\n\
+fdlibm (fdlibm license)\\n\
+Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\"");
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+
+/* clang-format off */
+/* origin: FreeBSD /usr/src/lib/msun/src/s_cos.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/* cos(x)
+ * Return cosine function of x.
+ *
+ * kernel function:
+ *      __sin           ... sine function on [-pi/4,pi/4]
+ *      __cos           ... cosine function on [-pi/4,pi/4]
+ *      __rem_pio2      ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *      in [-pi/4 , +pi/4], and let n = k mod 4.
+ *      We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *          0          S           C             T
+ *          1          C          -S            -1/T
+ *          2         -S          -C             T
+ *          3         -C           S            -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *      TRIG(x) returns trig(x) nearly rounded
+ */
+
+#define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i
+#define gethighw(hi,d) (hi) = asuint64(d) >> 32
+
+double cos(double x)
+{
+	double y[2];
+	uint32_t ix;
+	unsigned n;
+
+	gethighw(ix, x);
+	ix &= 0x7fffffff;
+
+	/* |x| ~< pi/4 */
+	if (ix <= 0x3fe921fb) {
+		if (ix < 0x3e46a09e) {  /* |x| < 2**-27 * sqrt(2) */
+			/* raise inexact if x!=0 */
+			feval(x + 0x1p120f);
+			return 1.0;
+		}
+		return __cos(x, 0);
+	}
+
+	/* cos(Inf or NaN) is NaN */
+	if (ix >= 0x7ff00000)
+		return x-x;
+
+	/* argument reduction */
+	n = __rem_pio2(x, y);
+	switch (n&3) {
+	case 0: return  __cos(y[0], y[1]);
+	case 1: return -__sin(y[0], y[1], 1);
+	case 2: return -__cos(y[0], y[1]);
+	default:
+		return  __sin(y[0], y[1], 1);
+	}
+}

libc/tinymath/cosh.c

@@ -1,26 +1,78 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/expo.internal.h"
+#include "libc/tinymath/feval.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns hyperbolic cosine of 𝑥.
+ *
+ *     cosh(x) = (exp(x) + 1/exp(x))/2
+ *             = 1 + 0.5*(exp(x)-1)*(exp(x)-1)/exp(x)
+ *             = 1 + x*x/2 + o(x^4)
  */
-double cosh(double x) {
+double cosh(double x)
-  return (exp(x) + exp(-x)) / 2;
+{
+	union {double f; uint64_t i;} u = {.f = x};
+	uint32_t w;
+	double t;
+
+	/* |x| */
+	u.i &= (uint64_t)-1/2;
+	x = u.f;
+	w = u.i >> 32;
+
+	/* |x| < log(2) */
+	if (w < 0x3fe62e42) {
+		if (w < 0x3ff00000 - (26<<20)) {
+			/* raise inexact if x!=0 */
+			feval(x + 0x1p120f);
+			return 1;
+		}
+		t = expm1(x);
+		return 1 + t*t/(2*(1+t));
+	}
+
+	/* |x| < log(DBL_MAX) */
+	if (w < 0x40862e42) {
+		t = exp(x);
+		/* note: if x>log(0x1p26) then the 1/t is not needed */
+		return 0.5*(t + 1/t);
+	}
+
+	/* |x| > log(DBL_MAX) or nan */
+	/* note: the result is stored to handle overflow */
+	t = __expo2(x, 1.0);
+	return t;
 }

libc/tinymath/coshf.c

@@ -1,26 +1,75 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/expo.internal.h"
+#include "libc/tinymath/feval.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns hyperbolic cosine of 𝑥.
+ *
+ *     cosh(x) = (exp(x) + 1/exp(x))/2
+ *             = 1 + 0.5*(exp(x)-1)*(exp(x)-1)/exp(x)
+ *             = 1 + x*x/2 + o(x^4)
  */
-float coshf(float x) {
+float coshf(float x)
-  return (expf(x) + expf(-x)) / 2;
+{
+	union {float f; uint32_t i;} u = {.f = x};
+	uint32_t w;
+	float t;
+
+	/* |x| */
+	u.i &= 0x7fffffff;
+	x = u.f;
+	w = u.i;
+
+	/* |x| < log(2) */
+	if (w < 0x3f317217) {
+		if (w < 0x3f800000 - (12<<23)) {
+			feval(x + 0x1p120f);
+			return 1;
+		}
+		t = expm1f(x);
+		return 1 + t*t/(2*(1+t));
+	}
+
+	/* |x| < log(FLT_MAX) */
+	if (w < 0x42b17217) {
+		t = expf(x);
+		return 0.5f*(t + 1/t);
+	}
+
+	/* |x| > log(FLT_MAX) or nan */
+	t = __expo2f(x, 1.0f);
+	return t;
 }

libc/tinymath/coshl.c

@@ -1,26 +1,77 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/expo.internal.h"
+#include "libc/tinymath/feval.internal.h"
+#include "libc/tinymath/ldshape.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns hyperbolic cosine of 𝑥.
+ *
+ *     cosh(x) = (exp(x) + 1/exp(x))/2
+ *             = 1 + 0.5*(exp(x)-1)*(exp(x)-1)/exp(x)
+ *             = 1 + x*x/2 + o(x^4)
  */
-long double coshl(long double x) {
+long double coshl(long double x)
-  return (expl(x) + expl(-x)) / 2;
+{
+	union ldshape u = {x};
+	unsigned ex = u.i.se & 0x7fff;
+	uint32_t w;
+	long double t;
+
+	/* |x| */
+	u.i.se = ex;
+	x = u.f;
+	w = u.i.m >> 32;
+
+	/* |x| < log(2) */
+	if (ex < 0x3fff-1 || (ex == 0x3fff-1 && w < 0xb17217f7)) {
+		if (ex < 0x3fff-32) {
+			fevalf(x + 0x1p120f);
+			return 1;
+		}
+		t = expm1l(x);
+		return 1 + t*t/(2*(1+t));
+	}
+
+	/* |x| < log(LDBL_MAX) */
+	if (ex < 0x3fff+13 || (ex == 0x3fff+13 && w < 0xb17217f7)) {
+		t = expl(x);
+		return 0.5*(t + 1/t);
+	}
+
+	/* |x| > log(LDBL_MAX) or nan */
+	t = expl(0.5*x);
+	return 0.5*t*t;
 }

libc/tinymath/erf.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -34,8 +34,8 @@ asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2014 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
+
 /* origin: FreeBSD /usr/src/lib/msun/src/s_erf.c */
 /*
  * ====================================================

libc/tinymath/erff.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -25,7 +25,6 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
 
 asm(".ident\t\"\\n\\n\
@@ -35,8 +34,8 @@ asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2020 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
+
 /* origin: FreeBSD /usr/src/lib/msun/src/s_erff.c */
 /*
  * ====================================================

libc/tinymath/expo.internal.h

@@ -0,0 +1,11 @@
+#ifndef COSMOPOLITAN_LIBC_TINYMATH_EXPO_H_
+#define COSMOPOLITAN_LIBC_TINYMATH_EXPO_H_
+#if !(__ASSEMBLER__ + __LINKER__ + 0)
+COSMOPOLITAN_C_START_
+
+double __expo2(double, double);
+float __expo2f(float, float);
+
+COSMOPOLITAN_C_END_
+#endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */
+#endif /* COSMOPOLITAN_LIBC_TINYMATH_EXPO_H_ */

libc/tinymath/expo2.c

@@ -0,0 +1,56 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│  Musl Libc                                                                   │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/math.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
+
+#define asdouble(i) ((union{uint64_t _i; double _f;}){i})._f
+#define INSERT_WORDS(d,hi,lo)                     \
+do {                                              \
+  (d) = asdouble(((uint64_t)(hi)<<32) | (uint32_t)(lo)); \
+} while (0)
+
+/* k is such that k*ln2 has minimal relative error and x - kln2 > log(DBL_MIN) */
+static const int k = 2043;
+static const double kln2 = 0x1.62066151add8bp+10;
+
+/* exp(x)/2 for x >= log(DBL_MAX), slightly better than 0.5*exp(x/2)*exp(x/2) */
+double __expo2(double x, double sign)
+{
+	double scale;
+
+	/* note that k is odd and scale*scale overflows */
+	INSERT_WORDS(scale, (uint32_t)(0x3ff + k/2) << 20, 0);
+	/* exp(x - k ln2) * 2**(k-1) */
+	/* in directed rounding correct sign before rounding or overflow is important */
+	return exp(x - kln2) * (sign * scale) * scale;
+}

libc/tinymath/expo2f.c

@@ -0,0 +1,56 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│  Musl Libc                                                                   │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/math.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
+
+#define asfloat(i) ((union{uint32_t _i; float _f;}){i})._f
+#define SET_FLOAT_WORD(d,w)                       \
+do {                                              \
+  (d) = asfloat(w);                               \
+} while (0)
+
+/* k is such that k*ln2 has minimal relative error and x - kln2 > log(FLT_MIN) */
+static const int k = 235;
+static const float kln2 = 0x1.45c778p+7f;
+
+/* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */
+float __expo2f(float x, float sign)
+{
+	float scale;
+
+	/* note that k is odd and scale*scale overflows */
+	SET_FLOAT_WORD(scale, (uint32_t)(0x7f + k/2) << 23);
+	/* exp(x - k ln2) * 2**(k-1) */
+	/* in directed rounding correct sign before rounding or overflow is important */
+	return expf(x - kln2) * (sign * scale) * scale;
+}

libc/tinymath/fdim.c

@@ -1,45 +1,27 @@
 /*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
 │vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
-╚──────────────────────────────────────────────────────────────────────────────╝
+╞══════════════════════════════════════════════════════════════════════════════╡
+│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│  Musl Libc                                                                   │
+│ Permission to use, copy, modify, and/or distribute this software for         │
-│  Copyright © 2005-2020 Rich Felker, et al.                                   │
+│ any purpose with or without fee is hereby granted, provided that the         │
-│                                                                              │
+│ above copyright notice and this permission notice appear in all copies.      │
-│  Permission is hereby granted, free of charge, to any person obtaining       │
-│  a copy of this software and associated documentation files (the             │
-│  "Software"), to deal in the Software without restriction, including         │
-│  without limitation the rights to use, copy, modify, merge, publish,         │
-│  distribute, sublicense, and/or sell copies of the Software, and to          │
-│  permit persons to whom the Software is furnished to do so, subject to       │
-│  the following conditions:                                                   │
-│                                                                              │
-│  The above copyright notice and this permission notice shall be              │
-│  included in all copies or substantial portions of the Software.             │
-│                                                                              │
-│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
-│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
-│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
-│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
-│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
-│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
-│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
+│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
 
-asm(".ident\t\"\\n\\n\
+/**
-Musl libc (MIT License)\\n\
+ * Returns positive difference.
-Copyright 2005-2020 Rich Felker, et. al.\"");
+ */
-asm(".include \"libc/disclaimer.inc\"");
+double fdim(double x, double y) {
-
+  if (isnan(x) || isnan(y)) return NAN;
-/* clang-format off */
+  return x > y ? x - y : 0;
-
-double fdim(double x, double y)
-{
-	if (isnan(x))
-		return x;
-	if (isnan(y))
-		return y;
-	return x > y ? x - y : 0;
 }

libc/tinymath/fdimf.c

@@ -1,45 +1,27 @@
 /*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
 │vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
-╚──────────────────────────────────────────────────────────────────────────────╝
+╞══════════════════════════════════════════════════════════════════════════════╡
+│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│  Musl Libc                                                                   │
+│ Permission to use, copy, modify, and/or distribute this software for         │
-│  Copyright © 2005-2020 Rich Felker, et al.                                   │
+│ any purpose with or without fee is hereby granted, provided that the         │
-│                                                                              │
+│ above copyright notice and this permission notice appear in all copies.      │
-│  Permission is hereby granted, free of charge, to any person obtaining       │
-│  a copy of this software and associated documentation files (the             │
-│  "Software"), to deal in the Software without restriction, including         │
-│  without limitation the rights to use, copy, modify, merge, publish,         │
-│  distribute, sublicense, and/or sell copies of the Software, and to          │
-│  permit persons to whom the Software is furnished to do so, subject to       │
-│  the following conditions:                                                   │
-│                                                                              │
-│  The above copyright notice and this permission notice shall be              │
-│  included in all copies or substantial portions of the Software.             │
-│                                                                              │
-│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
-│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
-│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
-│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
-│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
-│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
-│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
+│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
 
-asm(".ident\t\"\\n\\n\
+/**
-Musl libc (MIT License)\\n\
+ * Returns positive difference.
-Copyright 2005-2020 Rich Felker, et. al.\"");
+ */
-asm(".include \"libc/disclaimer.inc\"");
+float fdimf(float x, float y) {
-
+  if (isnan(x) || isnan(y)) return NAN;
-/* clang-format off */
+  return x > y ? x - y : 0;
-
-float fdimf(float x, float y)
-{
-	if (isnan(x))
-		return x;
-	if (isnan(y))
-		return y;
-	return x > y ? x - y : 0;
 }

libc/tinymath/fdiml.c

@@ -1,45 +1,27 @@
 /*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
 │vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
-╚──────────────────────────────────────────────────────────────────────────────╝
+╞══════════════════════════════════════════════════════════════════════════════╡
+│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│  Musl Libc                                                                   │
+│ Permission to use, copy, modify, and/or distribute this software for         │
-│  Copyright © 2005-2020 Rich Felker, et al.                                   │
+│ any purpose with or without fee is hereby granted, provided that the         │
-│                                                                              │
+│ above copyright notice and this permission notice appear in all copies.      │
-│  Permission is hereby granted, free of charge, to any person obtaining       │
-│  a copy of this software and associated documentation files (the             │
-│  "Software"), to deal in the Software without restriction, including         │
-│  without limitation the rights to use, copy, modify, merge, publish,         │
-│  distribute, sublicense, and/or sell copies of the Software, and to          │
-│  permit persons to whom the Software is furnished to do so, subject to       │
-│  the following conditions:                                                   │
-│                                                                              │
-│  The above copyright notice and this permission notice shall be              │
-│  included in all copies or substantial portions of the Software.             │
-│                                                                              │
-│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
-│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
-│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
-│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
-│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
-│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
-│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
+│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
 
-asm(".ident\t\"\\n\\n\
+/**
-Musl libc (MIT License)\\n\
+ * Returns positive difference.
-Copyright 2005-2020 Rich Felker, et. al.\"");
+ */
-asm(".include \"libc/disclaimer.inc\"");
+long double fdiml(long double x, long double y) {
-
+  if (isnan(x) || isnan(y)) return NAN;
-/* clang-format off */
+  return x > y ? x - y : 0;
-
-long double fdiml(long double x, long double y)
-{
-	if (isnan(x))
-		return x;
-	if (isnan(y))
-		return y;
-	return x > y ? x - y : 0;
 }

libc/tinymath/feval.internal.h

@@ -0,0 +1,20 @@
+#ifndef COSMOPOLITAN_LIBC_TINYMATH_FEVAL_INTERNAL_H_
+#define COSMOPOLITAN_LIBC_TINYMATH_FEVAL_INTERNAL_H_
+#if !(__ASSEMBLER__ + __LINKER__ + 0)
+COSMOPOLITAN_C_START_
+
+static inline void fevalf(float x) {
+  volatile float y = x;
+}
+
+static inline void feval(double x) {
+  volatile double y = x;
+}
+
+static inline void fevall(long double x) {
+  volatile long double y = x;
+}
+
+COSMOPOLITAN_C_END_
+#endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */
+#endif /* COSMOPOLITAN_LIBC_TINYMATH_FEVAL_INTERNAL_H_ */

libc/tinymath/fma.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -25,14 +25,12 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2020 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
 #define ASUINT64(x) ((union {double f; uint64_t i;}){x}).i

libc/tinymath/frexp.S

@@ -1,58 +0,0 @@
-/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8     -*-│
-│vi: set et ft=asm ts=8 tw=8 fenc=utf-8                                     :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
-│ Copyright 2020 Justine Alexandra Roberts Tunney                              │
-│                                                                              │
-│ Permission to use, copy, modify, and/or 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.                                                │
-╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/macros.internal.h"
-.source	__FILE__
-
-frexp:	.leafprologue
-	push	%rbx
-	push	%rdx
-	mov	%rdi,%rbx
-	movq	%xmm0,%rax
-	shr	$52,%rax
-	and	$0x7ff,%eax
-	jne	3f
-	xorps	%xmm1,%xmm1
-	ucomisd	%xmm1,%xmm0
-	jp	1f
-	je	2f
-1:	mulsd	6f(%rip),%xmm0
-	mov	%rbx,%rdi
-	call	frexp
-	subl	$64,(%rbx)
-	jmp	5f
-2:	movl	$0,(%rdi)
-	jmp	5f
-3:	cmp	$0x7ff,%eax
-	je	5f
-	movq	%xmm0,%rdx
-	sub	$0x3fe,%eax
-	mov	%eax,(%rdi)
-	movabs	$-9218868437227405313,%rax
-	and	%rax,%rdx
-	mov	$511,%eax
-	sal	$53,%rax
-	or	%rax,%rdx
-	movq	%rdx,%xmm0
-5:	pop	%rax
-	pop	%rbx
-	.leafepilogue
-	.endfn	frexp,globl
-
-	.rodata.cst8
-6:	.long	0,0x43f00000

libc/tinymath/frexp.c

@@ -0,0 +1,58 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│  Musl Libc                                                                   │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/math.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2020 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
+
+/**
+ * Splits number normalized fraction and exponent.
+ */
+double frexp(double x, int *e)
+{
+	union { double d; uint64_t i; } y = { x };
+	int ee = y.i>>52 & 0x7ff;
+
+	if (!ee) {
+		if (x) {
+			x = frexp(x*0x1p64, e);
+			*e -= 64;
+		} else *e = 0;
+		return x;
+	} else if (ee == 0x7ff) {
+		return x;
+	}
+
+	*e = ee - 0x3fe;
+	y.i &= 0x800fffffffffffffull;
+	y.i |= 0x3fe0000000000000ull;
+	return y.d;
+}

libc/tinymath/frexpf.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -25,16 +25,17 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2020 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
+/**
+ * Splits number normalized fraction and exponent.
+ */
 float frexpf(float x, int *e)
 {
 	union { float f; uint32_t i; } y = { x };

libc/tinymath/frexpl.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -25,24 +25,18 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
+#include "libc/tinymath/ldshape.internal.h"
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2020 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
-union ldshape {
+/**
-	long double f;
+ * Splits number normalized fraction and exponent.
-	struct {
+ */
-		uint64_t m;
-		uint16_t se;
-	} i;
-};
-
 long double frexpl(long double x, int *e)
 {
 	union ldshape u = {x};

libc/tinymath/gamma.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │

libc/tinymath/kcos.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -26,6 +26,7 @@
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/kernel.internal.h"
 
 asm(".ident\t\"\\n\\n\
 fdlibm (fdlibm license)\\n\

libc/tinymath/kernel.internal.h

@@ -0,0 +1,14 @@
+#ifndef COSMOPOLITAN_LIBC_TINYMATH_KERNEL_INTERNAL_H_
+#define COSMOPOLITAN_LIBC_TINYMATH_KERNEL_INTERNAL_H_
+#if !(__ASSEMBLER__ + __LINKER__ + 0)
+COSMOPOLITAN_C_START_
+
+double __sin(double, double, int);
+double __tan(double, double, int);
+double __cos(double, double);
+int __rem_pio2(double, double *);
+int __rem_pio2_large(double *, double *, int, int, int);
+
+COSMOPOLITAN_C_END_
+#endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */
+#endif /* COSMOPOLITAN_LIBC_TINYMATH_KERNEL_INTERNAL_H_ */

libc/tinymath/ksin.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -26,6 +26,7 @@
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/kernel.internal.h"
 
 asm(".ident\t\"\\n\\n\
 fdlibm (fdlibm license)\\n\

libc/tinymath/ktan.c

@@ -0,0 +1,154 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│  Musl Libc                                                                   │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/math.h"
+#include "libc/tinymath/kernel.internal.h"
+
+asm(".ident\t\"\\n\\n\
+fdlibm (fdlibm license)\\n\
+Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\"");
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+
+/* clang-format off */
+/* origin: FreeBSD /usr/src/lib/msun/src/k_tan.c */
+/*
+ * ====================================================
+ * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/* __tan( x, y, k )
+ * kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
+ * Input x is assumed to be bounded by ~pi/4 in magnitude.
+ * Input y is the tail of x.
+ * Input odd indicates whether tan (if odd = 0) or -1/tan (if odd = 1) is returned.
+ *
+ * Algorithm
+ *      1. Since tan(-x) = -tan(x), we need only to consider positive x.
+ *      2. Callers must return tan(-0) = -0 without calling here since our
+ *         odd polynomial is not evaluated in a way that preserves -0.
+ *         Callers may do the optimization tan(x) ~ x for tiny x.
+ *      3. tan(x) is approximated by a odd polynomial of degree 27 on
+ *         [0,0.67434]
+ *                               3             27
+ *              tan(x) ~ x + T1*x + ... + T13*x
+ *         where
+ *
+ *              |tan(x)         2     4            26   |     -59.2
+ *              |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
+ *              |  x                                    |
+ *
+ *         Note: tan(x+y) = tan(x) + tan'(x)*y
+ *                        ~ tan(x) + (1+x*x)*y
+ *         Therefore, for better accuracy in computing tan(x+y), let
+ *                   3      2      2       2       2
+ *              r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
+ *         then
+ *                                  3    2
+ *              tan(x+y) = x + (T1*x + (x *(r+y)+y))
+ *
+ *      4. For x in [0.67434,pi/4],  let y = pi/4 - x, then
+ *              tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
+ *                     = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
+ */
+
+#define ASUINT64(f)         ((union{double _f; uint64_t _i;}){f})._i
+#define ASDOUBLE(i)         ((union{uint64_t _i; double _f;}){i})._f
+#define GET_HIGH_WORD(hi,d) (hi) = ASUINT64(d) >> 32
+#define GET_LOW_WORD(lo,d)  (lo) = (uint32_t)ASUINT64(d)
+#define SET_LOW_WORD(d,lo)  INSERT(d, ASUINT64(d)>>32, lo)
+#define INSERT(d,hi,lo)     (d)=ASDOUBLE((uint64_t)(hi)<<32|(uint32_t)(lo))
+
+static const double T[] = {
+             3.33333333333334091986e-01, /* 3FD55555, 55555563 */
+             1.33333333333201242699e-01, /* 3FC11111, 1110FE7A */
+             5.39682539762260521377e-02, /* 3FABA1BA, 1BB341FE */
+             2.18694882948595424599e-02, /* 3F9664F4, 8406D637 */
+             8.86323982359930005737e-03, /* 3F8226E3, E96E8493 */
+             3.59207910759131235356e-03, /* 3F6D6D22, C9560328 */
+             1.45620945432529025516e-03, /* 3F57DBC8, FEE08315 */
+             5.88041240820264096874e-04, /* 3F4344D8, F2F26501 */
+             2.46463134818469906812e-04, /* 3F3026F7, 1A8D1068 */
+             7.81794442939557092300e-05, /* 3F147E88, A03792A6 */
+             7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */
+            -1.85586374855275456654e-05, /* BEF375CB, DB605373 */
+             2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */
+},
+pio4 =       7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */
+pio4lo =     3.06161699786838301793e-17; /* 3C81A626, 33145C07 */
+
+double __tan(double x, double y, int odd)
+{
+	double_t z, r, v, w, s, a;
+	double w0, a0;
+	uint32_t hx;
+	int big, sign;
+
+	GET_HIGH_WORD(hx,x);
+	big = (hx&0x7fffffff) >= 0x3FE59428; /* |x| >= 0.6744 */
+	if (big) {
+		sign = hx>>31;
+		if (sign) {
+			x = -x;
+			y = -y;
+		}
+		x = (pio4 - x) + (pio4lo - y);
+		y = 0.0;
+	}
+	z = x * x;
+	w = z * z;
+	/*
+	 * Break x^5*(T[1]+x^2*T[2]+...) into
+	 * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
+	 * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
+	 */
+	r = T[1] + w*(T[3] + w*(T[5] + w*(T[7] + w*(T[9] + w*T[11]))));
+	v = z*(T[2] + w*(T[4] + w*(T[6] + w*(T[8] + w*(T[10] + w*T[12])))));
+	s = z * x;
+	r = y + z*(s*(r + v) + y) + s*T[0];
+	w = x + r;
+	if (big) {
+		s = 1 - 2*odd;
+		v = s - 2.0 * (x + (r - w*w/(w + s)));
+		return sign ? -v : v;
+	}
+	if (!odd)
+		return w;
+	/* -1.0/(x+r) has up to 2ulp error, so compute it accurately */
+	w0 = w;
+	SET_LOW_WORD(w0, 0);
+	v = r - (w0 - x);       /* w0+v = r+x */
+	a0 = a = -1.0 / w;
+	SET_LOW_WORD(a0, 0);
+	return a0 + a*(1.0 + a0*w0 + a0*v);
+}

libc/tinymath/ldshape.internal.h

@@ -0,0 +1,16 @@
+#ifndef COSMOPOLITAN_LIBC_TINYMATH_LDSHAPE_INTERNAL_H_
+#define COSMOPOLITAN_LIBC_TINYMATH_LDSHAPE_INTERNAL_H_
+#if !(__ASSEMBLER__ + __LINKER__ + 0)
+COSMOPOLITAN_C_START_
+
+union ldshape {
+  long double f;
+  struct {
+    uint64_t m;
+    uint16_t se;
+  } i;
+};
+
+COSMOPOLITAN_C_END_
+#endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */
+#endif /* COSMOPOLITAN_LIBC_TINYMATH_LDSHAPE_INTERNAL_H_ */

libc/tinymath/modf.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │

libc/tinymath/modff.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │

libc/tinymath/modfl.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │

libc/tinymath/nextafter.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -31,7 +31,6 @@ asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2014 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
 static inline void force_eval(double x)

libc/tinymath/nextafterf.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -26,20 +26,14 @@
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2014 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
-static inline void force_eval(float x)
-{
-	volatile float y;
-	y = x;
-}
-
 float nextafterf(float x, float y)
 {
 	union {float f; uint32_t i;} ux={x}, uy={y};
@@ -62,9 +56,9 @@ float nextafterf(float x, float y)
 	e = ux.i & 0x7f800000;
 	/* raise overflow if ux.f is infinite and x is finite */
 	if (e == 0x7f800000)
-		force_eval(x+x);
+		feval(x+x);
 	/* raise underflow if ux.f is subnormal or zero */
 	if (e == 0)
-		force_eval(x*x + ux.f*ux.f);
+		feval(x*x + ux.f*ux.f);
 	return ux.f;
 }

libc/tinymath/nextafterl.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -26,28 +26,15 @@
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
+#include "libc/tinymath/ldshape.internal.h"
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2014 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
-union ldshape {
-	long double f;
-	struct {
-		uint64_t m;
-		uint16_t se;
-	} i;
-};
-
-static inline void force_eval(long double x)
-{
-	volatile long double y;
-	y = x;
-}
-
 long double nextafterl(long double x, long double y)
 {
 	union ldshape ux, uy;
@@ -80,6 +67,6 @@ long double nextafterl(long double x, long double y)
 		return x + x;
 	/* raise underflow if ux is subnormal or zero */
 	if ((ux.i.se & 0x7fff) == 0)
-		force_eval(x*x + ux.f*ux.f);
+		fevall(x*x + ux.f*ux.f);
 	return ux.f;
 }

libc/tinymath/nexttoward.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -25,14 +25,12 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2020 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
 static inline void fp_force_eval(double x)

libc/tinymath/nexttowardf.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -25,22 +25,15 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2020 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
-static inline void fp_force_evalf(float x)
-{
-	volatile float y;
-	y = x;
-}
-
 float nexttowardf(float x, long double y)
 {
 	union {float f; uint32_t i;} ux = {x};
@@ -68,9 +61,9 @@ float nexttowardf(float x, long double y)
 	e = ux.i & 0x7f800000;
 	/* raise overflow if ux.f is infinite and x is finite */
 	if (e == 0x7f800000)
-		fp_force_evalf(x+x);
+		feval(x+x);
 	/* raise underflow if ux.f is subnormal or zero */
 	if (e == 0)
-		fp_force_evalf(x*x + ux.f*ux.f);
+		feval(x*x + ux.f*ux.f);
 	return ux.f;
 }

libc/tinymath/nexttowardl.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -25,14 +25,12 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2020 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
 long double nexttowardl(long double x, long double y)

libc/tinymath/rempio2.c

@@ -0,0 +1,226 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│  Musl Libc                                                                   │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/bits/likely.h"
+#include "libc/math.h"
+#include "libc/tinymath/kernel.internal.h"
+
+asm(".ident\t\"\\n\\n\
+fdlibm (fdlibm license)\\n\
+Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\"");
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+
+/* clang-format off */
+/* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+/* __rem_pio2(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __rem_pio2_large() for large x
+ */
+
+#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
+#define EPS DBL_EPSILON
+#elif FLT_EVAL_METHOD==2
+#define EPS LDBL_EPSILON
+#endif
+
+/*
+ * invpio2:  53 bits of 2/pi
+ * pio2_1:   first  33 bit of pi/2
+ * pio2_1t:  pi/2 - pio2_1
+ * pio2_2:   second 33 bit of pi/2
+ * pio2_2t:  pi/2 - (pio2_1+pio2_2)
+ * pio2_3:   third  33 bit of pi/2
+ * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+#define toint (1.5/EPS)
+#define pio4 0x1.921fb54442d18p-1
+#define invpio2 6.36619772367581382433e-01
+#define pio2_1 1.57079632673412561417e+00
+#define pio2_1t 6.07710050650619224932e-11
+#define pio2_2 6.07710050630396597660e-11
+#define pio2_2t 2.02226624879595063154e-21
+#define pio2_3 2.02226624871116645580e-21
+#define pio2_3t 8.47842766036889956997e-32
+
+/* caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
+int __rem_pio2(double x, double *y)
+{
+	union {double f; uint64_t i;} u = {x};
+	double_t z,w,t,r,fn;
+	double tx[3],ty[2];
+	uint32_t ix;
+	int sign, n, ex, ey, i;
+	sign = u.i>>63;
+	ix = u.i>>32 & 0x7fffffff;
+	if (ix <= 0x400f6a7a) {  /* |x| ~<= 5pi/4 */
+		if ((ix & 0xfffff) == 0x921fb)  /* |x| ~= pi/2 or 2pi/2 */
+			goto medium;  /* cancellation -- use medium case */
+		if (ix <= 0x4002d97c) {  /* |x| ~<= 3pi/4 */
+			if (!sign) {
+				z = x - pio2_1;  /* one round good to 85 bits */
+				y[0] = z - pio2_1t;
+				y[1] = (z-y[0]) - pio2_1t;
+				return 1;
+			} else {
+				z = x + pio2_1;
+				y[0] = z + pio2_1t;
+				y[1] = (z-y[0]) + pio2_1t;
+				return -1;
+			}
+		} else {
+			if (!sign) {
+				z = x - 2*pio2_1;
+				y[0] = z - 2*pio2_1t;
+				y[1] = (z-y[0]) - 2*pio2_1t;
+				return 2;
+			} else {
+				z = x + 2*pio2_1;
+				y[0] = z + 2*pio2_1t;
+				y[1] = (z-y[0]) + 2*pio2_1t;
+				return -2;
+			}
+		}
+	}
+	if (ix <= 0x401c463b) {  /* |x| ~<= 9pi/4 */
+		if (ix <= 0x4015fdbc) {  /* |x| ~<= 7pi/4 */
+			if (ix == 0x4012d97c)  /* |x| ~= 3pi/2 */
+				goto medium;
+			if (!sign) {
+				z = x - 3*pio2_1;
+				y[0] = z - 3*pio2_1t;
+				y[1] = (z-y[0]) - 3*pio2_1t;
+				return 3;
+			} else {
+				z = x + 3*pio2_1;
+				y[0] = z + 3*pio2_1t;
+				y[1] = (z-y[0]) + 3*pio2_1t;
+				return -3;
+			}
+		} else {
+			if (ix == 0x401921fb)  /* |x| ~= 4pi/2 */
+				goto medium;
+			if (!sign) {
+				z = x - 4*pio2_1;
+				y[0] = z - 4*pio2_1t;
+				y[1] = (z-y[0]) - 4*pio2_1t;
+				return 4;
+			} else {
+				z = x + 4*pio2_1;
+				y[0] = z + 4*pio2_1t;
+				y[1] = (z-y[0]) + 4*pio2_1t;
+				return -4;
+			}
+		}
+	}
+	if (ix < 0x413921fb) {  /* |x| ~< 2^20*(pi/2), medium size */
+medium:
+		/* rint(x/(pi/2)) */
+		fn = (double_t)x*invpio2 + toint - toint;
+		n = (int32_t)fn;
+		r = x - fn*pio2_1;
+		w = fn*pio2_1t;  /* 1st round, good to 85 bits */
+		/* Matters with directed rounding. */
+		if (UNLIKELY(r - w < -pio4)) {
+			n--;
+			fn--;
+			r = x - fn*pio2_1;
+			w = fn*pio2_1t;
+		} else if (UNLIKELY(r - w > pio4)) {
+			n++;
+			fn++;
+			r = x - fn*pio2_1;
+			w = fn*pio2_1t;
+		}
+		y[0] = r - w;
+		u.f = y[0];
+		ey = u.i>>52 & 0x7ff;
+		ex = ix>>20;
+		if (ex - ey > 16) { /* 2nd round, good to 118 bits */
+			t = r;
+			w = fn*pio2_2;
+			r = t - w;
+			w = fn*pio2_2t - ((t-r)-w);
+			y[0] = r - w;
+			u.f = y[0];
+			ey = u.i>>52 & 0x7ff;
+			if (ex - ey > 49) {  /* 3rd round, good to 151 bits, covers all cases */
+				t = r;
+				w = fn*pio2_3;
+				r = t - w;
+				w = fn*pio2_3t - ((t-r)-w);
+				y[0] = r - w;
+			}
+		}
+		y[1] = (r - y[0]) - w;
+		return n;
+	}
+	/*
+	 * all other (large) arguments
+	 */
+	if (ix >= 0x7ff00000) {  /* x is inf or NaN */
+		y[0] = y[1] = x - x;
+		return 0;
+	}
+	/* set z = scalbn(|x|,-ilogb(x)+23) */
+	u.f = x;
+	u.i &= (uint64_t)-1>>12;
+	u.i |= (uint64_t)(0x3ff + 23)<<52;
+	z = u.f;
+	for (i=0; i < 2; i++) {
+		tx[i] = (double)(int32_t)z;
+		z     = (z-tx[i])*0x1p24;
+	}
+	tx[i] = z;
+	/* skip zero terms, first term is non-zero */
+	while (tx[i] == 0.0)
+		i--;
+	n = __rem_pio2_large(tx,ty,(int)(ix>>20)-(0x3ff+23),i+1,1);
+	if (sign) {
+		y[0] = -ty[0];
+		y[1] = -ty[1];
+		return -n;
+	}
+	y[0] = ty[0];
+	y[1] = ty[1];
+	return n;
+}

libc/tinymath/rempio2large.c

@@ -0,0 +1,478 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│  Musl Libc                                                                   │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/math.h"
+#include "libc/tinymath/kernel.internal.h"
+
+asm(".ident\t\"\\n\\n\
+fdlibm (fdlibm license)\\n\
+Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\"");
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+
+/* clang-format off */
+/* origin: FreeBSD /usr/src/lib/msun/src/k_rem_pio2.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/*
+ * __rem_pio2_large(x,y,e0,nx,prec)
+ * double x[],y[]; int e0,nx,prec;
+ *
+ * __rem_pio2_large return the last three digits of N with
+ *              y = x - N*pi/2
+ * so that |y| < pi/2.
+ *
+ * The method is to compute the integer (mod 8) and fraction parts of
+ * (2/pi)*x without doing the full multiplication. In general we
+ * skip the part of the product that are known to be a huge integer (
+ * more accurately, = 0 mod 8 ). Thus the number of operations are
+ * independent of the exponent of the input.
+ *
+ * (2/pi) is represented by an array of 24-bit integers in ipio2[].
+ *
+ * Input parameters:
+ *      x[]     The input value (must be positive) is broken into nx
+ *              pieces of 24-bit integers in double precision format.
+ *              x[i] will be the i-th 24 bit of x. The scaled exponent
+ *              of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
+ *              match x's up to 24 bits.
+ *
+ *              Example of breaking a double positive z into x[0]+x[1]+x[2]:
+ *                      e0 = ilogb(z)-23
+ *                      z  = scalbn(z,-e0)
+ *              for i = 0,1,2
+ *                      x[i] = floor(z)
+ *                      z    = (z-x[i])*2**24
+ *
+ *
+ *      y[]     ouput result in an array of double precision numbers.
+ *              The dimension of y[] is:
+ *                      24-bit  precision       1
+ *                      53-bit  precision       2
+ *                      64-bit  precision       2
+ *                      113-bit precision       3
+ *              The actual value is the sum of them. Thus for 113-bit
+ *              precison, one may have to do something like:
+ *
+ *              long double t,w,r_head, r_tail;
+ *              t = (long double)y[2] + (long double)y[1];
+ *              w = (long double)y[0];
+ *              r_head = t+w;
+ *              r_tail = w - (r_head - t);
+ *
+ *      e0      The exponent of x[0]. Must be <= 16360 or you need to
+ *              expand the ipio2 table.
+ *
+ *      nx      dimension of x[]
+ *
+ *      prec    an integer indicating the precision:
+ *                      0       24  bits (single)
+ *                      1       53  bits (double)
+ *                      2       64  bits (extended)
+ *                      3       113 bits (quad)
+ *
+ * External function:
+ *      double scalbn(), floor();
+ *
+ *
+ * Here is the description of some local variables:
+ *
+ *      jk      jk+1 is the initial number of terms of ipio2[] needed
+ *              in the computation. The minimum and recommended value
+ *              for jk is 3,4,4,6 for single, double, extended, and quad.
+ *              jk+1 must be 2 larger than you might expect so that our
+ *              recomputation test works. (Up to 24 bits in the integer
+ *              part (the 24 bits of it that we compute) and 23 bits in
+ *              the fraction part may be lost to cancelation before we
+ *              recompute.)
+ *
+ *      jz      local integer variable indicating the number of
+ *              terms of ipio2[] used.
+ *
+ *      jx      nx - 1
+ *
+ *      jv      index for pointing to the suitable ipio2[] for the
+ *              computation. In general, we want
+ *                      ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
+ *              is an integer. Thus
+ *                      e0-3-24*jv >= 0 or (e0-3)/24 >= jv
+ *              Hence jv = max(0,(e0-3)/24).
+ *
+ *      jp      jp+1 is the number of terms in PIo2[] needed, jp = jk.
+ *
+ *      q[]     double array with integral value, representing the
+ *              24-bits chunk of the product of x and 2/pi.
+ *
+ *      q0      the corresponding exponent of q[0]. Note that the
+ *              exponent for q[i] would be q0-24*i.
+ *
+ *      PIo2[]  double precision array, obtained by cutting pi/2
+ *              into 24 bits chunks.
+ *
+ *      f[]     ipio2[] in floating point
+ *
+ *      iq[]    integer array by breaking up q[] in 24-bits chunk.
+ *
+ *      fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]
+ *
+ *      ih      integer. If >0 it indicates q[] is >= 0.5, hence
+ *              it also indicates the *sign* of the result.
+ *
+ */
+/*
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+static const int init_jk[] = {3,4,4,6}; /* initial value for jk */
+
+/*
+ * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
+ *
+ *              integer array, contains the (24*i)-th to (24*i+23)-th
+ *              bit of 2/pi after binary point. The corresponding
+ *              floating value is
+ *
+ *                      ipio2[i] * 2^(-24(i+1)).
+ *
+ * NB: This table must have at least (e0-3)/24 + jk terms.
+ *     For quad precision (e0 <= 16360, jk = 6), this is 686.
+ */
+static const int32_t ipio2[] = {
+0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
+0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
+0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
+0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
+0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
+0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
+0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
+0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
+0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
+0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
+0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
+#if LDBL_MAX_EXP > 1024
+0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
+0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
+0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
+0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
+0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
+0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
+0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
+0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
+0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
+0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
+0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
+0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
+0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
+0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
+0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
+0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
+0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
+0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
+0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
+0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
+0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
+0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
+0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
+0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
+0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
+0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
+0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
+0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
+0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
+0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
+0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
+0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
+0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
+0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
+0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
+0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
+0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
+0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
+0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
+0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
+0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
+0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
+0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
+0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
+0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
+0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
+0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
+0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
+0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
+0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
+0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
+0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
+0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
+0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
+0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
+0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
+0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
+0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
+0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
+0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
+0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
+0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
+0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
+0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
+0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
+0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
+0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
+0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
+0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
+0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
+0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
+0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
+0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
+0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
+0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
+0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
+0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
+0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
+0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
+0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
+0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
+0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
+0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
+0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
+0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
+0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
+0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
+0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
+0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
+0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
+0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
+0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
+0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
+0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
+0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
+0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
+0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
+0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
+0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
+0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
+0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
+0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
+0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
+0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
+#endif
+};
+
+static const double PIo2[] = {
+  1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
+  7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
+  5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
+  3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
+  1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
+  1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
+  2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
+  2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
+};
+
+int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec)
+{
+	int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
+	double z,fw,f[20],fq[20],q[20];
+
+	/* initialize jk*/
+	jk = init_jk[prec];
+	jp = jk;
+
+	/* determine jx,jv,q0, note that 3>q0 */
+	jx = nx-1;
+	jv = (e0-3)/24;  if(jv<0) jv=0;
+	q0 = e0-24*(jv+1);
+
+	/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
+	j = jv-jx; m = jx+jk;
+	for (i=0; i<=m; i++,j++)
+		f[i] = j<0 ? 0.0 : (double)ipio2[j];
+
+	/* compute q[0],q[1],...q[jk] */
+	for (i=0; i<=jk; i++) {
+		for (j=0,fw=0.0; j<=jx; j++)
+			fw += x[j]*f[jx+i-j];
+		q[i] = fw;
+	}
+
+	jz = jk;
+recompute:
+	/* distill q[] into iq[] reversingly */
+	for (i=0,j=jz,z=q[jz]; j>0; i++,j--) {
+		fw    = (double)(int32_t)(0x1p-24*z);
+		iq[i] = (int32_t)(z - 0x1p24*fw);
+		z     = q[j-1]+fw;
+	}
+
+	/* compute n */
+	z  = scalbn(z,q0);       /* actual value of z */
+	z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
+	n  = (int32_t)z;
+	z -= (double)n;
+	ih = 0;
+	if (q0 > 0) {  /* need iq[jz-1] to determine n */
+		i  = iq[jz-1]>>(24-q0); n += i;
+		iq[jz-1] -= i<<(24-q0);
+		ih = iq[jz-1]>>(23-q0);
+	}
+	else if (q0 == 0) ih = iq[jz-1]>>23;
+	else if (z >= 0.5) ih = 2;
+
+	if (ih > 0) {  /* q > 0.5 */
+		n += 1; carry = 0;
+		for (i=0; i<jz; i++) {  /* compute 1-q */
+			j = iq[i];
+			if (carry == 0) {
+				if (j != 0) {
+					carry = 1;
+					iq[i] = 0x1000000 - j;
+				}
+			} else
+				iq[i] = 0xffffff - j;
+		}
+		if (q0 > 0) {  /* rare case: chance is 1 in 12 */
+			switch(q0) {
+			case 1:
+				iq[jz-1] &= 0x7fffff; break;
+			case 2:
+				iq[jz-1] &= 0x3fffff; break;
+			}
+		}
+		if (ih == 2) {
+			z = 1.0 - z;
+			if (carry != 0)
+				z -= scalbn(1.0,q0);
+		}
+	}
+
+	/* check if recomputation is needed */
+	if (z == 0.0) {
+		j = 0;
+		for (i=jz-1; i>=jk; i--) j |= iq[i];
+		if (j == 0) {  /* need recomputation */
+			for (k=1; iq[jk-k]==0; k++);  /* k = no. of terms needed */
+
+			for (i=jz+1; i<=jz+k; i++) {  /* add q[jz+1] to q[jz+k] */
+				f[jx+i] = (double)ipio2[jv+i];
+				for (j=0,fw=0.0; j<=jx; j++)
+					fw += x[j]*f[jx+i-j];
+				q[i] = fw;
+			}
+			jz += k;
+			goto recompute;
+		}
+	}
+
+	/* chop off zero terms */
+	if (z == 0.0) {
+		jz -= 1;
+		q0 -= 24;
+		while (iq[jz] == 0) {
+			jz--;
+			q0 -= 24;
+		}
+	} else { /* break z into 24-bit if necessary */
+		z = scalbn(z,-q0);
+		if (z >= 0x1p24) {
+			fw = (double)(int32_t)(0x1p-24*z);
+			iq[jz] = (int32_t)(z - 0x1p24*fw);
+			jz += 1;
+			q0 += 24;
+			iq[jz] = (int32_t)fw;
+		} else
+			iq[jz] = (int32_t)z;
+	}
+
+	/* convert integer "bit" chunk to floating-point value */
+	fw = scalbn(1.0,q0);
+	for (i=jz; i>=0; i--) {
+		q[i] = fw*(double)iq[i];
+		fw *= 0x1p-24;
+	}
+
+	/* compute PIo2[0,...,jp]*q[jz,...,0] */
+	for(i=jz; i>=0; i--) {
+		for (fw=0.0,k=0; k<=jp && k<=jz-i; k++)
+			fw += PIo2[k]*q[i+k];
+		fq[jz-i] = fw;
+	}
+
+	/* compress fq[] into y[] */
+	switch(prec) {
+	case 0:
+		fw = 0.0;
+		for (i=jz; i>=0; i--)
+			fw += fq[i];
+		y[0] = ih==0 ? fw : -fw;
+		break;
+	case 1:
+	case 2:
+		fw = 0.0;
+		for (i=jz; i>=0; i--)
+			fw += fq[i];
+		// TODO: drop excess precision here once double_t is used
+		fw = (double)fw;
+		y[0] = ih==0 ? fw : -fw;
+		fw = fq[0]-fw;
+		for (i=1; i<=jz; i++)
+			fw += fq[i];
+		y[1] = ih==0 ? fw : -fw;
+		break;
+	case 3:  /* painful */
+		for (i=jz; i>0; i--) {
+			fw      = fq[i-1]+fq[i];
+			fq[i]  += fq[i-1]-fw;
+			fq[i-1] = fw;
+		}
+		for (i=jz; i>1; i--) {
+			fw      = fq[i-1]+fq[i];
+			fq[i]  += fq[i-1]-fw;
+			fq[i-1] = fw;
+		}
+		for (fw=0.0,i=jz; i>=2; i--)
+			fw += fq[i];
+		if (ih==0) {
+			y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
+		} else {
+			y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
+		}
+	}
+	return n&7;
+}

libc/tinymath/remquo.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -25,16 +25,17 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2020 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
+/**
+ * Computes remainder and part of quotient.
+ */
 double remquo(double x, double y, int *quo)
 {
 	union {double f; uint64_t i;} ux = {x}, uy = {y};

libc/tinymath/remquof.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -25,16 +25,17 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2020 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
+/**
+ * Computes remainder and part of quotient.
+ */
 float remquof(float x, float y, int *quo)
 {
 	union {float f; uint32_t i;} ux = {x}, uy = {y};

libc/tinymath/remquol.c

@@ -1,5 +1,5 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Musl Libc                                                                   │
@@ -25,24 +25,18 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-
 #include "libc/math.h"
-
+#include "libc/tinymath/ldshape.internal.h"
-union ldshape {
-	long double f;
-	struct {
-		uint64_t m;
-		uint16_t se;
-	} i;
-};
 
 asm(".ident\t\"\\n\\n\
 Musl libc (MIT License)\\n\
 Copyright 2005-2020 Rich Felker, et. al.\"");
 asm(".include \"libc/disclaimer.inc\"");
-
 /* clang-format off */
 
+/**
+ * Computes remainder and part of quotient.
+ */
 long double remquol(long double x, long double y, int *quo)
 {
 	union ldshape ux = {x}, uy = {y};

libc/tinymath/sin.S

@@ -1,28 +0,0 @@
-/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8     -*-│
-│vi: set et ft=asm ts=8 tw=8 fenc=utf-8                                     :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
-│ Copyright 2020 Justine Alexandra Roberts Tunney                              │
-│                                                                              │
-│ Permission to use, copy, modify, and/or 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.                                                │
-╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/macros.internal.h"
-
-//	Returns sine of 𝑥.
-//
-//	@param	𝑥 is double scalar in low half of %xmm0
-//	@return	double scalar in low half of %xmm0
-//	@domain	-(3π/8) < 𝑥 < 3π/8 for best accuracy
-sin:	ezlea	sinl,ax
-	jmp	_d2ld2
-	.endfn	sin,globl

libc/tinymath/sin.c

@@ -0,0 +1,119 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│  Musl Libc                                                                   │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/math.h"
+#include "libc/tinymath/feval.internal.h"
+#include "libc/tinymath/kernel.internal.h"
+
+asm(".ident\t\"\\n\\n\
+fdlibm (fdlibm license)\\n\
+Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\"");
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+
+/* clang-format off */
+/* origin: FreeBSD /usr/src/lib/msun/src/s_sin.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/* sin(x)
+ * Return sine function of x.
+ *
+ * kernel function:
+ *      __sin            ... sine function on [-pi/4,pi/4]
+ *      __cos            ... cose function on [-pi/4,pi/4]
+ *      __rem_pio2       ... argument reduction routine
+ *
+ * Method.
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ *      in [-pi/4 , +pi/4], and let n = k mod 4.
+ *      We have
+ *
+ *          n        sin(x)      cos(x)        tan(x)
+ *     ----------------------------------------------------------
+ *          0          S           C             T
+ *          1          C          -S            -1/T
+ *          2         -S          -C             T
+ *          3         -C           S            -1/T
+ *     ----------------------------------------------------------
+ *
+ * Special cases:
+ *      Let trig be any of sin, cos, or tan.
+ *      trig(+-INF)  is NaN, with signals;
+ *      trig(NaN)    is that NaN;
+ *
+ * Accuracy:
+ *      TRIG(x) returns trig(x) nearly rounded
+ */
+
+#define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i
+#define gethighw(hi,d) (hi) = asuint64(d) >> 32
+
+double sin(double x)
+{
+	double y[2];
+	uint32_t ix;
+	unsigned n;
+
+	/* High word of x. */
+	gethighw(ix, x);
+	ix &= 0x7fffffff;
+
+	/* |x| ~< pi/4 */
+	if (ix <= 0x3fe921fb) {
+		if (ix < 0x3e500000) {  /* |x| < 2**-26 */
+			/* raise inexact if x != 0 and underflow if subnormal*/
+			feval(ix < 0x00100000 ? x/0x1p120f : x+0x1p120f);
+			return x;
+		}
+		return __sin(x, 0.0, 0);
+	}
+
+	/* sin(Inf or NaN) is NaN */
+	if (ix >= 0x7ff00000)
+		return x - x;
+
+	/* argument reduction needed */
+	n = __rem_pio2(x, y);
+	switch (n&3) {
+	case 0: return  __sin(y[0], y[1], 1);
+	case 1: return  __cos(y[0], y[1]);
+	case 2: return -__sin(y[0], y[1], 1);
+	default:
+		return -__cos(y[0], y[1]);
+	}
+}

libc/tinymath/sincos.S

@@ -1,39 +0,0 @@
-/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8     -*-│
-│vi: set et ft=asm ts=8 tw=8 fenc=utf-8                                     :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
-│ Copyright 2020 Justine Alexandra Roberts Tunney                              │
-│                                                                              │
-│ Permission to use, copy, modify, and/or 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.                                                │
-╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/macros.internal.h"
-
-//	Returns sine and cosine of 𝑥.
-//
-//	@param	𝑥 is double scalar in low half of %xmm0
-//	@param	%rdi is double *out_sin
-//	@param	%rsi is double *out_cos
-//	@domain	-(3π/8) < 𝑥 < 3π/8 for best accuracy
-sincos:	push	%rbp
-	mov	%rsp,%rbp
-	.profilable
-	push	%rax
-	movsd	%xmm0,-8(%rbp)
-	fldl	-8(%rbp)
-	fsincos
-	fxch
-	fstpl	(%rdi)
-	fstpl	(%rsi)
-	leave
-	ret
-	.endfn	sincos,globl

libc/tinymath/sincos.c

@@ -0,0 +1,110 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│  Musl Libc                                                                   │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/math.h"
+#include "libc/runtime/runtime.h"
+#include "libc/tinymath/feval.internal.h"
+#include "libc/tinymath/kernel.internal.h"
+
+asm(".ident\t\"\\n\\n\
+fdlibm (fdlibm license)\\n\
+Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\"");
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+
+/* clang-format off */
+/* origin: FreeBSD /usr/src/lib/msun/src/s_sin.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i
+#define gethighw(hi,d) (hi) = asuint64(d) >> 32
+
+void sincos(double x, double *sin, double *cos)
+{
+	double y[2], s, c;
+	uint32_t ix;
+	unsigned n;
+
+	gethighw(ix, x);
+	ix &= 0x7fffffff;
+
+	/* |x| ~< pi/4 */
+	if (ix <= 0x3fe921fb) {
+		/* if |x| < 2**-27 * sqrt(2) */
+		if (ix < 0x3e46a09e) {
+			/* raise inexact if x!=0 and underflow if subnormal */
+			feval(ix < 0x00100000 ? x/0x1p120f : x+0x1p120f);
+			*sin = x;
+			*cos = 1.0;
+			return;
+		}
+		*sin = __sin(x, 0.0, 0);
+		*cos = __cos(x, 0.0);
+		return;
+	}
+
+	/* sincos(Inf or NaN) is NaN */
+	if (ix >= 0x7ff00000) {
+		*sin = *cos = x - x;
+		return;
+	}
+
+	/* argument reduction needed */
+	n = __rem_pio2(x, y);
+	s = __sin(y[0], y[1], 1);
+	c = __cos(y[0], y[1]);
+	switch (n&3) {
+	case 0:
+		*sin = s;
+		*cos = c;
+		break;
+	case 1:
+		*sin = c;
+		*cos = -s;
+		break;
+	case 2:
+		*sin = -s;
+		*cos = -c;
+		break;
+	case 3:
+	default:
+		*sin = -c;
+		*cos = s;
+		break;
+	}
+}

libc/tinymath/sinh.c

@@ -1,27 +1,76 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/expo.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns hyperbolic sine of 𝑥.
+ * 
+ *     sinh(x) = (exp(x) - 1/exp(x))/2
+ *             = (exp(x)-1 + (exp(x)-1)/exp(x))/2
+ *             = x + x^3/6 + o(x^5)
  */
-double sinh(double x) {
+double sinh(double x)
-  if (!x) return x;
+{
-  return (exp(x) - exp(-x)) / 2;
+	union {double f; uint64_t i;} u = {.f = x};
+	uint32_t w;
+	double t, h, absx;
+
+	h = 0.5;
+	if (u.i >> 63)
+		h = -h;
+	/* |x| */
+	u.i &= (uint64_t)-1/2;
+	absx = u.f;
+	w = u.i >> 32;
+
+	/* |x| < log(DBL_MAX) */
+	if (w < 0x40862e42) {
+		t = expm1(absx);
+		if (w < 0x3ff00000) {
+			if (w < 0x3ff00000 - (26<<20))
+				/* note: inexact and underflow are raised by expm1 */
+				/* note: this branch avoids spurious underflow */
+				return x;
+			return h*(2*t - t*t/(t+1));
+		}
+		/* note: |x|>log(0x1p26)+eps could be just h*exp(x) */
+		return h*(t + t/(t+1));
+	}
+
+	/* |x| > log(DBL_MAX) or nan */
+	/* note: the result is stored to handle overflow */
+	t = __expo2(absx, 2*h);
+	return t;
 }

libc/tinymath/sinhf.c

@@ -1,27 +1,72 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8                                :vi│
+│vi: set et ft=c ts=8 tw=8 fenc=utf-8                                       :vi│
-╞══════════════════════════════════════════════════════════════════════════════╡
+╚──────────────────────────────────────────────────────────────────────────────╝
-│ Copyright 2021 Justine Alexandra Roberts Tunney                              │
 │                                                                              │
-│ Permission to use, copy, modify, and/or distribute this software for         │
+│  Musl Libc                                                                   │
-│ any purpose with or without fee is hereby granted, provided that the         │
+│  Copyright © 2005-2014 Rich Felker, et al.                                   │
-│ above copyright notice and this permission notice appear in all copies.      │
+│                                                                              │
+│  Permission is hereby granted, free of charge, to any person obtaining       │
+│  a copy of this software and associated documentation files (the             │
+│  "Software"), to deal in the Software without restriction, including         │
+│  without limitation the rights to use, copy, modify, merge, publish,         │
+│  distribute, sublicense, and/or sell copies of the Software, and to          │
+│  permit persons to whom the Software is furnished to do so, subject to       │
+│  the following conditions:                                                   │
+│                                                                              │
+│  The above copyright notice and this permission notice shall be              │
+│  included in all copies or substantial portions of the Software.             │
+│                                                                              │
+│  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,             │
+│  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF          │
+│  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.      │
+│  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY        │
+│  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,        │
+│  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE           │
+│  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
-│ 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.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/math.h"
+#include "libc/tinymath/expo.internal.h"
+
+asm(".ident\t\"\\n\\n\
+Musl libc (MIT License)\\n\
+Copyright 2005-2014 Rich Felker, et. al.\"");
+asm(".include \"libc/disclaimer.inc\"");
+/* clang-format off */
 
 /**
  * Returns hyperbolic sine of 𝑥.
+ * 
+ *     sinh(x) = (exp(x) - 1/exp(x))/2
+ *             = (exp(x)-1 + (exp(x)-1)/exp(x))/2
+ *             = x + x^3/6 + o(x^5)
  */
-float sinhf(float x) {
+float sinhf(float x)
-  if (!x) return x;
+{
-  return (expf(x) - expf(-x)) / 2;
+	union {float f; uint32_t i;} u = {.f = x};
+	uint32_t w;
+	float t, h, absx;
+
+	h = 0.5;
+	if (u.i >> 31)
+		h = -h;
+	/* |x| */
+	u.i &= 0x7fffffff;
+	absx = u.f;
+	w = u.i;
+
+	/* |x| < log(FLT_MAX) */
+	if (w < 0x42b17217) {
+		t = expm1f(absx);
+		if (w < 0x3f800000) {
+			if (w < 0x3f800000 - (12<<23))
+				return x;
+			return h*(2*t - t*t/(t+1));
+		}
+		return h*(t + t/(t+1));
+	}
+
+	/* |x| > logf(FLT_MAX) or nan */
+	t = __expo2f(absx, 2*h);
+	return t;
 }