Windows: Pull in some math and formatting functions from Musl

Source/Windows/Common/CRT/CMakeLists.txt

@@ -1 +1,2 @@
 target_sources(CommonWindows PRIVATE Alloc.cpp IO.cpp Math.cpp String.cpp Misc.cpp CRT.cpp)
+add_subdirectory(Musl)

Source/Windows/Common/CRT/Musl/CMakeLists.txt

@@ -0,0 +1 @@
+target_sources(CommonWindows PRIVATE exp2.c log2_data.c remainder.c strtoimax.c strtoull.c exp_data.c fmod.c log2.c remquo.c strtoll.c strtoumax.c __math_uflow.c __math_oflow.c __math_xflow.c __math_invalid.c __math_divzero.c)

Source/Windows/Common/CRT/Musl/__math_divzero.c

@@ -0,0 +1,8 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2020 Rich Felker, et al.
+
+#include "libm.h"
+
+double __math_divzero(uint32_t sign) {
+  return fp_barrier(sign ? -1.0 : 1.0) / 0.0;
+}

Source/Windows/Common/CRT/Musl/__math_invalid.c

@@ -0,0 +1,8 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2020 Rich Felker, et al.
+
+#include "libm.h"
+
+double __math_invalid(double x) {
+  return (x - x) / (x - x);
+}

Source/Windows/Common/CRT/Musl/__math_oflow.c

@@ -0,0 +1,8 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2020 Rich Felker, et al.
+
+#include "libm.h"
+
+double __math_oflow(uint32_t sign) {
+  return __math_xflow(sign, 0x1p769);
+}

Source/Windows/Common/CRT/Musl/__math_uflow.c

@@ -0,0 +1,8 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2020 Rich Felker, et al.
+
+#include "libm.h"
+
+double __math_uflow(uint32_t sign) {
+  return __math_xflow(sign, 0x1p-767);
+}

Source/Windows/Common/CRT/Musl/__math_xflow.c

@@ -0,0 +1,8 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2020 Rich Felker, et al.
+
+#include "libm.h"
+
+double __math_xflow(uint32_t sign, double y) {
+  return eval_as_double(fp_barrier(sign ? -y : y) * y);
+}

Source/Windows/Common/CRT/Musl/exp2.c

@@ -0,0 +1,121 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright (c) 2018, Arm Limited.
+
+#include <math.h>
+#include <stdint.h>
+#include "libm.h"
+#include "exp_data.h"
+
+#define N (1 << EXP_TABLE_BITS)
+#define Shift __exp_data.exp2_shift
+#define T __exp_data.tab
+#define C1 __exp_data.exp2_poly[0]
+#define C2 __exp_data.exp2_poly[1]
+#define C3 __exp_data.exp2_poly[2]
+#define C4 __exp_data.exp2_poly[3]
+#define C5 __exp_data.exp2_poly[4]
+
+/* Handle cases that may overflow or underflow when computing the result that
+   is scale*(1+TMP) without intermediate rounding.  The bit representation of
+   scale is in SBITS, however it has a computed exponent that may have
+   overflown into the sign bit so that needs to be adjusted before using it as
+   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) {
+  double_t scale, y;
+
+  if ((ki & 0x80000000) == 0) {
+    /* k > 0, the exponent of scale might have overflowed by 1.  */
+    sbits -= 1ull << 52;
+    scale = asdouble(sbits);
+    y = 2 * (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);
+}
+
+/* Top 12 bits of a double (sign and exponent bits).  */
+static inline uint32_t top12(double x) {
+  return asuint64(x) >> 52;
+}
+
+double exp2(double x) {
+  uint32_t abstop;
+  uint64_t ki, idx, top, sbits;
+  double_t kd, r, r2, scale, tail, tmp;
+
+  abstop = top12(x) & 0x7ff;
+  if (predict_false(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_oflow(0);
+      } else if (asuint64(x) >= asuint64(-1075.0)) {
+        return __math_uflow(0);
+      }
+    }
+    if (2 * asuint64(x) > 2 * asuint64(928.0)) {
+      /* Large x is special cased below.  */
+      abstop = 0;
+    }
+  }
+
+  /* exp2(x) = 2^(k/N) * 2^r, with 2^r in [2^(-1/2N),2^(1/2N)].  */
+  /* x = k/N + r, with int k and r in [-1/2N, 1/2N].  */
+  kd = eval_as_double(x + Shift);
+  ki = asuint64(kd); /* k.  */
+  kd -= Shift;       /* k/N for int k.  */
+  r = x - kd;
+  /* 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;
+  /* exp2(x) = 2^(k/N) * 2^r ~= scale + scale * (tail + 2^r - 1).  */
+  /* Evaluation is optimized assuming superscalar pipelined execution.  */
+  r2 = r * r;
+  /* Without fma the worst case error is 0.5/N ulp larger.  */
+  /* Worst case error is less than 0.5+0.86/N+(abs poly error * 2^53) ulp.  */
+  tmp = tail + r * C1 + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
+  if (predict_false(abstop == 0)) {
+    return specialcase(tmp, sbits, ki);
+  }
+  scale = asdouble(sbits);
+  /* Note: tmp == 0 or |tmp| > 2^-65 and scale > 2^-928, so there
+     is no spurious underflow here even without fma.  */
+  return eval_as_double(scale + scale * tmp);
+}

Source/Windows/Common/CRT/Musl/exp_data.c

@@ -0,0 +1,309 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright (c) 2018, Arm Limited.
+
+#include "exp_data.h"
+
+#define N (1 << EXP_TABLE_BITS)
+
+const struct exp_data __exp_data = {
+  // N/ln2
+  .invln2N = 0x1.71547652b82fep0 * N,
+  // -ln2/N
+  .negln2hiN = -0x1.62e42fefa0000p-8,
+  .negln2loN = -0x1.cf79abc9e3b3ap-47,
+// Used for rounding when !TOINT_INTRINSICS
+#if EXP_USE_TOINT_NARROW
+  .shift = 0x1800000000.8p0,
+#else
+  .shift = 0x1.8p52,
+#endif
+  // exp polynomial coefficients.
+  .poly =
+    {
+      // abs error: 1.555*2^-66
+      // ulp error: 0.509 (0.511 without fma)
+      // if |x| < ln2/256+eps
+      // abs error if |x| < ln2/256+0x1p-15: 1.09*2^-65
+      // abs error if |x| < ln2/128: 1.7145*2^-56
+      0x1.ffffffffffdbdp-2,
+      0x1.555555555543cp-3,
+      0x1.55555cf172b91p-5,
+      0x1.1111167a4d017p-7,
+    },
+  .exp2_shift = 0x1.8p52 / N,
+  // exp2 polynomial coefficients.
+  .exp2_poly =
+    {
+      // abs error: 1.2195*2^-65
+      // ulp error: 0.507 (0.511 without fma)
+      // if |x| < 1/256
+      // abs error if |x| < 1/128: 1.9941*2^-56
+      0x1.62e42fefa39efp-1,
+      0x1.ebfbdff82c424p-3,
+      0x1.c6b08d70cf4b5p-5,
+      0x1.3b2abd24650ccp-7,
+      0x1.5d7e09b4e3a84p-10,
+    },
+  // 2^(k/N) ~= H[k]*(1 + T[k]) for int k in [0,N)
+  // tab[2*k] = asuint64(T[k])
+  // tab[2*k+1] = asuint64(H[k]) - (k << 52)/N
+  .tab =
+    {
+      0x0,
+      0x3ff0000000000000,
+      0x3c9b3b4f1a88bf6e,
+      0x3feff63da9fb3335,
+      0xbc7160139cd8dc5d,
+      0x3fefec9a3e778061,
+      0xbc905e7a108766d1,
+      0x3fefe315e86e7f85,
+      0x3c8cd2523567f613,
+      0x3fefd9b0d3158574,
+      0xbc8bce8023f98efa,
+      0x3fefd06b29ddf6de,
+      0x3c60f74e61e6c861,
+      0x3fefc74518759bc8,
+      0x3c90a3e45b33d399,
+      0x3fefbe3ecac6f383,
+      0x3c979aa65d837b6d,
+      0x3fefb5586cf9890f,
+      0x3c8eb51a92fdeffc,
+      0x3fefac922b7247f7,
+      0x3c3ebe3d702f9cd1,
+      0x3fefa3ec32d3d1a2,
+      0xbc6a033489906e0b,
+      0x3fef9b66affed31b,
+      0xbc9556522a2fbd0e,
+      0x3fef9301d0125b51,
+      0xbc5080ef8c4eea55,
+      0x3fef8abdc06c31cc,
+      0xbc91c923b9d5f416,
+      0x3fef829aaea92de0,
+      0x3c80d3e3e95c55af,
+      0x3fef7a98c8a58e51,
+      0xbc801b15eaa59348,
+      0x3fef72b83c7d517b,
+      0xbc8f1ff055de323d,
+      0x3fef6af9388c8dea,
+      0x3c8b898c3f1353bf,
+      0x3fef635beb6fcb75,
+      0xbc96d99c7611eb26,
+      0x3fef5be084045cd4,
+      0x3c9aecf73e3a2f60,
+      0x3fef54873168b9aa,
+      0xbc8fe782cb86389d,
+      0x3fef4d5022fcd91d,
+      0x3c8a6f4144a6c38d,
+      0x3fef463b88628cd6,
+      0x3c807a05b0e4047d,
+      0x3fef3f49917ddc96,
+      0x3c968efde3a8a894,
+      0x3fef387a6e756238,
+      0x3c875e18f274487d,
+      0x3fef31ce4fb2a63f,
+      0x3c80472b981fe7f2,
+      0x3fef2b4565e27cdd,
+      0xbc96b87b3f71085e,
+      0x3fef24dfe1f56381,
+      0x3c82f7e16d09ab31,
+      0x3fef1e9df51fdee1,
+      0xbc3d219b1a6fbffa,
+      0x3fef187fd0dad990,
+      0x3c8b3782720c0ab4,
+      0x3fef1285a6e4030b,
+      0x3c6e149289cecb8f,
+      0x3fef0cafa93e2f56,
+      0x3c834d754db0abb6,
+      0x3fef06fe0a31b715,
+      0x3c864201e2ac744c,
+      0x3fef0170fc4cd831,
+      0x3c8fdd395dd3f84a,
+      0x3feefc08b26416ff,
+      0xbc86a3803b8e5b04,
+      0x3feef6c55f929ff1,
+      0xbc924aedcc4b5068,
+      0x3feef1a7373aa9cb,
+      0xbc9907f81b512d8e,
+      0x3feeecae6d05d866,
+      0xbc71d1e83e9436d2,
+      0x3feee7db34e59ff7,
+      0xbc991919b3ce1b15,
+      0x3feee32dc313a8e5,
+      0x3c859f48a72a4c6d,
+      0x3feedea64c123422,
+      0xbc9312607a28698a,
+      0x3feeda4504ac801c,
+      0xbc58a78f4817895b,
+      0x3feed60a21f72e2a,
+      0xbc7c2c9b67499a1b,
+      0x3feed1f5d950a897,
+      0x3c4363ed60c2ac11,
+      0x3feece086061892d,
+      0x3c9666093b0664ef,
+      0x3feeca41ed1d0057,
+      0x3c6ecce1daa10379,
+      0x3feec6a2b5c13cd0,
+      0x3c93ff8e3f0f1230,
+      0x3feec32af0d7d3de,
+      0x3c7690cebb7aafb0,
+      0x3feebfdad5362a27,
+      0x3c931dbdeb54e077,
+      0x3feebcb299fddd0d,
+      0xbc8f94340071a38e,
+      0x3feeb9b2769d2ca7,
+      0xbc87deccdc93a349,
+      0x3feeb6daa2cf6642,
+      0xbc78dec6bd0f385f,
+      0x3feeb42b569d4f82,
+      0xbc861246ec7b5cf6,
+      0x3feeb1a4ca5d920f,
+      0x3c93350518fdd78e,
+      0x3feeaf4736b527da,
+      0x3c7b98b72f8a9b05,
+      0x3feead12d497c7fd,
+      0x3c9063e1e21c5409,
+      0x3feeab07dd485429,
+      0x3c34c7855019c6ea,
+      0x3feea9268a5946b7,
+      0x3c9432e62b64c035,
+      0x3feea76f15ad2148,
+      0xbc8ce44a6199769f,
+      0x3feea5e1b976dc09,
+      0xbc8c33c53bef4da8,
+      0x3feea47eb03a5585,
+      0xbc845378892be9ae,
+      0x3feea34634ccc320,
+      0xbc93cedd78565858,
+      0x3feea23882552225,
+      0x3c5710aa807e1964,
+      0x3feea155d44ca973,
+      0xbc93b3efbf5e2228,
+      0x3feea09e667f3bcd,
+      0xbc6a12ad8734b982,
+      0x3feea012750bdabf,
+      0xbc6367efb86da9ee,
+      0x3fee9fb23c651a2f,
+      0xbc80dc3d54e08851,
+      0x3fee9f7df9519484,
+      0xbc781f647e5a3ecf,
+      0x3fee9f75e8ec5f74,
+      0xbc86ee4ac08b7db0,
+      0x3fee9f9a48a58174,
+      0xbc8619321e55e68a,
+      0x3fee9feb564267c9,
+      0x3c909ccb5e09d4d3,
+      0x3feea0694fde5d3f,
+      0xbc7b32dcb94da51d,
+      0x3feea11473eb0187,
+      0x3c94ecfd5467c06b,
+      0x3feea1ed0130c132,
+      0x3c65ebe1abd66c55,
+      0x3feea2f336cf4e62,
+      0xbc88a1c52fb3cf42,
+      0x3feea427543e1a12,
+      0xbc9369b6f13b3734,
+      0x3feea589994cce13,
+      0xbc805e843a19ff1e,
+      0x3feea71a4623c7ad,
+      0xbc94d450d872576e,
+      0x3feea8d99b4492ed,
+      0x3c90ad675b0e8a00,
+      0x3feeaac7d98a6699,
+      0x3c8db72fc1f0eab4,
+      0x3feeace5422aa0db,
+      0xbc65b6609cc5e7ff,
+      0x3feeaf3216b5448c,
+      0x3c7bf68359f35f44,
+      0x3feeb1ae99157736,
+      0xbc93091fa71e3d83,
+      0x3feeb45b0b91ffc6,
+      0xbc5da9b88b6c1e29,
+      0x3feeb737b0cdc5e5,
+      0xbc6c23f97c90b959,
+      0x3feeba44cbc8520f,
+      0xbc92434322f4f9aa,
+      0x3feebd829fde4e50,
+      0xbc85ca6cd7668e4b,
+      0x3feec0f170ca07ba,
+      0x3c71affc2b91ce27,
+      0x3feec49182a3f090,
+      0x3c6dd235e10a73bb,
+      0x3feec86319e32323,
+      0xbc87c50422622263,
+      0x3feecc667b5de565,
+      0x3c8b1c86e3e231d5,
+      0x3feed09bec4a2d33,
+      0xbc91bbd1d3bcbb15,
+      0x3feed503b23e255d,
+      0x3c90cc319cee31d2,
+      0x3feed99e1330b358,
+      0x3c8469846e735ab3,
+      0x3feede6b5579fdbf,
+      0xbc82dfcd978e9db4,
+      0x3feee36bbfd3f37a,
+      0x3c8c1a7792cb3387,
+      0x3feee89f995ad3ad,
+      0xbc907b8f4ad1d9fa,
+      0x3feeee07298db666,
+      0xbc55c3d956dcaeba,
+      0x3feef3a2b84f15fb,
+      0xbc90a40e3da6f640,
+      0x3feef9728de5593a,
+      0xbc68d6f438ad9334,
+      0x3feeff76f2fb5e47,
+      0xbc91eee26b588a35,
+      0x3fef05b030a1064a,
+      0x3c74ffd70a5fddcd,
+      0x3fef0c1e904bc1d2,
+      0xbc91bdfbfa9298ac,
+      0x3fef12c25bd71e09,
+      0x3c736eae30af0cb3,
+      0x3fef199bdd85529c,
+      0x3c8ee3325c9ffd94,
+      0x3fef20ab5fffd07a,
+      0x3c84e08fd10959ac,
+      0x3fef27f12e57d14b,
+      0x3c63cdaf384e1a67,
+      0x3fef2f6d9406e7b5,
+      0x3c676b2c6c921968,
+      0x3fef3720dcef9069,
+      0xbc808a1883ccb5d2,
+      0x3fef3f0b555dc3fa,
+      0xbc8fad5d3ffffa6f,
+      0x3fef472d4a07897c,
+      0xbc900dae3875a949,
+      0x3fef4f87080d89f2,
+      0x3c74a385a63d07a7,
+      0x3fef5818dcfba487,
+      0xbc82919e2040220f,
+      0x3fef60e316c98398,
+      0x3c8e5a50d5c192ac,
+      0x3fef69e603db3285,
+      0x3c843a59ac016b4b,
+      0x3fef7321f301b460,
+      0xbc82d52107b43e1f,
+      0x3fef7c97337b9b5f,
+      0xbc892ab93b470dc9,
+      0x3fef864614f5a129,
+      0x3c74b604603a88d3,
+      0x3fef902ee78b3ff6,
+      0x3c83c5ec519d7271,
+      0x3fef9a51fbc74c83,
+      0xbc8ff7128fd391f0,
+      0x3fefa4afa2a490da,
+      0xbc8dae98e223747d,
+      0x3fefaf482d8e67f1,
+      0x3c8ec3bc41aa2008,
+      0x3fefba1bee615a27,
+      0x3c842b94c3a9eb32,
+      0x3fefc52b376bba97,
+      0x3c8a64a931d185ee,
+      0x3fefd0765b6e4540,
+      0xbc8e37bae43be3ed,
+      0x3fefdbfdad9cbe14,
+      0x3c77893b4d91cd9d,
+      0x3fefe7c1819e90d8,
+      0x3c5305c14160cc89,
+      0x3feff3c22b8f71f1,
+    },
+};

Source/Windows/Common/CRT/Musl/exp_data.h

@@ -0,0 +1,24 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright (c) 2018, Arm Limited.
+
+#ifndef _EXP_DATA_H
+#define _EXP_DATA_H
+
+#include <stdint.h>
+
+#define EXP_TABLE_BITS 7
+#define EXP_POLY_ORDER 5
+#define EXP_USE_TOINT_NARROW 0
+#define EXP2_POLY_ORDER 5
+extern const struct exp_data {
+  double invln2N;
+  double shift;
+  double negln2hiN;
+  double negln2loN;
+  double poly[4]; /* Last four coefficients.  */
+  double exp2_shift;
+  double exp2_poly[EXP2_POLY_ORDER];
+  uint64_t tab[2 * (1 << EXP_TABLE_BITS)];
+} __exp_data;
+
+#endif

Source/Windows/Common/CRT/Musl/fmod.c

@@ -0,0 +1,80 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2020 Rich Felker, et al.
+
+#include <math.h>
+#include <stdint.h>
+
+double fmod(double x, double y) {
+  union {
+    double f;
+    uint64_t i;
+  } ux = {x}, uy = {y};
+  int ex = ux.i >> 52 & 0x7ff;
+  int ey = uy.i >> 52 & 0x7ff;
+  int sx = ux.i >> 63;
+  uint64_t i;
+
+  /* in the followings uxi should be ux.i, but then gcc wrongly adds */
+  /* float load/store to inner loops ruining performance and code size */
+  uint64_t uxi = ux.i;
+
+  if (uy.i << 1 == 0 || isnan(y) || ex == 0x7ff) {
+    return (x * y) / (x * y);
+  }
+  if (uxi << 1 <= uy.i << 1) {
+    if (uxi << 1 == uy.i << 1) {
+      return 0 * x;
+    }
+    return x;
+  }
+
+  /* normalize x and y */
+  if (!ex) {
+    for (i = uxi << 12; i >> 63 == 0; ex--, i <<= 1)
+      ;
+    uxi <<= -ex + 1;
+  } else {
+    uxi &= -1ULL >> 12;
+    uxi |= 1ULL << 52;
+  }
+  if (!ey) {
+    for (i = uy.i << 12; i >> 63 == 0; ey--, i <<= 1)
+      ;
+    uy.i <<= -ey + 1;
+  } else {
+    uy.i &= -1ULL >> 12;
+    uy.i |= 1ULL << 52;
+  }
+
+  /* x mod y */
+  for (; ex > ey; ex--) {
+    i = uxi - uy.i;
+    if (i >> 63 == 0) {
+      if (i == 0) {
+        return 0 * x;
+      }
+      uxi = i;
+    }
+    uxi <<= 1;
+  }
+  i = uxi - uy.i;
+  if (i >> 63 == 0) {
+    if (i == 0) {
+      return 0 * x;
+    }
+    uxi = i;
+  }
+  for (; uxi >> 52 == 0; uxi <<= 1, ex--)
+    ;
+
+  /* scale result */
+  if (ex > 0) {
+    uxi -= 1ULL << 52;
+    uxi |= (uint64_t)ex << 52;
+  } else {
+    uxi >>= -ex + 1;
+  }
+  uxi |= (uint64_t)sx << 63;
+  ux.i = uxi;
+  return ux.f;
+}

Source/Windows/Common/CRT/Musl/libm.h

@@ -0,0 +1,288 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2020 Rich Felker, et al.
+
+#ifndef _LIBM_H
+#define _LIBM_H
+
+#include <stdint.h>
+#include <float.h>
+#include <math.h>
+
+#define hidden
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN
+union ldshape {
+  long double f;
+  struct {
+    uint64_t m;
+    uint16_t se;
+  } i;
+};
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN
+/* This is the m68k variant of 80-bit long double, and this definition only works
+ * on archs where the alignment requirement of uint64_t is <= 4. */
+union ldshape {
+  long double f;
+  struct {
+    uint16_t se;
+    uint16_t pad;
+    uint64_t m;
+  } i;
+};
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN
+union ldshape {
+  long double f;
+  struct {
+    uint64_t lo;
+    uint32_t mid;
+    uint16_t top;
+    uint16_t se;
+  } i;
+  struct {
+    uint64_t lo;
+    uint64_t hi;
+  } i2;
+};
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __BIG_ENDIAN
+union ldshape {
+  long double f;
+  struct {
+    uint16_t se;
+    uint16_t top;
+    uint32_t mid;
+    uint64_t lo;
+  } i;
+  struct {
+    uint64_t hi;
+    uint64_t lo;
+  } i2;
+};
+#else
+#error Unsupported long double representation
+#endif
+
+/* Support non-nearest rounding mode.  */
+#define WANT_ROUNDING 1
+/* Support signaling NaNs.  */
+#define WANT_SNAN 0
+
+#if WANT_SNAN
+#error SNaN is unsupported
+#else
+#define issignalingf_inline(x) 0
+#define issignaling_inline(x) 0
+#endif
+
+#ifndef TOINT_INTRINSICS
+#define TOINT_INTRINSICS 0
+#endif
+
+#if TOINT_INTRINSICS
+/* 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 double_t roundtoint(double_t);
+
+/* 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 int32_t converttoint(double_t);
+#endif
+
+/* Helps static branch prediction so hot path can be better optimized.  */
+#ifdef __GNUC__
+#define predict_true(x) __builtin_expect(!!(x), 1)
+#define predict_false(x) __builtin_expect(x, 0)
+#else
+#define predict_true(x) (x)
+#define predict_false(x) (x)
+#endif
+
+/* Evaluate an expression as the specified type. With standard excess
+   precision handling a type cast or assignment is enough (with
+   -ffloat-store an assignment is required, in old compilers argument
+   passing and return statement may not drop excess precision).  */
+
+static inline float eval_as_float(float x) {
+  float y = x;
+  return y;
+}
+
+static inline double eval_as_double(double x) {
+  double y = x;
+  return y;
+}
+
+/* fp_barrier returns its input, but limits code transformations
+   as if it had a side-effect (e.g. observable io) and returned
+   an arbitrary value.  */
+
+#ifndef fp_barrierf
+#define fp_barrierf fp_barrierf
+static inline float fp_barrierf(float x) {
+  volatile float y = x;
+  return y;
+}
+#endif
+
+#ifndef fp_barrier
+#define fp_barrier fp_barrier
+static inline double fp_barrier(double x) {
+  volatile double y = x;
+  return y;
+}
+#endif
+
+#ifndef fp_barrierl
+#define fp_barrierl fp_barrierl
+static inline long double fp_barrierl(long double x) {
+  volatile long double y = x;
+  return y;
+}
+#endif
+
+/* fp_force_eval ensures that the input value is computed when that's
+   otherwise unused.  To prevent the constant folding of the input
+   expression, an additional fp_barrier may be needed or a compilation
+   mode that does so (e.g. -frounding-math in gcc). Then it can be
+   used to evaluate an expression for its fenv side-effects only.   */
+
+#ifndef fp_force_evalf
+#define fp_force_evalf fp_force_evalf
+static inline void fp_force_evalf(float x) {
+  volatile float y;
+  y = x;
+}
+#endif
+
+#ifndef fp_force_eval
+#define fp_force_eval fp_force_eval
+static inline void fp_force_eval(double x) {
+  volatile double y;
+  y = x;
+}
+#endif
+
+#ifndef fp_force_evall
+#define fp_force_evall fp_force_evall
+static inline void fp_force_evall(long double x) {
+  volatile long double y;
+  y = x;
+}
+#endif
+
+#define FORCE_EVAL(x)                         \
+  do {                                        \
+    if (sizeof(x) == sizeof(float)) {         \
+      fp_force_evalf(x);                      \
+    } else if (sizeof(x) == sizeof(double)) { \
+      fp_force_eval(x);                       \
+    } else {                                  \
+      fp_force_evall(x);                      \
+    }                                         \
+  } while (0)
+
+#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 EXTRACT_WORDS(hi, lo, d) \
+  do {                           \
+    uint64_t __u = asuint64(d);  \
+    (hi) = __u >> 32;            \
+    (lo) = (uint32_t)__u;        \
+  } 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 INSERT_WORDS(d, hi, lo)                              \
+  do {                                                       \
+    (d) = asdouble(((uint64_t)(hi) << 32) | (uint32_t)(lo)); \
+  } 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)
+
+#define GET_FLOAT_WORD(w, d) \
+  do {                       \
+    (w) = asuint(d);         \
+  } while (0)
+
+#define SET_FLOAT_WORD(d, w) \
+  do {                       \
+    (d) = asfloat(w);        \
+  } while (0)
+
+hidden int __rem_pio2_large(double*, double*, int, int, int);
+
+hidden int __rem_pio2(double, double*);
+hidden double __sin(double, double, int);
+hidden double __cos(double, double);
+hidden double __tan(double, double, int);
+hidden double __expo2(double, double);
+
+hidden int __rem_pio2f(float, double*);
+hidden float __sindf(double);
+hidden float __cosdf(double);
+hidden float __tandf(double, int);
+hidden float __expo2f(float, float);
+
+hidden int __rem_pio2l(long double, long double*);
+hidden long double __sinl(long double, long double, int);
+hidden long double __cosl(long double, long double);
+hidden long double __tanl(long double, long double, int);
+
+hidden long double __polevll(long double, const long double*, int);
+hidden long double __p1evll(long double, const long double*, int);
+
+extern int __signgam;
+hidden double __lgamma_r(double, int*);
+hidden float __lgammaf_r(float, int*);
+
+/* error handling functions */
+hidden float __math_xflowf(uint32_t, float);
+hidden float __math_uflowf(uint32_t);
+hidden float __math_oflowf(uint32_t);
+hidden float __math_divzerof(uint32_t);
+hidden float __math_invalidf(float);
+hidden double __math_xflow(uint32_t, double);
+hidden double __math_uflow(uint32_t);
+hidden double __math_oflow(uint32_t);
+hidden double __math_divzero(uint32_t);
+hidden double __math_invalid(double);
+#if LDBL_MANT_DIG != DBL_MANT_DIG
+hidden long double __math_invalidl(long double);
+#endif
+
+#endif

Source/Windows/Common/CRT/Musl/log2.c

@@ -0,0 +1,119 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright (c) 2018, Arm Limited.
+
+#include <math.h>
+#include <stdint.h>
+#include "libm.h"
+#include "log2_data.h"
+
+#define T __log2_data.tab
+#define T2 __log2_data.tab2
+#define B __log2_data.poly1
+#define A __log2_data.poly
+#define InvLn2hi __log2_data.invln2hi
+#define InvLn2lo __log2_data.invln2lo
+#define N (1 << LOG2_TABLE_BITS)
+#define OFF 0x3fe6000000000000
+
+/* Top 16 bits of a double.  */
+static inline uint32_t top16(double x) {
+  return asuint64(x) >> 48;
+}
+
+double log2(double x) {
+  double_t z, r, r2, r4, y, invc, logc, kd, hi, lo, t1, t2, t3, p;
+  uint64_t ix, iz, tmp;
+  uint32_t top;
+  int k, i;
+
+  ix = asuint64(x);
+  top = top16(x);
+#define LO asuint64(1.0 - 0x1.5b51p-5)
+#define HI asuint64(1.0 + 0x1.6ab2p-5)
+  if (predict_false(ix - LO < HI - LO)) {
+    /* Handle close to 1.0 inputs separately.  */
+    /* Fix sign of zero with downward rounding when x==1.  */
+    if (WANT_ROUNDING && predict_false(ix == asuint64(1.0))) {
+      return 0;
+    }
+    r = x - 1.0;
+#if __FP_FAST_FMA
+    hi = r * InvLn2hi;
+    lo = r * InvLn2lo + __builtin_fma(r, InvLn2hi, -hi);
+#else
+    double_t rhi, rlo;
+    rhi = asdouble(asuint64(r) & -1ULL << 32);
+    rlo = r - rhi;
+    hi = rhi * InvLn2hi;
+    lo = rlo * InvLn2hi + r * InvLn2lo;
+#endif
+    r2 = r * r; /* rounding error: 0x1p-62.  */
+    r4 = r2 * r2;
+    /* Worst-case error is less than 0.54 ULP (0.55 ULP without fma).  */
+    p = r2 * (B[0] + r * B[1]);
+    y = hi + p;
+    lo += hi - y + p;
+    lo += r4 * (B[2] + r * B[3] + r2 * (B[4] + r * B[5]) + r4 * (B[6] + r * B[7] + r2 * (B[8] + r * B[9])));
+    y += lo;
+    return eval_as_double(y);
+  }
+  if (predict_false(top - 0x0010 >= 0x7ff0 - 0x0010)) {
+    /* x < 0x1p-1022 or inf or nan.  */
+    if (ix * 2 == 0) {
+      return __math_divzero(1);
+    }
+    if (ix == asuint64(INFINITY)) { /* log(inf) == inf.  */
+      return x;
+    }
+    if ((top & 0x8000) || (top & 0x7ff0) == 0x7ff0) {
+      return __math_invalid(x);
+    }
+    /* x is subnormal, normalize it.  */
+    ix = asuint64(x * 0x1p52);
+    ix -= 52ULL << 52;
+  }
+
+  /* x = 2^k z; where z is in range [OFF,2*OFF) and exact.
+     The range is split into N subintervals.
+     The ith subinterval contains z and c is near its center.  */
+  tmp = ix - OFF;
+  i = (tmp >> (52 - LOG2_TABLE_BITS)) % N;
+  k = (int64_t)tmp >> 52; /* arithmetic shift */
+  iz = ix - (tmp & 0xfffULL << 52);
+  invc = T[i].invc;
+  logc = T[i].logc;
+  z = asdouble(iz);
+  kd = (double_t)k;
+
+  /* log2(x) = log2(z/c) + log2(c) + k.  */
+  /* r ~= z/c - 1, |r| < 1/(2*N).  */
+#if __FP_FAST_FMA
+  /* rounding error: 0x1p-55/N.  */
+  r = __builtin_fma(z, invc, -1.0);
+  t1 = r * InvLn2hi;
+  t2 = r * InvLn2lo + __builtin_fma(r, InvLn2hi, -t1);
+#else
+  double_t rhi, rlo;
+  /* rounding error: 0x1p-55/N + 0x1p-65.  */
+  r = (z - T2[i].chi - T2[i].clo) * invc;
+  rhi = asdouble(asuint64(r) & -1ULL << 32);
+  rlo = r - rhi;
+  t1 = rhi * InvLn2hi;
+  t2 = rlo * InvLn2hi + r * InvLn2lo;
+#endif
+
+  /* hi + lo = r/ln2 + log2(c) + k.  */
+  t3 = kd + logc;
+  hi = t3 + t1;
+  lo = t3 - hi + t1 + t2;
+
+  /* log2(r+1) = r/ln2 + r^2*poly(r).  */
+  /* Evaluation is optimized assuming superscalar pipelined execution.  */
+  r2 = r * r; /* rounding error: 0x1p-54/N^2.  */
+  r4 = r2 * r2;
+  /* Worst-case error if |y| > 0x1p-4: 0.547 ULP (0.550 ULP without fma).
+     ~ 0.5 + 2/N/ln2 + abs-poly-error*0x1p56 ULP (+ 0.003 ULP without fma).  */
+  p = A[0] + r * A[1] + r2 * (A[2] + r * A[3]) + r4 * (A[4] + r * A[5]);
+  y = lo + r2 * p + hi;
+  return eval_as_double(y);
+}

Source/Windows/Common/CRT/Musl/log2_data.c

@@ -0,0 +1,137 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright (c) 2018, Arm Limited.
+
+#include "log2_data.h"
+
+#define N (1 << LOG2_TABLE_BITS)
+
+const struct log2_data __log2_data = {
+  // First coefficient: 0x1.71547652b82fe1777d0ffda0d24p0
+  .invln2hi = 0x1.7154765200000p+0,
+  .invln2lo = 0x1.705fc2eefa200p-33,
+  .poly1 =
+    {
+      // relative error: 0x1.2fad8188p-63
+      // in -0x1.5b51p-5 0x1.6ab2p-5
+      -0x1.71547652b82fep-1,
+      0x1.ec709dc3a03f7p-2,
+      -0x1.71547652b7c3fp-2,
+      0x1.2776c50f05be4p-2,
+      -0x1.ec709dd768fe5p-3,
+      0x1.a61761ec4e736p-3,
+      -0x1.7153fbc64a79bp-3,
+      0x1.484d154f01b4ap-3,
+      -0x1.289e4a72c383cp-3,
+      0x1.0b32f285aee66p-3,
+    },
+  .poly =
+    {
+      // relative error: 0x1.a72c2bf8p-58
+      // abs error: 0x1.67a552c8p-66
+      // in -0x1.f45p-8 0x1.f45p-8
+      -0x1.71547652b8339p-1,
+      0x1.ec709dc3a04bep-2,
+      -0x1.7154764702ffbp-2,
+      0x1.2776c50034c48p-2,
+      -0x1.ec7b328ea92bcp-3,
+      0x1.a6225e117f92ep-3,
+    },
+  /* Algorithm:
+
+          x = 2^k z
+          log2(x) = k + log2(c) + log2(z/c)
+          log2(z/c) = poly(z/c - 1)
+
+  where z is in [1.6p-1; 1.6p0] which is split into N subintervals and z falls
+  into the ith one, then table entries are computed as
+
+          tab[i].invc = 1/c
+          tab[i].logc = (double)log2(c)
+          tab2[i].chi = (double)c
+          tab2[i].clo = (double)(c - (double)c)
+
+  where c is near the center of the subinterval and is chosen by trying +-2^29
+  floating point invc candidates around 1/center and selecting one for which
+
+          1) the rounding error in 0x1.8p10 + logc is 0,
+          2) the rounding error in z - chi - clo is < 0x1p-64 and
+          3) the rounding error in (double)log2(c) is minimized (< 0x1p-68).
+
+  Note: 1) ensures that k + logc can be computed without rounding error, 2)
+  ensures that z/c - 1 can be computed as (z - chi - clo)*invc with close to a
+  single rounding error when there is no fast fma for z*invc - 1, 3) ensures
+  that logc + poly(z/c - 1) has small error, however near x == 1 when
+  |log2(x)| < 0x1p-4, this is not enough so that is special cased.  */
+  .tab =
+    {
+      {0x1.724286bb1acf8p+0, -0x1.1095feecdb000p-1}, {0x1.6e1f766d2cca1p+0, -0x1.08494bd76d000p-1},
+      {0x1.6a13d0e30d48ap+0, -0x1.00143aee8f800p-1}, {0x1.661ec32d06c85p+0, -0x1.efec5360b4000p-2},
+      {0x1.623fa951198f8p+0, -0x1.dfdd91ab7e000p-2}, {0x1.5e75ba4cf026cp+0, -0x1.cffae0cc79000p-2},
+      {0x1.5ac055a214fb8p+0, -0x1.c043811fda000p-2}, {0x1.571ed0f166e1ep+0, -0x1.b0b67323ae000p-2},
+      {0x1.53909590bf835p+0, -0x1.a152f5a2db000p-2}, {0x1.5014fed61adddp+0, -0x1.9217f5af86000p-2},
+      {0x1.4cab88e487bd0p+0, -0x1.8304db0719000p-2}, {0x1.49539b4334feep+0, -0x1.74189f9a9e000p-2},
+      {0x1.460cbdfafd569p+0, -0x1.6552bb5199000p-2}, {0x1.42d664ee4b953p+0, -0x1.56b23a29b1000p-2},
+      {0x1.3fb01111dd8a6p+0, -0x1.483650f5fa000p-2}, {0x1.3c995b70c5836p+0, -0x1.39de937f6a000p-2},
+      {0x1.3991c4ab6fd4ap+0, -0x1.2baa1538d6000p-2}, {0x1.3698e0ce099b5p+0, -0x1.1d98340ca4000p-2},
+      {0x1.33ae48213e7b2p+0, -0x1.0fa853a40e000p-2}, {0x1.30d191985bdb1p+0, -0x1.01d9c32e73000p-2},
+      {0x1.2e025cab271d7p+0, -0x1.e857da2fa6000p-3}, {0x1.2b404cf13cd82p+0, -0x1.cd3c8633d8000p-3},
+      {0x1.288b02c7ccb50p+0, -0x1.b26034c14a000p-3}, {0x1.25e2263944de5p+0, -0x1.97c1c2f4fe000p-3},
+      {0x1.234563d8615b1p+0, -0x1.7d6023f800000p-3}, {0x1.20b46e33eaf38p+0, -0x1.633a71a05e000p-3},
+      {0x1.1e2eefdcda3ddp+0, -0x1.494f5e9570000p-3}, {0x1.1bb4a580b3930p+0, -0x1.2f9e424e0a000p-3},
+      {0x1.19453847f2200p+0, -0x1.162595afdc000p-3}, {0x1.16e06c0d5d73cp+0, -0x1.f9c9a75bd8000p-4},
+      {0x1.1485f47b7e4c2p+0, -0x1.c7b575bf9c000p-4}, {0x1.12358ad0085d1p+0, -0x1.960c60ff48000p-4},
+      {0x1.0fef00f532227p+0, -0x1.64ce247b60000p-4}, {0x1.0db2077d03a8fp+0, -0x1.33f78b2014000p-4},
+      {0x1.0b7e6d65980d9p+0, -0x1.0387d1a42c000p-4}, {0x1.0953efe7b408dp+0, -0x1.a6f9208b50000p-5},
+      {0x1.07325cac53b83p+0, -0x1.47a954f770000p-5}, {0x1.05197e40d1b5cp+0, -0x1.d23a8c50c0000p-6},
+      {0x1.03091c1208ea2p+0, -0x1.16a2629780000p-6}, {0x1.0101025b37e21p+0, -0x1.720f8d8e80000p-8},
+      {0x1.fc07ef9caa76bp-1, 0x1.6fe53b1500000p-7},  {0x1.f4465d3f6f184p-1, 0x1.11ccce10f8000p-5},
+      {0x1.ecc079f84107fp-1, 0x1.c4dfc8c8b8000p-5},  {0x1.e573a99975ae8p-1, 0x1.3aa321e574000p-4},
+      {0x1.de5d6f0bd3de6p-1, 0x1.918a0d08b8000p-4},  {0x1.d77b681ff38b3p-1, 0x1.e72e9da044000p-4},
+      {0x1.d0cb5724de943p-1, 0x1.1dcd2507f6000p-3},  {0x1.ca4b2dc0e7563p-1, 0x1.476ab03dea000p-3},
+      {0x1.c3f8ee8d6cb51p-1, 0x1.7074377e22000p-3},  {0x1.bdd2b4f020c4cp-1, 0x1.98ede8ba94000p-3},
+      {0x1.b7d6c006015cap-1, 0x1.c0db86ad2e000p-3},  {0x1.b20366e2e338fp-1, 0x1.e840aafcee000p-3},
+      {0x1.ac57026295039p-1, 0x1.0790ab4678000p-2},  {0x1.a6d01bc2731ddp-1, 0x1.1ac056801c000p-2},
+      {0x1.a16d3bc3ff18bp-1, 0x1.2db11d4fee000p-2},  {0x1.9c2d14967feadp-1, 0x1.406464ec58000p-2},
+      {0x1.970e4f47c9902p-1, 0x1.52dbe093af000p-2},  {0x1.920fb3982bcf2p-1, 0x1.651902050d000p-2},
+      {0x1.8d30187f759f1p-1, 0x1.771d2cdeaf000p-2},  {0x1.886e5ebb9f66dp-1, 0x1.88e9c857d9000p-2},
+      {0x1.83c97b658b994p-1, 0x1.9a80155e16000p-2},  {0x1.7f405ffc61022p-1, 0x1.abe186ed3d000p-2},
+      {0x1.7ad22181415cap-1, 0x1.bd0f2aea0e000p-2},  {0x1.767dcf99eff8cp-1, 0x1.ce0a43dbf4000p-2},
+    },
+#if !__FP_FAST_FMA
+  .tab2 =
+    {
+      {0x1.6200012b90a8ep-1, 0x1.904ab0644b605p-55},  {0x1.66000045734a6p-1, 0x1.1ff9bea62f7a9p-57},
+      {0x1.69fffc325f2c5p-1, 0x1.27ecfcb3c90bap-55},  {0x1.6e00038b95a04p-1, 0x1.8ff8856739326p-55},
+      {0x1.71fffe09994e3p-1, 0x1.afd40275f82b1p-55},  {0x1.7600015590e1p-1, -0x1.2fd75b4238341p-56},
+      {0x1.7a00012655bd5p-1, 0x1.808e67c242b76p-56},  {0x1.7e0003259e9a6p-1, -0x1.208e426f622b7p-57},
+      {0x1.81fffedb4b2d2p-1, -0x1.402461ea5c92fp-55}, {0x1.860002dfafcc3p-1, 0x1.df7f4a2f29a1fp-57},
+      {0x1.89ffff78c6b5p-1, -0x1.e0453094995fdp-55},  {0x1.8e00039671566p-1, -0x1.a04f3bec77b45p-55},
+      {0x1.91fffe2bf1745p-1, -0x1.7fa34400e203cp-56}, {0x1.95fffcc5c9fd1p-1, -0x1.6ff8005a0695dp-56},
+      {0x1.9a0003bba4767p-1, 0x1.0f8c4c4ec7e03p-56},  {0x1.9dfffe7b92da5p-1, 0x1.e7fd9478c4602p-55},
+      {0x1.a1fffd72efdafp-1, -0x1.a0c554dcdae7ep-57}, {0x1.a5fffde04ff95p-1, 0x1.67da98ce9b26bp-55},
+      {0x1.a9fffca5e8d2bp-1, -0x1.284c9b54c13dep-55}, {0x1.adfffddad03eap-1, 0x1.812c8ea602e3cp-58},
+      {0x1.b1ffff10d3d4dp-1, -0x1.efaddad27789cp-55}, {0x1.b5fffce21165ap-1, 0x1.3cb1719c61237p-58},
+      {0x1.b9fffd950e674p-1, 0x1.3f7d94194cep-56},    {0x1.be000139ca8afp-1, 0x1.50ac4215d9bcp-56},
+      {0x1.c20005b46df99p-1, 0x1.beea653e9c1c9p-57},  {0x1.c600040b9f7aep-1, -0x1.c079f274a70d6p-56},
+      {0x1.ca0006255fd8ap-1, -0x1.a0b4076e84c1fp-56}, {0x1.cdfffd94c095dp-1, 0x1.8f933f99ab5d7p-55},
+      {0x1.d1ffff975d6cfp-1, -0x1.82c08665fe1bep-58}, {0x1.d5fffa2561c93p-1, -0x1.b04289bd295f3p-56},
+      {0x1.d9fff9d228b0cp-1, 0x1.70251340fa236p-55},  {0x1.de00065bc7e16p-1, -0x1.5011e16a4d80cp-56},
+      {0x1.e200002f64791p-1, 0x1.9802f09ef62ep-55},   {0x1.e600057d7a6d8p-1, -0x1.e0b75580cf7fap-56},
+      {0x1.ea00027edc00cp-1, -0x1.c848309459811p-55}, {0x1.ee0006cf5cb7cp-1, -0x1.f8027951576f4p-55},
+      {0x1.f2000782b7dccp-1, -0x1.f81d97274538fp-55}, {0x1.f6000260c450ap-1, -0x1.071002727ffdcp-59},
+      {0x1.f9fffe88cd533p-1, -0x1.81bdce1fda8bp-58},  {0x1.fdfffd50f8689p-1, 0x1.7f91acb918e6ep-55},
+      {0x1.0200004292367p+0, 0x1.b7ff365324681p-54},  {0x1.05fffe3e3d668p+0, 0x1.6fa08ddae957bp-55},
+      {0x1.0a0000a85a757p+0, -0x1.7e2de80d3fb91p-58}, {0x1.0e0001a5f3fccp+0, -0x1.1823305c5f014p-54},
+      {0x1.11ffff8afbaf5p+0, -0x1.bfabb6680bac2p-55}, {0x1.15fffe54d91adp+0, -0x1.d7f121737e7efp-54},
+      {0x1.1a00011ac36e1p+0, 0x1.c000a0516f5ffp-54},  {0x1.1e00019c84248p+0, -0x1.082fbe4da5dap-54},
+      {0x1.220000ffe5e6ep+0, -0x1.8fdd04c9cfb43p-55}, {0x1.26000269fd891p+0, 0x1.cfe2a7994d182p-55},
+      {0x1.2a00029a6e6dap+0, -0x1.00273715e8bc5p-56}, {0x1.2dfffe0293e39p+0, 0x1.b7c39dab2a6f9p-54},
+      {0x1.31ffff7dcf082p+0, 0x1.df1336edc5254p-56},  {0x1.35ffff05a8b6p+0, -0x1.e03564ccd31ebp-54},
+      {0x1.3a0002e0eaeccp+0, 0x1.5f0e74bd3a477p-56},  {0x1.3e000043bb236p+0, 0x1.c7dcb149d8833p-54},
+      {0x1.4200002d187ffp+0, 0x1.e08afcf2d3d28p-56},  {0x1.460000d387cb1p+0, 0x1.20837856599a6p-55},
+      {0x1.4a00004569f89p+0, -0x1.9fa5c904fbcd2p-55}, {0x1.4e000043543f3p+0, -0x1.81125ed175329p-56},
+      {0x1.51fffcc027f0fp+0, 0x1.883d8847754dcp-54},  {0x1.55ffffd87b36fp+0, -0x1.709e731d02807p-55},
+      {0x1.59ffff21df7bap+0, 0x1.7f79f68727b02p-55},  {0x1.5dfffebfc3481p+0, -0x1.180902e30e93ep-54},
+    },
+#endif
+};

Source/Windows/Common/CRT/Musl/log2_data.h

@@ -0,0 +1,25 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright (c) 2018, Arm Limited.
+
+#ifndef _LOG2_DATA_H
+#define _LOG2_DATA_H
+
+#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 !__FP_FAST_FMA
+  struct {
+    double chi, clo;
+  } tab2[1 << LOG2_TABLE_BITS];
+#endif
+} __log2_data;
+
+#endif

Source/Windows/Common/CRT/Musl/remainder.c

@@ -0,0 +1,9 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2020 Rich Felker, et al.
+
+#include <math.h>
+
+double remainder(double x, double y) {
+  int q;
+  return remquo(x, y, &q);
+}

Source/Windows/Common/CRT/Musl/remquo.c

@@ -0,0 +1,95 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2020 Rich Felker, et al.
+
+#include <math.h>
+#include <stdint.h>
+
+double remquo(double x, double y, int* quo) {
+  union {
+    double f;
+    uint64_t i;
+  } ux = {x}, uy = {y};
+  int ex = ux.i >> 52 & 0x7ff;
+  int ey = uy.i >> 52 & 0x7ff;
+  int sx = ux.i >> 63;
+  int sy = uy.i >> 63;
+  uint32_t q;
+  uint64_t i;
+  uint64_t uxi = ux.i;
+
+  *quo = 0;
+  if (uy.i << 1 == 0 || isnan(y) || ex == 0x7ff) {
+    return (x * y) / (x * y);
+  }
+  if (ux.i << 1 == 0) {
+    return x;
+  }
+
+  /* normalize x and y */
+  if (!ex) {
+    for (i = uxi << 12; i >> 63 == 0; ex--, i <<= 1)
+      ;
+    uxi <<= -ex + 1;
+  } else {
+    uxi &= -1ULL >> 12;
+    uxi |= 1ULL << 52;
+  }
+  if (!ey) {
+    for (i = uy.i << 12; i >> 63 == 0; ey--, i <<= 1)
+      ;
+    uy.i <<= -ey + 1;
+  } else {
+    uy.i &= -1ULL >> 12;
+    uy.i |= 1ULL << 52;
+  }
+
+  q = 0;
+  if (ex < ey) {
+    if (ex + 1 == ey) {
+      goto end;
+    }
+    return x;
+  }
+
+  /* x mod y */
+  for (; ex > ey; ex--) {
+    i = uxi - uy.i;
+    if (i >> 63 == 0) {
+      uxi = i;
+      q++;
+    }
+    uxi <<= 1;
+    q <<= 1;
+  }
+  i = uxi - uy.i;
+  if (i >> 63 == 0) {
+    uxi = i;
+    q++;
+  }
+  if (uxi == 0) {
+    ex = -60;
+  } else {
+    for (; uxi >> 52 == 0; uxi <<= 1, ex--)
+      ;
+  }
+end:
+  /* scale result and decide between |x| and |x|-|y| */
+  if (ex > 0) {
+    uxi -= 1ULL << 52;
+    uxi |= (uint64_t)ex << 52;
+  } else {
+    uxi >>= -ex + 1;
+  }
+  ux.i = uxi;
+  x = ux.f;
+  if (sy) {
+    y = -y;
+  }
+  if (ex == ey || (ex + 1 == ey && (2 * x > y || (2 * x == y && q % 2)))) {
+    x -= y;
+    q++;
+  }
+  q &= 0x7fffffff;
+  *quo = sx ^ sy ? -(int)q : (int)q;
+  return sx ? -x : x;
+}

Source/Windows/Common/CRT/Musl/strtoimax.c

@@ -0,0 +1,33 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2011 Rich Felker, et al.
+// NOTE: From an older musl release that avoids stdio usage
+
+#include <inttypes.h>
+#include <errno.h>
+#include <ctype.h>
+
+intmax_t strtoimax(const char* s1, char** p, int base) {
+  const unsigned char* s = s1;
+  int sign = 0;
+  uintmax_t x;
+
+  /* Initial whitespace */
+  for (; isspace(*s); s++)
+    ;
+
+  /* Optional sign */
+  if (*s == '-') {
+    sign = *s++;
+  } else if (*s == '+') {
+    s++;
+  }
+
+  x = strtoumax(s, p, base);
+  if (x > INTMAX_MAX) {
+    if (!sign || -x != INTMAX_MIN) {
+      errno = ERANGE;
+    }
+    return sign ? INTMAX_MIN : INTMAX_MAX;
+  }
+  return sign ? -x : x;
+}

Source/Windows/Common/CRT/Musl/strtoll.c

@@ -0,0 +1,20 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2011 Rich Felker, et al.
+// NOTE: From an older musl release that avoids stdio usage
+
+#include <stdlib.h>
+#include <inttypes.h>
+#include <errno.h>
+#include <limits.h>
+
+long long strtoll(const char* s, char** p, int base) {
+  intmax_t x = strtoimax(s, p, base);
+  if (x > LLONG_MAX) {
+    errno = ERANGE;
+    return LLONG_MAX;
+  } else if (x < LLONG_MIN) {
+    errno = ERANGE;
+    return LLONG_MIN;
+  }
+  return x;
+}

Source/Windows/Common/CRT/Musl/strtoull.c

@@ -0,0 +1,17 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2011 Rich Felker, et al.
+// NOTE: From an older musl release that avoids stdio usage
+
+#include <stdlib.h>
+#include <inttypes.h>
+#include <errno.h>
+#include <limits.h>
+
+unsigned long long strtoull(const char* s, char** p, int base) {
+  uintmax_t x = strtoumax(s, p, base);
+  if (x > ULLONG_MAX) {
+    errno = ERANGE;
+    return ULLONG_MAX;
+  }
+  return x;
+}

Source/Windows/Common/CRT/Musl/strtoumax.c

@@ -0,0 +1,140 @@
+// SPDX-License-Identifier: MIT
+// SPDX-FileCopyrightText: Copyright © 2005-2011 Rich Felker, et al.
+// NOTE: From an older musl release that avoids stdio usage
+
+#include <inttypes.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <ctype.h>
+#include <stdio.h>
+
+/* Lookup table for digit values. -1==255>=36 -> invalid */
+static const unsigned char digits[] = {
+  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0,  1,  2,  3,  4,  5,  6,  7,  8,  9,  -1, -1, -1, -1, -1, -1,
+  -1, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, -1, -1, -1, -1, -1,
+  -1, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, -1, -1, -1, -1, -1,
+  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
+};
+
+uintmax_t strtoumax(const char* s1, char** p, int base) {
+  const unsigned char* s = s1;
+  size_t x1, z1;
+  uintmax_t x, z = 0;
+  int sign = 0;
+  int shift;
+
+  if (!p) {
+    p = (char**)&s1;
+  }
+
+  /* Initial whitespace */
+  for (; isspace(*s); s++)
+    ;
+
+  /* Optional sign */
+  if (*s == '-') {
+    sign = *s++;
+  } else if (*s == '+') {
+    s++;
+  }
+
+  /* Default base 8, 10, or 16 depending on prefix */
+  if (base == 0) {
+    if (s[0] == '0') {
+      if ((s[1] | 32) == 'x') {
+        base = 16;
+      } else {
+        base = 8;
+      }
+    } else {
+      base = 10;
+    }
+  }
+
+  if ((unsigned)base - 2 > 36 - 2 || digits[*s] >= base) {
+    *p = (char*)s1;
+    errno = EINVAL;
+    return 0;
+  }
+
+  /* Main loops. Only use big types if we have to. */
+  if (base == 10) {
+    for (x1 = 0; isdigit(*s) && x1 <= SIZE_MAX / 10 - 10; s++) {
+      x1 = 10 * x1 + *s - '0';
+    }
+    for (x = x1; isdigit(*s) && x <= UINTMAX_MAX / 10 - 10; s++) {
+      x = 10 * x + *s - '0';
+    }
+    if (isdigit(*s)) {
+      if (isdigit(s[1]) || 10 * x > UINTMAX_MAX - (*s - '0')) {
+        goto overflow;
+      }
+      x = 10 * x + *s - '0';
+    }
+  } else if (!(base & base / 2)) {
+    if (base == 16) {
+      if (s[0] == '0' && (s[1] | 32) == 'x' && digits[s[2]] < 16) {
+        s += 2;
+      }
+      shift = 4;
+      z1 = SIZE_MAX / 16;
+      z = UINTMAX_MAX / 16;
+    } else if (base == 8) {
+      shift = 3;
+      z1 = SIZE_MAX / 8;
+      z = UINTMAX_MAX / 8;
+    } else if (base == 2) {
+      shift = 1;
+      z1 = SIZE_MAX / 2;
+      z = UINTMAX_MAX / 2;
+    } else if (base == 4) {
+      shift = 2;
+      z1 = SIZE_MAX / 4;
+      z = UINTMAX_MAX / 4;
+    } else /* if (base == 32) */ {
+      shift = 5;
+      z1 = SIZE_MAX / 32;
+      z = UINTMAX_MAX / 32;
+    }
+    for (x1 = 0; digits[*s] < base && x1 <= z1; s++) {
+      x1 = (x1 << shift) + digits[*s];
+    }
+    for (x = x1; digits[*s] < base && x <= z; s++) {
+      x = (x << shift) + digits[*s];
+    }
+    if (digits[*s] < base) {
+      goto overflow;
+    }
+  } else {
+    z1 = SIZE_MAX / base - base;
+    for (x1 = 0; digits[*s] < base && x1 <= z1; s++) {
+      x1 = x1 * base + digits[*s];
+    }
+    if (digits[*s] < base) {
+      z = UINTMAX_MAX / base - base;
+    }
+    for (x = x1; digits[*s] < base && x <= z; s++) {
+      x = x * base + digits[*s];
+    }
+    if (digits[*s] < base) {
+      if (digits[s[1]] < base || x * base > UINTMAX_MAX - digits[*s]) {
+        goto overflow;
+      }
+      x = x * base + digits[*s];
+    }
+  }
+
+  *p = (char*)s;
+  return sign ? -x : x;
+
+overflow:
+  for (; digits[*s] < base; s++)
+    ;
+  *p = (char*)s;
+  errno = ERANGE;
+  return UINTMAX_MAX;
+}