Release Cosmopolitan v3.3.1

- Write some more unit tests
- memcpy() on ARM is now faster
- Address the Musl complex math FIXME comments
- Some libm funcs like pow() now support setting errno
- Import the latest and greatest math functions from ARM
- Use more accurate atan2f() and log1pf() implementations
- atoi() and atol() will no longer saturate or clobber errno

CONTRIBUTING.md

@@ -16,10 +16,25 @@ contributors who prefer to remain anonymous to the public.
 
 The first time you send a pull request, you need to send an email to
 Justine Tunney <jtunney@gmail.com> stating that you intend to assign her
-the copyright to the changes you contribute to Cosmopolitan. This only
-applies to the code you *choose* to contribute. It only has to happen
-once. The email should be sent from an email address associated with
-your identity. Your email should link to your pull request.
+the copyright to the changes you contribute to Cosmopolitan. It only
+needs to happen once. This only applies to the code you *choose* to
+contribute. The email should be sent from an email address associated
+with your identity. Your email should link to your pull request.
+
+To make things easy, here's an example of a good email you can use:
+
+> **From**: YOUR NAME ([email protected])  
+> **To**: Justine Tunney ([email protected])  
+> **Subject**: Cosmopolitan Copyright Assignment for YOUR NAME
+>
+> Hi Justine,
+>
+> I made my first contribution to Cosmopolitan in
+> https://github.com/jart/cosmopolitan/pull/XXXX could you please take a
+> look? I intend to assign you the copyright to the changes I contribute
+> to Cosmopolitan.
+>
+> Thanks!
 
 Please note that in order to give Justine the copyright, it has to be
 yours to give in the first place. If you're employed, then you should
@@ -38,6 +53,27 @@ owners and the code should go in the `third_party/` folder. Every third
 party project should have a `README.cosmo` file that documents its
 provenance as well as any local changes you've made.
 
+## Copyright Policy Exceptions
+
+### Tests
+
+You're encoraged to claim ownership of your test code. If you add a new
+file under the `test/` directory, then you should put your name in the
+ISC license header at the top of the file. If you add new test cases to
+an existing unit test file, then you're encouraged to append a line with
+your name to the existing copyright header of that file.
+
+### Exceptional Features
+
+Let's say you discovered a faster better way to implement `log10()` and
+you want to give it to Cosmopolitan. In cases like this, it really isn't
+appropriate for Justine to own your code. What you could do instead, is
+write your own new and improved `log10.c` from scratch, put your name on
+the top with the ISC license, and then add a `__notice()` directive so
+that your name will be embedded inside every executable that links the
+`log10()` function. This will help you get your name out there. Please
+note you need get approval from Justine each time you want to do this.
+
 ## Style Guide
 
 You can use clang-format to automatically format your files:
@@ -47,34 +83,4 @@ clang-format -i -style=file tool/net/redbean.c
 ```
 
 If you use Emacs this can be automated on save for Cosmopolitan using
-[tool/emacs/cosmo-format.el]([tool/emacs/cosmo-format.el]).
-
-### Source Files
-
-- Must use include paths relative to the root of the repository
-- Must have comment at top of file documenting copyright and license
-- Must have notice embedding if not owned by Justine (exception: tests)
-- May use language extensions that are supported by both GCC and Clang
-- Should use Google indentation (otherwise use `/* clang-format off */`)
-- Should use asm() instead of compiler APIs (exception: ctz, clz, memcpy)
-
-### Header Files
-
-- Must not have copyright or license comments
-- Must have once guards (otherwise change `.h` to `.inc`)
-- Must be ANSI C89 compatible to be included in the amalgamation header
-- Must include its dependencies (exception: libc/integral/normalize.inc)
-- Must not define objects (i.e. `cc -c -xc foo.h` will produce empty `.o`)
-- Should not use typedefs
-- Should not use forward declarations
-- Should not include documentation comments
-- Should not include parameter names in prototypes
-- Should not pose problems if included by C++ or Assembly sources
-- Should not declare non-ANSI code, at all, when the user requests ANSI
-
-### Build Config
-
-- Must not write files outside `o/`
-- Must not communicate with Internet
-- Must not depend on system libraries
-- Must not depend on system commands (exception: sh, make, gzip, zip)
+[tool/emacs/cosmo-format.el](tool/emacs/cosmo-format.el).

Makefile

@@ -149,9 +149,9 @@ export MODE
 export SOURCE_DATE_EPOCH
 export TMPDIR
 
-COSMOCC = .cosmocc/3.3
+COSMOCC = .cosmocc/3.3.1
 TOOLCHAIN = $(COSMOCC)/bin/$(ARCH)-linux-cosmo-
-DOWNLOAD := $(shell build/download-cosmocc.sh $(COSMOCC) 3.3 d26ec8f4e48f6db004fc6a9677c7ff3b50c3b21e936e9393158aa2ed51b0b549)
+DOWNLOAD := $(shell build/download-cosmocc.sh $(COSMOCC) 3.3.1 079622f9772c19e24db4bfa71c1ed90b7d436e9a4b6f4c579b525df5a48fe620)
 
 AS = $(TOOLCHAIN)as
 CC = $(TOOLCHAIN)gcc
@@ -206,7 +206,7 @@ endif
 .UNVEIL +=					\
 	libc/integral				\
 	libc/stdbool.h				\
-		rwc:/dev/shm				\
+	rwc:/dev/shm				\
 	rx:.cosmocc				\
 	rx:build/bootstrap			\
 	r:build/portcosmo.h			\
@@ -297,6 +297,7 @@ include third_party/nsync/testing/BUILD.mk
 include libc/testlib/BUILD.mk
 include tool/viz/lib/BUILD.mk
 include tool/args/BUILD.mk
+include test/math/BUILD.mk
 include test/posix/BUILD.mk
 include test/libcxx/BUILD.mk
 include test/tool/args/BUILD.mk
@@ -555,9 +556,9 @@ o/cosmopolitan.html: private .UNSANDBOXED = 1
 o/cosmopolitan.html:							\
 		o/$(MODE)/third_party/chibicc/chibicc.com.dbg		\
 		$(filter-out %.s,$(foreach x,$(COSMOPOLITAN_OBJECTS),$($(x)_SRCS)))	\
-		$(filter-out %.cc,$(SRCS))				\
+		$(filter-out %.cpp,$(filter-out %.cc,$(SRCS)))				\
 		$(HDRS)
-	$(file >$(TMPDIR)/$(subst /,_,$@),$(filter-out %.cc,$(filter-out %.s,$(foreach x,$(COSMOPOLITAN_OBJECTS),$($(x)_SRCS)))))
+	$(file >$(TMPDIR)/$(subst /,_,$@),$(filter-out %.cpp,$(filter-out %.cc,$(filter-out %.s,$(foreach x,$(COSMOPOLITAN_OBJECTS),$($(x)_SRCS))))))
 	o/$(MODE)/third_party/chibicc/chibicc.com.dbg -J		\
 		-fno-common -include libc/integral/normalize.inc -o $@	\
 		-DCOSMO @$(TMPDIR)/$(subst /,_,$@)

ape/ape.lds

@@ -528,6 +528,7 @@ SECTIONS {
     *(.piro.data.sort.iat.*)
 #endif
     *(__patchable_function_entries)
+    *(.note.gnu.property)
     *(__mcount_loc)
     *(.rela.dyn)
     *(.discard)

build/definitions.mk

@@ -95,7 +95,6 @@ DEFAULT_CCFLAGS +=							\
 DEFAULT_COPTS ?=							\
 	-fno-ident							\
 	-fno-common							\
-	-fno-math-errno							\
 	-fno-gnu-unique							\
 	-fstrict-aliasing						\
 	-fstrict-overflow						\

libc/calls/calls.h

@@ -208,6 +208,7 @@ int execvpe(const char *, char *const[], char *const[]) libcesque;
 int euidaccess(const char *, int) libcesque;
 int eaccess(const char *, int) libcesque;
 int madvise(void *, uint64_t, int) libcesque;
+int getcpu(unsigned *, unsigned *) libcesque;
 #endif
 
 #ifdef _COSMO_SOURCE

libc/calls/clock_nanosleep-cosmo.c

@@ -0,0 +1,84 @@
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8                               :vi │
+╞══════════════════════════════════════════════════════════════════════════════╡
+│ Copyright 2024 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/calls/struct/timespec.h"
+#include "libc/errno.h"
+#include "libc/runtime/clktck.h"
+#include "libc/runtime/runtime.h"
+#include "libc/sysv/consts/clock.h"
+#include "libc/sysv/consts/timer.h"
+
+/**
+ * Sleeps with higher accuracy at the cost of cpu.
+ */
+int cosmo_clock_nanosleep(int clock, int flags, const struct timespec *req,
+                          struct timespec *rem) {
+
+  // pick clocks
+  int time_clock;
+  int sleep_clock;
+  if (clock == CLOCK_REALTIME ||  //
+      clock == CLOCK_REALTIME_PRECISE) {
+    time_clock = clock;
+    sleep_clock = CLOCK_REALTIME;
+  } else if (clock == CLOCK_MONOTONIC ||  //
+             clock == CLOCK_MONOTONIC_PRECISE) {
+    time_clock = clock;
+    sleep_clock = CLOCK_MONOTONIC;
+  } else if (clock == CLOCK_REALTIME_COARSE ||  //
+             clock == CLOCK_REALTIME_FAST) {
+    return sys_clock_nanosleep(CLOCK_REALTIME, flags, req, rem);
+  } else if (clock == CLOCK_MONOTONIC_COARSE ||  //
+             clock == CLOCK_MONOTONIC_FAST) {
+    return sys_clock_nanosleep(CLOCK_MONOTONIC, flags, req, rem);
+  } else {
+    return sys_clock_nanosleep(clock, flags, req, rem);
+  }
+
+  // sleep bulk of time in kernel
+  struct timespec start, deadline, remain, waitfor, now;
+  struct timespec quantum = timespec_fromnanos(1000000000 / CLK_TCK);
+  clock_gettime(time_clock, &start);
+  deadline = flags & TIMER_ABSTIME ? *req : timespec_add(start, *req);
+  if (timespec_cmp(start, deadline) >= 0) return 0;
+  remain = timespec_sub(deadline, start);
+  if (timespec_cmp(remain, quantum) > 0) {
+    waitfor = timespec_sub(remain, quantum);
+    if (sys_clock_nanosleep(sleep_clock, 0, &waitfor, rem) == -1) {
+      if (!flags && rem && errno == EINTR) {
+        *rem = timespec_add(*rem, quantum);
+      }
+      return -1;
+    }
+  }
+
+  // spin through final scheduling quantum
+  int rc = 0;
+  ftrace_enabled(-1);
+  do {
+    if (_check_cancel()) {
+      rc = -1;
+      break;
+    }
+    clock_gettime(time_clock, &now);
+  } while (timespec_cmp(now, deadline) < 0);
+  ftrace_enabled(+1);
+  return rc;
+}

libc/calls/clock_nanosleep-sys.c

@@ -0,0 +1,61 @@
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8                               :vi │
+╞══════════════════════════════════════════════════════════════════════════════╡
+│ Copyright 2024 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/cp.internal.h"
+#include "libc/calls/struct/timespec.h"
+#include "libc/calls/struct/timespec.internal.h"
+#include "libc/dce.h"
+#include "libc/errno.h"
+#include "libc/intrin/describeflags.internal.h"
+#include "libc/intrin/strace.internal.h"
+#include "libc/intrin/weaken.h"
+#include "libc/sysv/errfuns.h"
+#include "libc/thread/thread.h"
+
+int sys_clock_nanosleep(int clock, int flags,  //
+                        const struct timespec *req, struct timespec *rem) {
+  int rc;
+  BEGIN_CANCELATION_POINT;
+  if (IsLinux() || IsFreebsd() || IsNetbsd()) {
+    rc = __sys_clock_nanosleep(clock, flags, req, rem);
+  } else if (IsXnu()) {
+    rc = sys_clock_nanosleep_xnu(clock, flags, req, rem);
+  } else if (IsOpenbsd()) {
+    rc = sys_clock_nanosleep_openbsd(clock, flags, req, rem);
+  } else if (IsWindows()) {
+    rc = sys_clock_nanosleep_nt(clock, flags, req, rem);
+  } else {
+    rc = enosys();
+  }
+  if (rc > 0) {
+    errno = rc;
+    rc = -1;
+  }
+  // system call support might not detect cancelation on bsds
+  if (rc == -1 && errno == EINTR &&      //
+      _weaken(pthread_testcancel_np) &&  //
+      _weaken(pthread_testcancel_np)()) {
+    rc = ecanceled();
+  }
+  END_CANCELATION_POINT;
+  STRACE("sys_clock_nanosleep(%s, %s, %s, [%s]) → %d% m",
+         DescribeClockName(clock), DescribeSleepFlags(flags),
+         DescribeTimespec(0, req), DescribeTimespec(rc, rem), rc);
+  return rc;
+}

libc/calls/clock_nanosleep.c

@@ -16,112 +16,10 @@
 │ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR             │
 │ PERFORMANCE OF THIS SOFTWARE.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/assert.h"
-#include "libc/calls/cp.internal.h"
-#include "libc/calls/internal.h"
 #include "libc/calls/struct/timespec.h"
-#include "libc/calls/struct/timespec.internal.h"
 #include "libc/dce.h"
 #include "libc/errno.h"
-#include "libc/intrin/describeflags.internal.h"
-#include "libc/intrin/kprintf.h"
-#include "libc/intrin/strace.internal.h"
-#include "libc/intrin/weaken.h"
-#include "libc/runtime/clktck.h"
-#include "libc/runtime/runtime.h"
-#include "libc/sysv/consts/clock.h"
 #include "libc/sysv/consts/timer.h"
-#include "libc/sysv/errfuns.h"
-#include "libc/thread/thread.h"
-
-static int sys_clock_nanosleep(int clock, int flags,  //
-                               const struct timespec *req,
-                               struct timespec *rem) {
-  int rc;
-  BEGIN_CANCELATION_POINT;
-  if (IsLinux() || IsFreebsd() || IsNetbsd()) {
-    rc = __sys_clock_nanosleep(clock, flags, req, rem);
-  } else if (IsXnu()) {
-    rc = sys_clock_nanosleep_xnu(clock, flags, req, rem);
-  } else if (IsOpenbsd()) {
-    rc = sys_clock_nanosleep_openbsd(clock, flags, req, rem);
-  } else if (IsWindows()) {
-    rc = sys_clock_nanosleep_nt(clock, flags, req, rem);
-  } else {
-    rc = enosys();
-  }
-  if (rc > 0) {
-    errno = rc;
-    rc = -1;
-  }
-  // system call support might not detect cancelation on bsds
-  if (rc == -1 && errno == EINTR &&      //
-      _weaken(pthread_testcancel_np) &&  //
-      _weaken(pthread_testcancel_np)()) {
-    rc = ecanceled();
-  }
-  END_CANCELATION_POINT;
-  STRACE("sys_clock_nanosleep(%s, %s, %s, [%s]) → %d% m",
-         DescribeClockName(clock), DescribeSleepFlags(flags),
-         DescribeTimespec(0, req), DescribeTimespec(rc, rem), rc);
-  return rc;
-}
-
-static int cosmo_clock_nanosleep(int clock, int flags,
-                                 const struct timespec *req,
-                                 struct timespec *rem) {
-
-  // pick clocks
-  int time_clock;
-  int sleep_clock;
-  if (clock == CLOCK_REALTIME ||  //
-      clock == CLOCK_REALTIME_PRECISE) {
-    time_clock = clock;
-    sleep_clock = CLOCK_REALTIME;
-  } else if (clock == CLOCK_MONOTONIC ||  //
-             clock == CLOCK_MONOTONIC_PRECISE) {
-    time_clock = clock;
-    sleep_clock = CLOCK_MONOTONIC;
-  } else if (clock == CLOCK_REALTIME_COARSE ||  //
-             clock == CLOCK_REALTIME_FAST) {
-    return sys_clock_nanosleep(CLOCK_REALTIME, flags, req, rem);
-  } else if (clock == CLOCK_MONOTONIC_COARSE ||  //
-             clock == CLOCK_MONOTONIC_FAST) {
-    return sys_clock_nanosleep(CLOCK_MONOTONIC, flags, req, rem);
-  } else {
-    return sys_clock_nanosleep(clock, flags, req, rem);
-  }
-
-  // sleep bulk of time in kernel
-  struct timespec start, deadline, remain, waitfor, now;
-  struct timespec quantum = timespec_fromnanos(1000000000 / CLK_TCK);
-  unassert(!clock_gettime(time_clock, &start));
-  deadline = flags & TIMER_ABSTIME ? *req : timespec_add(start, *req);
-  if (timespec_cmp(start, deadline) >= 0) return 0;
-  remain = timespec_sub(deadline, start);
-  if (timespec_cmp(remain, quantum) > 0) {
-    waitfor = timespec_sub(remain, quantum);
-    if (sys_clock_nanosleep(sleep_clock, 0, &waitfor, rem) == -1) {
-      if (!flags && rem && errno == EINTR) {
-        *rem = timespec_add(*rem, quantum);
-      }
-      return -1;
-    }
-  }
-
-  // spin through final scheduling quantum
-  int rc = 0;
-  ftrace_enabled(-1);
-  do {
-    if (_check_cancel()) {
-      rc = -1;
-      break;
-    }
-    unassert(!clock_gettime(time_clock, &now));
-  } while (timespec_cmp(now, deadline) < 0);
-  ftrace_enabled(+1);
-  return rc;
-}
 
 /**
  * Sleeps for particular amount of time.
@@ -157,10 +55,8 @@ static int cosmo_clock_nanosleep(int clock, int flags,
  * on OpenBSD it's good; on XNU it's bad; and on Windows it's ugly.
  *
  * @param clock may be
- *     - `CLOCK_REALTIME` to have nanosecond-accurate wall time sleeps
- *     - `CLOCK_REALTIME_COARSE` to not spin through scheduler quantum
- *     - `CLOCK_MONOTONIC` to base the sleep off the monotinic clock
- *     - `CLOCK_MONOTONIC_COARSE` to once again not do userspace spin
+ *     - `CLOCK_REALTIME`
+ *     - `CLOCK_MONOTONIC`
  * @param flags can be 0 for relative and `TIMER_ABSTIME` for absolute
  * @param req can be a relative or absolute time, depending on `flags`
  * @param rem shall be updated with the remainder of unslept time when
@@ -193,7 +89,7 @@ errno_t clock_nanosleep(int clock, int flags,        //
     return EINVAL;
   }
   errno_t old = errno;
-  int rc = cosmo_clock_nanosleep(clock, flags, req, rem);
+  int rc = sys_clock_nanosleep(clock, flags, req, rem);
   errno_t err = !rc ? 0 : errno;
   errno = old;
   return err;

libc/calls/getcpu.c

@@ -0,0 +1,76 @@
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8                               :vi │
+╞══════════════════════════════════════════════════════════════════════════════╡
+│ Copyright 2024 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/struct/cpuset.h"
+#include "libc/calls/syscall_support-nt.internal.h"
+#include "libc/dce.h"
+#include "libc/errno.h"
+#include "libc/nexgen32e/rdtscp.h"
+#include "libc/nexgen32e/x86feature.h"
+#include "libc/nt/struct/processornumber.h"
+#include "libc/nt/synchronization.h"
+#include "libc/runtime/syslib.internal.h"
+#include "libc/sysv/errfuns.h"
+
+int sys_getcpu(unsigned *opt_cpu, unsigned *opt_node, void *tcache);
+
+int getcpu(unsigned *out_opt_cpu, unsigned *out_opt_node) {
+  unsigned cpu;
+  unsigned node;
+  if (X86_HAVE(RDTSCP)) {
+    unsigned tsc_aux;
+    rdtscp(&tsc_aux);
+    cpu = TSC_AUX_CORE(tsc_aux);
+    node = TSC_AUX_NODE(tsc_aux);
+  } else if (IsWindows()) {
+    struct NtProcessorNumber pn;
+    GetCurrentProcessorNumberEx(&pn);
+    cpu = 64 * pn.Group + pn.Number;
+    unsigned short node16;
+    if (GetNumaProcessorNodeEx(&pn, &node16)) {
+      node = node16;
+    } else {
+      return __winerr();
+    }
+  } else if (IsXnuSilicon()) {
+    if (__syslib->__version >= 9) {
+      size_t cpu64;
+      errno_t err = __syslib->__pthread_cpu_number_np(&cpu64);
+      if (!err) {
+        cpu = cpu64;
+        node = 0;
+      } else {
+        errno = err;
+        return -1;
+      }
+    } else {
+      return enosys();
+    }
+  } else {
+    int rc = sys_getcpu(&cpu, &node, 0);
+    if (rc == -1) return -1;
+  }
+  if (out_opt_cpu) {
+    *out_opt_cpu = cpu;
+  }
+  if (out_opt_node) {
+    *out_opt_node = node;
+  }
+  return 0;
+}

libc/calls/struct/timespec.h

@@ -19,8 +19,8 @@ int timespec_getres(struct timespec *, int) libcesque;
 int timespec_get(struct timespec *, int) libcesque;
 
 #ifdef _COSMO_SOURCE
-/* cosmopolitan libc's non-posix timespec library
-   removed by default due to emacs codebase clash */
+int sys_clock_nanosleep(int, int, const struct timespec *, struct timespec *);
+int cosmo_clock_nanosleep(int, int, const struct timespec *, struct timespec *);
 #define timespec_zero ((struct timespec){0})
 #define timespec_max  ((struct timespec){0x7fffffffffffffff, 999999999})
 libcesque int timespec_cmp(struct timespec, struct timespec) pureconst;

libc/dce.h

@@ -1,6 +1,6 @@
+#ifdef _COSMO_SOURCE
 #ifndef COSMOPOLITAN_LIBC_DCE_H_
 #define COSMOPOLITAN_LIBC_DCE_H_
-#ifdef _COSMO_SOURCE
 /*─────────────────────────────────────────────────────────────────────────────╗
 │ cosmopolitan § autotune » dead code elimination                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
@@ -125,5 +125,5 @@ int IsQemuUser(void);
 
 COSMOPOLITAN_C_END_
 #endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */
-#endif /* _COSMO_SOURCE */
 #endif /* COSMOPOLITAN_LIBC_DCE_H_ */
+#endif /* _COSMO_SOURCE */

libc/dlopen/dlfcn.h

@@ -16,13 +16,11 @@ void *dlopen(const char *, int) libcesque;
 void *dlsym(void *, const char *) libcesque;
 int dlclose(void *) libcesque;
 
-#ifdef _COSMO_SOURCE
 char *cosmo_dlerror(void) libcesque;
 void *cosmo_dlopen(const char *, int) libcesque;
 void *cosmo_dlsym(void *, const char *) libcesque;
 void *cosmo_dltramp(void *) libcesque;
 int cosmo_dlclose(void *) libcesque;
-#endif
 
 COSMOPOLITAN_C_END_
 #endif /* COSMOPOLITAN_LIBC_DLFCN_H_ */

libc/fmt/atoi.c

@@ -16,44 +16,28 @@
 │ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR             │
 │ PERFORMANCE OF THIS SOFTWARE.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/errno.h"
 #include "libc/fmt/conv.h"
-#include "libc/limits.h"
-#include "libc/stdckdint.h"
 #include "libc/str/str.h"
 
 /**
- * Decodes decimal integer from ASCII string.
+ * Turns string into int.
  *
- *     atoi 10⁸              22𝑐         7𝑛𝑠
- *     strtol 10⁸            37𝑐        12𝑛𝑠
- *     strtoul 10⁸           35𝑐        11𝑛𝑠
- *     wcstol 10⁸            30𝑐        10𝑛𝑠
- *     wcstoul 10⁸           30𝑐        10𝑛𝑠
- *     strtoimax 10⁸         80𝑐        26𝑛𝑠
- *     strtoumax 10⁸         78𝑐        25𝑛𝑠
- *     wcstoimax 10⁸         77𝑐        25𝑛𝑠
- *     wcstoumax 10⁸         76𝑐        25𝑛𝑠
+ * Decimal is the only radix supported. Leading whitespace (as specified
+ * by the isspace() function) is skipped over. Unlike strtol(), the atoi
+ * function has undefined behavior on error and it never changes `errno`
  *
- * @param s is a non-null nul-terminated string
+ * @param nptr is a non-null nul-terminated string
  * @return the decoded signed saturated integer
- * @raise ERANGE on overflow
  */
-int atoi(const char *s) {
+int atoi(const char *nptr) {
   int x, c, d;
-  do c = *s++;
-  while (c == ' ' || c == '\t');
+  do c = *nptr++;
+  while (isspace(c));
   d = c == '-' ? -1 : 1;
-  if (c == '-' || c == '+') c = *s++;
-  for (x = 0; isdigit(c); c = *s++) {
-    if (ckd_mul(&x, x, 10) || ckd_add(&x, x, (c - '0') * d)) {
-      errno = ERANGE;
-      if (d > 0) {
-        return INT_MAX;
-      } else {
-        return INT_MIN;
-      }
-    }
+  if (c == '-' || c == '+') c = *nptr++;
+  for (x = 0; isdigit(c); c = *nptr++) {
+    x *= 10;
+    x += (c - '0') * d;
   }
   return x;
 }

libc/fmt/atol.c

@@ -16,34 +16,29 @@
 │ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR             │
 │ PERFORMANCE OF THIS SOFTWARE.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/errno.h"
 #include "libc/fmt/conv.h"
-#include "libc/limits.h"
-#include "libc/stdckdint.h"
 #include "libc/str/str.h"
 
 /**
- * Decodes decimal integer from ASCII string.
+ * Turns string into long.
  *
- * @param s is a non-null nul-terminated string
+ * Decimal is the only radix supported. Leading whitespace (as specified
+ * by the isspace() function) is skipped over. Unlike strtol(), the atoi
+ * function has undefined behavior on error and it never changes `errno`
+ *
+ * @param nptr is a non-null nul-terminated string
  * @return the decoded signed saturated integer
  */
-long atol(const char *s) {
+long atol(const char *nptr) {
   long x;
   int c, d;
-  do c = *s++;
-  while (c == ' ' || c == '\t');
+  do c = *nptr++;
+  while (isspace(c));
   d = c == '-' ? -1 : 1;
-  if (c == '-' || c == '+') c = *s++;
-  for (x = 0; isdigit(c); c = *s++) {
-    if (ckd_mul(&x, x, 10) || ckd_add(&x, x, (c - '0') * d)) {
-      errno = ERANGE;
-      if (d > 0) {
-        return LONG_MAX;
-      } else {
-        return LONG_MIN;
-      }
-    }
+  if (c == '-' || c == '+') c = *nptr++;
+  for (x = 0; isdigit(c); c = *nptr++) {
+    x *= 10;
+    x += (c - '0') * d;
   }
   return x;
 }

libc/fmt/itoa.h

@@ -1,6 +1,6 @@
+#ifdef _COSMO_SOURCE
 #ifndef COSMOPOLITAN_LIBC_FMT_ITOA_H_
 #define COSMOPOLITAN_LIBC_FMT_ITOA_H_
-#ifdef _COSMO_SOURCE
 COSMOPOLITAN_C_START_
 
 #define LengthInt64           __LengthInt64
@@ -42,5 +42,5 @@ libcesque size_t uint64toarray_fixed16(uint64_t, char[hasatleast 17], uint8_t);
 libcesque size_t uint64toarray_radix8(uint64_t, char[hasatleast 24]);
 
 COSMOPOLITAN_C_END_
-#endif /* _COSMO_SOURCE */
 #endif /* COSMOPOLITAN_LIBC_FMT_ITOA_H_ */
+#endif /* _COSMO_SOURCE */

libc/fmt/leb128.h

@@ -1,6 +1,6 @@
+#ifdef _COSMO_SOURCE
 #ifndef COSMOPOLITAN_LIBC_FMT_LEB128_H_
 #define COSMOPOLITAN_LIBC_FMT_LEB128_H_
-#ifdef _COSMO_SOURCE
 COSMOPOLITAN_C_START_
 
 #define sleb64   __sleb64
@@ -16,5 +16,5 @@ int unzleb64(const char *, size_t, int64_t *) libcesque;
 int unuleb64(const char *, size_t, uint64_t *) libcesque;
 
 COSMOPOLITAN_C_END_
-#endif /* _COSMO_SOURCE */
 #endif /* COSMOPOLITAN_LIBC_FMT_LEB128_H_ */
+#endif /* _COSMO_SOURCE */

libc/integral/normalize.inc

@@ -4,7 +4,7 @@
 
 #define __COSMOPOLITAN_MAJOR__ 3
 #define __COSMOPOLITAN_MINOR__ 3
-#define __COSMOPOLITAN_PATCH__ 0
+#define __COSMOPOLITAN_PATCH__ 1
 #define __COSMOPOLITAN__                                                   \
   (100000000 * __COSMOPOLITAN_MAJOR__ + 1000000 * __COSMOPOLITAN_MINOR__ + \
    __COSMOPOLITAN_PATCH__)

libc/intrin/aarch64/memchr.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/memcpy.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │
@@ -80,11 +80,12 @@ ENTRY (__memcpy_aarch64_simd)
 	PTR_ARG (1)
 	SIZE_ARG (2)
 	add	srcend, src, count
-	add	dstend, dstin, count
 	cmp	count, 128
 	b.hi	L(copy_long)
+	add	dstend, dstin, count
 	cmp	count, 32
 	b.hi	L(copy32_128)
+	nop
 
 	/* Small copies: 0..32 bytes.  */
 	cmp	count, 16
@@ -95,6 +96,18 @@ ENTRY (__memcpy_aarch64_simd)
 	str	B_q, [dstend, -16]
 	ret
 
+	.p2align 4
+	/* Medium copies: 33..128 bytes.  */
+L(copy32_128):
+	ldp	A_q, B_q, [src]
+	ldp	C_q, D_q, [srcend, -32]
+	cmp	count, 64
+	b.hi	L(copy128)
+	stp	A_q, B_q, [dstin]
+	stp	C_q, D_q, [dstend, -32]
+	ret
+
+	.p2align 4
 	/* Copy 8-15 bytes.  */
 L(copy16):
 	tbz	count, 3, L(copy8)
@@ -104,7 +117,6 @@ L(copy16):
 	str	A_h, [dstend, -8]
 	ret
 
-	.p2align 3
 	/* Copy 4-7 bytes.  */
 L(copy8):
 	tbz	count, 2, L(copy4)
@@ -114,6 +126,19 @@ L(copy8):
 	str	B_lw, [dstend, -4]
 	ret
 
+	/* Copy 65..128 bytes.  */
+L(copy128):
+	ldp	E_q, F_q, [src, 32]
+	cmp	count, 96
+	b.ls	L(copy96)
+	ldp	G_q, H_q, [srcend, -64]
+	stp	G_q, H_q, [dstend, -64]
+L(copy96):
+	stp	A_q, B_q, [dstin]
+	stp	E_q, F_q, [dstin, 32]
+	stp	C_q, D_q, [dstend, -32]
+	ret
+
 	/* Copy 0..3 bytes using a branchless sequence.  */
 L(copy4):
 	cbz	count, L(copy0)
@@ -127,33 +152,11 @@ L(copy4):
 L(copy0):
 	ret
 
-	.p2align 4
-	/* Medium copies: 33..128 bytes.  */
-L(copy32_128):
-	ldp	A_q, B_q, [src]
-	ldp	C_q, D_q, [srcend, -32]
-	cmp	count, 64
-	b.hi	L(copy128)
-	stp	A_q, B_q, [dstin]
-	stp	C_q, D_q, [dstend, -32]
-	ret
-
-	.p2align 4
-	/* Copy 65..128 bytes.  */
-L(copy128):
-	ldp	E_q, F_q, [src, 32]
-	cmp	count, 96
-	b.ls	L(copy96)
-	ldp	G_q, H_q, [srcend, -64]
-	stp	G_q, H_q, [dstend, -64]
-L(copy96):
-	stp	A_q, B_q, [dstin]
-	stp	E_q, F_q, [dstin, 32]
-	stp	C_q, D_q, [dstend, -32]
-	ret
-
+	.p2align 3
 	/* Copy more than 128 bytes.  */
 L(copy_long):
+	add	dstend, dstin, count
+
 	/* Use backwards copy if there is an overlap.  */
 	sub	tmp1, dstin, src
 	cmp	tmp1, count
@@ -190,6 +193,9 @@ L(copy64_from_end):
 	stp	A_q, B_q, [dstend, -32]
 	ret
 
+	.p2align 4
+	nop
+
 	/* Large backwards copy for overlapping copies.
 	   Copy 16 bytes and then align srcend to 16-byte alignment.  */
 L(copy_long_backwards):

libc/intrin/aarch64/memrchr.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/memset.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/stpcpy.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/strchr.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/strchrnul.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/strcmp.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/strcpy.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/strlen.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/strncmp.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/strnlen.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/aarch64/strrchr.S

@@ -3,7 +3,7 @@
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │

libc/intrin/describebacktrace.c

@@ -17,6 +17,8 @@
 │ PERFORMANCE OF THIS SOFTWARE.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
 #include "libc/intrin/describebacktrace.internal.h"
+#include "libc/intrin/kprintf.h"
+#include "libc/intrin/weaken.h"
 #include "libc/log/libfatal.internal.h"
 #include "libc/nexgen32e/stackframe.h"
 
@@ -28,6 +30,9 @@ dontinstrument const char *(DescribeBacktrace)(char buf[N],
   char *pe = p + N;
   bool gotsome = false;
   while (fr) {
+    if (_weaken(kisdangerous) && _weaken(kisdangerous)(fr)) {
+      break;
+    }
     if (p + 16 + 1 + 1 <= pe) {
       if (gotsome) {
         *p++ = ' ';

libc/intrin/dll.h

@@ -1,6 +1,6 @@
+#ifdef _COSMO_SOURCE
 #ifndef COSMOPOLITAN_LIBC_INTRIN_DLL_H_
 #define COSMOPOLITAN_LIBC_INTRIN_DLL_H_
-#ifdef _COSMO_SOURCE
 #define dll_make_first   __dll_make_first
 #define dll_make_last    __dll_make_last
 #define dll_remove       __dll_remove
@@ -55,5 +55,5 @@ void dll_make_first(struct Dll **, struct Dll *) paramsnonnull((1)) libcesque;
 void dll_splice_after(struct Dll *, struct Dll *) paramsnonnull((1)) libcesque;
 
 COSMOPOLITAN_C_END_
-#endif /* _COSMO_SOURCE */
 #endif /* COSMOPOLITAN_LIBC_INTRIN_DLL_H_ */
+#endif /* _COSMO_SOURCE */

libc/intrin/fbclibm.c

@@ -0,0 +1,2 @@
+__notice(freebsd_complex_notice, "FreeBSD Complex Math (BSD-2 License)\n\
+Copyright (c) 2012 Stephen Montgomery-Smith <[email protected]>");

libc/intrin/kprintf.h

@@ -1,6 +1,6 @@
+#ifdef _COSMO_SOURCE
 #ifndef COSMOPOLITAN_LIBC_INTRIN_KPRINTF_H_
 #define COSMOPOLITAN_LIBC_INTRIN_KPRINTF_H_
-#ifdef _COSMO_SOURCE
 
 #define klog         __klog
 #define kprintf      __kprintf
@@ -54,5 +54,5 @@ void uvprintf(const char *, va_list) libcesque;
 
 COSMOPOLITAN_C_END_
 #endif /* !(__ASSEMBLER__ + __LINKER__ + 0) */
-#endif /* _COSMO_SOURCE */
 #endif /* COSMOPOLITAN_LIBC_INTRIN_KPRINTF_H_ */
+#endif /* _COSMO_SOURCE */

libc/isystem/stdlib.h

@@ -14,6 +14,5 @@
 #include "libc/str/str.h"
 #include "libc/sysv/consts/exit.h"
 #include "libc/temp.h"
-#include "third_party/musl/crypt.h"
 #include "third_party/musl/rand48.h"
 #endif /* _STDLIB_H */

libc/isystem/unistd.h

@@ -10,12 +10,8 @@
 #include "libc/sysv/consts/o.h"
 #include "libc/sysv/consts/ok.h"
 #include "libc/time/time.h"
+#include "libc/unistd.h"
 #include "third_party/getopt/long1.h"
 #include "third_party/musl/crypt.h"
 #include "third_party/musl/lockf.h"
-
-#ifndef _CS_PATH
-#define _CS_PATH 0
-#endif
-
 #endif /* _UNISTD_H */

libc/log/log.h

@@ -1,7 +1,7 @@
+#ifdef _COSMO_SOURCE
 #ifndef COSMOPOLITAN_LIBC_LOG_LOG_H_
 #define COSMOPOLITAN_LIBC_LOG_LOG_H_
 #include "libc/stdio/stdio.h"
-#ifdef _COSMO_SOURCE
 
 #define kLogFatal   0
 #define kLogError   1
@@ -224,5 +224,5 @@ void vffatalf(ARGS, va_list) asm("vflogf") ATTRV relegated wontreturn libcesque;
 
 #endif /* __STRICT_ANSI__ */
 COSMOPOLITAN_C_END_
-#endif /* _COSMO_SOURCE */
 #endif /* COSMOPOLITAN_LIBC_LOG_LOG_H_ */
+#endif /* _COSMO_SOURCE */

libc/log/logfile.c

@@ -21,6 +21,6 @@
 
 FILE *__log_file;
 
-__attribute__((__constructor__(60))) static void loginit(void) {
+__attribute__((__constructor__(61))) static void loginit(void) {
   __log_file = stderr;
 }

libc/math.h

@@ -9,15 +9,47 @@
 #define M_LOG10E   0.43429448190325182765 /* log₁₀𝑒 */
 #define M_LN2      0.69314718055994530942 /* logₑ2 */
 #define M_LN10     2.30258509299404568402 /* logₑ10 */
-#define M_PI       3.14159265358979323846 /* pi */
-#define M_PI_2     1.57079632679489661923 /* pi/2 */
-#define M_PI_4     0.78539816339744830962 /* pi/4 */
-#define M_1_PI     0.31830988618379067154 /* 1/pi */
-#define M_2_PI     0.63661977236758134308 /* 2/pi */
-#define M_2_SQRTPI 1.12837916709551257390 /* 2/sqrt(pi) */
+#define M_PI       3.14159265358979323846 /* 𝜋 */
+#define M_PI_2     1.57079632679489661923 /* 𝜋/2 */
+#define M_PI_4     0.78539816339744830962 /* 𝜋/4 */
+#define M_1_PI     0.31830988618379067154 /* 1/𝜋 */
+#define M_2_PI     0.63661977236758134308 /* 2/𝜋 */
+#define M_2_SQRTPI 1.12837916709551257390 /* 2/sqrt(𝜋) */
 #define M_SQRT2    1.41421356237309504880 /* sqrt(2) */
 #define M_SQRT1_2  0.70710678118654752440 /* 1/sqrt(2) */
 
+#if defined(_GNU_SOURCE) || defined(_COSMO_SOURCE)
+#define M_Ef        2.7182818284590452354f  /* 𝑒 */
+#define M_LOG2Ef    1.4426950408889634074f  /* log₂𝑒 */
+#define M_LOG10Ef   0.43429448190325182765f /* log₁₀𝑒 */
+#define M_LN2f      0.69314718055994530942f /* logₑ2 */
+#define M_LN10f     2.30258509299404568402f /* logₑ10 */
+#define M_PIf       3.14159265358979323846f /* 𝜋 */
+#define M_PI_2f     1.57079632679489661923f /* 𝜋/2 */
+#define M_PI_4f     0.78539816339744830962f /* 𝜋/4 */
+#define M_1_PIf     0.31830988618379067154f /* 1/𝜋 */
+#define M_2_PIf     0.63661977236758134308f /* 2/𝜋 */
+#define M_2_SQRTPIf 1.12837916709551257390f /* 2/sqrt(𝜋) */
+#define M_SQRT2f    1.41421356237309504880f /* sqrt(2) */
+#define M_SQRT1_2f  0.70710678118654752440f /* 1/sqrt(2) */
+#endif
+
+#if defined(_GNU_SOURCE) || defined(_COSMO_SOURCE)
+#define M_El        2.718281828459045235360287471352662498L /* 𝑒 */
+#define M_LOG2El    1.442695040888963407359924681001892137L /* log₂𝑒 */
+#define M_LOG10El   0.434294481903251827651128918916605082L /* log₁₀𝑒 */
+#define M_LN2l      0.693147180559945309417232121458176568L /* logₑ2 */
+#define M_LN10l     2.302585092994045684017991454684364208L /* logₑ10 */
+#define M_PIl       3.141592653589793238462643383279502884L /* 𝜋 */
+#define M_PI_2l     1.570796326794896619231321691639751442L /* 𝜋/2 */
+#define M_PI_4l     0.785398163397448309615660845819875721L /* 𝜋/4 */
+#define M_1_PIl     0.318309886183790671537767526745028724L /* 1/𝜋 */
+#define M_2_PIl     0.636619772367581343075535053490057448L /* 2/𝜋 */
+#define M_2_SQRTPIl 1.128379167095512573896158903121545172L /* 2/sqrt(𝜋) */
+#define M_SQRT2l    1.414213562373095048801688724209698079L /* sqrt(2) */
+#define M_SQRT1_2l  0.707106781186547524400844362104849039L /* 1/sqrt(2) */
+#endif
+
 #define DBL_DECIMAL_DIG   __DBL_DECIMAL_DIG__
 #define DBL_DIG           __DBL_DIG__
 #define DBL_EPSILON       __DBL_EPSILON__
@@ -76,6 +108,27 @@
 #define FP_ILOGB0    (-2147483647 - 1)
 #define FP_ILOGBNAN  (-2147483647 - 1)
 
+#define MATH_ERRNO     1
+#define MATH_ERREXCEPT 2
+
+#ifdef __FAST_MATH__
+#define math_errhandling 0
+#elif defined(__NO_MATH_ERRNO__)
+#define math_errhandling (MATH_ERREXCEPT)
+#else
+#define math_errhandling (MATH_ERRNO | MATH_ERREXCEPT)
+#endif
+
+#ifdef __FP_FAST_FMA
+#define FP_FAST_FMA 1
+#endif
+#ifdef __FP_FAST_FMAF
+#define FP_FAST_FMAF 1
+#endif
+#ifdef __FP_FAST_FMAL
+#define FP_FAST_FMAL 1
+#endif
+
 COSMOPOLITAN_C_START_
 
 #define NAN       __builtin_nanf("")

libc/mem/gc.h

@@ -1,6 +1,6 @@
+#ifdef _COSMO_SOURCE
 #ifndef COSMOPOLITAN_LIBC_MEM_GC_H_
 #define COSMOPOLITAN_LIBC_MEM_GC_H_
-#ifdef _COSMO_SOURCE
 COSMOPOLITAN_C_START_
 
 libcesque void _gc_free(void *);
@@ -19,5 +19,5 @@ libcesque void gclongjmp(void *, int) wontreturn;
   })
 
 COSMOPOLITAN_C_END_
-#endif /* _COSMO_SOURCE */
 #endif /* COSMOPOLITAN_LIBC_MEM_GC_H_ */
+#endif /* _COSMO_SOURCE */

libc/nexgen32e/stackframe.h

@@ -1,6 +1,6 @@
+#ifdef _COSMO_SOURCE
 #ifndef COSMOPOLITAN_LIBC_NEXGEN32E_STACKFRAME_H_
 #define COSMOPOLITAN_LIBC_NEXGEN32E_STACKFRAME_H_
-#ifdef _COSMO_SOURCE
 COSMOPOLITAN_C_START_
 
 struct StackFrame {
@@ -9,5 +9,5 @@ struct StackFrame {
 };
 
 COSMOPOLITAN_C_END_
-#endif /* _COSMO_SOURCE */
 #endif /* COSMOPOLITAN_LIBC_NEXGEN32E_STACKFRAME_H_ */
+#endif /* _COSMO_SOURCE */

libc/nt/dll.h

@@ -27,6 +27,7 @@ COSMOPOLITAN_C_START_
 ╚────────────────────────────────────────────────────────────────────────────│*/
 
 int64_t LoadLibrary(const char16_t *lpLibFileName);
+int64_t LoadLibraryA(const char *lpLibFileName);
 int64_t LoadLibraryEx(const char16_t *lpLibFileName, int64_t hFile,
                       uint32_t dwFlags);
 uint32_t GetModuleFileNameA(int64_t hModule, char *lpFilename, uint32_t nSize);

libc/nt/events.h

@@ -57,6 +57,9 @@ int32_t DeregisterEventSource(uint64_t handle);
 int64_t CreateEvent(const struct NtSecurityAttributes *opt_lpEventAttributes,
                     bool32 bManualReset, bool32 bInitialState,
                     const char16_t *opt_lpName);
+int64_t CreateEventA(const struct NtSecurityAttributes *opt_lpEventAttributes,
+                     bool32 bManualReset, bool32 bInitialState,
+                     const char *opt_lpName);
 int64_t CreateEventEx(const struct NtSecurityAttributes *lpEventAttributes,
                       const char16_t *lpName, uint32_t dwFlags,
                       uint32_t dwDesiredAccess);

libc/nt/kernel32/CreateEventA.S

@@ -0,0 +1,18 @@
+#include "libc/nt/codegen.h"
+.imp	kernel32,__imp_CreateEventA,CreateEventA
+
+	.text.windows
+	.ftrace1
+CreateEventA:
+	.ftrace2
+#ifdef __x86_64__
+	push	%rbp
+	mov	%rsp,%rbp
+	mov	__imp_CreateEventA(%rip),%rax
+	jmp	__sysv2nt
+#elif defined(__aarch64__)
+	mov	x0,#0
+	ret
+#endif
+	.endfn	CreateEventA,globl
+	.previous

libc/nt/kernel32/GetNumaProcessorNodeEx.S

@@ -0,0 +1,18 @@
+#include "libc/nt/codegen.h"
+.imp	kernel32,__imp_GetNumaProcessorNodeEx,GetNumaProcessorNodeEx
+
+	.text.windows
+	.ftrace1
+GetNumaProcessorNodeEx:
+	.ftrace2
+#ifdef __x86_64__
+	push	%rbp
+	mov	%rsp,%rbp
+	mov	__imp_GetNumaProcessorNodeEx(%rip),%rax
+	jmp	__sysv2nt
+#elif defined(__aarch64__)
+	mov	x0,#0
+	ret
+#endif
+	.endfn	GetNumaProcessorNodeEx,globl
+	.previous

libc/nt/kernel32/LoadLibraryA.S

@@ -0,0 +1,20 @@
+#include "libc/nt/codegen.h"
+.imp	kernel32,__imp_LoadLibraryA,LoadLibraryA
+
+	.text.windows
+	.ftrace1
+LoadLibraryA:
+	.ftrace2
+#ifdef __x86_64__
+	push	%rbp
+	mov	%rsp,%rbp
+	mov	%rdi,%rcx
+	sub	$32,%rsp
+	call	*__imp_LoadLibraryA(%rip)
+	leave
+#elif defined(__aarch64__)
+	mov	x0,#0
+#endif
+	ret
+	.endfn	LoadLibraryA,globl
+	.previous

libc/nt/master.sh

@@ -63,6 +63,7 @@ imp	'ConnectNamedPipe'					ConnectNamedPipe					kernel32	2
 imp	'ContinueDebugEvent'					ContinueDebugEvent					kernel32	3
 imp	'CopyFile'						CopyFileW						kernel32	3
 imp	'CreateEvent'						CreateEventW						kernel32	4
+imp	'CreateEventA'						CreateEventA						kernel32	4
 imp	'CreateEventEx'						CreateEventExW						kernel32	4
 imp	'CreateHardLink'					CreateHardLinkW						kernel32	3
 imp	'CreateIoCompletionPort'				CreateIoCompletionPort					kernel32	4
@@ -134,6 +135,7 @@ imp	'GetModuleHandle'					GetModuleHandleA					kernel32	1
 imp	'GetModuleHandleEx'					GetModuleHandleExW					kernel32	3
 imp	'GetModuleHandleW'					GetModuleHandleW					kernel32	1
 imp	'GetNamedPipeInfo'					GetNamedPipeInfo					kernel32	5
+imp	'GetNumaProcessorNodeEx'				GetNumaProcessorNodeEx					kernel32	2
 imp	'GetNumberOfConsoleInputEvents'				GetNumberOfConsoleInputEvents				kernel32	2
 imp	'GetNumberOfConsoleMouseButtons'			GetNumberOfConsoleMouseButtons				kernel32	1
 imp	'GetOverlappedResult'					GetOverlappedResult					kernel32	4
@@ -195,6 +197,7 @@ imp	'InitializeProcThreadAttributeList'			InitializeProcThreadAttributeList			ke
 imp	'InitializeSRWLock'					InitializeSRWLock					kernel32	1
 imp	'LeaveCriticalSection'					LeaveCriticalSection					kernel32	1
 imp	'LoadLibrary'						LoadLibraryW						kernel32	1
+imp	'LoadLibraryA'						LoadLibraryA						kernel32	1
 imp	'LoadLibraryEx'						LoadLibraryExW						kernel32	3
 imp	'LoadResource'						LoadResource						kernel32	2
 imp	'LocalFree'						LocalFree						kernel32	1

libc/nt/synchronization.h

@@ -118,6 +118,9 @@ bool32 GetSystemTimeAdjustment(uint32_t *lpTimeAdjustment,
 
 void GetCurrentProcessorNumberEx(struct NtProcessorNumber *out_ProcNumber);
 
+bool32 GetNumaProcessorNodeEx(const struct NtProcessorNumber *Processor,
+                              unsigned short *out_NodeNumber);
+
 #if ShouldUseMsabiAttribute()
 #include "libc/nt/thunk/synchronization.inc"
 #endif /* ShouldUseMsabiAttribute() */

libc/runtime/ftracer.c

@@ -51,7 +51,7 @@ static struct CosmoFtrace g_ftrace;
 
 __funline int GetNestingLevelImpl(struct StackFrame *frame) {
   int nesting = -2;
-  while (frame) {
+  while (frame && !kisdangerous(frame)) {
     ++nesting;
     frame = frame->next;
   }

libc/runtime/stack.h

@@ -1,6 +1,6 @@
+#ifdef _COSMO_SOURCE
 #ifndef COSMOPOLITAN_LIBC_RUNTIME_STACK_H_
 #define COSMOPOLITAN_LIBC_RUNTIME_STACK_H_
-#ifdef _COSMO_SOURCE
 
 /**
  * Returns preferred size and alignment of thread stack.
@@ -156,5 +156,5 @@ int FreeCosmoStack(void *) libcesque;
 
 COSMOPOLITAN_C_END_
 #endif /* GNU ELF */
-#endif /* _COSMO_SOURCE */
 #endif /* COSMOPOLITAN_LIBC_RUNTIME_STACK_H_ */
+#endif /* _COSMO_SOURCE */

libc/stdio/append.h

@@ -1,6 +1,6 @@
+#ifdef _COSMO_SOURCE
 #ifndef COSMOPOLITAN_LIBC_STDIO_APPEND_H_
 #define COSMOPOLITAN_LIBC_STDIO_APPEND_H_
-#ifdef _COSMO_SOURCE
 
 #define APPEND_COOKIE 21578
 
@@ -32,5 +32,5 @@ ssize_t kappendf(char **, const char *, ...) libcesque;
 ssize_t kvappendf(char **, const char *, va_list) libcesque;
 
 COSMOPOLITAN_C_END_
-#endif /* _COSMO_SOURCE */
 #endif /* COSMOPOLITAN_LIBC_STDIO_APPEND_H_ */
+#endif /* _COSMO_SOURCE */

libc/stdio/confstr.c

@@ -20,6 +20,7 @@
 #include "libc/stdio/stdio.h"
 #include "libc/str/str.h"
 #include "libc/sysv/errfuns.h"
+#include "libc/unistd.h"
 
 size_t confstr(int name, char *buf, size_t len) {
   if (name == _CS_PATH) {

libc/stdio/stdio.h

@@ -1,12 +1,11 @@
 #ifndef COSMOPOLITAN_LIBC_STDIO_H_
 #define COSMOPOLITAN_LIBC_STDIO_H_
 
-#define EOF      -1  /* end of file */
-#define WEOF     -1u /* end of file (multibyte) */
-#define _IOFBF   0   /* fully buffered */
-#define _IOLBF   1   /* line buffered */
-#define _IONBF   2   /* no buffering */
-#define _CS_PATH 0
+#define EOF    -1  /* end of file */
+#define WEOF   -1u /* end of file (multibyte) */
+#define _IOFBF 0   /* fully buffered */
+#define _IOLBF 1   /* line buffered */
+#define _IONBF 2   /* no buffering */
 
 #define L_tmpnam     20
 #define L_ctermid    20

libc/stdio/syscall.c

@@ -43,8 +43,15 @@ long syscall(long number, ...) {
       size_t buflen = va_arg(va, size_t);
       unsigned flags = va_arg(va, unsigned);
       va_end(va);
-      ssize_t rc = getrandom(buf, buflen, flags);
-      return rc;
+      return getrandom(buf, buflen, flags);
+    }
+    case SYS_getcpu: {
+      va_list va;
+      va_start(va, number);
+      unsigned *cpu = va_arg(va, unsigned *);
+      unsigned *node = va_arg(va, unsigned *);
+      va_end(va);
+      return getcpu(cpu, node);
     }
   }
 }

libc/stdio/syscall.h

@@ -2,8 +2,9 @@
 #define COSMOPOLITAN_LIBC_STDIO_SYSCALL_H_
 COSMOPOLITAN_C_START_
 
-#define SYS_gettid    186
-#define SYS_getrandom 318
+#define SYS_gettid    1
+#define SYS_getrandom 2
+#define SYS_getcpu    3
 
 long syscall(long, ...) libcesque;
 

libc/stdio/vcscanf.c

@@ -50,6 +50,12 @@
     }                                         \
     c;                                        \
   })
+#define UNBUFFER                \
+  ({                            \
+    if (c != -1) {              \
+      fpbuf[--fpbufcur] = '\0'; \
+    }                           \
+  })
 
 /**
  * String / file / stream decoder.
@@ -369,10 +375,11 @@ int __vcscanf(int callback(void *),    //
                   }
                 } while ((c = BUFFER) != -1 && c != ')');
                 if (c == ')') {
-                  c = BUFFER;
+                  c = READ;
                 }
                 goto GotFloatingPointNumber;
               } else {
+                UNBUFFER;
                 goto GotFloatingPointNumber;
               }
             } else {
@@ -410,9 +417,7 @@ int __vcscanf(int callback(void *),    //
                   goto Done;
                 }
               } else {
-                if (c != -1 && unget) {
-                  unget(c, arg);
-                }
+                UNBUFFER;
                 goto GotFloatingPointNumber;
               }
             } else {
@@ -465,13 +470,24 @@ int __vcscanf(int callback(void *),    //
         Continue:
           continue;
         Break:
-          if (c != -1 && unget) {
-            unget(c, arg);
-          }
+          UNBUFFER;
           break;
         } while ((c = BUFFER) != -1);
       GotFloatingPointNumber:
-        fp = strtod((char *)fpbuf, NULL);
+        /* An empty buffer can't be a valid float; don't even bother parsing. */
+        bool valid = fpbufcur > 0;
+        if (valid) {
+          char *ep;
+          fp = strtod((char *)fpbuf, &ep);
+          /* We should have parsed the whole buffer. */
+          valid = ep == (char *)fpbuf + fpbufcur;
+        }
+        free(fpbuf);
+        fpbuf = NULL;
+        fpbufcur = fpbufsize = 0;
+        if (!valid) {
+          goto Done;
+        }
         if (!discard) {
           ++items;
           void *out = va_arg(va, void *);
@@ -481,9 +497,6 @@ int __vcscanf(int callback(void *),    //
             *(double *)out = (double)fp;
           }
         }
-        free(fpbuf);
-        fpbuf = NULL;
-        fpbufcur = fpbufsize = 0;
         continue;
       ReportConsumed:
         n_ptr = va_arg(va, int *);
@@ -537,6 +550,11 @@ int __vcscanf(int callback(void *),    //
               if (!j && c == -1 && !items) {
                 items = -1;
                 goto Done;
+              } else if (rawmode && j != width) {
+                /* The C standard says that %c "matches a sequence of characters of
+                 * **exactly** the number specified by the field width". If we have
+                 * fewer characters, what we've just read is invalid. */
+                goto Done;
               } else if (!rawmode && j < bufsize) {
                 if (charbytes == sizeof(char)) {
                   buf[j] = '\0';

libc/stdlib.h

@@ -6,5 +6,10 @@ char *fcvt(double, int, int *, int *) libcesque;
 char *ecvt(double, int, int *, int *) libcesque;
 char *gcvt(double, int, char *) libcesque;
 
+#if defined(_XOPEN_SOURCE) || defined(_GNU_SOURCE) || defined(_BSD_SOURCE) || \
+    defined(_COSMO_SOURCE)
+void setkey(const char *) libcesque;
+#endif
+
 COSMOPOLITAN_C_END_
 #endif /* COSMOPOLITAN_LIBC_STDLIB_H_ */

libc/tinymath/BUILD.mk

@@ -54,6 +54,7 @@ o/$(MODE)/libc/tinymath/loglq.o: private		\
 
 $(LIBC_TINYMATH_A_OBJS): private			\
 		CFLAGS +=				\
+			-fmath-errno			\
 			-fsigned-zeros			\
 			-ftrapping-math			\
 			-frounding-math			\

libc/tinymath/LICENSE.optimized-routines

@@ -5,7 +5,7 @@ MIT OR Apache-2.0 WITH LLVM-exception
 MIT License
 -----------
 
-Copyright (c) 1999-2022, Arm Limited.
+Copyright (c) 2018-2024, Arm Limited.
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

libc/tinymath/__math_divzero.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/tinymath/internal.h"
-
-double __math_divzero(uint32_t sign)
-{
-	return fp_barrier(sign ? -1.0 : 1.0) / 0.0;
-}

libc/tinymath/__math_divzerof.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/tinymath/internal.h"
-
-float __math_divzerof(uint32_t sign)
-{
-	return fp_barrierf(sign ? -1.0f : 1.0f) / 0.0f;
-}

libc/tinymath/__math_invalid.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/tinymath/internal.h"
-
-double __math_invalid(double x)
-{
-	return (x - x) / (x - x);
-}

libc/tinymath/__math_invalidf.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/tinymath/internal.h"
-
-float __math_invalidf(float x)
-{
-	return (x - x) / (x - x);
-}

libc/tinymath/__math_oflow.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/tinymath/internal.h"
-
-dontinstrument double __math_oflow(uint32_t sign)
-{
-	return __math_xflow(sign, 0x1p769);
-}

libc/tinymath/__math_oflowf.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/tinymath/internal.h"
-
-dontinstrument float __math_oflowf(uint32_t sign)
-{
-	return __math_xflowf(sign, 0x1p97f);
-}

libc/tinymath/__math_uflow.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/tinymath/internal.h"
-
-dontinstrument double __math_uflow(uint32_t sign)
-{
-	return __math_xflow(sign, 0x1p-767);
-}

libc/tinymath/__math_uflowf.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/tinymath/internal.h"
-
-dontinstrument float __math_uflowf(uint32_t sign)
-{
-	return __math_xflowf(sign, 0x1p-95f);
-}

libc/tinymath/__math_xflow.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/tinymath/internal.h"
-
-dontinstrument double __math_xflow(uint32_t sign, double y)
-{
-	return eval_as_double(fp_barrier(sign ? -y : y) * y);
-}

libc/tinymath/__math_xflowf.c

@@ -1,33 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/tinymath/internal.h"
-
-dontinstrument float __math_xflowf(uint32_t sign, float y)
-{
-	return eval_as_float(fp_barrierf(sign ? -y : y) * y);
-}

libc/tinymath/arm.internal.h

@@ -0,0 +1,509 @@
+#ifndef COSMOPOLITAN_LIBC_TINYMATH_ARM_H_
+#define COSMOPOLITAN_LIBC_TINYMATH_ARM_H_
+#include "libc/math.h"
+COSMOPOLITAN_C_START_
+
+#define USE_GLIBC_ABI 1
+
+/* If defined to 1, return correct results for special cases in non-nearest
+   rounding modes (logf (1.0f) returns 0.0f with FE_DOWNWARD rather than -0.0f).
+   This may be set to 0 if there is no fenv support or if math functions only
+   get called in round to nearest mode.  */
+#ifdef __ROUNDING_MATH__
+#define WANT_ROUNDING 1
+#else
+#define WANT_ROUNDING 0
+#endif
+
+/* If defined to 1, set errno in math functions according to ISO C.  Many math
+   libraries do not set errno, so this is 0 by default.  It may need to be
+   set to 1 if math.h has (math_errhandling & MATH_ERRNO) != 0.  */
+#ifdef __NO_MATH_ERRNO__
+#define WANT_ERRNO 0
+#else
+#define WANT_ERRNO 1
+#endif
+
+/*------------------------------------------------------------------------------*/
+/* optimized-routines/math/math_config.h */
+
+#ifndef WANT_ROUNDING
+/* If defined to 1, return correct results for special cases in non-nearest
+   rounding modes (logf (1.0f) returns 0.0f with FE_DOWNWARD rather than -0.0f).
+   This may be set to 0 if there is no fenv support or if math functions only
+   get called in round to nearest mode.  */
+# define WANT_ROUNDING 1
+#endif
+#ifndef WANT_ERRNO
+/* If defined to 1, set errno in math functions according to ISO C.  Many math
+   libraries do not set errno, so this is 0 by default.  It may need to be
+   set to 1 if math.h has (math_errhandling & MATH_ERRNO) != 0.  */
+# define WANT_ERRNO 0
+#endif
+#ifndef WANT_ERRNO_UFLOW
+/* Set errno to ERANGE if result underflows to 0 (in all rounding modes).  */
+# define WANT_ERRNO_UFLOW (WANT_ROUNDING && WANT_ERRNO)
+#endif
+
+/* Compiler can inline round as a single instruction.  */
+#ifndef HAVE_FAST_ROUND
+# if __aarch64__
+#   define HAVE_FAST_ROUND 1
+# else
+#   define HAVE_FAST_ROUND 0
+# endif
+#endif
+
+/* Compiler can inline lround, but not (long)round(x).  */
+#ifndef HAVE_FAST_LROUND
+# if __aarch64__ && (100*__GNUC__ + __GNUC_MINOR__) >= 408 && __NO_MATH_ERRNO__
+#   define HAVE_FAST_LROUND 1
+# else
+#   define HAVE_FAST_LROUND 0
+# endif
+#endif
+
+/* Compiler can inline fma as a single instruction.  */
+#ifndef HAVE_FAST_FMA
+# if defined FP_FAST_FMA || __aarch64__
+#   define HAVE_FAST_FMA 1
+# else
+#   define HAVE_FAST_FMA 0
+# endif
+#endif
+
+/* Provide *_finite symbols and some of the glibc hidden symbols
+   so libmathlib can be used with binaries compiled against glibc
+   to interpose math functions with both static and dynamic linking.  */
+#ifndef USE_GLIBC_ABI
+# if __GNUC__
+#   define USE_GLIBC_ABI 1
+# else
+#   define USE_GLIBC_ABI 0
+# endif
+#endif
+
+/* Optionally used extensions.  */
+#ifdef __GNUC__
+# define HIDDEN __attribute__ ((__visibility__ ("hidden")))
+# define NOINLINE __attribute__ ((noinline))
+# define UNUSED __attribute__ ((unused))
+# define likely(x) __builtin_expect (!!(x), 1)
+# define unlikely(x) __builtin_expect (x, 0)
+# if __GNUC__ >= 9
+#   define attribute_copy(f) __attribute__ ((copy (f)))
+# else
+#   define attribute_copy(f)
+# endif
+# define strong_alias(f, a) \
+  extern __typeof (f) a __attribute__ ((alias (#f))) attribute_copy (f);
+# define hidden_alias(f, a) \
+  extern __typeof (f) a __attribute__ ((alias (#f), visibility ("hidden"))) \
+  attribute_copy (f);
+#else
+# define HIDDEN
+# define NOINLINE
+# define UNUSED
+# define likely(x) (x)
+# define unlikely(x) (x)
+#endif
+
+/* Return ptr but hide its value from the compiler so accesses through it
+   cannot be optimized based on the contents.  */
+#define ptr_barrier(ptr)                                                      \
+  ({                                                                          \
+    __typeof (ptr) __ptr = (ptr);                                             \
+    __asm("" : "+r"(__ptr));                                                  \
+    __ptr;                                                                    \
+  })
+
+#if HAVE_FAST_ROUND
+/* When set, the roundtoint and converttoint functions are provided with
+   the semantics documented below.  */
+# define TOINT_INTRINSICS 1
+
+/* Round x to nearest int in all rounding modes, ties have to be rounded
+   consistently with converttoint so the results match.  If the result
+   would be outside of [-2^31, 2^31-1] then the semantics is unspecified.  */
+static inline double_t
+roundtoint (double_t x)
+{
+  return round (x);
+}
+
+/* Convert x to nearest int in all rounding modes, ties have to be rounded
+   consistently with roundtoint.  If the result is not representible in an
+   int32_t then the semantics is unspecified.  */
+static inline int32_t
+converttoint (double_t x)
+{
+# if HAVE_FAST_LROUND
+  return lround (x);
+# else
+  return (long) round (x);
+# endif
+}
+#endif
+
+static inline uint32_t
+asuint (float f)
+{
+  union
+  {
+    float f;
+    uint32_t i;
+  } u = {f};
+  return u.i;
+}
+
+static inline float
+asfloat (uint32_t i)
+{
+  union
+  {
+    uint32_t i;
+    float f;
+  } u = {i};
+  return u.f;
+}
+
+static inline uint64_t
+asuint64 (double f)
+{
+  union
+  {
+    double f;
+    uint64_t i;
+  } u = {f};
+  return u.i;
+}
+
+static inline double
+asdouble (uint64_t i)
+{
+  union
+  {
+    uint64_t i;
+    double f;
+  } u = {i};
+  return u.f;
+}
+
+#ifndef IEEE_754_2008_SNAN
+# define IEEE_754_2008_SNAN 1
+#endif
+static inline int
+issignalingf_inline (float x)
+{
+  uint32_t ix = asuint (x);
+  if (!IEEE_754_2008_SNAN)
+    return (ix & 0x7fc00000) == 0x7fc00000;
+  return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000;
+}
+
+static inline int
+issignaling_inline (double x)
+{
+  uint64_t ix = asuint64 (x);
+  if (!IEEE_754_2008_SNAN)
+    return (ix & 0x7ff8000000000000) == 0x7ff8000000000000;
+  return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL;
+}
+
+#if __aarch64__ && __GNUC__
+/* Prevent the optimization of a floating-point expression.  */
+static inline float
+opt_barrier_float (float x)
+{
+  __asm__ __volatile__ ("" : "+w" (x));
+  return x;
+}
+static inline double
+opt_barrier_double (double x)
+{
+  __asm__ __volatile__ ("" : "+w" (x));
+  return x;
+}
+/* Force the evaluation of a floating-point expression for its side-effect.  */
+static inline void
+force_eval_float (float x)
+{
+  __asm__ __volatile__ ("" : "+w" (x));
+}
+static inline void
+force_eval_double (double x)
+{
+  __asm__ __volatile__ ("" : "+w" (x));
+}
+#else
+static inline float
+opt_barrier_float (float x)
+{
+  volatile float y = x;
+  return y;
+}
+static inline double
+opt_barrier_double (double x)
+{
+  volatile double y = x;
+  return y;
+}
+static inline void
+force_eval_float (float x)
+{
+  volatile float y UNUSED = x;
+}
+static inline void
+force_eval_double (double x)
+{
+  volatile double y UNUSED = x;
+}
+#endif
+
+/* Evaluate an expression as the specified type, normally a type
+   cast should be enough, but compilers implement non-standard
+   excess-precision handling, so when FLT_EVAL_METHOD != 0 then
+   these functions may need to be customized.  */
+static inline float
+eval_as_float (float x)
+{
+  return x;
+}
+static inline double
+eval_as_double (double x)
+{
+  return x;
+}
+
+/* Error handling tail calls for special cases, with a sign argument.
+   The sign of the return value is set if the argument is non-zero.  */
+
+/* The result overflows.  */
+HIDDEN float __math_oflowf (uint32_t);
+/* The result underflows to 0 in nearest rounding mode.  */
+HIDDEN float __math_uflowf (uint32_t);
+/* The result underflows to 0 in some directed rounding mode only.  */
+HIDDEN float __math_may_uflowf (uint32_t);
+/* Division by zero.  */
+HIDDEN float __math_divzerof (uint32_t);
+/* The result overflows.  */
+HIDDEN double __math_oflow (uint32_t);
+/* The result underflows to 0 in nearest rounding mode.  */
+HIDDEN double __math_uflow (uint32_t);
+/* The result underflows to 0 in some directed rounding mode only.  */
+HIDDEN double __math_may_uflow (uint32_t);
+/* Division by zero.  */
+HIDDEN double __math_divzero (uint32_t);
+
+/* Error handling using input checking.  */
+
+/* Invalid input unless it is a quiet NaN.  */
+HIDDEN float __math_invalidf (float);
+/* Invalid input unless it is a quiet NaN.  */
+HIDDEN double __math_invalid (double);
+/* Invalid input unless it is a quiet NaN.  */
+HIDDEN long double __math_invalidl (long double);
+
+/* Error handling using output checking, only for errno setting.  */
+
+/* Check if the result overflowed to infinity.  */
+HIDDEN double __math_check_oflow (double);
+/* Check if the result underflowed to 0.  */
+HIDDEN double __math_check_uflow (double);
+
+/* Check if the result overflowed to infinity.  */
+static inline double
+check_oflow (double x)
+{
+  return WANT_ERRNO ? __math_check_oflow (x) : x;
+}
+
+/* Check if the result underflowed to 0.  */
+static inline double
+check_uflow (double x)
+{
+  return WANT_ERRNO ? __math_check_uflow (x) : x;
+}
+
+/* Check if the result overflowed to infinity.  */
+HIDDEN float __math_check_oflowf (float);
+/* Check if the result underflowed to 0.  */
+HIDDEN float __math_check_uflowf (float);
+
+/* Check if the result overflowed to infinity.  */
+static inline float
+check_oflowf (float x)
+{
+  return WANT_ERRNO ? __math_check_oflowf (x) : x;
+}
+
+/* Check if the result underflowed to 0.  */
+static inline float
+check_uflowf (float x)
+{
+  return WANT_ERRNO ? __math_check_uflowf (x) : x;
+}
+
+/* Shared between expf, exp2f and powf.  */
+#define EXP2F_TABLE_BITS 5
+#define EXP2F_POLY_ORDER 3
+extern const struct exp2f_data
+{
+  uint64_t tab[1 << EXP2F_TABLE_BITS];
+  double shift_scaled;
+  double poly[EXP2F_POLY_ORDER];
+  double shift;
+  double invln2_scaled;
+  double poly_scaled[EXP2F_POLY_ORDER];
+} __exp2f_data HIDDEN;
+
+#define LOGF_TABLE_BITS 4
+#define LOGF_POLY_ORDER 4
+extern const struct logf_data
+{
+  struct
+  {
+    double invc, logc;
+  } tab[1 << LOGF_TABLE_BITS];
+  double ln2;
+  double poly[LOGF_POLY_ORDER - 1]; /* First order coefficient is 1.  */
+} __logf_data HIDDEN;
+
+#define LOG2F_TABLE_BITS 4
+#define LOG2F_POLY_ORDER 4
+extern const struct log2f_data
+{
+  struct
+  {
+    double invc, logc;
+  } tab[1 << LOG2F_TABLE_BITS];
+  double poly[LOG2F_POLY_ORDER];
+} __log2f_data HIDDEN;
+
+#define POWF_LOG2_TABLE_BITS 4
+#define POWF_LOG2_POLY_ORDER 5
+#if TOINT_INTRINSICS
+# define POWF_SCALE_BITS EXP2F_TABLE_BITS
+#else
+# define POWF_SCALE_BITS 0
+#endif
+#define POWF_SCALE ((double) (1 << POWF_SCALE_BITS))
+extern const struct powf_log2_data
+{
+  struct
+  {
+    double invc, logc;
+  } tab[1 << POWF_LOG2_TABLE_BITS];
+  double poly[POWF_LOG2_POLY_ORDER];
+} __powf_log2_data HIDDEN;
+
+
+#define EXP_TABLE_BITS 7
+#define EXP_POLY_ORDER 5
+/* Use polynomial that is optimized for a wider input range.  This may be
+   needed for good precision in non-nearest rounding and !TOINT_INTRINSICS.  */
+#define EXP_POLY_WIDE 0
+/* Use close to nearest rounding toint when !TOINT_INTRINSICS.  This may be
+   needed for good precision in non-nearest rouning and !EXP_POLY_WIDE.  */
+#define EXP_USE_TOINT_NARROW 0
+#define EXP2_POLY_ORDER 5
+#define EXP2_POLY_WIDE 0
+/* Wider exp10 polynomial necessary for good precision in non-nearest rounding
+   and !TOINT_INTRINSICS.  */
+#define EXP10_POLY_WIDE 0
+extern const struct exp_data
+{
+  double invln2N;
+  double invlog10_2N;
+  double shift;
+  double negln2hiN;
+  double negln2loN;
+  double neglog10_2hiN;
+  double neglog10_2loN;
+  double poly[4]; /* Last four coefficients.  */
+  double exp2_shift;
+  double exp2_poly[EXP2_POLY_ORDER];
+  double exp10_poly[5];
+  uint64_t tab[2*(1 << EXP_TABLE_BITS)];
+} __exp_data HIDDEN;
+
+#define LOG_TABLE_BITS 7
+#define LOG_POLY_ORDER 6
+#define LOG_POLY1_ORDER 12
+extern const struct log_data
+{
+  double ln2hi;
+  double ln2lo;
+  double poly[LOG_POLY_ORDER - 1]; /* First coefficient is 1.  */
+  double poly1[LOG_POLY1_ORDER - 1];
+  struct {double invc, logc;} tab[1 << LOG_TABLE_BITS];
+#if !HAVE_FAST_FMA
+  struct {double chi, clo;} tab2[1 << LOG_TABLE_BITS];
+#endif
+} __log_data HIDDEN;
+
+#define LOG2_TABLE_BITS 6
+#define LOG2_POLY_ORDER 7
+#define LOG2_POLY1_ORDER 11
+extern const struct log2_data
+{
+  double invln2hi;
+  double invln2lo;
+  double poly[LOG2_POLY_ORDER - 1];
+  double poly1[LOG2_POLY1_ORDER - 1];
+  struct {double invc, logc;} tab[1 << LOG2_TABLE_BITS];
+#if !HAVE_FAST_FMA
+  struct {double chi, clo;} tab2[1 << LOG2_TABLE_BITS];
+#endif
+} __log2_data HIDDEN;
+
+#define POW_LOG_TABLE_BITS 7
+#define POW_LOG_POLY_ORDER 8
+extern const struct pow_log_data
+{
+  double ln2hi;
+  double ln2lo;
+  double poly[POW_LOG_POLY_ORDER - 1]; /* First coefficient is 1.  */
+  /* Note: the pad field is unused, but allows slightly faster indexing.  */
+  struct {double invc, pad, logc, logctail;} tab[1 << POW_LOG_TABLE_BITS];
+} __pow_log_data HIDDEN;
+
+extern const struct erff_data
+{
+  float erff_poly_A[6];
+  float erff_poly_B[7];
+} __erff_data HIDDEN;
+
+#define ERF_POLY_A_ORDER 19
+#define ERF_POLY_A_NCOEFFS 10
+#define ERFC_POLY_C_NCOEFFS 16
+#define ERFC_POLY_D_NCOEFFS 18
+#define ERFC_POLY_E_NCOEFFS 14
+#define ERFC_POLY_F_NCOEFFS 17
+extern const struct erf_data
+{
+  double erf_poly_A[ERF_POLY_A_NCOEFFS];
+  double erf_ratio_N_A[5];
+  double erf_ratio_D_A[5];
+  double erf_ratio_N_B[7];
+  double erf_ratio_D_B[6];
+  double erfc_poly_C[ERFC_POLY_C_NCOEFFS];
+  double erfc_poly_D[ERFC_POLY_D_NCOEFFS];
+  double erfc_poly_E[ERFC_POLY_E_NCOEFFS];
+  double erfc_poly_F[ERFC_POLY_F_NCOEFFS];
+} __erf_data HIDDEN;
+
+#define V_EXP_TABLE_BITS 7
+extern const uint64_t __v_exp_data[1 << V_EXP_TABLE_BITS] HIDDEN;
+
+#define V_LOG_TABLE_BITS 7
+extern const struct v_log_data
+{
+  struct
+  {
+    double invc, logc;
+  } table[1 << V_LOG_TABLE_BITS];
+} __v_log_data HIDDEN;
+
+COSMOPOLITAN_C_END_
+#endif /* COSMOPOLITAN_LIBC_TINYMATH_ARM_H_ */

libc/tinymath/asinhl.c

@@ -5,13 +5,6 @@
 │ FreeBSD lib/msun/src/s_asinhl.c                                              │
 │ Converted to ldbl by David Schultz <[email protected]> and Bruce D. Evans.     │
 │                                                                              │
-│ 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.                                                                │
-│                                                                              │
 │ Copyright (c) 1992-2023 The FreeBSD Project.                                 │
 │                                                                              │
 │ Redistribution and use in source and binary forms, with or without           │
@@ -35,12 +28,17 @@
 │ OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF       │
 │ SUCH DAMAGE.                                                                 │
 │                                                                              │
+│ Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.            │
+│                                                                              │
+│ Developed at SunPro, a Sun Microsystems, Inc. business.                      │
+│ Permission to use, copy, modify, and distribute this                         │
+│ software is freely granted, provided that this notice                        │
+│ is preserved.                                                                │
+│                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/math.h"
 #include "libc/tinymath/freebsd.internal.h"
-
-__static_yoink("fdlibm_notice");
 __static_yoink("freebsd_libm_notice");
+__static_yoink("fdlibm_notice");
 
 #if !(LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024)
 

libc/tinymath/atan2f.c

@@ -1,177 +1,120 @@
-/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
-│ vi: set et ft=c ts=8 sts=2 sw=2 fenc=utf-8                               :vi │
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│ vi: set noet ft=c ts=8 sw=8 fenc=utf-8                                   :vi │
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
-│  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│ FreeBSD lib/msun/src/s_asinhl.c                                              │
+│ Converted to ldbl by David Schultz <[email protected]> and Bruce D. Evans.     │
 │                                                                              │
-│  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:                                                   │
+│ Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.            │
 │                                                                              │
-│  The above copyright notice and this permission notice shall be              │
-│  included in all copies or substantial portions of the Software.             │
+│ 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.                                                                │
 │                                                                              │
-│  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.                      │
+│ Copyright (c) 1992-2023 The FreeBSD Project.                                 │
+│                                                                              │
+│ Redistribution and use in source and binary forms, with or without           │
+│ modification, are permitted provided that the following conditions           │
+│ are met:                                                                     │
+│ 1. Redistributions of source code must retain the above copyright            │
+│    notice, this list of conditions and the following disclaimer.             │
+│ 2. Redistributions in binary form must reproduce the above copyright         │
+│    notice, this list of conditions and the following disclaimer in the       │
+│    documentation and/or other materials provided with the distribution.      │
+│                                                                              │
+│ THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND       │
+│ ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE        │
+│ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE   │
+│ ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE      │
+│ FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL   │
+│ DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS      │
+│ OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)        │
+│ HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT   │
+│ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY    │
+│ OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF       │
+│ SUCH DAMAGE.                                                                 │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/intrin/likely.h"
-#include "libc/math.h"
-#include "libc/tinymath/atanf_common.internal.h"
-#include "libc/tinymath/internal.h"
-__static_yoink("arm_optimized_routines_notice");
+#include "libc/tinymath/freebsd.internal.h"
+__static_yoink("freebsd_libm_notice");
+__static_yoink("fdlibm_notice");
 
-#define Pi (0x1.921fb6p+1f)
-#define PiOver2 (0x1.921fb6p+0f)
-#define PiOver4 (0x1.921fb6p-1f)
-#define SignMask (0x80000000)
+static volatile float
+tiny  = 1.0e-30;
+static const float
+zero  = 0.0,
+pi_o_4  = 7.8539818525e-01, /* 0x3f490fdb */
+pi_o_2  = 1.5707963705e+00, /* 0x3fc90fdb */
+pi      = 3.1415927410e+00; /* 0x40490fdb */
+static volatile float
+pi_lo   = -8.7422776573e-08; /* 0xb3bbbd2e */
 
-/* We calculate atan2f by P(n/d), where n and d are similar to the input
-   arguments, and P is a polynomial. The polynomial may underflow.
-   POLY_UFLOW_BOUND is the lower bound of the difference in exponents of n and d
-   for which P underflows, and is used to special-case such inputs.  */
-#define POLY_UFLOW_BOUND 24
-
-static inline int32_t
-biased_exponent (float f)
-{
-  uint32_t fi = asuint (f);
-  int32_t ex = (int32_t) ((fi & 0x7f800000) >> 23);
-  if (UNLIKELY (ex == 0))
-    {
-      /* Subnormal case - we still need to get the exponent right for subnormal
-	 numbers as division may take us back inside the normal range.  */
-      return ex - __builtin_clz (fi << 9);
-    }
-  return ex;
-}
-
-/* Fast implementation of scalar atan2f. Largest observed error is
-   2.88ulps in [99.0, 101.0] x [99.0, 101.0]:
-   atan2f(0x1.9332d8p+6, 0x1.8cb6c4p+6) got 0x1.964646p-1
-				       want 0x1.964640p-1.  */
+/**
+ * Returns arc tangent of 𝑦/𝑥.
+ */
 float
-atan2f (float y, float x)
+atan2f(float y, float x)
 {
-  uint32_t ix = asuint (x);
-  uint32_t iy = asuint (y);
+	float z;
+	int32_t k,m,hx,hy,ix,iy;
 
-  uint32_t sign_x = ix & SignMask;
-  uint32_t sign_y = iy & SignMask;
+	GET_FLOAT_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	GET_FLOAT_WORD(hy,y);
+	iy = hy&0x7fffffff;
+	if((ix>0x7f800000)||
+	   (iy>0x7f800000))	/* x or y is NaN */
+	    return nan_mix(x, y);
+	if(hx==0x3f800000) return atanf(y);   /* x=1.0 */
+	m = ((hy>>31)&1)|((hx>>30)&2);	/* 2*sign(x)+sign(y) */
 
-  uint32_t iax = ix & ~SignMask;
-  uint32_t iay = iy & ~SignMask;
-
-  /* x or y is NaN.  */
-  if ((iax > 0x7f800000) || (iay > 0x7f800000))
-    return x + y;
-
-  /* m = 2 * sign(x) + sign(y).  */
-  uint32_t m = ((iy >> 31) & 1) | ((ix >> 30) & 2);
-
-  /* The following follows glibc ieee754 implementation, except
-     that we do not use +-tiny shifts (non-nearest rounding mode).  */
-
-  int32_t exp_diff = biased_exponent (x) - biased_exponent (y);
-
-  /* Special case for (x, y) either on or very close to the x axis. Either y =
-     0, or y is tiny and x is huge (difference in exponents >=
-     POLY_UFLOW_BOUND). In the second case, we only want to use this special
-     case when x is negative (i.e. quadrants 2 or 3).  */
-  if (UNLIKELY (iay == 0 || (exp_diff >= POLY_UFLOW_BOUND && m >= 2)))
-    {
-      switch (m)
-	{
-	case 0:
-	case 1:
-	  return y; /* atan(+-0,+anything)=+-0.  */
-	case 2:
-	  return Pi; /* atan(+0,-anything) = pi.  */
-	case 3:
-	  return -Pi; /* atan(-0,-anything) =-pi.  */
-	}
-    }
-  /* Special case for (x, y) either on or very close to the y axis. Either x =
-     0, or x is tiny and y is huge (difference in exponents >=
-     POLY_UFLOW_BOUND).  */
-  if (UNLIKELY (iax == 0 || exp_diff <= -POLY_UFLOW_BOUND))
-    return sign_y ? -PiOver2 : PiOver2;
-
-  /* x is INF.  */
-  if (iax == 0x7f800000)
-    {
-      if (iay == 0x7f800000)
-	{
-	  switch (m)
-	    {
-	    case 0:
-	      return PiOver4; /* atan(+INF,+INF).  */
-	    case 1:
-	      return -PiOver4; /* atan(-INF,+INF).  */
-	    case 2:
-	      return 3.0f * PiOver4; /* atan(+INF,-INF).  */
-	    case 3:
-	      return -3.0f * PiOver4; /* atan(-INF,-INF).  */
+    /* when y = 0 */
+	if(iy==0) {
+	    switch(m) {
+		case 0:
+		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
+		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
+		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
 	    }
 	}
-      else
-	{
-	  switch (m)
-	    {
-	    case 0:
-	      return 0.0f; /* atan(+...,+INF).  */
-	    case 1:
-	      return -0.0f; /* atan(-...,+INF).  */
-	    case 2:
-	      return Pi; /* atan(+...,-INF).  */
-	    case 3:
-	      return -Pi; /* atan(-...,-INF).  */
+    /* when x = 0 */
+	if(ix==0) return (hy<0)?  -pi_o_2-tiny: pi_o_2+tiny;
+
+    /* when x is INF */
+	if(ix==0x7f800000) {
+	    if(iy==0x7f800000) {
+		switch(m) {
+		    case 0: return  pi_o_4+tiny;/* atan(+INF,+INF) */
+		    case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
+		    case 2: return  (float)3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
+		    case 3: return (float)-3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
+		}
+	    } else {
+		switch(m) {
+		    case 0: return  zero  ;	/* atan(+...,+INF) */
+		    case 1: return -zero  ;	/* atan(-...,+INF) */
+		    case 2: return  pi+tiny  ;	/* atan(+...,-INF) */
+		    case 3: return -pi-tiny  ;	/* atan(-...,-INF) */
+		}
 	    }
 	}
-    }
-  /* y is INF.  */
-  if (iay == 0x7f800000)
-    return sign_y ? -PiOver2 : PiOver2;
+    /* when y is INF */
+	if(iy==0x7f800000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
 
-  uint32_t sign_xy = sign_x ^ sign_y;
-
-  float ax = asfloat (iax);
-  float ay = asfloat (iay);
-
-  bool pred_aygtax = (ay > ax);
-
-  /* Set up z for call to atanf.  */
-  float n = pred_aygtax ? -ax : ay;
-  float d = pred_aygtax ? ay : ax;
-  float z = n / d;
-
-  float ret;
-  if (UNLIKELY (m < 2 && exp_diff >= POLY_UFLOW_BOUND))
-    {
-      /* If (x, y) is very close to x axis and x is positive, the polynomial
-	 will underflow and evaluate to z.  */
-      ret = z;
-    }
-  else
-    {
-      /* Work out the correct shift.  */
-      float shift = sign_x ? -2.0f : 0.0f;
-      shift = pred_aygtax ? shift + 1.0f : shift;
-      shift *= PiOver2;
-
-      ret = eval_poly (z, z, shift);
-    }
-
-  /* Account for the sign of x and y.  */
-  return asfloat (asuint (ret) ^ sign_xy);
+    /* compute y/x */
+	k = (iy-ix)>>23;
+	if(k > 26) {			/* |y/x| >  2**26 */
+	    z=pi_o_2+(float)0.5*pi_lo;
+	    m&=1;
+	}
+	else if(k<-26&&hx<0) z=0.0; 	/* 0 > |y|/x > -2**-26 */
+	else z=atanf(fabsf(y/x));	/* safe to do y/x */
+	switch (m) {
+	    case 0: return       z  ;	/* atan(+,+) */
+	    case 1: return      -z  ;	/* atan(-,+) */
+	    case 2: return  pi-(z-pi_lo);/* atan(+,-) */
+	    default: /* case 3 */
+	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
+	}
 }

libc/tinymath/atan_common.internal.h

@@ -1,54 +0,0 @@
-#ifndef COSMOPOLITAN_LIBC_TINYMATH_ATAN_COMMON_H_
-#define COSMOPOLITAN_LIBC_TINYMATH_ATAN_COMMON_H_
-#include "libc/tinymath/atan_data.internal.h"
-#include "libc/tinymath/estrin_wrap.internal.h"
-#include "libc/tinymath/horner.internal.h"
-COSMOPOLITAN_C_START_
-
-/*
- * Double-precision polynomial evaluation function for scalar and vector atan(x)
- * and atan2(y,x).
- *
- * Copyright (c) 2021-2023, Arm Limited.
- * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
- */
-
-#if WANT_VMATH
-
-#define DBL_T float64x2_t
-#define P(i) v_f64 (__atan_poly_data.poly[i])
-
-#else
-
-#define DBL_T double
-#define P(i) __atan_poly_data.poly[i]
-
-#endif
-
-/* Polynomial used in fast atan(x) and atan2(y,x) implementations
-   The order 19 polynomial P approximates (atan(sqrt(x))-sqrt(x))/x^(3/2).  */
-static inline DBL_T
-eval_poly (DBL_T z, DBL_T az, DBL_T shift)
-{
-  /* Use split Estrin scheme for P(z^2) with deg(P)=19. Use split instead of
-     full scheme to avoid underflow in x^16.  */
-  DBL_T z2 = z * z;
-  DBL_T x2 = z2 * z2;
-  DBL_T x4 = x2 * x2;
-  DBL_T x8 = x4 * x4;
-  DBL_T y
-    = FMA (ESTRIN_11_ (z2, x2, x4, x8, P, 8), x8, ESTRIN_7 (z2, x2, x4, P));
-
-  /* Finalize. y = shift + z + z^3 * P(z^2).  */
-  y = FMA (y, z2 * az, az);
-  y = y + shift;
-
-  return y;
-}
-
-#undef DBL_T
-#undef FMA
-#undef P
-
-COSMOPOLITAN_C_END_
-#endif /* COSMOPOLITAN_LIBC_TINYMATH_ATAN_COMMON_H_ */

libc/tinymath/atan_data.internal.h

@@ -1,11 +0,0 @@
-#ifndef COSMOPOLITAN_LIBC_TINYMATH_ATAN_DATA_H_
-#define COSMOPOLITAN_LIBC_TINYMATH_ATAN_DATA_H_
-COSMOPOLITAN_C_START_
-
-#define ATAN_POLY_NCOEFFS 20
-extern const struct atan_poly_data {
-  double poly[ATAN_POLY_NCOEFFS];
-} __atan_poly_data;
-
-COSMOPOLITAN_C_END_
-#endif /* COSMOPOLITAN_LIBC_TINYMATH_ATAN_DATA_H_ */

libc/tinymath/atanf_common.internal.h

@@ -1,43 +0,0 @@
-#ifndef COSMOPOLITAN_LIBC_TINYMATH_ATANF_COMMON_H_
-#define COSMOPOLITAN_LIBC_TINYMATH_ATANF_COMMON_H_
-#include "libc/tinymath/atanf_data.internal.h"
-#include "libc/tinymath/estrin_wrap.internal.h"
-#include "libc/tinymath/hornerf.internal.h"
-COSMOPOLITAN_C_START_
-
-#if WANT_VMATH
-
-#define FLT_T float32x4_t
-#define P(i) v_f32 (__atanf_poly_data.poly[i])
-
-#else
-
-#define FLT_T float
-#define P(i) __atanf_poly_data.poly[i]
-
-#endif
-
-/* Polynomial used in fast atanf(x) and atan2f(y,x) implementations
-   The order 7 polynomial P approximates (atan(sqrt(x))-sqrt(x))/x^(3/2).  */
-static inline FLT_T
-eval_poly (FLT_T z, FLT_T az, FLT_T shift)
-{
-  /* Use 2-level Estrin scheme for P(z^2) with deg(P)=7. However,
-     a standard implementation using z8 creates spurious underflow
-     in the very last fma (when z^8 is small enough).
-     Therefore, we split the last fma into a mul and and an fma.
-     Horner and single-level Estrin have higher errors that exceed
-     threshold.  */
-  FLT_T z2 = z * z;
-  FLT_T z4 = z2 * z2;
-
-  /* Then assemble polynomial.  */
-  FLT_T y = FMA (z4, z4 * ESTRIN_3_ (z2, z4, P, 4), ESTRIN_3 (z2, z4, P));
-
-  /* Finalize:
-     y = shift + z * P(z^2).  */
-  return FMA (y, z2 * az, az) + shift;
-}
-
-COSMOPOLITAN_C_END_
-#endif /* COSMOPOLITAN_LIBC_TINYMATH_ATANF_COMMON_H_ */

libc/tinymath/atanf_data.internal.h

@@ -1,11 +0,0 @@
-#ifndef COSMOPOLITAN_LIBC_TINYMATH_ATANF_DATA_H_
-#define COSMOPOLITAN_LIBC_TINYMATH_ATANF_DATA_H_
-COSMOPOLITAN_C_START_
-
-#define ATANF_POLY_NCOEFFS 8
-extern const struct atanf_poly_data {
-  float poly[ATANF_POLY_NCOEFFS];
-} __atanf_poly_data;
-
-COSMOPOLITAN_C_END_
-#endif /* COSMOPOLITAN_LIBC_TINYMATH_ATANF_DATA_H_ */

libc/tinymath/cacos.c

@@ -1,41 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-
-// FIXME: Hull et al. "Implementing the complex arcsine and arccosine functions using exception handling" 1997
-
-/* acos(z) = pi/2 - asin(z) */
-
-double complex cacos(double complex z)
-{
-	z = casin(z);
-	return CMPLX(M_PI_2 - creal(z), -cimag(z));
-}

libc/tinymath/cacosf.c

@@ -1,39 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-
-// FIXME
-
-float complex cacosf(float complex z)
-{
-	z = casinf(z);
-	return CMPLXF((float)M_PI_2 - crealf(z), -cimagf(z));
-}

libc/tinymath/cacosh.c

@@ -1,44 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-
-/* acosh(z) = i acos(z) */
-double complex cacosh(double complex z)
-{
-	int zineg = signbit(cimag(z));
-	z = cacos(z);
-	if (zineg) return CMPLX(cimag(z), -creal(z));
-	else       return CMPLX(-cimag(z), creal(z));
-}
-
-#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
-__weak_reference(cacosh, cacoshl);
-#endif

libc/tinymath/cacoshf.c

@@ -1,39 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-
-float complex cacoshf(float complex z)
-{
-	int zineg = signbit(cimagf(z));
-	z = cacosf(z);
-	if (zineg) return CMPLXF(cimagf(z), -crealf(z));
-	else       return CMPLXF(-cimagf(z), crealf(z));
-}

libc/tinymath/casinf.c

@@ -1,45 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-
-// FIXME
-
-float complex casinf(float complex z)
-{
-	float complex w;
-	float x, y;
-
-	x = crealf(z);
-	y = cimagf(z);
-	w = CMPLXF(1.0 - (x - y)*(x + y), -2.0*x*y);
-	float complex r = clogf(CMPLXF(-y, x) + csqrtf(w));
-	return CMPLXF(cimagf(r), -crealf(r));
-}

libc/tinymath/casinh.c

@@ -1,39 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-
-/* asinh(z) = -i asin(i z) */
-
-double complex casinh(double complex z)
-{
-	z = casin(CMPLX(-cimag(z), creal(z)));
-	return CMPLX(cimag(z), -creal(z));
-}

libc/tinymath/casinhf.c

@@ -1,37 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-
-float complex casinhf(float complex z)
-{
-	z = casinf(CMPLXF(-cimagf(z), crealf(z)));
-	return CMPLXF(cimagf(z), -crealf(z));
-}

libc/tinymath/catan.c

@@ -1,142 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-__static_yoink("openbsd_libm_notice");
-
-/* origin: OpenBSD /usr/src/lib/libm/src/s_catan.c */
-/*
- * Copyright (c) 2008 Stephen L. Moshier <[email protected]>
- *
- * Permission to use, copy, modify, and distribute this software for any
- * purpose with or without fee is hereby granted, provided that the above
- * copyright notice and this permission notice appear in all copies.
- *
- * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
- * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
- * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
- * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
- * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
- * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
- * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
- */
-/*
- *      Complex circular arc tangent
- *
- *
- * SYNOPSIS:
- *
- * double complex catan();
- * double complex z, w;
- *
- * w = catan (z);
- *
- *
- * DESCRIPTION:
- *
- * If
- *     z = x + iy,
- *
- * then
- *          1       (    2x     )
- * Re w  =  - arctan(-----------)  +  k PI
- *          2       (     2    2)
- *                  (1 - x  - y )
- *
- *               ( 2         2)
- *          1    (x  +  (y+1) )
- * Im w  =  - log(------------)
- *          4    ( 2         2)
- *               (x  +  (y-1) )
- *
- * Where k is an arbitrary integer.
- *
- * catan(z) = -i catanh(iz).
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    DEC       -10,+10      5900       1.3e-16     7.8e-18
- *    IEEE      -10,+10     30000       2.3e-15     8.5e-17
- * The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,
- * had peak relative error 1.5e-16, rms relative error
- * 2.9e-17.  See also clog().
- */
-
-#define MAXNUM 1.0e308
-
-static const double DP1 = 3.14159265160560607910E0;
-static const double DP2 = 1.98418714791870343106E-9;
-static const double DP3 = 1.14423774522196636802E-17;
-
-static double _redupi(double x)
-{
-	double t;
-	long i;
-
-	t = x/M_PI;
-	if (t >= 0.0)
-		t += 0.5;
-	else
-		t -= 0.5;
-
-	i = t;  /* the multiple */
-	t = i;
-	t = ((x - t * DP1) - t * DP2) - t * DP3;
-	return t;
-}
-
-double complex catan(double complex z)
-{
-	double complex w;
-	double a, t, x, x2, y;
-
-	x = creal(z);
-	y = cimag(z);
-
-	x2 = x * x;
-	a = 1.0 - x2 - (y * y);
-
-	t = 0.5 * atan2(2.0 * x, a);
-	w = _redupi(t);
-
-	t = y - 1.0;
-	a = x2 + (t * t);
-
-	t = y + 1.0;
-	a = (x2 + t * t)/a;
-	w = CMPLX(w, 0.25 * log(a));
-	return w;
-}
-
-#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
-__weak_reference(catan, catanl);
-#endif

libc/tinymath/catanf.c

@@ -1,135 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-__static_yoink("openbsd_libm_notice");
-
-/* origin: OpenBSD /usr/src/lib/libm/src/s_catanf.c */
-/*
- * Copyright (c) 2008 Stephen L. Moshier <[email protected]>
- *
- * Permission to use, copy, modify, and distribute this software for any
- * purpose with or without fee is hereby granted, provided that the above
- * copyright notice and this permission notice appear in all copies.
- *
- * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
- * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
- * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
- * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
- * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
- * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
- * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
- */
-/*
- *      Complex circular arc tangent
- *
- *
- * SYNOPSIS:
- *
- * float complex catanf();
- * float complex z, w;
- *
- * w = catanf( z );
- *
- *
- * DESCRIPTION:
- *
- * If
- *     z = x + iy,
- *
- * then
- *          1       (    2x     )
- * Re w  =  - arctan(-----------)  +  k PI
- *          2       (     2    2)
- *                  (1 - x  - y )
- *
- *               ( 2         2)
- *          1    (x  +  (y+1) )
- * Im w  =  - log(------------)
- *          4    ( 2         2)
- *               (x  +  (y-1) )
- *
- * Where k is an arbitrary integer.
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      -10,+10     30000        2.3e-6      5.2e-8
- */
-
-
-#define MAXNUMF 1.0e38F
-
-static const double DP1 = 3.140625;
-static const double DP2 = 9.67502593994140625E-4;
-static const double DP3 = 1.509957990978376432E-7;
-
-static float _redupif(float xx)
-{
-	float x, t;
-	long i;
-
-	x = xx;
-	t = x/(float)M_PI;
-	if (t >= 0.0f)
-		t += 0.5f;
-	else
-		t -= 0.5f;
-
-	i = t;  /* the multiple */
-	t = i;
-	t = ((x - t * DP1) - t * DP2) - t * DP3;
-	return t;
-}
-
-float complex catanf(float complex z)
-{
-	float complex w;
-	float a, t, x, x2, y;
-
-	x = crealf(z);
-	y = cimagf(z);
-
-	x2 = x * x;
-	a = 1.0f - x2 - (y * y);
-
-	t = 0.5f * atan2f(2.0f * x, a);
-	w = _redupif(t);
-
-	t = y - 1.0f;
-	a = x2 + (t * t);
-
-	t = y + 1.0f;
-	a = (x2 + (t * t))/a;
-	w = CMPLXF(w, 0.25f * logf(a));
-	return w;
-}

libc/tinymath/catanh.c

@@ -1,39 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-
-/* atanh = -i atan(i z) */
-
-double complex catanh(double complex z)
-{
-	z = catan(CMPLX(-cimag(z), creal(z)));
-	return CMPLX(cimag(z), -creal(z));
-}

libc/tinymath/catanhf.c

@@ -1,37 +0,0 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=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/complex.h"
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-__static_yoink("musl_libc_notice");
-
-float complex catanhf(float complex z)
-{
-	z = catanf(CMPLXF(-cimagf(z), crealf(z)));
-	return CMPLXF(cimagf(z), -crealf(z));
-}

libc/tinymath/catrig.c

@@ -0,0 +1,651 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│ vi: set noet ft=c ts=8 sw=8 fenc=utf-8                                   :vi │
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│ FreeBSD lib/msun/src/catrig.c                                                │
+│ Copyright (c) 2012 Stephen Montgomery-Smith <[email protected]>            │
+│                                                                              │
+│ Redistribution and use in source and binary forms, with or without           │
+│ modification, are permitted provided that the following conditions           │
+│ are met:                                                                     │
+│ 1. Redistributions of source code must retain the above copyright            │
+│    notice, this list of conditions and the following disclaimer.             │
+│ 2. Redistributions in binary form must reproduce the above copyright         │
+│    notice, this list of conditions and the following disclaimer in the       │
+│    documentation and/or other materials provided with the distribution.      │
+│                                                                              │
+│ THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND       │
+│ ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE        │
+│ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE   │
+│ ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE      │
+│ FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL   │
+│ DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS      │
+│ OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)        │
+│ HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT   │
+│ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY    │
+│ OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF       │
+│ SUCH DAMAGE.                                                                 │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/tinymath/freebsd.internal.h"
+__static_yoink("freebsd_complex_notice");
+
+#undef isinf
+#define isinf(x)	(fabs(x) == INFINITY)
+#undef isnan
+#define isnan(x)	((x) != (x))
+#define	raise_inexact()	do { volatile float _j = 1 + tiny; (void)_j; } while(0)
+#undef signbit
+#define signbit(x)	(__builtin_signbit(x))
+
+/* We need that DBL_EPSILON^2/128 is larger than FOUR_SQRT_MIN. */
+static const double
+A_crossover =		10, /* Hull et al suggest 1.5, but 10 works better */
+B_crossover =		0.6417,			/* suggested by Hull et al */
+FOUR_SQRT_MIN =		0x1p-509,		/* >= 4 * sqrt(DBL_MIN) */
+QUARTER_SQRT_MAX =	0x1p509,		/* <= sqrt(DBL_MAX) / 4 */
+m_e =			2.7182818284590452e0,	/*  0x15bf0a8b145769.0p-51 */
+m_ln2 =			6.9314718055994531e-1,	/*  0x162e42fefa39ef.0p-53 */
+pio2_hi =		1.5707963267948966e0,	/*  0x1921fb54442d18.0p-52 */
+RECIP_EPSILON =		1 / DBL_EPSILON,
+SQRT_3_EPSILON =	2.5809568279517849e-8,	/*  0x1bb67ae8584caa.0p-78 */
+SQRT_6_EPSILON =	3.6500241499888571e-8,	/*  0x13988e1409212e.0p-77 */
+SQRT_MIN =		0x1p-511;		/* >= sqrt(DBL_MIN) */
+
+static const volatile double
+pio2_lo =		6.1232339957367659e-17;	/*  0x11a62633145c07.0p-106 */
+static const volatile float
+tiny =			0x1p-100;
+
+static double complex clog_for_large_values(double complex z);
+
+/*
+ * Testing indicates that all these functions are accurate up to 4 ULP.
+ * The functions casin(h) and cacos(h) are about 2.5 times slower than asinh.
+ * The functions catan(h) are a little under 2 times slower than atanh.
+ *
+ * The code for casinh, casin, cacos, and cacosh comes first.  The code is
+ * rather complicated, and the four functions are highly interdependent.
+ *
+ * The code for catanh and catan comes at the end.  It is much simpler than
+ * the other functions, and the code for these can be disconnected from the
+ * rest of the code.
+ */
+
+/*
+ *			================================
+ *			| casinh, casin, cacos, cacosh |
+ *			================================
+ */
+
+/*
+ * The algorithm is very close to that in "Implementing the complex arcsine
+ * and arccosine functions using exception handling" by T. E. Hull, Thomas F.
+ * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
+ * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
+ * http://dl.acm.org/citation.cfm?id=275324.
+ *
+ * Throughout we use the convention z = x + I*y.
+ *
+ * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B)
+ * where
+ * A = (|z+I| + |z-I|) / 2
+ * B = (|z+I| - |z-I|) / 2 = y/A
+ *
+ * These formulas become numerically unstable:
+ *   (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that
+ *       is, Re(casinh(z)) is close to 0);
+ *   (b) for Im(casinh(z)) when z is close to either of the intervals
+ *       [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is
+ *       close to PI/2).
+ *
+ * These numerical problems are overcome by defining
+ * f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2
+ * Then if A < A_crossover, we use
+ *   log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1)))
+ *   A-1 = f(x, 1+y) + f(x, 1-y)
+ * and if B > B_crossover, we use
+ *   asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y)))
+ *   A-y = f(x, y+1) + f(x, y-1)
+ * where without loss of generality we have assumed that x and y are
+ * non-negative.
+ *
+ * Much of the difficulty comes because the intermediate computations may
+ * produce overflows or underflows.  This is dealt with in the paper by Hull
+ * et al by using exception handling.  We do this by detecting when
+ * computations risk underflow or overflow.  The hardest part is handling the
+ * underflows when computing f(a, b).
+ *
+ * Note that the function f(a, b) does not appear explicitly in the paper by
+ * Hull et al, but the idea may be found on pages 308 and 309.  Introducing the
+ * function f(a, b) allows us to concentrate many of the clever tricks in this
+ * paper into one function.
+ */
+
+/*
+ * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2.
+ * Pass hypot(a, b) as the third argument.
+ */
+static inline double
+f(double a, double b, double hypot_a_b)
+{
+	if (b < 0)
+		return ((hypot_a_b - b) / 2);
+	if (b == 0)
+		return (a / 2);
+	return (a * a / (hypot_a_b + b) / 2);
+}
+
+/*
+ * All the hard work is contained in this function.
+ * x and y are assumed positive or zero, and less than RECIP_EPSILON.
+ * Upon return:
+ * rx = Re(casinh(z)) = -Im(cacos(y + I*x)).
+ * B_is_usable is set to 1 if the value of B is usable.
+ * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y.
+ * If returning sqrt_A2my2 has potential to result in an underflow, it is
+ * rescaled, and new_y is similarly rescaled.
+ */
+static inline void
+do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
+    double *sqrt_A2my2, double *new_y)
+{
+	double R, S, A; /* A, B, R, and S are as in Hull et al. */
+	double Am1, Amy; /* A-1, A-y. */
+
+	R = hypot(x, y + 1);		/* |z+I| */
+	S = hypot(x, y - 1);		/* |z-I| */
+
+	/* A = (|z+I| + |z-I|) / 2 */
+	A = (R + S) / 2;
+	/*
+	 * Mathematically A >= 1.  There is a small chance that this will not
+	 * be so because of rounding errors.  So we will make certain it is
+	 * so.
+	 */
+	if (A < 1)
+		A = 1;
+
+	if (A < A_crossover) {
+		/*
+		 * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).
+		 * rx = log1p(Am1 + sqrt(Am1*(A+1)))
+		 */
+		if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) {
+			/*
+			 * fp is of order x^2, and fm = x/2.
+			 * A = 1 (inexactly).
+			 */
+			*rx = sqrt(x);
+		} else if (x >= DBL_EPSILON * fabs(y - 1)) {
+			/*
+			 * Underflow will not occur because
+			 * x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN
+			 */
+			Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
+			*rx = log1p(Am1 + sqrt(Am1 * (A + 1)));
+		} else if (y < 1) {
+			/*
+			 * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and
+			 * A = 1 (inexactly).
+			 */
+			*rx = x / sqrt((1 - y) * (1 + y));
+		} else {		/* if (y > 1) */
+			/*
+			 * A-1 = y-1 (inexactly).
+			 */
+			*rx = log1p((y - 1) + sqrt((y - 1) * (y + 1)));
+		}
+	} else {
+		*rx = log(A + sqrt(A * A - 1));
+	}
+
+	*new_y = y;
+
+	if (y < FOUR_SQRT_MIN) {
+		/*
+		 * Avoid a possible underflow caused by y/A.  For casinh this
+		 * would be legitimate, but will be picked up by invoking atan2
+		 * later on.  For cacos this would not be legitimate.
+		 */
+		*B_is_usable = 0;
+		*sqrt_A2my2 = A * (2 / DBL_EPSILON);
+		*new_y = y * (2 / DBL_EPSILON);
+		return;
+	}
+
+	/* B = (|z+I| - |z-I|) / 2 = y/A */
+	*B = y / A;
+	*B_is_usable = 1;
+
+	if (*B > B_crossover) {
+		*B_is_usable = 0;
+		/*
+		 * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1).
+		 * sqrt_A2my2 = sqrt(Amy*(A+y))
+		 */
+		if (y == 1 && x < DBL_EPSILON / 128) {
+			/*
+			 * fp is of order x^2, and fm = x/2.
+			 * A = 1 (inexactly).
+			 */
+			*sqrt_A2my2 = sqrt(x) * sqrt((A + y) / 2);
+		} else if (x >= DBL_EPSILON * fabs(y - 1)) {
+			/*
+			 * Underflow will not occur because
+			 * x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN
+			 * and
+			 * x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN
+			 */
+			Amy = f(x, y + 1, R) + f(x, y - 1, S);
+			*sqrt_A2my2 = sqrt(Amy * (A + y));
+		} else if (y > 1) {
+			/*
+			 * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and
+			 * A = y (inexactly).
+			 *
+			 * y < RECIP_EPSILON.  So the following
+			 * scaling should avoid any underflow problems.
+			 */
+			*sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y /
+			    sqrt((y + 1) * (y - 1));
+			*new_y = y * (4 / DBL_EPSILON / DBL_EPSILON);
+		} else {		/* if (y < 1) */
+			/*
+			 * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and
+			 * A = 1 (inexactly).
+			 */
+			*sqrt_A2my2 = sqrt((1 - y) * (1 + y));
+		}
+	}
+}
+
+/*
+ * casinh(z) = z + O(z^3)   as z -> 0
+ *
+ * casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2)   as z -> infinity
+ * The above formula works for the imaginary part as well, because
+ * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3)
+ *    as z -> infinity, uniformly in y
+ */
+double complex
+casinh(double complex z)
+{
+	double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
+	int B_is_usable;
+	double complex w;
+
+	x = creal(z);
+	y = cimag(z);
+	ax = fabs(x);
+	ay = fabs(y);
+
+	if (isnan(x) || isnan(y)) {
+		/* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */
+		if (isinf(x))
+			return (CMPLX(x, y + y));
+		/* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
+		if (isinf(y))
+			return (CMPLX(y, x + x));
+		/* casinh(NaN + I*0) = NaN + I*0 */
+		if (y == 0)
+			return (CMPLX(x + x, y));
+		/*
+		 * All other cases involving NaN return NaN + I*NaN.
+		 * C99 leaves it optional whether to raise invalid if one of
+		 * the arguments is not NaN, so we opt not to raise it.
+		 */
+		return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+		/* clog...() will raise inexact unless x or y is infinite. */
+		if (signbit(x) == 0)
+			w = clog_for_large_values(z) + m_ln2;
+		else
+			w = clog_for_large_values(-z) + m_ln2;
+		return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y)));
+	}
+
+	/* Avoid spuriously raising inexact for z = 0. */
+	if (x == 0 && y == 0)
+		return (z);
+
+	/* All remaining cases are inexact. */
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+		return (z);
+
+	do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
+	if (B_is_usable)
+		ry = asin(B);
+	else
+		ry = atan2(new_y, sqrt_A2my2);
+	return (CMPLX(copysign(rx, x), copysign(ry, y)));
+}
+
+/*
+ * casin(z) = reverse(casinh(reverse(z)))
+ * where reverse(x + I*y) = y + I*x = I*conj(z).
+ */
+double complex
+casin(double complex z)
+{
+	double complex w = casinh(CMPLX(cimag(z), creal(z)));
+
+	return (CMPLX(cimag(w), creal(w)));
+}
+
+/*
+ * cacos(z) = PI/2 - casin(z)
+ * but do the computation carefully so cacos(z) is accurate when z is
+ * close to 1.
+ *
+ * cacos(z) = PI/2 - z + O(z^3)   as z -> 0
+ *
+ * cacos(z) = -sign(y)*I*clog(z) + O(1/z^2)   as z -> infinity
+ * The above formula works for the real part as well, because
+ * Re(cacos(z)) = atan2(fabs(y), x) + O(y/z^3)
+ *    as z -> infinity, uniformly in y
+ */
+double complex
+cacos(double complex z)
+{
+	double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
+	int sx, sy;
+	int B_is_usable;
+	double complex w;
+
+	x = creal(z);
+	y = cimag(z);
+	sx = signbit(x);
+	sy = signbit(y);
+	ax = fabs(x);
+	ay = fabs(y);
+
+	if (isnan(x) || isnan(y)) {
+		/* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
+		if (isinf(x))
+			return (CMPLX(y + y, -INFINITY));
+		/* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */
+		if (isinf(y))
+			return (CMPLX(x + x, -y));
+		/* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */
+		if (x == 0)
+			return (CMPLX(pio2_hi + pio2_lo, y + y));
+		/*
+		 * All other cases involving NaN return NaN + I*NaN.
+		 * C99 leaves it optional whether to raise invalid if one of
+		 * the arguments is not NaN, so we opt not to raise it.
+		 */
+		return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+		/* clog...() will raise inexact unless x or y is infinite. */
+		w = clog_for_large_values(z);
+		rx = fabs(cimag(w));
+		ry = creal(w) + m_ln2;
+		if (sy == 0)
+			ry = -ry;
+		return (CMPLX(rx, ry));
+	}
+
+	/* Avoid spuriously raising inexact for z = 1. */
+	if (x == 1 && y == 0)
+		return (CMPLX(0, -y));
+
+	/* All remaining cases are inexact. */
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+		return (CMPLX(pio2_hi - (x - pio2_lo), -y));
+
+	do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
+	if (B_is_usable) {
+		if (sx == 0)
+			rx = acos(B);
+		else
+			rx = acos(-B);
+	} else {
+		if (sx == 0)
+			rx = atan2(sqrt_A2mx2, new_x);
+		else
+			rx = atan2(sqrt_A2mx2, -new_x);
+	}
+	if (sy == 0)
+		ry = -ry;
+	return (CMPLX(rx, ry));
+}
+
+/*
+ * cacosh(z) = I*cacos(z) or -I*cacos(z)
+ * where the sign is chosen so Re(cacosh(z)) >= 0.
+ */
+double complex
+cacosh(double complex z)
+{
+	double complex w;
+	double rx, ry;
+
+	w = cacos(z);
+	rx = creal(w);
+	ry = cimag(w);
+	/* cacosh(NaN + I*NaN) = NaN + I*NaN */
+	if (isnan(rx) && isnan(ry))
+		return (CMPLX(ry, rx));
+	/* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */
+	/* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */
+	if (isnan(rx))
+		return (CMPLX(fabs(ry), rx));
+	/* cacosh(0 + I*NaN) = NaN + I*NaN */
+	if (isnan(ry))
+		return (CMPLX(ry, ry));
+	return (CMPLX(fabs(ry), copysign(rx, cimag(z))));
+}
+
+/*
+ * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON.
+ */
+static double complex
+clog_for_large_values(double complex z)
+{
+	double x, y;
+	double ax, ay, t;
+
+	x = creal(z);
+	y = cimag(z);
+	ax = fabs(x);
+	ay = fabs(y);
+	if (ax < ay) {
+		t = ax;
+		ax = ay;
+		ay = t;
+	}
+
+	/*
+	 * Avoid overflow in hypot() when x and y are both very large.
+	 * Divide x and y by E, and then add 1 to the logarithm.  This
+	 * depends on E being larger than sqrt(2), since the return value of
+	 * hypot cannot overflow if neither argument is greater in magnitude
+	 * than 1/sqrt(2) of the maximum value of the return type.  Likewise
+	 * this determines the necessary threshold for using this method
+	 * (however, actually use 1/2 instead as it is simpler).
+	 *
+	 * Dividing by E causes an insignificant loss of accuracy; however
+	 * this method is still poor since it is uneccessarily slow.
+	 */
+	if (ax > DBL_MAX / 2)
+		return (CMPLX(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x)));
+
+	/*
+	 * Avoid overflow when x or y is large.  Avoid underflow when x or
+	 * y is small.
+	 */
+	if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
+		return (CMPLX(log(hypot(x, y)), atan2(y, x)));
+
+	return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x)));
+}
+
+/*
+ *				=================
+ *				| catanh, catan |
+ *				=================
+ */
+
+/*
+ * sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow).
+ * Assumes x*x and y*y will not overflow.
+ * Assumes x and y are finite.
+ * Assumes y is non-negative.
+ * Assumes fabs(x) >= DBL_EPSILON.
+ */
+static inline double
+sum_squares(double x, double y)
+{
+
+	/* Avoid underflow when y is small. */
+	if (y < SQRT_MIN)
+		return (x * x);
+
+	return (x * x + y * y);
+}
+
+/*
+ * real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y).
+ * Assumes x and y are not NaN, and one of x and y is larger than
+ * RECIP_EPSILON.  We avoid unwarranted underflow.  It is important to not use
+ * the code creal(1/z), because the imaginary part may produce an unwanted
+ * underflow.
+ * This is only called in a context where inexact is always raised before
+ * the call, so no effort is made to avoid or force inexact.
+ */
+static inline double
+real_part_reciprocal(double x, double y)
+{
+	double scale;
+	uint32_t hx, hy;
+	int32_t ix, iy;
+
+	/*
+	 * This code is inspired by the C99 document n1124.pdf, Section G.5.1,
+	 * example 2.
+	 */
+	GET_HIGH_WORD(hx, x);
+	ix = hx & 0x7ff00000;
+	GET_HIGH_WORD(hy, y);
+	iy = hy & 0x7ff00000;
+#undef BIAS
+#define	BIAS	(DBL_MAX_EXP - 1)
+/* XXX more guard digits are useful iff there is extra precision. */
+#define	CUTOFF	(DBL_MANT_DIG / 2 + 1)	/* just half or 1 guard digit */
+	if (ix - iy >= CUTOFF << 20 || isinf(x))
+		return (1 / x);		/* +-Inf -> +-0 is special */
+	if (iy - ix >= CUTOFF << 20)
+		return (x / y / y);	/* should avoid double div, but hard */
+	if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20)
+		return (x / (x * x + y * y));
+	scale = 1;
+	SET_HIGH_WORD(scale, 0x7ff00000 - ix);	/* 2**(1-ilogb(x)) */
+	x *= scale;
+	y *= scale;
+	return (x / (x * x + y * y) * scale);
+}
+
+/*
+ * catanh(z) = log((1+z)/(1-z)) / 2
+ *           = log1p(4*x / |z-1|^2) / 4
+ *             + I * atan2(2*y, (1-x)*(1+x)-y*y) / 2
+ *
+ * catanh(z) = z + O(z^3)   as z -> 0
+ *
+ * catanh(z) = 1/z + sign(y)*I*PI/2 + O(1/z^3)   as z -> infinity
+ * The above formula works for the real part as well, because
+ * Re(catanh(z)) = x/|z|^2 + O(x/z^4)
+ *    as z -> infinity, uniformly in x
+ */
+double complex
+catanh(double complex z)
+{
+	double x, y, ax, ay, rx, ry;
+
+	x = creal(z);
+	y = cimag(z);
+	ax = fabs(x);
+	ay = fabs(y);
+
+	/* This helps handle many cases. */
+	if (y == 0 && ax <= 1)
+		return (CMPLX(atanh(x), y));
+
+	/* To ensure the same accuracy as atan(), and to filter out z = 0. */
+	if (x == 0)
+		return (CMPLX(x, atan(y)));
+
+	if (isnan(x) || isnan(y)) {
+		/* catanh(+-Inf + I*NaN) = +-0 + I*NaN */
+		if (isinf(x))
+			return (CMPLX(copysign(0, x), y + y));
+		/* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
+		if (isinf(y))
+			return (CMPLX(copysign(0, x),
+			    copysign(pio2_hi + pio2_lo, y)));
+		/*
+		 * All other cases involving NaN return NaN + I*NaN.
+		 * C99 leaves it optional whether to raise invalid if one of
+		 * the arguments is not NaN, so we opt not to raise it.
+		 */
+		return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
+		return (CMPLX(real_part_reciprocal(x, y),
+		    copysign(pio2_hi + pio2_lo, y)));
+
+	if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
+		/*
+		 * z = 0 was filtered out above.  All other cases must raise
+		 * inexact, but this is the only case that needs to do it
+		 * explicitly.
+		 */
+		raise_inexact();
+		return (z);
+	}
+
+	if (ax == 1 && ay < DBL_EPSILON)
+		rx = (m_ln2 - log(ay)) / 2;
+	else
+		rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4;
+
+	if (ax == 1)
+		ry = atan2(2, -ay) / 2;
+	else if (ay < DBL_EPSILON)
+		ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2;
+	else
+		ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
+
+	return (CMPLX(copysign(rx, x), copysign(ry, y)));
+}
+
+/*
+ * catan(z) = reverse(catanh(reverse(z)))
+ * where reverse(x + I*y) = y + I*x = I*conj(z).
+ */
+double complex
+catan(double complex z)
+{
+	double complex w = catanh(CMPLX(cimag(z), creal(z)));
+
+	return (CMPLX(cimag(w), creal(w)));
+}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(cacosh, cacoshl);
+__weak_reference(cacos, cacosl);
+__weak_reference(casinh, casinhl);
+__weak_reference(casin, casinl);
+__weak_reference(catanh, catanhl);
+__weak_reference(catan, catanl);
+#endif

libc/tinymath/catrigf.c

@@ -0,0 +1,377 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│ vi: set noet ft=c ts=8 sw=8 fenc=utf-8                                   :vi │
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│ FreeBSD lib/msun/src/catrigf.c                                               │
+│ Copyright (c) 2012 Stephen Montgomery-Smith <[email protected]>            │
+│                                                                              │
+│ Redistribution and use in source and binary forms, with or without           │
+│ modification, are permitted provided that the following conditions           │
+│ are met:                                                                     │
+│ 1. Redistributions of source code must retain the above copyright            │
+│    notice, this list of conditions and the following disclaimer.             │
+│ 2. Redistributions in binary form must reproduce the above copyright         │
+│    notice, this list of conditions and the following disclaimer in the       │
+│    documentation and/or other materials provided with the distribution.      │
+│                                                                              │
+│ THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND       │
+│ ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE        │
+│ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE   │
+│ ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE      │
+│ FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL   │
+│ DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS      │
+│ OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)        │
+│ HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT   │
+│ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY    │
+│ OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF       │
+│ SUCH DAMAGE.                                                                 │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/tinymath/freebsd.internal.h"
+__static_yoink("freebsd_complex_notice");
+
+#undef isinf
+#define isinf(x)	(fabsf(x) == INFINITY)
+#undef isnan
+#define isnan(x)	((x) != (x))
+#define	raise_inexact()	do { volatile float _j = 1 + tiny; (void)_j; } while(0)
+#undef signbit
+#define signbit(x)	(__builtin_signbitf(x))
+
+static const float
+A_crossover =		10,
+B_crossover =		0.6417,
+FOUR_SQRT_MIN =		0x1p-61,
+QUARTER_SQRT_MAX =	0x1p61,
+m_e =			2.7182818285e0,		/*  0xadf854.0p-22 */
+m_ln2 =			6.9314718056e-1,	/*  0xb17218.0p-24 */
+pio2_hi =		1.5707962513e0,		/*  0xc90fda.0p-23 */
+RECIP_EPSILON =		1 / FLT_EPSILON,
+SQRT_3_EPSILON =	5.9801995673e-4,	/*  0x9cc471.0p-34 */
+SQRT_6_EPSILON =	8.4572793338e-4,	/*  0xddb3d7.0p-34 */
+SQRT_MIN =		0x1p-63;
+
+static const volatile float
+pio2_lo =		7.5497899549e-8,	/*  0xa22169.0p-47 */
+tiny =			0x1p-100;
+
+static float complex clog_for_large_values(float complex z);
+
+static inline float
+f(float a, float b, float hypot_a_b)
+{
+	if (b < 0)
+		return ((hypot_a_b - b) / 2);
+	if (b == 0)
+		return (a / 2);
+	return (a * a / (hypot_a_b + b) / 2);
+}
+
+static inline void
+do_hard_work(float x, float y, float *rx, int *B_is_usable, float *B,
+    float *sqrt_A2my2, float *new_y)
+{
+	float R, S, A;
+	float Am1, Amy;
+
+	R = hypotf(x, y + 1);
+	S = hypotf(x, y - 1);
+
+	A = (R + S) / 2;
+	if (A < 1)
+		A = 1;
+
+	if (A < A_crossover) {
+		if (y == 1 && x < FLT_EPSILON * FLT_EPSILON / 128) {
+			*rx = sqrtf(x);
+		} else if (x >= FLT_EPSILON * fabsf(y - 1)) {
+			Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
+			*rx = log1pf(Am1 + sqrtf(Am1 * (A + 1)));
+		} else if (y < 1) {
+			*rx = x / sqrtf((1 - y) * (1 + y));
+		} else {
+			*rx = log1pf((y - 1) + sqrtf((y - 1) * (y + 1)));
+		}
+	} else {
+		*rx = logf(A + sqrtf(A * A - 1));
+	}
+
+	*new_y = y;
+
+	if (y < FOUR_SQRT_MIN) {
+		*B_is_usable = 0;
+		*sqrt_A2my2 = A * (2 / FLT_EPSILON);
+		*new_y = y * (2 / FLT_EPSILON);
+		return;
+	}
+
+	*B = y / A;
+	*B_is_usable = 1;
+
+	if (*B > B_crossover) {
+		*B_is_usable = 0;
+		if (y == 1 && x < FLT_EPSILON / 128) {
+			*sqrt_A2my2 = sqrtf(x) * sqrtf((A + y) / 2);
+		} else if (x >= FLT_EPSILON * fabsf(y - 1)) {
+			Amy = f(x, y + 1, R) + f(x, y - 1, S);
+			*sqrt_A2my2 = sqrtf(Amy * (A + y));
+		} else if (y > 1) {
+			*sqrt_A2my2 = x * (4 / FLT_EPSILON / FLT_EPSILON) * y /
+			    sqrtf((y + 1) * (y - 1));
+			*new_y = y * (4 / FLT_EPSILON / FLT_EPSILON);
+		} else {
+			*sqrt_A2my2 = sqrtf((1 - y) * (1 + y));
+		}
+	}
+}
+
+float complex
+casinhf(float complex z)
+{
+	float x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
+	int B_is_usable;
+	float complex w;
+
+	x = crealf(z);
+	y = cimagf(z);
+	ax = fabsf(x);
+	ay = fabsf(y);
+
+	if (isnan(x) || isnan(y)) {
+		if (isinf(x))
+			return (CMPLXF(x, y + y));
+		if (isinf(y))
+			return (CMPLXF(y, x + x));
+		if (y == 0)
+			return (CMPLXF(x + x, y));
+		return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+		if (signbit(x) == 0)
+			w = clog_for_large_values(z) + m_ln2;
+		else
+			w = clog_for_large_values(-z) + m_ln2;
+		return (CMPLXF(copysignf(crealf(w), x),
+		    copysignf(cimagf(w), y)));
+	}
+
+	if (x == 0 && y == 0)
+		return (z);
+
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+		return (z);
+
+	do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
+	if (B_is_usable)
+		ry = asinf(B);
+	else
+		ry = atan2f(new_y, sqrt_A2my2);
+	return (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
+}
+
+float complex
+casinf(float complex z)
+{
+	float complex w = casinhf(CMPLXF(cimagf(z), crealf(z)));
+
+	return (CMPLXF(cimagf(w), crealf(w)));
+}
+
+float complex
+cacosf(float complex z)
+{
+	float x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
+	int sx, sy;
+	int B_is_usable;
+	float complex w;
+
+	x = crealf(z);
+	y = cimagf(z);
+	sx = signbit(x);
+	sy = signbit(y);
+	ax = fabsf(x);
+	ay = fabsf(y);
+
+	if (isnan(x) || isnan(y)) {
+		if (isinf(x))
+			return (CMPLXF(y + y, -INFINITY));
+		if (isinf(y))
+			return (CMPLXF(x + x, -y));
+		if (x == 0)
+			return (CMPLXF(pio2_hi + pio2_lo, y + y));
+		return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+		w = clog_for_large_values(z);
+		rx = fabsf(cimagf(w));
+		ry = crealf(w) + m_ln2;
+		if (sy == 0)
+			ry = -ry;
+		return (CMPLXF(rx, ry));
+	}
+
+	if (x == 1 && y == 0)
+		return (CMPLXF(0, -y));
+
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+		return (CMPLXF(pio2_hi - (x - pio2_lo), -y));
+
+	do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
+	if (B_is_usable) {
+		if (sx == 0)
+			rx = acosf(B);
+		else
+			rx = acosf(-B);
+	} else {
+		if (sx == 0)
+			rx = atan2f(sqrt_A2mx2, new_x);
+		else
+			rx = atan2f(sqrt_A2mx2, -new_x);
+	}
+	if (sy == 0)
+		ry = -ry;
+	return (CMPLXF(rx, ry));
+}
+
+float complex
+cacoshf(float complex z)
+{
+	float complex w;
+	float rx, ry;
+
+	w = cacosf(z);
+	rx = crealf(w);
+	ry = cimagf(w);
+	if (isnan(rx) && isnan(ry))
+		return (CMPLXF(ry, rx));
+	if (isnan(rx))
+		return (CMPLXF(fabsf(ry), rx));
+	if (isnan(ry))
+		return (CMPLXF(ry, ry));
+	return (CMPLXF(fabsf(ry), copysignf(rx, cimagf(z))));
+}
+
+static float complex
+clog_for_large_values(float complex z)
+{
+	float x, y;
+	float ax, ay, t;
+
+	x = crealf(z);
+	y = cimagf(z);
+	ax = fabsf(x);
+	ay = fabsf(y);
+	if (ax < ay) {
+		t = ax;
+		ax = ay;
+		ay = t;
+	}
+
+	if (ax > FLT_MAX / 2)
+		return (CMPLXF(logf(hypotf(x / m_e, y / m_e)) + 1,
+		    atan2f(y, x)));
+
+	if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
+		return (CMPLXF(logf(hypotf(x, y)), atan2f(y, x)));
+
+	return (CMPLXF(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
+}
+
+static inline float
+sum_squares(float x, float y)
+{
+
+	if (y < SQRT_MIN)
+		return (x * x);
+
+	return (x * x + y * y);
+}
+
+static inline float
+real_part_reciprocal(float x, float y)
+{
+	float scale;
+	uint32_t hx, hy;
+	int32_t ix, iy;
+
+	GET_FLOAT_WORD(hx, x);
+	ix = hx & 0x7f800000;
+	GET_FLOAT_WORD(hy, y);
+	iy = hy & 0x7f800000;
+#undef BIAS
+#define	BIAS	(FLT_MAX_EXP - 1)
+#define	CUTOFF	(FLT_MANT_DIG / 2 + 1)
+	if (ix - iy >= CUTOFF << 23 || isinf(x))
+		return (1 / x);
+	if (iy - ix >= CUTOFF << 23)
+		return (x / y / y);
+	if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23)
+		return (x / (x * x + y * y));
+	SET_FLOAT_WORD(scale, 0x7f800000 - ix);
+	x *= scale;
+	y *= scale;
+	return (x / (x * x + y * y) * scale);
+}
+
+float complex
+catanhf(float complex z)
+{
+	float x, y, ax, ay, rx, ry;
+
+	x = crealf(z);
+	y = cimagf(z);
+	ax = fabsf(x);
+	ay = fabsf(y);
+
+	if (y == 0 && ax <= 1)
+		return (CMPLXF(atanhf(x), y));
+
+	if (x == 0)
+		return (CMPLXF(x, atanf(y)));
+
+	if (isnan(x) || isnan(y)) {
+		if (isinf(x))
+			return (CMPLXF(copysignf(0, x), y + y));
+		if (isinf(y))
+			return (CMPLXF(copysignf(0, x),
+			    copysignf(pio2_hi + pio2_lo, y)));
+		return (CMPLXF(nan_mix(x, y), nan_mix(x, y)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
+		return (CMPLXF(real_part_reciprocal(x, y),
+		    copysignf(pio2_hi + pio2_lo, y)));
+
+	if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
+		raise_inexact();
+		return (z);
+	}
+
+	if (ax == 1 && ay < FLT_EPSILON)
+		rx = (m_ln2 - logf(ay)) / 2;
+	else
+		rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4;
+
+	if (ax == 1)
+		ry = atan2f(2, -ay) / 2;
+	else if (ay < FLT_EPSILON)
+		ry = atan2f(2 * ay, (1 - ax) * (1 + ax)) / 2;
+	else
+		ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
+
+	return (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
+}
+
+float complex
+catanf(float complex z)
+{
+	float complex w = catanhf(CMPLXF(cimagf(z), crealf(z)));
+
+	return (CMPLXF(cimagf(w), crealf(w)));
+}

libc/tinymath/cos.c

@@ -78,7 +78,6 @@ __static_yoink("fdlibm_notice");
 
 /**
  * Returns cosine of 𝑥.
- * @note should take ~5ns
  */
 double cos(double x)
 {

libc/tinymath/cosf.c

@@ -1,9 +1,9 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=8 fenc=utf-8                                   :vi │
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8                               :vi │
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
-│  Musl Libc                                                                   │
-│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│  Optimized Routines                                                          │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │
@@ -25,92 +25,63 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/math.h"
-#include "libc/tinymath/complex.internal.h"
-#include "libc/tinymath/feval.internal.h"
-#include "libc/tinymath/kernel.internal.h"
-__static_yoink("freebsd_libm_notice");
-__static_yoink("musl_libc_notice");
-__static_yoink("fdlibm_notice");
-
-
-/* origin: FreeBSD /usr/src/lib/msun/src/s_cosf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
- * Optimized by Bruce D. Evans.
- */
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
-
-/* Small multiples of pi/2 rounded to double precision. */
-static const double
-c1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
-c2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
-c3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
-c4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */
+#include "libc/tinymath/sincosf.internal.h"
+__static_yoink("arm_optimized_routines_notice");
 
 /**
- * Returns cosine of 𝑥.
- * @note should take about 5ns
+ * Returns cosine of y.
+ *
+ * This is a fast cosf implementation. The worst-case ULP is 0.5607, and
+ * the maximum relative error is 0.5303 * 2^-23. A single-step range
+ * reduction is used for small values. Large inputs have their range
+ * reduced using fast integer arithmetic.
+ *
+ * @raise EDOM and FE_INVALID if y is an infinity
  */
-float cosf(float x)
+float
+cosf (float y)
 {
-	double y;
-	uint32_t ix;
-	unsigned n, sign;
+  double x = y;
+  double s;
+  int n;
+  const sincos_t *p = &__sincosf_table[0];
 
-	GET_FLOAT_WORD(ix, x);
-	sign = ix >> 31;
-	ix &= 0x7fffffff;
+  if (abstop12 (y) < abstop12 (pio4f))
+    {
+      double x2 = x * x;
 
-	if (ix <= 0x3f490fda) {  /* |x| ~<= pi/4 */
-		if (ix < 0x39800000) {  /* |x| < 2**-12 */
-			/* raise inexact if x != 0 */
-			FORCE_EVAL(x + 0x1p120f);
-			return 1.0f;
-		}
-		return __cosdf(x);
-	}
-	if (ix <= 0x407b53d1) {  /* |x| ~<= 5*pi/4 */
-		if (ix > 0x4016cbe3)  /* |x|  ~> 3*pi/4 */
-			return -__cosdf(sign ? x+c2pio2 : x-c2pio2);
-		else {
-			if (sign)
-				return __sindf(x + c1pio2);
-			else
-				return __sindf(c1pio2 - x);
-		}
-	}
-	if (ix <= 0x40e231d5) {  /* |x| ~<= 9*pi/4 */
-		if (ix > 0x40afeddf)  /* |x| ~> 7*pi/4 */
-			return __cosdf(sign ? x+c4pio2 : x-c4pio2);
-		else {
-			if (sign)
-				return __sindf(-x - c3pio2);
-			else
-				return __sindf(x - c3pio2);
-		}
-	}
+      if (unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
+	return 1.0f;
 
-	/* cos(Inf or NaN) is NaN */
-	if (ix >= 0x7f800000)
-		return x-x;
+      return sinf_poly (x, x2, p, 1);
+    }
+  else if (likely (abstop12 (y) < abstop12 (120.0f)))
+    {
+      x = reduce_fast (x, p, &n);
 
-	/* general argument reduction needed */
-	n = __rem_pio2f(x,&y);
-	switch (n&3) {
-	case 0: return  __cosdf(y);
-	case 1: return  __sindf(-y);
-	case 2: return -__cosdf(y);
-	default:
-		return  __sindf(y);
-	}
+      /* Setup the signs for sin and cos.  */
+      s = p->sign[n & 3];
+
+      if (n & 2)
+	p = &__sincosf_table[1];
+
+      return sinf_poly (x * s, x * x, p, n ^ 1);
+    }
+  else if (abstop12 (y) < abstop12 (INFINITY))
+    {
+      uint32_t xi = asuint (y);
+      int sign = xi >> 31;
+
+      x = reduce_large (xi, &n);
+
+      /* Setup signs for sin and cos - include original sign.  */
+      s = p->sign[(n + sign) & 3];
+
+      if ((n + sign) & 2)
+	p = &__sincosf_table[1];
+
+      return sinf_poly (x * s, x * x, p, n ^ 1);
+    }
+  else
+    return __math_invalidf (y);
 }

libc/tinymath/erf.c

@@ -1,9 +1,9 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=8 fenc=utf-8                                   :vi │
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8                               :vi │
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
-│  Musl Libc                                                                   │
-│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│  Optimized Routines                                                          │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │
@@ -25,310 +25,247 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/math.h"
-__static_yoink("musl_libc_notice");
-__static_yoink("fdlibm_notice");
+#include "libc/tinymath/arm.internal.h"
+__static_yoink("arm_optimized_routines_notice");
 
-/* origin: FreeBSD /usr/src/lib/msun/src/s_erf.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* double erf(double x)
- * double erfc(double x)
- *                           x
- *                    2      |\
- *     erf(x)  =  ---------  | exp(-t*t)dt
- *                 sqrt(pi) \|
- *                           0
- *
- *     erfc(x) =  1-erf(x)
- *  Note that
- *              erf(-x) = -erf(x)
- *              erfc(-x) = 2 - erfc(x)
- *
- * Method:
- *      1. For |x| in [0, 0.84375]
- *          erf(x)  = x + x*R(x^2)
- *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
- *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
- *         where R = P/Q where P is an odd poly of degree 8 and
- *         Q is an odd poly of degree 10.
- *                                               -57.90
- *                      | R - (erf(x)-x)/x | <= 2
- *
- *
- *         Remark. The formula is derived by noting
- *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
- *         and that
- *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
- *         is close to one. The interval is chosen because the fix
- *         point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
- *         near 0.6174), and by some experiment, 0.84375 is chosen to
- *         guarantee the error is less than one ulp for erf.
- *
- *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
- *         c = 0.84506291151 rounded to single (24 bits)
- *              erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
- *              erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
- *                        1+(c+P1(s)/Q1(s))    if x < 0
- *              |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
- *         Remark: here we use the taylor series expansion at x=1.
- *              erf(1+s) = erf(1) + s*Poly(s)
- *                       = 0.845.. + P1(s)/Q1(s)
- *         That is, we use rational approximation to approximate
- *                      erf(1+s) - (c = (single)0.84506291151)
- *         Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
- *         where
- *              P1(s) = degree 6 poly in s
- *              Q1(s) = degree 6 poly in s
- *
- *      3. For x in [1.25,1/0.35(~2.857143)],
- *              erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
- *              erf(x)  = 1 - erfc(x)
- *         where
- *              R1(z) = degree 7 poly in z, (z=1/x^2)
- *              S1(z) = degree 8 poly in z
- *
- *      4. For x in [1/0.35,28]
- *              erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
- *                      = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
- *                      = 2.0 - tiny            (if x <= -6)
- *              erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6, else
- *              erf(x)  = sign(x)*(1.0 - tiny)
- *         where
- *              R2(z) = degree 6 poly in z, (z=1/x^2)
- *              S2(z) = degree 7 poly in z
- *
- *      Note1:
- *         To compute exp(-x*x-0.5625+R/S), let s be a single
- *         precision number and s := x; then
- *              -x*x = -s*s + (s-x)*(s+x)
- *              exp(-x*x-0.5626+R/S) =
- *                      exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
- *      Note2:
- *         Here 4 and 5 make use of the asymptotic series
- *                        exp(-x*x)
- *              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
- *                        x*sqrt(pi)
- *         We use rational approximation to approximate
- *              g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
- *         Here is the error bound for R1/S1 and R2/S2
- *              |R1/S1 - f(x)|  < 2**(-62.57)
- *              |R2/S2 - f(x)|  < 2**(-61.52)
- *
- *      5. For inf > x >= 28
- *              erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
- *              erfc(x) = tiny*tiny (raise underflow) if x > 0
- *                      = 2 - tiny if x<0
- *
- *      7. Special case:
- *              erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
- *              erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
- *              erfc/erf(NaN) is NaN
- */
+#define TwoOverSqrtPiMinusOne 0x1.06eba8214db69p-3
+#define C 0x1.b0ac16p-1
+#define PA __erf_data.erf_poly_A
+#define NA __erf_data.erf_ratio_N_A
+#define DA __erf_data.erf_ratio_D_A
+#define NB __erf_data.erf_ratio_N_B
+#define DB __erf_data.erf_ratio_D_B
+#define PC __erf_data.erfc_poly_C
+#define PD __erf_data.erfc_poly_D
+#define PE __erf_data.erfc_poly_E
+#define PF __erf_data.erfc_poly_F
 
-static const double
-erx  = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
-/*
- * Coefficients for approximation to  erf on [0,0.84375]
- */
-efx8 =  1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
-pp0  =  1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
-pp1  = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
-pp2  = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
-pp3  = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
-pp4  = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
-qq1  =  3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
-qq2  =  6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
-qq3  =  5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
-qq4  =  1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
-qq5  = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
-/*
- * Coefficients for approximation to  erf  in [0.84375,1.25]
- */
-pa0  = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
-pa1  =  4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
-pa2  = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
-pa3  =  3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
-pa4  = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
-pa5  =  3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
-pa6  = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
-qa1  =  1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
-qa2  =  5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
-qa3  =  7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
-qa4  =  1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
-qa5  =  1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
-qa6  =  1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
-/*
- * Coefficients for approximation to  erfc in [1.25,1/0.35]
- */
-ra0  = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
-ra1  = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
-ra2  = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
-ra3  = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
-ra4  = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
-ra5  = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
-ra6  = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
-ra7  = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
-sa1  =  1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
-sa2  =  1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
-sa3  =  4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
-sa4  =  6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
-sa5  =  4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
-sa6  =  1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
-sa7  =  6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
-sa8  = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
-/*
- * Coefficients for approximation to  erfc in [1/.35,28]
- */
-rb0  = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
-rb1  = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
-rb2  = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
-rb3  = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
-rb4  = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
-rb5  = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
-rb6  = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
-sb1  =  3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
-sb2  =  3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
-sb3  =  1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
-sb4  =  3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
-sb5  =  2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
-sb6  =  4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
-sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
-
-#define asuint(f) ((union{float _f; uint32_t _i;}){f})._i
-#define asfloat(i) ((union{uint32_t _i; float _f;}){i})._f
-#define asuint64(f) ((union{double _f; uint64_t _i;}){f})._i
-#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)
-#define GET_HIGH_WORD(hi,d)                       \
-do {                                              \
-  (hi) = asuint64(d) >> 32;                       \
-} while (0)
-#define GET_LOW_WORD(lo,d)                        \
-do {                                              \
-  (lo) = (uint32_t)asuint64(d);                   \
-} while (0)
-#define SET_HIGH_WORD(d,hi)                       \
-  INSERT_WORDS(d, hi, (uint32_t)asuint64(d))
-#define SET_LOW_WORD(d,lo)                        \
-  INSERT_WORDS(d, asuint64(d)>>32, lo)
-
-static double erfc1(double x)
+/* Top 32 bits of a double.  */
+static inline uint32_t
+top32 (double x)
 {
-	double_t s,P,Q;
-
-	s = fabs(x) - 1;
-	P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
-	Q = 1+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
-	return 1 - erx - P/Q;
-}
-
-static double erfc2(uint32_t ix, double x)
-{
-	double_t s,R,S;
-	double z;
-
-	if (ix < 0x3ff40000)  /* |x| < 1.25 */
-		return erfc1(x);
-
-	x = fabs(x);
-	s = 1/(x*x);
-	if (ix < 0x4006db6d) {  /* |x| < 1/.35 ~ 2.85714 */
-		R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
-		     ra5+s*(ra6+s*ra7))))));
-		S = 1.0+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
-		     sa5+s*(sa6+s*(sa7+s*sa8)))))));
-	} else {                /* |x| > 1/.35 */
-		R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
-		     rb5+s*rb6)))));
-		S = 1.0+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
-		     sb5+s*(sb6+s*sb7))))));
-	}
-	z = x;
-	SET_LOW_WORD(z,0);
-	return exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S)/x;
+  return asuint64 (x) >> 32;
 }
 
 /**
- * Returns error function of 𝑥.
+ * Returns error function of x.
+ *
+ * Highest measured error is 1.01 ULPs at 0x1.39956ac43382fp+0.
+ *
+ * @raise ERANGE on underflow
  */
-double erf(double x)
+double
+erf (double x)
 {
-	double r,s,z,y;
-	uint32_t ix;
-	int sign;
+  /* Get top word and sign.  */
+  uint32_t ix = top32 (x);
+  uint32_t ia = ix & 0x7fffffff;
+  uint32_t sign = ix >> 31;
 
-	GET_HIGH_WORD(ix, x);
-	sign = ix>>31;
-	ix &= 0x7fffffff;
-	if (ix >= 0x7ff00000) {
-		/* erf(nan)=nan, erf(+-inf)=+-1 */
-		return 1-2*sign + 1/x;
+  /* Normalized and subnormal cases */
+  if (ia < 0x3feb0000)
+    { /* a = |x| < 0.84375.  */
+
+      if (ia < 0x3e300000)
+	{ /* a < 2^(-28).  */
+	  if (ia < 0x00800000)
+	    { /* a < 2^(-1015).  */
+	      double y =  fma (TwoOverSqrtPiMinusOne, x, x);
+	      return check_uflow (y);
+	    }
+	  return x + TwoOverSqrtPiMinusOne * x;
 	}
-	if (ix < 0x3feb0000) {  /* |x| < 0.84375 */
-		if (ix < 0x3e300000) {  /* |x| < 2**-28 */
-			/* avoid underflow */
-			return 0.125*(8*x + efx8*x);
-		}
-		z = x*x;
-		r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
-		s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
-		y = r/s;
-		return x + x*y;
+
+      double x2 = x * x;
+
+      if (ia < 0x3fe00000)
+	{ /* a < 0.5  - Use polynomial approximation.  */
+	  double r1 = fma (x2, PA[1], PA[0]);
+	  double r2 = fma (x2, PA[3], PA[2]);
+	  double r3 = fma (x2, PA[5], PA[4]);
+	  double r4 = fma (x2, PA[7], PA[6]);
+	  double r5 = fma (x2, PA[9], PA[8]);
+	  double x4 = x2 * x2;
+	  double r = r5;
+	  r = fma (x4, r, r4);
+	  r = fma (x4, r, r3);
+	  r = fma (x4, r, r2);
+	  r = fma (x4, r, r1);
+	  return fma (r, x, x); /* This fma is crucial for accuracy.  */
 	}
-	if (ix < 0x40180000)  /* 0.84375 <= |x| < 6 */
-		y = 1 - erfc2(ix,x);
-	else
-		y = 1 - 0x1p-1022;
-	return sign ? -y : y;
+      else
+	{ /* 0.5 <= a < 0.84375 - Use rational approximation.  */
+	  double x4, x8, r1n, r2n, r1d, r2d, r3d;
+
+	  r1n = fma (x2, NA[1], NA[0]);
+	  x4 = x2 * x2;
+	  r2n = fma (x2, NA[3], NA[2]);
+	  x8 = x4 * x4;
+	  r1d = fma (x2, DA[0], 1.0);
+	  r2d = fma (x2, DA[2], DA[1]);
+	  r3d = fma (x2, DA[4], DA[3]);
+	  double P = r1n + x4 * r2n + x8 * NA[4];
+	  double Q = r1d + x4 * r2d + x8 * r3d;
+	  return fma (P / Q, x, x);
+	}
+    }
+  else if (ia < 0x3ff40000)
+    { /* 0.84375 <= |x| < 1.25.  */
+      double a2, a4, a6, r1n, r2n, r3n, r4n, r1d, r2d, r3d, r4d;
+      double a = fabs (x) - 1.0;
+      r1n = fma (a, NB[1], NB[0]);
+      a2 = a * a;
+      r1d = fma (a, DB[0], 1.0);
+      a4 = a2 * a2;
+      r2n = fma (a, NB[3], NB[2]);
+      a6 = a4 * a2;
+      r2d = fma (a, DB[2], DB[1]);
+      r3n = fma (a, NB[5], NB[4]);
+      r3d = fma (a, DB[4], DB[3]);
+      r4n = NB[6];
+      r4d = DB[5];
+      double P = r1n + a2 * r2n + a4 * r3n + a6 * r4n;
+      double Q = r1d + a2 * r2d + a4 * r3d + a6 * r4d;
+      if (sign)
+	return -C - P / Q;
+      else
+	return C + P / Q;
+    }
+  else if (ia < 0x40000000)
+    { /* 1.25 <= |x| < 2.0.  */
+      double a = fabs (x);
+      a = a - 1.25;
+
+      double r1 = fma (a, PC[1], PC[0]);
+      double r2 = fma (a, PC[3], PC[2]);
+      double r3 = fma (a, PC[5], PC[4]);
+      double r4 = fma (a, PC[7], PC[6]);
+      double r5 = fma (a, PC[9], PC[8]);
+      double r6 = fma (a, PC[11], PC[10]);
+      double r7 = fma (a, PC[13], PC[12]);
+      double r8 = fma (a, PC[15], PC[14]);
+
+      double a2 = a * a;
+
+      double r = r8;
+      r = fma (a2, r, r7);
+      r = fma (a2, r, r6);
+      r = fma (a2, r, r5);
+      r = fma (a2, r, r4);
+      r = fma (a2, r, r3);
+      r = fma (a2, r, r2);
+      r = fma (a2, r, r1);
+
+      if (sign)
+	return -1.0 + r;
+      else
+	return 1.0 - r;
+    }
+  else if (ia < 0x400a0000)
+    { /* 2 <= |x| < 3.25.  */
+      double a = fabs (x);
+      a = fma (0.5, a, -1.0);
+
+      double r1 = fma (a, PD[1], PD[0]);
+      double r2 = fma (a, PD[3], PD[2]);
+      double r3 = fma (a, PD[5], PD[4]);
+      double r4 = fma (a, PD[7], PD[6]);
+      double r5 = fma (a, PD[9], PD[8]);
+      double r6 = fma (a, PD[11], PD[10]);
+      double r7 = fma (a, PD[13], PD[12]);
+      double r8 = fma (a, PD[15], PD[14]);
+      double r9 = fma (a, PD[17], PD[16]);
+
+      double a2 = a * a;
+
+      double r = r9;
+      r = fma (a2, r, r8);
+      r = fma (a2, r, r7);
+      r = fma (a2, r, r6);
+      r = fma (a2, r, r5);
+      r = fma (a2, r, r4);
+      r = fma (a2, r, r3);
+      r = fma (a2, r, r2);
+      r = fma (a2, r, r1);
+
+      if (sign)
+	return -1.0 + r;
+      else
+	return 1.0 - r;
+    }
+  else if (ia < 0x40100000)
+    { /* 3.25 <= |x| < 4.0.  */
+      double a = fabs (x);
+      a = a - 3.25;
+
+      double r1 = fma (a, PE[1], PE[0]);
+      double r2 = fma (a, PE[3], PE[2]);
+      double r3 = fma (a, PE[5], PE[4]);
+      double r4 = fma (a, PE[7], PE[6]);
+      double r5 = fma (a, PE[9], PE[8]);
+      double r6 = fma (a, PE[11], PE[10]);
+      double r7 = fma (a, PE[13], PE[12]);
+
+      double a2 = a * a;
+
+      double r = r7;
+      r = fma (a2, r, r6);
+      r = fma (a2, r, r5);
+      r = fma (a2, r, r4);
+      r = fma (a2, r, r3);
+      r = fma (a2, r, r2);
+      r = fma (a2, r, r1);
+
+      if (sign)
+	return -1.0 + r;
+      else
+	return 1.0 - r;
+    }
+  else if (ia < 0x4017a000)
+    { /* 4 <= |x| < 5.90625.  */
+      double a = fabs (x);
+      a = fma (0.5, a, -2.0);
+
+      double r1 = fma (a, PF[1], PF[0]);
+      double r2 = fma (a, PF[3], PF[2]);
+      double r3 = fma (a, PF[5], PF[4]);
+      double r4 = fma (a, PF[7], PF[6]);
+      double r5 = fma (a, PF[9], PF[8]);
+      double r6 = fma (a, PF[11], PF[10]);
+      double r7 = fma (a, PF[13], PF[12]);
+      double r8 = fma (a, PF[15], PF[14]);
+      double r9 = PF[16];
+
+      double a2 = a * a;
+
+      double r = r9;
+      r = fma (a2, r, r8);
+      r = fma (a2, r, r7);
+      r = fma (a2, r, r6);
+      r = fma (a2, r, r5);
+      r = fma (a2, r, r4);
+      r = fma (a2, r, r3);
+      r = fma (a2, r, r2);
+      r = fma (a2, r, r1);
+
+      if (sign)
+	return -1.0 + r;
+      else
+	return 1.0 - r;
+    }
+  else
+    {
+      /* Special cases : erf(nan)=nan, erf(+inf)=+1 and erf(-inf)=-1.  */
+      if (unlikely (ia >= 0x7ff00000))
+	return (double) (1.0 - (sign << 1)) + 1.0 / x;
+
+      if (sign)
+	return -1.0;
+      else
+	return 1.0;
+    }
 }
 
-/**
- * Returns complementary error function of 𝑥.
- */
-double erfc(double x)
-{
-	double r,s,z,y;
-	uint32_t ix;
-	int sign;
-
-	GET_HIGH_WORD(ix, x);
-	sign = ix>>31;
-	ix &= 0x7fffffff;
-	if (ix >= 0x7ff00000) {
-		/* erfc(nan)=nan, erfc(+-inf)=0,2 */
-		return 2*sign + 1/x;
-	}
-	if (ix < 0x3feb0000) {  /* |x| < 0.84375 */
-		if (ix < 0x3c700000)  /* |x| < 2**-56 */
-			return 1.0 - x;
-		z = x*x;
-		r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
-		s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
-		y = r/s;
-		if (sign || ix < 0x3fd00000) {  /* x < 1/4 */
-			return 1.0 - (x+x*y);
-		}
-		return 0.5 - (x - 0.5 + x*y);
-	}
-	if (ix < 0x403c0000) {  /* 0.84375 <= |x| < 28 */
-		return sign ? 2 - erfc2(ix,x) : erfc2(ix,x);
-	}
-	return sign ? 2 - 0x1p-1022 : 0x1p-1022*0x1p-1022;
-}
-
-#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+#if LDBL_MANT_DIG == 53
 __weak_reference(erf, erfl);
-__weak_reference(erfc, erfcl);
 #endif

libc/tinymath/erf_data.c

@@ -0,0 +1,105 @@
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8                               :vi │
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│  Optimized Routines                                                          │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
+│                                                                              │
+│  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/tinymath/arm.internal.h"
+__static_yoink("arm_optimized_routines_notice");
+
+/*
+Minimax approximation of erf
+*/
+const struct erf_data __erf_data = {
+.erf_poly_A = {
+#if ERF_POLY_A_NCOEFFS == 10
+0x1.06eba8214db68p-3, -0x1.812746b037948p-2, 0x1.ce2f21a03872p-4,
+-0x1.b82ce30e6548p-6, 0x1.565bcc360a2f2p-8, -0x1.c02d812bc979ap-11,
+0x1.f99bddfc1ebe9p-14, -0x1.f42c457cee912p-17, 0x1.b0e414ec20ee9p-20,
+-0x1.18c47fd143c5ep-23
+#endif
+},
+/* Rational approximation on [0x1p-28, 0.84375] */
+.erf_ratio_N_A = {
+0x1.06eba8214db68p-3, -0x1.4cd7d691cb913p-2, -0x1.d2a51dbd7194fp-6,
+-0x1.7a291236668e4p-8, -0x1.8ead6120016acp-16
+},
+.erf_ratio_D_A = {
+0x1.97779cddadc09p-2, 0x1.0a54c5536cebap-4, 0x1.4d022c4d36b0fp-8,
+0x1.15dc9221c1a1p-13, -0x1.09c4342a2612p-18
+},
+/* Rational approximation on [0.84375, 1.25] */
+.erf_ratio_N_B = {
+-0x1.359b8bef77538p-9, 0x1.a8d00ad92b34dp-2, -0x1.7d240fbb8c3f1p-2,
+0x1.45fca805120e4p-2, -0x1.c63983d3e28ecp-4, 0x1.22a36599795ebp-5,
+-0x1.1bf380a96073fp-9
+},
+.erf_ratio_D_B = {
+0x1.b3e6618eee323p-4, 0x1.14af092eb6f33p-1, 0x1.2635cd99fe9a7p-4,
+0x1.02660e763351fp-3, 0x1.bedc26b51dd1cp-7, 0x1.88b545735151dp-7
+},
+.erfc_poly_C = {
+#if ERFC_POLY_C_NCOEFFS == 16
+/* Generated using Sollya::remez(f(c*x+d), deg, [(a-d)/c;(b-d)/c], 1, 1e-16), [|D ...|] with deg=15 a=1.25 b=2 c=1 d=1.25 */
+0x1.3bcd133aa0ffcp-4, -0x1.e4652fadcb702p-3, 0x1.2ebf3dcca0446p-2,
+-0x1.571d01c62d66p-3, 0x1.93a9a8f5b3413p-8, 0x1.8281cbcc2cd52p-5,
+-0x1.5cffd86b4de16p-6, -0x1.db4ccf595053ep-9, 0x1.757cbf8684edap-8,
+-0x1.ce7dfd2a9e56ap-11, -0x1.99ee3bc5a3263p-11, 0x1.3c57cf9213f5fp-12,
+0x1.60692996bf254p-14, -0x1.6e44cb7c1fa2ap-14, 0x1.9d4484ac482b2p-16,
+-0x1.578c9e375d37p-19
+#endif
+},
+.erfc_poly_D = {
+#if ERFC_POLY_D_NCOEFFS == 18
+/* Generated using Sollya::remez(f(c*x+d), deg, [(a-d)/c;(b-d)/c], 1, 1e-16), [|D ...|] with deg=17 a=2 b=3.25 c=2 d=2 */
+0x1.328f5ec350e5p-8, -0x1.529b9e8cf8e99p-5, 0x1.529b9e8cd9e71p-3,
+-0x1.8b0ae3a023bf2p-2, 0x1.1a2c592599d82p-1, -0x1.ace732477e494p-2,
+-0x1.e1a06a27920ffp-6, 0x1.bae92a6d27af6p-2, -0x1.a15470fcf5ce7p-2,
+0x1.bafe45d18e213p-6, 0x1.0d950680d199ap-2, -0x1.8c9481e8f22e3p-3,
+-0x1.158450ed5c899p-4, 0x1.c01f2973b44p-3, -0x1.73ed2827546a7p-3,
+0x1.47733687d1ff7p-4, -0x1.2dec70d00b8e1p-6, 0x1.a947ab83cd4fp-10
+#endif
+},
+.erfc_poly_E = {
+#if ERFC_POLY_E_NCOEFFS == 14
+/* Generated using Sollya::remez(f(c*x+d), deg, [(a-d)/c;(b-d)/c], 1, 1e-16), [|D ...|] with deg=13 a=3.25 b=4 c=1 d=3.25 */
+0x1.20c13035539e4p-18, -0x1.e9b5e8d16df7ep-16, 0x1.8de3cd4733bf9p-14,
+-0x1.9aa48beb8382fp-13, 0x1.2c7d713370a9fp-12, -0x1.490b12110b9e2p-12,
+0x1.1459c5d989d23p-12, -0x1.64b28e9f1269p-13, 0x1.57c76d9d05cf8p-14,
+-0x1.bf271d9951cf8p-16, 0x1.db7ea4d4535c9p-19, 0x1.91c2e102d5e49p-20,
+-0x1.e9f0826c2149ep-21, 0x1.60eebaea236e1p-23
+#endif
+},
+.erfc_poly_F = {
+#if ERFC_POLY_F_NCOEFFS == 17
+/* Generated using Sollya::remez(f(c*x+d), deg, [(a-d)/c;(b-d)/c], 1, 1e-16), [|D ...|] with deg=16 a=4 b=5.90625 c=2 d=4 */
+0x1.08ddd130d1fa6p-26, -0x1.10b146f59ff06p-22, 0x1.10b135328b7b2p-19,
+-0x1.6039988e7575fp-17, 0x1.497d365e19367p-15, -0x1.da48d9afac83ep-14,
+0x1.1024c9b1fbb48p-12, -0x1.fc962e7066272p-12, 0x1.87297282d4651p-11,
+-0x1.f057b255f8c59p-11, 0x1.0228d0eee063p-10, -0x1.b1b21b84ec41cp-11,
+0x1.1ead8ae9e1253p-11, -0x1.1e708fba37fccp-12, 0x1.9559363991edap-14,
+-0x1.68c827b783d9cp-16, 0x1.2ec4adeccf4a2p-19
+#endif
+}
+};

libc/tinymath/erfc.c

@@ -0,0 +1,279 @@
+/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
+│ vi: set noet ft=c ts=8 sw=8 fenc=utf-8                                   :vi │
+╚──────────────────────────────────────────────────────────────────────────────╝
+│                                                                              │
+│ Copyright (c) 1992-2024 The FreeBSD Project                                  │
+│                                                                              │
+│ Redistribution and use in source and binary forms, with or without           │
+│ modification, are permitted provided that the following conditions           │
+│ are met:                                                                     │
+│ 1. Redistributions of source code must retain the above copyright            │
+│    notice, this list of conditions and the following disclaimer.             │
+│ 2. Redistributions in binary form must reproduce the above copyright         │
+│    notice, this list of conditions and the following disclaimer in the       │
+│    documentation and/or other materials provided with the distribution.      │
+│                                                                              │
+│ THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND       │
+│ ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE        │
+│ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE   │
+│ ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE      │
+│ FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL   │
+│ DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS      │
+│ OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)        │
+│ HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT   │
+│ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY    │
+│ OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF       │
+│ SUCH DAMAGE.                                                                 │
+│                                                                              │
+│ Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.            │
+│                                                                              │
+│ Developed at SunPro, a Sun Microsystems, Inc. business.                      │
+│ Permission to use, copy, modify, and distribute this                         │
+│ software is freely granted, provided that this notice                        │
+│ is preserved.                                                                │
+│                                                                              │
+╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/tinymath/freebsd.internal.h"
+__static_yoink("freebsd_libm_notice");
+__static_yoink("fdlibm_notice");
+
+/* double erf(double x)
+ * double erfc(double x)
+ *			     x
+ *		      2      |\
+ *     erf(x)  =  ---------  | exp(-t*t)dt
+ *	 	   sqrt(pi) \|
+ *			     0
+ *
+ *     erfc(x) =  1-erf(x)
+ *  Note that
+ *		erf(-x) = -erf(x)
+ *		erfc(-x) = 2 - erfc(x)
+ *
+ * Method:
+ *	1. For |x| in [0, 0.84375]
+ *	    erf(x)  = x + x*R(x^2)
+ *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
+ *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
+ *	   where R = P/Q where P is an odd poly of degree 8 and
+ *	   Q is an odd poly of degree 10.
+ *						 -57.90
+ *			| R - (erf(x)-x)/x | <= 2
+ *
+ *
+ *	   Remark. The formula is derived by noting
+ *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
+ *	   and that
+ *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
+ *	   is close to one. The interval is chosen because the fix
+ *	   point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
+ *	   near 0.6174), and by some experiment, 0.84375 is chosen to
+ * 	   guarantee the error is less than one ulp for erf.
+ *
+ *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
+ *         c = 0.84506291151 rounded to single (24 bits)
+ *         	erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
+ *         	erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
+ *			  1+(c+P1(s)/Q1(s))    if x < 0
+ *         	|P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
+ *	   Remark: here we use the taylor series expansion at x=1.
+ *		erf(1+s) = erf(1) + s*Poly(s)
+ *			 = 0.845.. + P1(s)/Q1(s)
+ *	   That is, we use rational approximation to approximate
+ *			erf(1+s) - (c = (single)0.84506291151)
+ *	   Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
+ *	   where
+ *		P1(s) = degree 6 poly in s
+ *		Q1(s) = degree 6 poly in s
+ *
+ *      3. For x in [1.25,1/0.35(~2.857143)],
+ *         	erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
+ *         	erf(x)  = 1 - erfc(x)
+ *	   where
+ *		R1(z) = degree 7 poly in z, (z=1/x^2)
+ *		S1(z) = degree 8 poly in z
+ *
+ *      4. For x in [1/0.35,28]
+ *         	erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
+ *			= 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
+ *			= 2.0 - tiny		(if x <= -6)
+ *         	erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6, else
+ *         	erf(x)  = sign(x)*(1.0 - tiny)
+ *	   where
+ *		R2(z) = degree 6 poly in z, (z=1/x^2)
+ *		S2(z) = degree 7 poly in z
+ *
+ *      Note1:
+ *	   To compute exp(-x*x-0.5625+R/S), let s be a single
+ *	   precision number and s := x; then
+ *		-x*x = -s*s + (s-x)*(s+x)
+ *	        exp(-x*x-0.5626+R/S) =
+ *			exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
+ *      Note2:
+ *	   Here 4 and 5 make use of the asymptotic series
+ *			  exp(-x*x)
+ *		erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
+ *			  x*sqrt(pi)
+ *	   We use rational approximation to approximate
+ *      	g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
+ *	   Here is the error bound for R1/S1 and R2/S2
+ *      	|R1/S1 - f(x)|  < 2**(-62.57)
+ *      	|R2/S2 - f(x)|  < 2**(-61.52)
+ *
+ *      5. For inf > x >= 28
+ *         	erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
+ *         	erfc(x) = tiny*tiny (raise underflow) if x > 0
+ *			= 2 - tiny if x<0
+ *
+ *      7. Special case:
+ *         	erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
+ *         	erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
+ *	   	erfc/erf(NaN) is NaN
+ */
+
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile double tiny= 1e-300;
+
+static const double
+half= 0.5,
+one = 1,
+two = 2,
+/* c = (float)0.84506291151 */
+erx =  8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
+/*
+ * In the domain [0, 2**-28], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-84.
+ */
+efx =  1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
+efx8=  1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
+/*
+ * Coefficients for approximation to erf on [0,0.84375]
+ */
+pp0  =  1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
+pp1  = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
+pp2  = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
+pp3  = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
+pp4  = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
+qq1  =  3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
+qq2  =  6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
+qq3  =  5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
+qq4  =  1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
+qq5  = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
+/*
+ * Coefficients for approximation to erf in [0.84375,1.25]
+ */
+pa0  = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
+pa1  =  4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
+pa2  = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
+pa3  =  3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
+pa4  = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
+pa5  =  3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
+pa6  = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
+qa1  =  1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
+qa2  =  5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
+qa3  =  7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
+qa4  =  1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
+qa5  =  1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
+qa6  =  1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
+/*
+ * Coefficients for approximation to erfc in [1.25,1/0.35]
+ */
+ra0  = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
+ra1  = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
+ra2  = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
+ra3  = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
+ra4  = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
+ra5  = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
+ra6  = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
+ra7  = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
+sa1  =  1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
+sa2  =  1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
+sa3  =  4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
+sa4  =  6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
+sa5  =  4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
+sa6  =  1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
+sa7  =  6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
+sa8  = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
+/*
+ * Coefficients for approximation to erfc in [1/.35,28]
+ */
+rb0  = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
+rb1  = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
+rb2  = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
+rb3  = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
+rb4  = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
+rb5  = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
+rb6  = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
+sb1  =  3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
+sb2  =  3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
+sb3  =  1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
+sb4  =  3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
+sb5  =  2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
+sb6  =  4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
+sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
+
+/**
+ * Returns complementary error function of x, i.e. 1.0 - erf(x).
+ */
+double
+erfc(double x)
+{
+	int32_t hx,ix;
+	double R,S,P,Q,s,y,z,r;
+	GET_HIGH_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if(ix>=0x7ff00000) {			/* erfc(nan)=nan */
+						/* erfc(+-inf)=0,2 */
+	    return (double)(((uint32_t)hx>>31)<<1)+one/x;
+	}
+
+	if(ix < 0x3feb0000) {		/* |x|<0.84375 */
+	    if(ix < 0x3c700000)  	/* |x|<2**-56 */
+		return one-x;
+	    z = x*x;
+	    r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
+	    s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
+	    y = r/s;
+	    if(hx < 0x3fd00000) {  	/* x<1/4 */
+		return one-(x+x*y);
+	    } else {
+		r = x*y;
+		r += (x-half);
+	        return half - r ;
+	    }
+	}
+	if(ix < 0x3ff40000) {		/* 0.84375 <= |x| < 1.25 */
+	    s = fabs(x)-one;
+	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
+	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
+	    if(hx>=0) {
+	        z  = one-erx; return z - P/Q;
+	    } else {
+		z = erx+P/Q; return one+z;
+	    }
+	}
+	if (ix < 0x403c0000) {		/* |x|<28 */
+	    x = fabs(x);
+ 	    s = one/(x*x);
+	    if(ix< 0x4006DB6D) {	/* |x| < 1/.35 ~ 2.857143*/
+		R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
+		S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		    s*sa8)))))));
+	    } else {			/* |x| >= 1/.35 ~ 2.857143 */
+		if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
+		R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
+		S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
+	    }
+	    z  = x;
+	    SET_LOW_WORD(z,0);
+	    r  =  exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
+	    if(hx>0) return r/x; else return two-r/x;
+	} else {
+	    if(hx>0) return tiny*tiny; else return two-tiny;
+	}
+}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(erfc, erfcl);
+#endif

libc/tinymath/erff.c

@@ -1,9 +1,9 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=8 fenc=utf-8                                   :vi │
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8                               :vi │
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
-│  Musl Libc                                                                   │
-│  Copyright © 2005-2020 Rich Felker, et al.                                   │
+│  Optimized Routines                                                          │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │
@@ -25,189 +25,99 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/math.h"
-__static_yoink("freebsd_libm_notice");
-__static_yoink("musl_libc_notice");
-__static_yoink("fdlibm_notice");
+#include "libc/tinymath/arm.internal.h"
+__static_yoink("arm_optimized_routines_notice");
 
-/* origin: FreeBSD /usr/src/lib/msun/src/s_erff.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 TwoOverSqrtPiMinusOne 0x1.06eba8p-3f
+#define A __erff_data.erff_poly_A
+#define B __erff_data.erff_poly_B
 
-#define asuint(f) ((union{float _f; uint32_t _i;}){f})._i
-#define asfloat(i) ((union{uint32_t _i; float _f;}){i})._f
-
-static const float
-erx  =  8.4506291151e-01, /* 0x3f58560b */
-/*
- * Coefficients for approximation to  erf on [0,0.84375]
- */
-efx8 =  1.0270333290e+00, /* 0x3f8375d4 */
-pp0  =  1.2837916613e-01, /* 0x3e0375d4 */
-pp1  = -3.2504209876e-01, /* 0xbea66beb */
-pp2  = -2.8481749818e-02, /* 0xbce9528f */
-pp3  = -5.7702702470e-03, /* 0xbbbd1489 */
-pp4  = -2.3763017452e-05, /* 0xb7c756b1 */
-qq1  =  3.9791721106e-01, /* 0x3ecbbbce */
-qq2  =  6.5022252500e-02, /* 0x3d852a63 */
-qq3  =  5.0813062117e-03, /* 0x3ba68116 */
-qq4  =  1.3249473704e-04, /* 0x390aee49 */
-qq5  = -3.9602282413e-06, /* 0xb684e21a */
-/*
- * Coefficients for approximation to  erf  in [0.84375,1.25]
- */
-pa0  = -2.3621185683e-03, /* 0xbb1acdc6 */
-pa1  =  4.1485610604e-01, /* 0x3ed46805 */
-pa2  = -3.7220788002e-01, /* 0xbebe9208 */
-pa3  =  3.1834661961e-01, /* 0x3ea2fe54 */
-pa4  = -1.1089469492e-01, /* 0xbde31cc2 */
-pa5  =  3.5478305072e-02, /* 0x3d1151b3 */
-pa6  = -2.1663755178e-03, /* 0xbb0df9c0 */
-qa1  =  1.0642088205e-01, /* 0x3dd9f331 */
-qa2  =  5.4039794207e-01, /* 0x3f0a5785 */
-qa3  =  7.1828655899e-02, /* 0x3d931ae7 */
-qa4  =  1.2617121637e-01, /* 0x3e013307 */
-qa5  =  1.3637083583e-02, /* 0x3c5f6e13 */
-qa6  =  1.1984500103e-02, /* 0x3c445aa3 */
-/*
- * Coefficients for approximation to  erfc in [1.25,1/0.35]
- */
-ra0  = -9.8649440333e-03, /* 0xbc21a093 */
-ra1  = -6.9385856390e-01, /* 0xbf31a0b7 */
-ra2  = -1.0558626175e+01, /* 0xc128f022 */
-ra3  = -6.2375331879e+01, /* 0xc2798057 */
-ra4  = -1.6239666748e+02, /* 0xc322658c */
-ra5  = -1.8460508728e+02, /* 0xc3389ae7 */
-ra6  = -8.1287437439e+01, /* 0xc2a2932b */
-ra7  = -9.8143291473e+00, /* 0xc11d077e */
-sa1  =  1.9651271820e+01, /* 0x419d35ce */
-sa2  =  1.3765776062e+02, /* 0x4309a863 */
-sa3  =  4.3456588745e+02, /* 0x43d9486f */
-sa4  =  6.4538726807e+02, /* 0x442158c9 */
-sa5  =  4.2900814819e+02, /* 0x43d6810b */
-sa6  =  1.0863500214e+02, /* 0x42d9451f */
-sa7  =  6.5702495575e+00, /* 0x40d23f7c */
-sa8  = -6.0424413532e-02, /* 0xbd777f97 */
-/*
- * Coefficients for approximation to  erfc in [1/.35,28]
- */
-rb0  = -9.8649431020e-03, /* 0xbc21a092 */
-rb1  = -7.9928326607e-01, /* 0xbf4c9dd4 */
-rb2  = -1.7757955551e+01, /* 0xc18e104b */
-rb3  = -1.6063638306e+02, /* 0xc320a2ea */
-rb4  = -6.3756646729e+02, /* 0xc41f6441 */
-rb5  = -1.0250950928e+03, /* 0xc480230b */
-rb6  = -4.8351919556e+02, /* 0xc3f1c275 */
-sb1  =  3.0338060379e+01, /* 0x41f2b459 */
-sb2  =  3.2579251099e+02, /* 0x43a2e571 */
-sb3  =  1.5367296143e+03, /* 0x44c01759 */
-sb4  =  3.1998581543e+03, /* 0x4547fdbb */
-sb5  =  2.5530502930e+03, /* 0x451f90ce */
-sb6  =  4.7452853394e+02, /* 0x43ed43a7 */
-sb7  = -2.2440952301e+01; /* 0xc1b38712 */
-
-static float erfc1(float x)
+/* Top 12 bits of a float.  */
+static inline uint32_t
+top12 (float x)
 {
-	float_t s,P,Q;
-
-	s = fabsf(x) - 1;
-	P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
-	Q = 1+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
-	return 1 - erx - P/Q;
+  return asuint (x) >> 20;
 }
 
-static float erfc2(uint32_t ix, float x)
+/* Efficient implementation of erff
+   using either a pure polynomial approximation or
+   the exponential of a polynomial.
+   Worst-case error is 1.09ulps at 0x1.c111acp-1.  */
+float
+erff (float x)
 {
-	float_t s,R,S;
-	float z;
+  float r, x2, u;
 
-	if (ix < 0x3fa00000)  /* |x| < 1.25 */
-		return erfc1(x);
+  /* Get top word.  */
+  uint32_t ix = asuint (x);
+  uint32_t sign = ix >> 31;
+  uint32_t ia12 = top12 (x) & 0x7ff;
 
-	x = fabsf(x);
-	s = 1/(x*x);
-	if (ix < 0x4036db6d) {   /* |x| < 1/0.35 */
-		R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
-		     ra5+s*(ra6+s*ra7))))));
-		S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
-		     sa5+s*(sa6+s*(sa7+s*sa8)))))));
-	} else {                 /* |x| >= 1/0.35 */
-		R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
-		     rb5+s*rb6)))));
-		S = 1.0f+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
-		     sb5+s*(sb6+s*sb7))))));
+  /* Limit of both intervals is 0.875 for performance reasons but coefficients
+     computed on [0.0, 0.921875] and [0.921875, 4.0], which brought accuracy
+     from 0.94 to 1.1ulps.  */
+  if (ia12 < 0x3f6)
+    { /* a = |x| < 0.875.  */
+
+      /* Tiny and subnormal cases.  */
+      if (unlikely (ia12 < 0x318))
+	{ /* |x| < 2^(-28).  */
+	  if (unlikely (ia12 < 0x040))
+	    { /* |x| < 2^(-119).  */
+	      float y = fmaf (TwoOverSqrtPiMinusOne, x, x);
+	      return check_uflowf (y);
+	    }
+	  return x + TwoOverSqrtPiMinusOne * x;
 	}
-	ix = asuint(x);
-	z = asfloat(ix&0xffffe000);
-	return expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S)/x;
-}
 
-float erff(float x)
-{
-	float r,s,z,y;
-	uint32_t ix;
-	int sign;
+      x2 = x * x;
 
-	ix = asuint(x);
-	sign = ix>>31;
-	ix &= 0x7fffffff;
-	if (ix >= 0x7f800000) {
-		/* erf(nan)=nan, erf(+-inf)=+-1 */
-		return 1-2*sign + 1/x;
-	}
-	if (ix < 0x3f580000) {  /* |x| < 0.84375 */
-		if (ix < 0x31800000) {  /* |x| < 2**-28 */
-			/*avoid underflow */
-			return 0.125f*(8*x + efx8*x);
-		}
-		z = x*x;
-		r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
-		s = 1+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
-		y = r/s;
-		return x + x*y;
-	}
-	if (ix < 0x40c00000)  /* |x| < 6 */
-		y = 1 - erfc2(ix,x);
-	else
-		y = 1 - 0x1p-120f;
-	return sign ? -y : y;
-}
+      /* Normalized cases (|x| < 0.921875). Use Horner scheme for x+x*P(x^2).  */
+      r = A[5];
+      r = fmaf (r, x2, A[4]);
+      r = fmaf (r, x2, A[3]);
+      r = fmaf (r, x2, A[2]);
+      r = fmaf (r, x2, A[1]);
+      r = fmaf (r, x2, A[0]);
+      r = fmaf (r, x, x);
+    }
+  else if (ia12 < 0x408)
+    { /* |x| < 4.0 - Use a custom Estrin scheme.  */
 
-float erfcf(float x)
-{
-	float r,s,z,y;
-	uint32_t ix;
-	int sign;
+      float a = fabsf (x);
+      /* Start with Estrin scheme on high order (small magnitude) coefficients.  */
+      r = fmaf (B[6], a, B[5]);
+      u = fmaf (B[4], a, B[3]);
+      x2 = x * x;
+      r = fmaf (r, x2, u);
+      /* Then switch to pure Horner scheme.  */
+      r = fmaf (r, a, B[2]);
+      r = fmaf (r, a, B[1]);
+      r = fmaf (r, a, B[0]);
+      r = fmaf (r, a, a);
+      /* Single precision exponential with ~0.5ulps,
+	 ensures erff has max. rel. error
+	 < 1ulp on [0.921875, 4.0],
+	 < 1.1ulps on [0.875, 4.0].  */
+      r = expf (-r);
+      /* Explicit copysign (calling copysignf increases latency).  */
+      if (sign)
+	r = -1.0f + r;
+      else
+	r = 1.0f - r;
+    }
+  else
+    { /* |x| >= 4.0.  */
 
-	ix = asuint(x);
-	sign = ix>>31;
-	ix &= 0x7fffffff;
-	if (ix >= 0x7f800000) {
-		/* erfc(nan)=nan, erfc(+-inf)=0,2 */
-		return 2*sign + 1/x;
-	}
+      /* Special cases : erff(nan)=nan, erff(+inf)=+1 and erff(-inf)=-1.  */
+      if (unlikely (ia12 >= 0x7f8))
+	return (1.f - (float) ((ix >> 31) << 1)) + 1.f / x;
 
-	if (ix < 0x3f580000) {  /* |x| < 0.84375 */
-		if (ix < 0x23800000)  /* |x| < 2**-56 */
-			return 1.0f - x;
-		z = x*x;
-		r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
-		s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
-		y = r/s;
-		if (sign || ix < 0x3e800000)  /* x < 1/4 */
-			return 1.0f - (x+x*y);
-		return 0.5f - (x - 0.5f + x*y);
-	}
-	if (ix < 0x41e00000) {  /* |x| < 28 */
-		return sign ? 2 - erfc2(ix,x) : erfc2(ix,x);
-	}
-	return sign ? 2 - 0x1p-120f : 0x1p-120f*0x1p-120f;
+      /* Explicit copysign (calling copysignf increases latency).  */
+      if (sign)
+	r = -1.0f;
+      else
+	r = 1.0f;
+    }
+  return r;
 }

libc/tinymath/erff_data.c

@@ -1,9 +1,9 @@
 /*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
-│ vi: set et ft=c ts=8 sts=2 sw=2 fenc=utf-8                               :vi │
+│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8                               :vi │
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │
@@ -25,12 +25,19 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/tinymath/log1pf_data.internal.h"
+#include "libc/tinymath/arm.internal.h"
 __static_yoink("arm_optimized_routines_notice");
 
-/* Polynomial coefficients generated using floating-point minimax
-   algorithm, see tools/log1pf.sollya for details.  */
-const struct log1pf_data __log1pf_data
-  = {.coeffs = {-0x1p-1f, 0x1.5555aap-2f, -0x1.000038p-2f, 0x1.99675cp-3f,
-		-0x1.54ef78p-3f, 0x1.28a1f4p-3f, -0x1.0da91p-3f, 0x1.abcb6p-4f,
-		-0x1.6f0d5ep-5f}};
+/* Minimax approximation of erff. */
+const struct erff_data __erff_data = {
+.erff_poly_A = {
+0x1.06eba6p-03f, -0x1.8126e0p-02f, 0x1.ce1a46p-04f,
+-0x1.b68bd2p-06f, 0x1.473f48p-08f, -0x1.3a1a82p-11f
+},
+.erff_poly_B = {
+0x1.079d0cp-3f, 0x1.450aa0p-1f, 0x1.b55cb0p-4f,
+-0x1.8d6300p-6f, 0x1.fd1336p-9f, -0x1.91d2ccp-12f,
+0x1.222900p-16f
+}
+};
+

libc/tinymath/exp.c

@@ -1,9 +1,9 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=8 fenc=utf-8                                   :vi │
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8                               :vi │
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
 │  Optimized Routines                                                          │
-│  Copyright (c) 1999-2022, Arm Limited.                                       │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │
@@ -25,19 +25,9 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/intrin/likely.h"
-#include "libc/math.h"
-#include "libc/tinymath/exp_data.internal.h"
-#include "libc/tinymath/internal.h"
+#include "libc/tinymath/arm.internal.h"
 __static_yoink("arm_optimized_routines_notice");
 
-/*
- * Double-precision e^x function.
- *
- * Copyright (c) 2018, Arm Limited.
- * SPDX-License-Identifier: MIT
- */
-
 #define N (1 << EXP_TABLE_BITS)
 #define InvLn2N __exp_data.invln2N
 #define NegLn2hiN __exp_data.negln2hiN
@@ -48,6 +38,7 @@ __static_yoink("arm_optimized_routines_notice");
 #define C3 __exp_data.poly[6 - EXP_POLY_ORDER]
 #define C4 __exp_data.poly[7 - EXP_POLY_ORDER]
 #define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
+#define C6 __exp_data.poly[9 - EXP_POLY_ORDER]
 
 /* Handle cases that may overflow or underflow when computing the result that
    is scale*(1+TMP) without intermediate rounding.  The bit representation of
@@ -56,114 +47,154 @@ __static_yoink("arm_optimized_routines_notice");
    a double.  (int32_t)KI is the k used in the argument reduction and exponent
    adjustment of scale, positive k here means the result may overflow and
    negative k means the result may underflow.  */
-static inline double specialcase(double_t tmp, uint64_t sbits, uint64_t ki)
+static inline double
+specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
 {
-	double_t scale, y;
+  double_t scale, y;
 
-	if ((ki & 0x80000000) == 0) {
-		/* k > 0, the exponent of scale might have overflowed by <= 460.  */
-		sbits -= 1009ull << 52;
-		scale = asdouble(sbits);
-		y = 0x1p1009 * (scale + scale * tmp);
-		return eval_as_double(y);
-	}
-	/* k < 0, need special care in the subnormal range.  */
-	sbits += 1022ull << 52;
-	scale = asdouble(sbits);
-	y = scale + scale * tmp;
-	if (y < 1.0) {
-		/* Round y to the right precision before scaling it into the subnormal
-		 range to avoid double rounding that can cause 0.5+E/2 ulp error where
-		 E is the worst-case ulp error outside the subnormal range.  So this
-		 is only useful if the goal is better than 1 ulp worst-case error.  */
-		double_t hi, lo;
-		lo = scale - y + scale * tmp;
-		hi = 1.0 + y;
-		lo = 1.0 - hi + y + lo;
-		y = eval_as_double(hi + lo) - 1.0;
-		/* Avoid -0.0 with downward rounding.  */
-		if (WANT_ROUNDING && y == 0.0)
-			y = 0.0;
-		/* The underflow exception needs to be signaled explicitly.  */
-		fp_force_eval(fp_barrier(0x1p-1022) * 0x1p-1022);
-	}
-	y = 0x1p-1022 * y;
-	return eval_as_double(y);
+  if ((ki & 0x80000000) == 0)
+    {
+      /* k > 0, the exponent of scale might have overflowed by <= 460.  */
+      sbits -= 1009ull << 52;
+      scale = asdouble (sbits);
+      y = 0x1p1009 * (scale + scale * tmp);
+      return check_oflow (eval_as_double (y));
+    }
+  /* k < 0, need special care in the subnormal range.  */
+  sbits += 1022ull << 52;
+  scale = asdouble (sbits);
+  y = scale + scale * tmp;
+  if (y < 1.0)
+    {
+      /* Round y to the right precision before scaling it into the subnormal
+	 range to avoid double rounding that can cause 0.5+E/2 ulp error where
+	 E is the worst-case ulp error outside the subnormal range.  So this
+	 is only useful if the goal is better than 1 ulp worst-case error.  */
+      double_t hi, lo;
+      lo = scale - y + scale * tmp;
+      hi = 1.0 + y;
+      lo = 1.0 - hi + y + lo;
+      y = eval_as_double (hi + lo) - 1.0;
+      /* Avoid -0.0 with downward rounding.  */
+      if (WANT_ROUNDING && y == 0.0)
+	y = 0.0;
+      /* The underflow exception needs to be signaled explicitly.  */
+      force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022);
+    }
+  y = 0x1p-1022 * y;
+  return check_uflow (eval_as_double (y));
 }
 
 /* Top 12 bits of a double (sign and exponent bits).  */
-static inline uint32_t top12(double x)
+static inline uint32_t
+top12 (double x)
 {
-	return asuint64(x) >> 52;
+  return asuint64 (x) >> 52;
+}
+
+/* Computes exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
+   If hastail is 0 then xtail is assumed to be 0 too.  */
+static inline double
+exp_inline (double x, double xtail, int hastail)
+{
+  uint32_t abstop;
+  uint64_t ki, idx, top, sbits;
+  /* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
+  double_t kd, z, r, r2, scale, tail, tmp;
+
+  abstop = top12 (x) & 0x7ff;
+  if (unlikely (abstop - top12 (0x1p-54) >= top12 (512.0) - top12 (0x1p-54)))
+    {
+      if (abstop - top12 (0x1p-54) >= 0x80000000)
+	/* Avoid spurious underflow for tiny x.  */
+	/* Note: 0 is common input.  */
+	return WANT_ROUNDING ? 1.0 + x : 1.0;
+      if (abstop >= top12 (1024.0))
+	{
+	  if (asuint64 (x) == asuint64 (-INFINITY))
+	    return 0.0;
+	  if (abstop >= top12 (INFINITY))
+	    return 1.0 + x;
+	  if (asuint64 (x) >> 63)
+	    return __math_uflow (0);
+	  else
+	    return __math_oflow (0);
+	}
+      /* Large x is special cased below.  */
+      abstop = 0;
+    }
+
+  /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)].  */
+  /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N].  */
+  z = InvLn2N * x;
+#if TOINT_INTRINSICS
+  kd = roundtoint (z);
+  ki = converttoint (z);
+#elif EXP_USE_TOINT_NARROW
+  /* z - kd is in [-0.5-2^-16, 0.5] in all rounding modes.  */
+  kd = eval_as_double (z + Shift);
+  ki = asuint64 (kd) >> 16;
+  kd = (double_t) (int32_t) ki;
+#else
+  /* z - kd is in [-1, 1] in non-nearest rounding modes.  */
+  kd = eval_as_double (z + Shift);
+  ki = asuint64 (kd);
+  kd -= Shift;
+#endif
+  r = x + kd * NegLn2hiN + kd * NegLn2loN;
+  /* The code assumes 2^-200 < |xtail| < 2^-8/N.  */
+  if (hastail)
+    r += xtail;
+  /* 2^(k/N) ~= scale * (1 + tail).  */
+  idx = 2 * (ki % N);
+  top = ki << (52 - EXP_TABLE_BITS);
+  tail = asdouble (T[idx]);
+  /* This is only a valid scale when -1023*N < k < 1024*N.  */
+  sbits = T[idx + 1] + top;
+  /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1).  */
+  /* Evaluation is optimized assuming superscalar pipelined execution.  */
+  r2 = r * r;
+  /* Without fma the worst case error is 0.25/N ulp larger.  */
+  /* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp.  */
+#if EXP_POLY_ORDER == 4
+  tmp = tail + r + r2 * C2 + r * r2 * (C3 + r * C4);
+#elif EXP_POLY_ORDER == 5
+  tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
+#elif EXP_POLY_ORDER == 6
+  tmp = tail + r + r2 * (0.5 + r * C3) + r2 * r2 * (C4 + r * C5 + r2 * C6);
+#endif
+  if (unlikely (abstop == 0))
+    return specialcase (tmp, sbits, ki);
+  scale = asdouble (sbits);
+  /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
+     is no spurious underflow here even without fma.  */
+  return eval_as_double (scale + scale * tmp);
 }
 
 /**
  * Returns 𝑒^x.
+ *
+ * @raise ERANGE on overflow or underflow
  */
-double exp(double x)
+double
+exp (double x)
 {
-	uint32_t abstop;
-	uint64_t ki, idx, top, sbits;
-	double_t kd, z, r, r2, scale, tail, tmp;
-
-	abstop = top12(x) & 0x7ff;
-	if (UNLIKELY(abstop - top12(0x1p-54) >= top12(512.0) - top12(0x1p-54))) {
-		if (abstop - top12(0x1p-54) >= 0x80000000)
-			/* Avoid spurious underflow for tiny x.  */
-			/* Note: 0 is common input.  */
-			return WANT_ROUNDING ? 1.0 + x : 1.0;
-		if (abstop >= top12(1024.0)) {
-			if (asuint64(x) == asuint64(-INFINITY))
-				return 0.0;
-			if (abstop >= top12(INFINITY))
-				return 1.0 + x;
-			if (asuint64(x) >> 63)
-				return __math_uflow(0);
-			else
-				return __math_oflow(0);
-		}
-		/* Large x is special cased below.  */
-		abstop = 0;
-	}
-
-	/* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)].  */
-	/* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N].  */
-	z = InvLn2N * x;
-#if TOINT_INTRINSICS
-	kd = roundtoint(z);
-	ki = converttoint(z);
-#elif EXP_USE_TOINT_NARROW
-	/* z - kd is in [-0.5-2^-16, 0.5] in all rounding modes.  */
-	kd = eval_as_double(z + Shift);
-	ki = asuint64(kd) >> 16;
-	kd = (double_t)(int32_t)ki;
-#else
-	/* z - kd is in [-1, 1] in non-nearest rounding modes.  */
-	kd = eval_as_double(z + Shift);
-	ki = asuint64(kd);
-	kd -= Shift;
-#endif
-	r = x + kd * NegLn2hiN + kd * NegLn2loN;
-	/* 2^(k/N) ~= scale * (1 + tail).  */
-	idx = 2 * (ki % N);
-	top = ki << (52 - EXP_TABLE_BITS);
-	tail = asdouble(T[idx]);
-	/* This is only a valid scale when -1023*N < k < 1024*N.  */
-	sbits = T[idx + 1] + top;
-	/* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1).  */
-	/* Evaluation is optimized assuming superscalar pipelined execution.  */
-	r2 = r * r;
-	/* Without fma the worst case error is 0.25/N ulp larger.  */
-	/* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp.  */
-	tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
-	if (UNLIKELY(abstop == 0))
-		return specialcase(tmp, sbits, ki);
-	scale = asdouble(sbits);
-	/* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
-	   is no spurious underflow here even without fma.  */
-	return eval_as_double(scale + scale * tmp);
+  return exp_inline (x, 0, 0);
 }
 
-#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
-__weak_reference(exp, expl);
+/* May be useful for implementing pow where more than double
+   precision input is needed.  */
+double
+__exp_dd (double x, double xtail)
+{
+  return exp_inline (x, xtail, 1);
+}
+
+#if USE_GLIBC_ABI
+strong_alias (exp, __exp_finite)
+hidden_alias (exp, __ieee754_exp)
+hidden_alias (__exp_dd, __exp1)
+# if LDBL_MANT_DIG == 53
+long double expl (long double x) { return exp (x); }
+# endif
 #endif

libc/tinymath/exp10.c

@@ -1,9 +1,9 @@
-/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:8;tab-width:8;coding:utf-8   -*-│
-│ vi: set noet ft=c ts=8 sw=8 fenc=utf-8                                   :vi │
+/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
+│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8                               :vi │
 ╚──────────────────────────────────────────────────────────────────────────────╝
 │                                                                              │
-│  Musl Libc                                                                   │
-│  Copyright © 2005-2014 Rich Felker, et al.                                   │
+│  Optimized Routines                                                          │
+│  Copyright (c) 2018-2024, Arm Limited.                                       │
 │                                                                              │
 │  Permission is hereby granted, free of charge, to any person obtaining       │
 │  a copy of this software and associated documentation files (the             │
@@ -25,33 +25,135 @@
 │  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                      │
 │                                                                              │
 ╚─────────────────────────────────────────────────────────────────────────────*/
-#include "libc/math.h"
-__static_yoink("musl_libc_notice");
+#include "libc/tinymath/arm.internal.h"
+__static_yoink("arm_optimized_routines_notice");
+
+#define N (1 << EXP_TABLE_BITS)
+#define IndexMask (N - 1)
+#define OFlowBound 0x1.34413509f79ffp8 /* log10(DBL_MAX).  */
+#define UFlowBound -0x1.5ep+8 /* -350.  */
+#define SmallTop 0x3c6 /* top12(0x1p-57).  */
+#define BigTop 0x407   /* top12(0x1p8).  */
+#define Thresh 0x41    /* BigTop - SmallTop.  */
+#define Shift __exp_data.shift
+#define C(i) __exp_data.exp10_poly[i]
+
+static double
+special_case (uint64_t sbits, double_t tmp, uint64_t ki)
+{
+  double_t scale, y;
+
+  if (ki - (1ull << 16) < 0x80000000)
+    {
+      /* The exponent of scale might have overflowed by 1.  */
+      sbits -= 1ull << 52;
+      scale = asdouble (sbits);
+      y = 2 * (scale + scale * tmp);
+      return check_oflow (eval_as_double (y));
+    }
+
+  /* n < 0, need special care in the subnormal range.  */
+  sbits += 1022ull << 52;
+  scale = asdouble (sbits);
+  y = scale + scale * tmp;
+
+  if (y < 1.0)
+    {
+      /* Round y to the right precision before scaling it into the subnormal
+	 range to avoid double rounding that can cause 0.5+E/2 ulp error where
+	 E is the worst-case ulp error outside the subnormal range.  So this
+	 is only useful if the goal is better than 1 ulp worst-case error.  */
+      double_t lo = scale - y + scale * tmp;
+      double_t hi = 1.0 + y;
+      lo = 1.0 - hi + y + lo;
+      y = eval_as_double (hi + lo) - 1.0;
+      /* Avoid -0.0 with downward rounding.  */
+      if (WANT_ROUNDING && y == 0.0)
+	y = 0.0;
+      /* The underflow exception needs to be signaled explicitly.  */
+      force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022);
+    }
+  y = 0x1p-1022 * y;
+
+  return check_uflow (y);
+}
 
 /**
  * Returns 10ˣ.
+ *
+ * The largest observed error is ~0.513 ULP.
  */
-double exp10(double x)
+double
+exp10 (double x)
 {
-	static const double p10[] = {
-		1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10,
-		1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
-		1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
-		1e10, 1e11, 1e12, 1e13, 1e14, 1e15
-	};
-	double n, y = modf(x, &n);
-	union {double f; uint64_t i;} u = {n};
-	/* fabs(n) < 16 without raising invalid on nan */
-	if ((u.i>>52 & 0x7ff) < 0x3ff+4) {
-		if (!y) return p10[(int)n+15];
-		y = exp2(3.32192809488736234787031942948939 * y);
-		return y * p10[(int)n+15];
-	}
-	return pow(10.0, x);
+  uint64_t ix = asuint64 (x);
+  uint32_t abstop = (ix >> 52) & 0x7ff;
+
+  if (unlikely (abstop - SmallTop >= Thresh))
+    {
+      if (abstop - SmallTop >= 0x80000000)
+	/* Avoid spurious underflow for tiny x.
+	   Note: 0 is common input.  */
+	return x + 1;
+      if (abstop == 0x7ff)
+	return ix == asuint64 (-INFINITY) ? 0.0 : x + 1.0;
+      if (x >= OFlowBound)
+	return __math_oflow (0);
+      if (x < UFlowBound)
+	return __math_uflow (0);
+
+      /* Large x is special-cased below.  */
+      abstop = 0;
+    }
+
+  /* Reduce x: z = x * N / log10(2), k = round(z).  */
+  double_t z = __exp_data.invlog10_2N * x;
+  double_t kd;
+  int64_t ki;
+#if TOINT_INTRINSICS
+  kd = roundtoint (z);
+  ki = converttoint (z);
+#else
+  kd = eval_as_double (z + Shift);
+  kd -= Shift;
+  ki = kd;
+#endif
+
+  /* r = x - k * log10(2), r in [-0.5, 0.5].  */
+  double_t r = x;
+  r = __exp_data.neglog10_2hiN * kd + r;
+  r = __exp_data.neglog10_2loN * kd + r;
+
+  /* exp10(x) = 2^(k/N) * 2^(r/N).
+     Approximate the two components separately.  */
+
+  /* s = 2^(k/N), using lookup table.  */
+  uint64_t e = ki << (52 - EXP_TABLE_BITS);
+  uint64_t i = (ki & IndexMask) * 2;
+  uint64_t u = __exp_data.tab[i + 1];
+  uint64_t sbits = u + e;
+
+  double_t tail = asdouble (__exp_data.tab[i]);
+
+  /* 2^(r/N) ~= 1 + r * Poly(r).  */
+  double_t r2 = r * r;
+  double_t p = C (0) + r * C (1);
+  double_t y = C (2) + r * C (3);
+  y = y + r2 * C (4);
+  y = p + r2 * y;
+  y = tail + y * r;
+
+  if (unlikely (abstop == 0))
+    return special_case (sbits, y, ki);
+
+  /* Assemble components:
+     y  = 2^(r/N) * 2^(k/N)
+       ~= (y + 1) * s.  */
+  double_t s = asdouble (sbits);
+  return eval_as_double (s * y + s);
 }
 
 __strong_reference(exp10, pow10);
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 __weak_reference(exp10, exp10l);
-__weak_reference(exp10, pow10l);
 #endif