Compute Math functions as V8 does and link libm only when a program needs it

packages/compiler/src/backend/link-plan-core.test.ts

@@ -0,0 +1,20 @@
+import { expect, test } from "vitest";
+import { systemLibraryLinkerArgs } from "./link-plan-core.js";
+import { MACOS_ARM64_TARGET, LINUX_ARM64_GNU_TARGET, LINUX_X64_GNU_TARGET } from "./targets.js";
+
+test("Linux links libm as needed and other system libraries as given", () => {
+  for (const target of [LINUX_X64_GNU_TARGET, LINUX_ARM64_GNU_TARGET]) {
+    expect(systemLibraryLinkerArgs(target, ["m"])).toEqual([
+      "-Wl,--as-needed",
+      "-lm",
+      "-Wl,--no-as-needed",
+    ]);
+    expect(systemLibraryLinkerArgs(target, ["z", "m"])).toEqual([
+      "-lz",
+      "-Wl,--as-needed",
+      "-lm",
+      "-Wl,--no-as-needed",
+    ]);
+  }
+  expect(systemLibraryLinkerArgs(MACOS_ARM64_TARGET, ["System"])).toEqual(["-lSystem"]);
+});

packages/compiler/src/backend/link-plan-core.ts

@@ -70,3 +70,19 @@ export function executableLinkInputs(options: {
     ],
   };
 }
+
+/** The linker arguments for an executable's system libraries. On Linux, libm
+ * links as needed: the runtime defines the libm functions most programs
+ * reach (scr_glibc_math.c) and Math's transcendentals (scr_ieee754.c), so a
+ * program that calls no other libm function, such as pow, starts without
+ * loading libm.so. */
+export function systemLibraryLinkerArgs(
+  target: Pick<NativeTargetSpec, "platform">,
+  names: readonly string[],
+): string[] {
+  return names.flatMap((name) =>
+    name === "m" && target.platform === "linux"
+      ? ["-Wl,--as-needed", "-lm", "-Wl,--no-as-needed"]
+      : [`-l${name}`],
+  );
+}

packages/compiler/src/backend/linker.ts

@@ -22,6 +22,7 @@ import { type NativeArtifactDependency } from "./native/contracts.js";
 import { RuntimePackError, stageRuntimePackArtifacts } from "./runtime-pack.js";
 import type { NativeLinkPlan } from "./link-plan.js";
 import { lldLinkerArgs } from "./lld.js";
+import { systemLibraryLinkerArgs } from "./link-plan-core.js";
 import {
   createDarwinDebugSymbols,
   installDarwinDebugSymbols,
@@ -248,7 +249,7 @@ export async function linkNativeExecutable(
       ...lld,
       ...plan.driverFlags,
       ...plan.inputs.map((input) => staged.replacements.get(input) ?? input),
-      ...plan.systemLibraries.map((name) => `-l${name}`),
+      ...systemLibraryLinkerArgs(plan.target, plan.systemLibraries),
       "-o",
       privateOut,
     ];

packages/compiler/src/backend/llvm/expr-primitives.ts

@@ -47,7 +47,7 @@ import {
   type ArrayAccess,
 } from "./dense-array-access.js";
 import { emitWrappingSum, integerOrNaNRange, markWrappingOperands } from "./wrapping-sums.js";
-import { emitPowerOfTwoQuotient, isPowerOfTwoCall } from "./lib-filesystem.js";
+import { emitPowF64, emitPowerOfTwoQuotient, isPowerOfTwoCall } from "./lib-filesystem.js";
 
 /** The signed integer predicates of the number compares. */
 const INTEGER_COMPARES: Record<string, string> = {
@@ -305,7 +305,6 @@ export function emitOperatorExpr(
         "===": "oeq",
         "!==": "une",
       };
-      const libm: Record<string, string> = { "**": "pow" };
       const bit: Record<string, "and" | "or" | "xor" | "shl" | "ashr" | "lshr"> = {
         "&": "and",
         "|": "or",
@@ -471,32 +470,8 @@ export function emitOperatorExpr(
         B.line(`${t} = ${e.op === ">>>" ? "uitofp" : "sitofp"} i32 ${result} to double`);
         return { name: t, type: e.type, uint32: result };
       } else {
-        const fn = libm[e.op];
+        if (e.op !== "**") throw new LlvmUnsupportedError(`bin:${e.op}`, e.loc);
-        if (fn === undefined) throw new LlvmUnsupportedError(`bin:${e.op}`, e.loc);
+        return { name: emitPowF64(host, l.name, r.name, e.right), type: e.type };
-        host.declare(`declare double @${fn}(double, double)`);
-        // libm's pow answers 1 for pow(±1, ±Infinity) and pow(1, NaN), where
-        // `**` answers NaN (and LLVM folds pow(1.0, x) to 1.0): an inline
-        // guard selects NaN for a base of ±1 with a non-finite exponent and
-        // keeps the direct pow call. A finite literal exponent, or a literal
-        // base other than ±1, cannot meet those cases and needs no guard.
-        if (fn === "pow" && !libmPowAgrees(e.left, e.right)) {
-          host.declare(`declare double @llvm.fabs.f64(double)`);
-          const power = B.tmp();
-          B.line(`${power} = call double @pow(double ${l.name}, double ${r.name})`);
-          const base = B.tmp();
-          const unit = B.tmp();
-          const magnitude = B.tmp();
-          const finite = B.tmp();
-          const reject = B.tmp();
-          B.line(`${base} = call double @llvm.fabs.f64(double ${l.name})`);
-          B.line(`${unit} = fcmp oeq double ${base}, 1.0`);
-          B.line(`${magnitude} = call double @llvm.fabs.f64(double ${r.name})`);
-          B.line(`${finite} = fcmp olt double ${magnitude}, 0x7FF0000000000000`);
-          B.line(`${reject} = select i1 ${finite}, i1 false, i1 ${unit}`);
-          B.line(`${t} = select i1 ${reject}, double 0x7FF8000000000000, double ${power}`);
-        } else {
-          B.line(`${t} = call double @${fn}(double ${l.name}, double ${r.name})`);
-        }
       }
       return { name: t, type: e.type };
     }
@@ -1473,19 +1448,3 @@ export function emitRecordLitFields(
     host.initField(ptr, type, inline, v);
   }
 }
-
-/** Whether libm pow matches `**` for these operands: every case except a
- * base of ±1 with a NaN or infinite exponent agrees. */
-function libmPowAgrees(base: IrExpr, exponent: IrExpr): boolean {
-  const power = numericLiteral(exponent);
-  if (power !== null && Number.isFinite(power)) return true;
-  const value = numericLiteral(base);
-  return value !== null && !Number.isNaN(value) && Math.abs(value) !== 1;
-}
-
-function numericLiteral(expr: IrExpr): number | null {
-  if (expr.kind === "numLit") return expr.value;
-  if (expr.kind === "unary" && expr.op === "-" && expr.operand.kind === "numLit")
-    return -expr.operand.value;
-  return null;
-}

packages/compiler/src/backend/llvm/lib-dispatch.ts

@@ -17,6 +17,7 @@ import type { LlvmEmitterContext, LibCallExpr, LibCallPrefix, LlValue } from "./
 import {
   ASSERT_VERDICT_LIB_FNS,
   LIB_FN_SYMS,
+  PURE_MATH_LIB_FNS,
   mockTimerSym,
   USES_TIMERS_LIB_FNS,
 } from "./lib-shared.js";
@@ -777,7 +778,8 @@ function emitSymbolLibCall(
     return { name: "", type: e.type };
   }
   const t = B.tmp();
-  B.line(`${t} = call ${retDecl} @${sym}(${argList})`);
+  const attrs = PURE_MATH_LIB_FNS.has(e.fn) ? " nounwind willreturn memory(none)" : "";
+  B.line(`${t} = call ${retDecl} @${sym}(${argList})${attrs}`);
   // The result joins its frame BEFORE the pending check so an unwind
   // releases the dummy (NULL for refcounted returns) harmlessly.
   const out = host.own({ name: t, type: e.type });

packages/compiler/src/backend/llvm/lib-filesystem.ts

@@ -1027,6 +1027,54 @@ function emitPowerOfTwo(host: LlvmEmitterContext, n: string): string {
   return out;
 }
 
+/** `base ** exponent` and `Math.pow(base, exponent)` as V8's math::pow
+ * answers them: a literal exponent 2 is `base * base`, a literal 0.5 is
+ * `sqrt(base + 0)` (+Infinity for an infinite base), and every other
+ * exponent calls scr_math_pow, which applies V8's other special cases
+ * before the platform pow. The platform pow is never called by its libm
+ * name, which LLVM would rewrite (pow(2, x) to exp2(x), pow(x, -1) to 1/x)
+ * into answers that can differ from V8's. */
+export function emitPowF64(
+  host: LlvmEmitterContext,
+  base: string,
+  exponent: string,
+  exponentExpr: IrExpr,
+): string {
+  const B = host.B;
+  const literal =
+    exponentExpr.kind === "numLit"
+      ? exponentExpr.value
+      : exponentExpr.kind === "unary" &&
+          exponentExpr.op === "-" &&
+          exponentExpr.operand.kind === "numLit"
+        ? -exponentExpr.operand.value
+        : null;
+  const out = B.tmp();
+  if (literal === 2) {
+    B.line(`${out} = fmul double ${base}, ${base}`);
+    return out;
+  }
+  if (literal === 0.5) {
+    host.declare(`declare double @llvm.sqrt.f64(double)`);
+    host.declare(`declare double @llvm.fabs.f64(double)`);
+    const plus = B.tmp(),
+      root = B.tmp(),
+      magnitude = B.tmp(),
+      infinite = B.tmp();
+    B.line(`${plus} = fadd double ${base}, ${f64Lit(0)}`);
+    B.line(`${root} = call double @llvm.sqrt.f64(double ${plus})`);
+    B.line(`${magnitude} = call double @llvm.fabs.f64(double ${base})`);
+    B.line(`${infinite} = fcmp oeq double ${magnitude}, 0x7FF0000000000000`);
+    B.line(`${out} = select i1 ${infinite}, double 0x7FF0000000000000, double ${root}`);
+    return out;
+  }
+  host.declare(`declare double @scr_math_pow(double, double)`);
+  B.line(
+    `${out} = call double @scr_math_pow(double ${base}, double ${exponent}) nounwind willreturn`,
+  );
+  return out;
+}
+
 /** `Math.pow(2, n)`, the divisor emitPowerOfTwoQuotient takes. */
 export function isPowerOfTwoCall(e: IrExpr): boolean {
   return (
@@ -1099,6 +1147,11 @@ export function emitPrimitiveLibCall(host: LlvmEmitterContext, e: LibCallExpr):
     const v = host.emitExpr(e.args[1]!);
     return { name: emitPowerOfTwo(host, v.name), type: e.type };
   }
+  if (e.fn === "math.pow") {
+    const base = host.emitExpr(e.args[0]!);
+    const exponent = host.emitExpr(e.args[1]!);
+    return { name: emitPowF64(host, base.name, exponent.name, e.args[1]!), type: e.type };
+  }
   if (e.fn === "math.floor" || e.fn === "math.ceil") {
     const intr = e.fn === "math.floor" ? "floor" : "ceil";
     const v = host.emitExpr(e.args[0]!);

packages/compiler/src/backend/llvm/lib-shared.ts

@@ -35,6 +35,31 @@ export const ASSERT_VERDICT_LIB_FNS: ReadonlyMap<string, number | null> = new Ma
   ["assert.partialDeepResult", null],
 ]);
 
+/** The Math functions that call scr_ieee754.c. They read and write no
+ * memory, so calls with the same argument fold and loops may hoist them. */
+export const PURE_MATH_LIB_FNS: ReadonlySet<string> = new Set([
+  "math.sin",
+  "math.sinh",
+  "math.cos",
+  "math.cosh",
+  "math.tan",
+  "math.tanh",
+  "math.asin",
+  "math.asinh",
+  "math.acos",
+  "math.acosh",
+  "math.atan",
+  "math.atanh",
+  "math.cbrt",
+  "math.exp",
+  "math.expm1",
+  "math.log",
+  "math.log1p",
+  "math.log2",
+  "math.log10",
+  "math.atan2",
+]);
+
 export const LIB_FN_SYMS: Record<string, string> = {
   "assert.expectsErrDyn": "scr_assert_expects_err_dyn",
   "assert.throwsValidate": "scr_assert_throws_validate",
@@ -118,27 +143,29 @@ export const LIB_FN_SYMS: Record<string, string> = {
   "str.codePoints": "scr_str_code_points",
   "num.fromString": "scr_string_to_number",
   "math.round": "scr_math_round",
-  "math.sin": "sin",
+  // Math's transcendental functions are V8's (scr_ieee754.c), so their
-  "math.sinh": "sinh",
+  // results are the ones Node prints.
-  "math.cos": "cos",
+  "math.sin": "scr_ieee754_sin",
-  "math.cosh": "cosh",
+  "math.sinh": "scr_ieee754_sinh",
-  "math.tan": "tan",
+  "math.cos": "scr_ieee754_cos",
-  "math.tanh": "tanh",
+  "math.cosh": "scr_ieee754_cosh",
-  "math.asin": "asin",
+  "math.tan": "scr_ieee754_tan",
-  "math.asinh": "asinh",
+  "math.tanh": "scr_ieee754_tanh",
-  "math.acos": "acos",
+  "math.asin": "scr_ieee754_asin",
-  "math.acosh": "acosh",
+  "math.asinh": "scr_ieee754_asinh",
-  "math.atan": "atan",
+  "math.acos": "scr_ieee754_acos",
-  "math.atanh": "atanh",
+  "math.acosh": "scr_ieee754_acosh",
-  "math.cbrt": "cbrt",
+  "math.atan": "scr_ieee754_atan",
-  "math.exp": "exp",
+  "math.atanh": "scr_ieee754_atanh",
-  "math.expm1": "expm1",
+  "math.cbrt": "scr_ieee754_cbrt",
+  "math.exp": "scr_ieee754_exp",
+  "math.expm1": "scr_ieee754_expm1",
   "math.sqrt": "sqrt",
-  "math.log": "log",
+  "math.log": "scr_ieee754_log",
-  "math.log1p": "log1p",
+  "math.log1p": "scr_ieee754_log1p",
-  "math.log2": "log2",
+  "math.log2": "scr_ieee754_log2",
-  "math.log10": "log10",
+  "math.log10": "scr_ieee754_log10",
-  "math.atan2": "atan2",
+  "math.atan2": "scr_ieee754_atan2",
   "math.pow": "scr_math_pow",
   // encodeURIComponent throws URIError on a lone surrogate
   // (MAY_THROW_LIB_FNS); decodeUriComponent is with the throwing tier below.

packages/compiler/src/backend/native/runtime-inputs.ts

@@ -74,6 +74,8 @@ export const EXECUTABLE_RUNTIME_SOURCES = [
   "scr_slots.c",
   "scr_alloc.c",
   "scr_cipher_info.c",
+  "scr_ieee754.c",
+  "scr_glibc_math.c",
 ] as const;
 
 export function runtimeSrcDir(): string {

packages/compiler/src/native/compiler.ts

@@ -21,7 +21,7 @@ import type {
   CompileRequestResult,
 } from "../compile-types.js";
 import { LlvmUnsupportedError } from "../backend/llvm/emitter.js";
-import { executableLinkInputs } from "../backend/link-plan-core.js";
+import { executableLinkInputs, systemLibraryLinkerArgs } from "../backend/link-plan-core.js";
 import { formatNativeLinkInfo } from "../backend/native-link-info-core.js";
 import {
   emitNativeObject,
@@ -527,7 +527,7 @@ export class NativeCompiler {
           ...lld,
           ...plan.driverFlags,
           ...plan.inputs,
-          ...plan.systemLibraries.map((name) => `-l${name}`),
+          ...systemLibraryLinkerArgs(target, plan.systemLibraries),
           "-o",
           stagedOutput,
         ];

packages/compiler/test/llvm-remainder.test.ts

@@ -90,7 +90,7 @@ console.log(remainder(5, 2.5), signed(-7.5), maybeZero(7.5, 2.5));
     expect(maybeZero).toContain("rem.float");
     expect(ll).not.toContain("@llvm.fptosi.sat");
     expect(ll).not.toContain("@fmod");
-    expect(ll).toContain("call double @pow(");
+    expect(ll).toContain("call double @scr_math_pow(");
 
     // Where truncation is an integer round trip, a dividend below 2^62 takes
     // the i64 remainder by the power of two and adds back its fraction; a

packages/compiler/test/ts7/baselines/order-parity/tests/corpus/7717-math-fdlibm.ts.json

@@ -0,0 +1,6 @@
+{
+ "order": [
+  "<repo>/tests/corpus/7717-math-fdlibm.ts"
+ ],
+ "diags": []
+}

packages/compiler/test/ts7/baselines/order-parity/tests/corpus/7718-math-fdlibm.cjs.json

@@ -0,0 +1,6 @@
+{
+ "order": [
+  "<repo>/tests/corpus/7718-math-fdlibm.cjs"
+ ],
+ "diags": []
+}

packages/runtime-darwin-arm64/scripts/build.mjs

@@ -348,6 +348,11 @@ async function build() {
         "artifacts/licenses/v8-number-radix.txt",
         "BSD-3-Clause",
       ],
+      [
+        join(vendorRoot, "v8-ieee754", "LICENSE"),
+        "artifacts/licenses/v8-ieee754.txt",
+        "BSD-3-Clause",
+      ],
       [join(vendorRoot, "unicode", "LICENSE"), "artifacts/licenses/unicode.txt", "Unicode-3.0"],
       [join(zlib, "LICENSE"), "artifacts/licenses/zlib.txt", "Zlib"],
       [join(vendorRoot, "brotli", "LICENSE"), "artifacts/licenses/brotli.txt", "MIT"],

packages/runtime-pack-common/runtime-pack-matrix.mjs

@@ -52,6 +52,10 @@ export const BASE_RUNTIME_SOURCES = [
   "scr_alloc.c",
   // crypto.getCipherInfo's table; it links only into programs that call it.
   "scr_cipher_info.c",
+  // Math's transcendental functions as V8 computes them, and on glibc the
+  // exact libm functions that let a program link without libm.so.
+  "scr_ieee754.c",
+  "scr_glibc_math.c",
 ];
 
 export const OPTIONAL_RUNTIME_SOURCES = [

packages/runtime-pack-common/scripts/build.mjs

@@ -447,6 +447,11 @@ async function build() {
         "artifacts/licenses/v8-number-radix.txt",
         "BSD-3-Clause",
       ],
+      [
+        join(vendorRoot, "v8-ieee754", "LICENSE"),
+        "artifacts/licenses/v8-ieee754.txt",
+        "BSD-3-Clause",
+      ],
       [join(zlib, "LICENSE"), "artifacts/licenses/zlib.txt", "Zlib"],
       [join(vendorRoot, "brotli", "LICENSE"), "artifacts/licenses/brotli.txt", "MIT"],
       [join(mbedtls, "LICENSE"), "artifacts/licenses/mbedtls.txt", "Apache-2.0"],

packages/runtime/src/scr_bytes.c

@@ -4141,8 +4141,8 @@ static bool scr_bytes_rw_check(const ScrBytes *b, double offset, size_t width) {
  * exactly Node — and writes zeros downstream. */
 static bool scr_bytes_check_int(const ScrBytes *b, double value, double offset,
                                 size_t width, bool sign) {
-  double max = sign ? exp2((double)(8 * width - 1)) - 1 : exp2((double)(8 * width)) - 1;
+  double max = sign ? ldexp(1.0, (int)(8 * width - 1)) - 1 : ldexp(1.0, (int)(8 * width)) - 1;
-  double min = sign ? -exp2((double)(8 * width - 1)) : 0;
+  double min = sign ? -ldexp(1.0, (int)(8 * width - 1)) : 0;
   if (value > max || value < min) {
     char recv[48];
     scr_num_received(value, recv);

packages/runtime/src/scr_bytes_dyn.c

@@ -44,8 +44,8 @@ static size_t scr_bytes_dyn_width(ScrBytesNumKind kind) {
 
 /* Node's checkInt range test and message, before the offset is looked at. */
 static bool scr_bytes_dyn_value_in_range(double value, size_t width, bool sign) {
-  double max = sign ? exp2((double)(8 * width - 1)) - 1 : exp2((double)(8 * width)) - 1;
+  double max = sign ? ldexp(1.0, (int)(8 * width - 1)) - 1 : ldexp(1.0, (int)(8 * width)) - 1;
-  double min = sign ? -exp2((double)(8 * width - 1)) : 0;
+  double min = sign ? -ldexp(1.0, (int)(8 * width - 1)) : 0;
   if (!(value > max || value < min)) return true;
   char recv[48];
   scr_num_received(value, recv);

packages/runtime/src/scr_glibc_math.c

@@ -0,0 +1,271 @@
+/* The libm functions that glibc executables reach without calling a Math
+ * transcendental: floor, ceil, trunc and round (LLVM's rounding intrinsics
+ * on x86-64 before SSE4.1, and the runtime's own calls), fmod (the %
+ * operator), sqrt (Math.sqrt's call for a negative argument), and ldexp,
+ * scalbn and frexp, which glibc exports from both libc and libm (the link
+ * would otherwise take them from libm). Each is exactly specified, so these
+ * definitions answer what glibc's do, the NaN payloads of NaN inputs
+ * included. Defining them in the executable lets a program that uses no
+ * other libm function link without libm.so (-lm is linked as needed),
+ * which saves loading it at startup. Library builds keep the embedder's
+ * libm. */
+
+#include <stdint.h>
+#include <string.h>
+
+#if defined(__GLIBC__) && !defined(SCR_LIB)
+
+/* The sanitizer runtime calls some of these functions while it initializes,
+ * and the dynamic loader calls the resolvers before that: none of the code
+ * here carries sanitizer instrumentation. no_sanitize("thread") alone keeps
+ * ThreadSanitizer's function entry and exit calls, which crash before its
+ * runtime initializes; disable_sanitizer_instrumentation removes them. */
+#pragma clang attribute push(__attribute__((no_sanitize("address", "undefined", "thread"), disable_sanitizer_instrumentation)), apply_to = function)
+
+static inline uint64_t sgm_bits(double d) {
+  uint64_t bits;
+  memcpy(&bits, &d, sizeof bits);
+  return bits;
+}
+
+static inline double sgm_double(uint64_t bits) {
+  double d;
+  memcpy(&d, &bits, sizeof d);
+  return d;
+}
+
+#define SGM_SIGN (UINT64_C(1) << 63)
+
+/* Hidden: the executable's own references bind to these definitions, which
+ * are not exported (libc.so and libm.so define the same names), so section
+ * GC drops the ones a program does not call. */
+#define SGM_DEFINE __attribute__((visibility("hidden")))
+
+/* Values of magnitude below 2^52 round through an integer conversion: the
+ * truncated value, stepped by one when the rounding direction needs it, with
+ * the value's sign so that a zero result keeps it (ceil(-0.5) is -0). */
+static inline int sgm_small(double x) { return __builtin_fabs(x) < 0x1p52; }
+
+/* The unbiased exponent; 1024 for the infinities and NaN. Values with an
+ * exponent of 52 or more are integers already. */
+static inline int sgm_exponent(uint64_t bits) { return (int)((bits >> 52) & 0x7FF) - 1023; }
+
+static double sgm_floor(double x) {
+  if (sgm_small(x)) {
+    double t = (double)(int64_t)x;
+    return __builtin_copysign(t > x ? t - 1.0 : t, x);
+  }
+  uint64_t bits = sgm_bits(x);
+  int e = sgm_exponent(bits);
+  if (e >= 52) return e == 1024 ? x + x : x;
+  if (e < 0) {
+    if ((bits << 1) == 0) return x;
+    return bits & SGM_SIGN ? -1.0 : 0.0;
+  }
+  uint64_t fraction = (UINT64_C(1) << (52 - e)) - 1;
+  if ((bits & fraction) == 0) return x;
+  if (bits & SGM_SIGN) bits += fraction;
+  return sgm_double(bits & ~fraction);
+}
+
+static double sgm_ceil(double x) {
+  if (sgm_small(x)) {
+    double t = (double)(int64_t)x;
+    return __builtin_copysign(t < x ? t + 1.0 : t, x);
+  }
+  uint64_t bits = sgm_bits(x);
+  int e = sgm_exponent(bits);
+  if (e >= 52) return e == 1024 ? x + x : x;
+  if (e < 0) {
+    if ((bits << 1) == 0) return x;
+    return bits & SGM_SIGN ? -0.0 : 1.0;
+  }
+  uint64_t fraction = (UINT64_C(1) << (52 - e)) - 1;
+  if ((bits & fraction) == 0) return x;
+  if (!(bits & SGM_SIGN)) bits += fraction;
+  return sgm_double(bits & ~fraction);
+}
+
+static double sgm_trunc(double x) {
+  if (sgm_small(x)) return __builtin_copysign((double)(int64_t)x, x);
+  uint64_t bits = sgm_bits(x);
+  int e = sgm_exponent(bits);
+  if (e >= 52) return e == 1024 ? x + x : x;
+  if (e < 0) return sgm_double(bits & SGM_SIGN);
+  uint64_t fraction = (UINT64_C(1) << (52 - e)) - 1;
+  return sgm_double(bits & ~fraction);
+}
+
+#if defined(__x86_64__)
+#include <cpuid.h>
+#include <smmintrin.h>
+
+/* glibc selects SSE4.1's rounding instruction for floor, ceil and trunc when
+ * the CPU has it, and so do these: the same results, with no call into the
+ * integer path. */
+__attribute__((target("sse4.1"))) static double sgm_floor_sse41(double x) {
+  return _mm_cvtsd_f64(_mm_floor_sd(_mm_set_sd(x), _mm_set_sd(x)));
+}
+__attribute__((target("sse4.1"))) static double sgm_ceil_sse41(double x) {
+  return _mm_cvtsd_f64(_mm_ceil_sd(_mm_set_sd(x), _mm_set_sd(x)));
+}
+__attribute__((target("sse4.1"))) static double sgm_trunc_sse41(double x) {
+  return _mm_cvtsd_f64(_mm_round_sd(_mm_set_sd(x), _mm_set_sd(x), _MM_FROUND_TO_ZERO | _MM_FROUND_NO_EXC));
+}
+
+static int sgm_has_sse41(void) {
+  unsigned eax, ebx, ecx, edx;
+  return __get_cpuid(1, &eax, &ebx, &ecx, &edx) && (ecx & bit_SSE4_1) != 0;
+}
+
+typedef double (*SgmUnary)(double);
+static SgmUnary sgm_resolve_floor(void) { return sgm_has_sse41() ? sgm_floor_sse41 : sgm_floor; }
+static SgmUnary sgm_resolve_ceil(void) { return sgm_has_sse41() ? sgm_ceil_sse41 : sgm_ceil; }
+static SgmUnary sgm_resolve_trunc(void) { return sgm_has_sse41() ? sgm_trunc_sse41 : sgm_trunc; }
+SGM_DEFINE double floor(double x) __attribute__((ifunc("sgm_resolve_floor")));
+SGM_DEFINE double ceil(double x) __attribute__((ifunc("sgm_resolve_ceil")));
+SGM_DEFINE double trunc(double x) __attribute__((ifunc("sgm_resolve_trunc")));
+#else
+SGM_DEFINE double floor(double x) { return sgm_floor(x); }
+SGM_DEFINE double ceil(double x) { return sgm_ceil(x); }
+SGM_DEFINE double trunc(double x) { return sgm_trunc(x); }
+#endif
+
+/* Math.sqrt is LLVM's sqrt libcall: the square root instruction inline, and
+ * a call to sqrt for a negative argument, where libm would set errno. This
+ * definition answers that call with the same instruction. */
+SGM_DEFINE double sqrt(double x) {
+#if defined(__x86_64__)
+  return _mm_cvtsd_f64(_mm_sqrt_sd(_mm_set_sd(x), _mm_set_sd(x)));
+#elif defined(__aarch64__)
+  double r;
+  __asm__("fsqrt %d0, %d1" : "=w"(r) : "w"(x));
+  return r;
+#else
+  return __builtin_sqrt(x);
+#endif
+}
+
+/* Halfway cases round away from zero. */
+SGM_DEFINE double round(double x) {
+  uint64_t bits = sgm_bits(x);
+  int e = sgm_exponent(bits);
+  if (e >= 52) return e == 1024 ? x + x : x;
+  if (e < 0) {
+    uint64_t sign = bits & SGM_SIGN;
+    return e == -1 ? sgm_double(sign | sgm_bits(1.0)) : sgm_double(sign);
+  }
+  uint64_t fraction = (UINT64_C(1) << (52 - e)) - 1;
+  if ((bits & fraction) == 0) return x;
+  bits += UINT64_C(1) << (51 - e);
+  return sgm_double(bits & ~fraction);
+}
+
+/* The significand of a finite nonzero double as an integer below 2^53, and
+ * its exponent so that the value's magnitude is significand * 2^exponent. */
+static inline uint64_t sgm_significand(uint64_t bits, int *exponent) {
+  int biased = (int)((bits >> 52) & 0x7FF);
+  uint64_t significand = bits & ((UINT64_C(1) << 52) - 1);
+  if (biased == 0) {
+    *exponent = -1074;
+    return significand;
+  }
+  *exponent = biased - 1075;
+  return significand | (UINT64_C(1) << 52);
+}
+
+/* x - n*y for the integer n = trunc(x/y), which is exact: the remainder
+ * comes from shift-and-subtract over the significands. NaN, an infinite x
+ * and a zero y answer (x*y)/(x*y), as glibc's fmod does. */
+SGM_DEFINE double fmod(double x, double y) {
+  uint64_t xb = sgm_bits(x), yb = sgm_bits(y);
+  uint64_t xmag = xb & ~SGM_SIGN, ymag = yb & ~SGM_SIGN;
+  uint64_t infinity = UINT64_C(0x7FF0000000000000);
+  if (ymag == 0 || xmag >= infinity || ymag > infinity) return (x * y) / (x * y);
+  if (xmag < ymag) return x;
+  if (xmag == ymag) return 0.0 * x;
+  int ex, ey;
+  uint64_t mx = sgm_significand(xb, &ex);
+  uint64_t my = sgm_significand(yb, &ey);
+  /* Normalize both to 53 significant bits, so the loop below shifts one
+   * exponent step at a time. */
+  while (mx < (UINT64_C(1) << 52)) {
+    mx <<= 1;
+    ex--;
+  }
+  while (my < (UINT64_C(1) << 52)) {
+    my <<= 1;
+    ey--;
+  }
+  for (; ex > ey; ex--) {
+    if (mx >= my) mx -= my;
+    if (mx == 0) return 0.0 * x;
+    mx <<= 1;
+  }
+  if (mx >= my) mx -= my;
+  if (mx == 0) return 0.0 * x;
+  /* The remainder is mx * 2^ey, below y: rebuild it as a double. */
+  while (mx < (UINT64_C(1) << 52)) {
+    mx <<= 1;
+    ey--;
+  }
+  uint64_t result;
+  if (ey + 1075 >= 1) {
+    result = ((uint64_t)(ey + 1075) << 52) | (mx & ((UINT64_C(1) << 52) - 1));
+  } else {
+    result = mx >> (1 - (ey + 1075));
+  }
+  return sgm_double(result | (xb & SGM_SIGN));
+}
+
+/* x * 2^n, rounded once. Exponents beyond one double's range scale in
+ * steps; the step toward the subnormals stays 53 bits above the smallest
+ * normal so that only the last multiplication rounds. */
+SGM_DEFINE double scalbn(double x, int n) {
+  if (n > 1023) {
+    x *= 0x1p1023;
+    n -= 1023;
+    if (n > 1023) {
+      x *= 0x1p1023;
+      n -= 1023;
+      if (n > 1023) n = 1023;
+    }
+  } else if (n < -1022) {
+    x *= 0x1p-1022 * 0x1p53;
+    n += 1022 - 53;
+    if (n < -1022) {
+      x *= 0x1p-1022 * 0x1p53;
+      n += 1022 - 53;
+      if (n < -1022) n = -1022;
+    }
+  }
+  return x * sgm_double((uint64_t)(0x3FF + n) << 52);
+}
+
+SGM_DEFINE double ldexp(double x, int n) { return scalbn(x, n); }
+
+/* The significand in [0.5, 1) and the power of two: x = significand * 2^e.
+ * Zero, the infinities and NaN answer themselves with e = 0. */
+SGM_DEFINE double frexp(double x, int *e) {
+  uint64_t bits = sgm_bits(x);
+  int biased = (int)((bits >> 52) & 0x7FF);
+  if (biased == 0) {
+    if ((bits << 1) == 0) {
+      *e = 0;
+      return x;
+    }
+    x = frexp(x * 0x1p64, e);
+    *e -= 64;
+    return x;
+  }
+  if (biased == 0x7FF) {
+    *e = 0;
+    return x;
+  }
+  *e = biased - 1022;
+  return sgm_double((bits & ~(UINT64_C(0x7FF) << 52)) | (UINT64_C(0x3FE) << 52));
+}
+
+#pragma clang attribute pop
+
+#endif

packages/runtime/src/scr_ieee754.h

@@ -0,0 +1,26 @@
+#ifndef SCR_IEEE754_H
+#define SCR_IEEE754_H
+
+/* Math's transcendental functions as V8 computes them (scr_ieee754.c). */
+double scr_ieee754_acos(double x);
+double scr_ieee754_acosh(double x);
+double scr_ieee754_asin(double x);
+double scr_ieee754_asinh(double x);
+double scr_ieee754_atan(double x);
+double scr_ieee754_atan2(double y, double x);
+double scr_ieee754_atanh(double x);
+double scr_ieee754_cbrt(double x);
+double scr_ieee754_cos(double x);
+double scr_ieee754_cosh(double x);
+double scr_ieee754_exp(double x);
+double scr_ieee754_expm1(double x);
+double scr_ieee754_log(double x);
+double scr_ieee754_log10(double x);
+double scr_ieee754_log1p(double x);
+double scr_ieee754_log2(double x);
+double scr_ieee754_sin(double x);
+double scr_ieee754_sinh(double x);
+double scr_ieee754_tan(double x);
+double scr_ieee754_tanh(double x);
+
+#endif

packages/runtime/src/scr_js_value.c

@@ -209,7 +209,7 @@ static double scr_js_number_arithmetic(double a, double b, const ScrStr *operati
   case '/': return a / b;
   case '%': return fmod(a, b);
   default:
-    if (operation->len == 2) return (fabs(a) == 1 && !isfinite(b)) ? NAN : pow(a, b);
+    if (operation->len == 2) return scr_math_pow(a, b);
     return a * b;
   }
 }

packages/runtime/src/scr_json.c

@@ -11160,9 +11160,8 @@ ScrDyn *scr_dyn_arithmetic(const ScrDyn *left, const ScrDyn *right, const ScrStr
     if (result) { out = scr_dyn_new_bigint(result); scr_bigint_release(result); }
   } else {
     double a = l->v.num, b = r->v.num;
-    // ECMAScript differs from libm for ±1 ** ±Infinity and 1 ** NaN.
     double result = op == '-' ? a - b : op == '/' ? a / b : op == '%' ? fmod(a, b)
-      : operation->len == 2 ? ((fabs(a) == 1 && !isfinite(b)) ? NAN : pow(a, b)) : a * b;
+      : operation->len == 2 ? scr_math_pow(a, b) : a * b;
     out = scr_dyn_new_num(result);
   }
   scr_dyn_release(l);

packages/runtime/src/scr_lib.c

@@ -6828,11 +6828,17 @@ double scr_math_round(double x) {
   return (r == 0.0 && x < 0.0) ? -0.0 : r;
 }
 
-/* IEEE pow(±1, ±Infinity) and pow(1, NaN) return 1; ECMAScript requires
+/* V8's math::pow: the platform pow, as Node's Math.pow and ** call it,
- * NaN. The other special cases, including signed zero and odd integer
+ * after the cases V8 answers itself. IEEE pow(±1, ±Infinity) and
- * powers, follow libm. */
+ * pow(1, NaN) return 1; ECMAScript requires NaN. A NaN base becomes the
+ * quiet NaN, and the exponents 2 and 0.5 answer base * base and
+ * sqrt(base + 0) (+Infinity for an infinite base), as V8's optimizing
+ * compilers do; pow can differ from those in the last bit. */
 double scr_math_pow(double base, double exponent) {
   if (isnan(exponent) || (isinf(exponent) && fabs(base) == 1.0)) return NAN;
+  if (isnan(base)) base = NAN;
+  if (exponent == 2.0) return base * base;
+  if (exponent == 0.5) return isinf(base) ? INFINITY : sqrt(base + 0.0);
   /* 2 ** n for an integer n in the double range is an exact power of two,
    * the value pow returns; its bits are the exponent alone. */
   if (base == 2.0 && exponent >= -1074.0 && exponent <= 1023.0) {
@@ -12815,8 +12821,15 @@ ScrStr *scr_num_to_string_radix_f64(double value, double requested) {
   bool negative = value < 0;
   value = fabs(value);
   double integer = floor(value), fraction = value - integer;
-  double tolerance = 0.5 * (nextafter(value, INFINITY) - value);
+  /* The next double above the positive finite value: its bits plus one
-  if (tolerance <= 0) tolerance = nextafter(0, INFINITY);
+   * (Infinity above the largest double). */
+  uint64_t value_bits;
+  memcpy(&value_bits, &value, sizeof value_bits);
+  value_bits += 1;
+  double next;
+  memcpy(&next, &value_bits, sizeof next);
+  double tolerance = 0.5 * (next - value);
+  if (tolerance <= 0) tolerance = 0x1p-1074;
   if (fraction >= tolerance) {
     buffer[end++] = '.';
     do {

packages/runtime/src/scr_runtime.h

@@ -31,6 +31,7 @@
 #include "scr_string_methods.h"
 #include "scr_qs.h"
 #include "scr_array.h"
+#include "scr_ieee754.h"
 #include "scr_regex.h"
 #include "scr_map.h"
 #include "scr_closure.h"

packages/runtime/vendor/README.md

@@ -16,6 +16,13 @@ The x86_64 and AArch64 register-save/restore and `makecontext` logic is adapted
 
 The double-to-shortest-decimal core of Ryū (d2s and its headers only — the float and fixed/exponent variants are not vendored). scr_number.c textually `#include`s d2s.c and calls its internal `d2d` digit generator; ECMA-262 digit *placement* stays in scr_number.c. The only modification to upstream files: `#include "ryu/x.h"` flattened to `#include "x.h"` in d2s.c and d2s_intrinsics.h so the directory is self-contained. To update, re-fetch the listed files at the new commit and re-apply that include flattening.
 
+## v8-ieee754/ (derived source)
+
+- Project: https://github.com/v8/v8, `src/base/ieee754.cc` as Node v24.22.0 bundles it (`deps/v8`, V8 13.6.233.17)
+- License: BSD-3-Clause (V8), with fdlibm's notice (see v8-ieee754/LICENSE)
+
+V8's adaptation of fdlibm, the code behind Node's Math.sin, Math.exp, Math.cbrt and the other transcendental functions on every platform. `src/scr_ieee754.c` is a C port of the file, so compiled Math members print the values Node prints. The function bodies are unchanged; the port renames the entry points (`scr_ieee754_*`), gives the kernels internal linkage, and replaces V8's C++ helpers (bit casts, wraparound arithmetic, NaN and infinity constants) with C. `legacy::pow` is not ported: Node's Math.pow calls the platform pow. To update, re-fetch `deps/v8/src/base/ieee754.cc` at the new Node version and repeat the same conversion.
+
 ## quickjs-ng/
 
 - Project: https://github.com/quickjs-ng/quickjs

packages/runtime/vendor/v8-ieee754/LICENSE

@@ -0,0 +1,75 @@
+This license applies to all parts of V8 that are not externally
+maintained libraries.  The externally maintained libraries used by V8
+are:
+
+  - PCRE test suite, located in
+    test/mjsunit/third_party/regexp-pcre/regexp-pcre.js.  This is based on the
+    test suite from PCRE-7.3, which is copyrighted by the University
+    of Cambridge and Google, Inc.  The copyright notice and license
+    are embedded in regexp-pcre.js.
+
+  - Layout tests, located in test/mjsunit/third_party/object-keys.  These are
+    based on layout tests from webkit.org which are copyrighted by
+    Apple Computer, Inc. and released under a 3-clause BSD license.
+
+  - Strongtalk assembler, the basis of the files assembler-arm-inl.h,
+    assembler-arm.cc, assembler-arm.h, assembler-ia32-inl.h,
+    assembler-ia32.cc, assembler-ia32.h, assembler-x64-inl.h,
+    assembler-x64.cc, assembler-x64.h, assembler.cc and assembler.h.
+    This code is copyrighted by Sun Microsystems Inc. and released
+    under a 3-clause BSD license.
+
+  - Valgrind client API header, located at third_party/valgrind/valgrind.h
+    This is released under the BSD license.
+
+  - The Wasm C/C++ API headers, located at third_party/wasm-api/wasm.{h,hh}
+    This is released under the Apache license. The API's upstream prototype
+    implementation also formed the basis of V8's implementation in
+    src/wasm/c-api.cc.
+
+These libraries have their own licenses; we recommend you read them,
+as their terms may differ from the terms below.
+
+Further license information can be found in LICENSE files located in
+sub-directories.
+
+Copyright 2014, the V8 project authors. All rights reserved.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * 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.
+    * Neither the name of Google Inc. nor the names of its
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 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 COPYRIGHT
+OWNER 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.
+
+---
+
+src/scr_ieee754.c is a port of V8's src/base/ieee754.cc, which V8 adapted
+from fdlibm. fdlibm's notice:
+
+====================================================
+Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+
+Developed at SunSoft, a Sun Microsystems, Inc. business.
+Permission to use, copy, modify, and distribute this
+software is freely granted, provided that this notice
+is preserved.
+====================================================

tests/corpus/1110-math-members.ts

@@ -1,20 +1,17 @@
-// Math members execute in the island engine; results exit as ordinary
+// Math members print the values Node prints, the transcendentals
-// static numbers. Members with exact IEEE results print directly. The
+// included: they are V8's fdlibm functions.
-// transcendentals round through toFixed(9): the engine uses the system
-// libm and Node uses V8's fdlibm port, which can differ in the last ulp
-// (tan(1), asin(0.5), acos(0.5) do on macOS).
 console.log(Math.floor(7.6), Math.ceil(7.2), Math.round(7.5), Math.round(-7.5), Math.trunc(-7.9));
 console.log(Math.abs(-3.5), Math.sign(-12), Math.sign(0), Math.sign(4.2));
 console.log(Math.sqrt(144), Math.sqrt(2), Math.pow(2, 10), Math.pow(2, 0.5));
 console.log(Math.hypot(3, 4), Math.min(2, -9), Math.max(2, -9));
 console.log(Math.log2(1024), Math.log10(100000), Math.exp(0), Math.log(1));
 console.log(Math.PI, Math.E);
-console.log(Math.cbrt(27).toFixed(9), Math.exp(1).toFixed(9), Math.log(10).toFixed(9));
+console.log(Math.cbrt(27), Math.exp(1), Math.log(10));
-console.log(Math.sin(1).toFixed(9), Math.cos(1).toFixed(9), Math.tan(1).toFixed(9));
+console.log(Math.sin(1), Math.cos(1), Math.tan(1));
-console.log(Math.asin(0.5).toFixed(9), Math.acos(0.5).toFixed(9));
+console.log(Math.asin(0.5), Math.acos(0.5));
-console.log(Math.atan(1).toFixed(9), Math.atan2(1, 1).toFixed(9), Math.atan2(-1, -1).toFixed(9));
+console.log(Math.atan(1), Math.atan2(1, 1), Math.atan2(-1, -1));
-// Island results chain into static arithmetic like any number.
+// Math results chain into static arithmetic like any number.
 const area = Math.PI * Math.pow(2, 2);
 console.log(area.toFixed(6), Math.floor(area) + 1);
-// Non-finite edges: strict exits accept every engine number.
+// Non-finite results.
 console.log(Math.log(0), Math.log(-1), Math.sqrt(-1), Math.pow(0, -1));

tests/corpus/7717-math-fdlibm.ts

@@ -0,0 +1,97 @@
+// Math's transcendental functions print the values Node prints: each line
+// digests the result bits of one function over a few thousand inputs, so
+// a difference in the last bit of any result changes the line. Sample
+// values print in full, and pow's special cases follow V8's.
+const bits = new Float64Array(1);
+const words = new Uint32Array(bits.buffer);
+
+function digest(values: number[]): string {
+  let h1 = 0x811c9dc5;
+  let h2 = 0x01000193;
+  for (const v of values) {
+    // NaN prints as NaN; its bits are not part of the result.
+    bits[0] = Number.isNaN(v) ? 0 : v;
+    h1 = Math.imul(h1 ^ words[0], 0x01000193) >>> 0;
+    h2 = Math.imul(h2 ^ words[1], 0x01000193) >>> 0;
+  }
+  return `${h1}:${h2}`;
+}
+
+let seed = 12345;
+function next(): number {
+  seed = (seed * 1103515245 + 12345) % 2147483648;
+  return seed / 2147483648;
+}
+
+const wide: number[] = [];
+const unit: number[] = [];
+const small: number[] = [];
+for (let i = 0; i < 3000; i++) {
+  wide.push((next() - 0.5) * Math.pow(10, Math.floor(next() * 12) - 4));
+  unit.push(next() * 2 - 1);
+  small.push((next() - 0.5) * 40);
+}
+const positive = wide.map((x) => Math.abs(x));
+const edges = [0, -0, 1, -1, 0.5, -0.5, 2, Math.PI, Math.E, 1e-300, 5e-324, 1e300, Infinity, -Infinity, NaN];
+
+const unary: [string, (x: number) => number, number[]][] = [
+  ["sin", Math.sin, wide],
+  ["cos", Math.cos, wide],
+  ["tan", Math.tan, wide],
+  ["asin", Math.asin, unit],
+  ["acos", Math.acos, unit],
+  ["atan", Math.atan, wide],
+  ["sinh", Math.sinh, small],
+  ["cosh", Math.cosh, small],
+  ["tanh", Math.tanh, small],
+  ["asinh", Math.asinh, wide],
+  ["acosh", Math.acosh, positive.map((x) => x + 1)],
+  ["atanh", Math.atanh, unit],
+  ["exp", Math.exp, small],
+  ["expm1", Math.expm1, small],
+  ["log", Math.log, positive],
+  ["log1p", Math.log1p, positive],
+  ["log2", Math.log2, positive],
+  ["log10", Math.log10, positive],
+  ["cbrt", Math.cbrt, wide],
+  ["sqrt", Math.sqrt, positive],
+];
+for (const [name, fn, inputs] of unary) {
+  console.log(name, digest(inputs.map(fn)), digest(edges.map(fn)), fn(0.3), fn(0.75));
+}
+
+// Math.sin(i) and Math.cos(i * 0.37) for i below 4096: main's libm calls
+// printed 284 of these 8192 values differently from Node.
+const sines: number[] = [];
+const cosines: number[] = [];
+for (let i = 0; i < 4096; i++) {
+  sines.push(Math.sin(i));
+  cosines.push(Math.cos(i * 0.37));
+}
+console.log("sin(i)", digest(sines), "cos(i*0.37)", digest(cosines), Math.cos(2.59));
+
+const atan2: number[] = [];
+const pow: number[] = [];
+for (let i = 0; i < wide.length; i++) {
+  atan2.push(Math.atan2(wide[i], small[i]));
+  pow.push(Math.pow(positive[i] % 50, small[i] / 4));
+}
+console.log("atan2", digest(atan2), Math.atan2(2.5, 3), Math.atan2(-0, -1));
+console.log("pow", digest(pow), Math.pow(1.1, 0.3), 1.0000001 ** 1e7);
+
+// V8 answers these itself: x * x, the square root, and NaN for ±1 ** ±Infinity.
+const squares = wide.map((x) => x ** 2);
+const roots = positive.map((x) => Math.pow(x, 0.5));
+console.log("pow 2", digest(squares), "pow 0.5", digest(roots));
+// The same exponents when they are only known at run time.
+const two = Number(process.argv[2] ?? "2");
+const half = two / 4;
+console.log("pow runtime 2", digest(wide.map((x) => x ** two)), "0.5", digest(positive.map((x) => Math.pow(x, half))));
+console.log(Math.pow(-Infinity, 0.5), Math.pow(-0, 0.5), (-8) ** 0.5, 1 ** Infinity, (-1) ** -Infinity, NaN ** 0);
+// Rounding is exact: Math.floor, Math.ceil and Math.trunc over the same
+// inputs and over halves, signed zeros included.
+const halves = wide.map((x, i) => Math.round(x * 4) / 4 - (i % 2 === 0 ? 0 : 0.5));
+console.log("floor", digest(wide.map(Math.floor)), digest(halves.map(Math.floor)), digest(edges.map(Math.floor)));
+console.log("ceil", digest(wide.map(Math.ceil)), digest(halves.map(Math.ceil)), digest(edges.map(Math.ceil)));
+console.log("trunc", digest(wide.map(Math.trunc)), digest(halves.map(Math.trunc)), Object.is(Math.ceil(-0.5), -0), Object.is(Math.floor(-0), -0));
+console.log(Math.hypot(3, 4), Math.hypot(1e200, 1e200), 7 % 2.5, -7.5 % 2, Math.floor(-0.5), Math.ceil(-0.5), Math.round(2.5), Math.trunc(-4.7));

tests/corpus/7718-math-fdlibm.cjs

@@ -0,0 +1,12 @@
+// Math functions and ** over untyped values print the values Node prints
+// (V8's fdlibm results, and V8's pow special cases).
+const values = JSON.parse("[0.3, 0.75, 2.5, -7.25, 1e-8, 123456.789, 1e21]");
+const rows = [];
+for (const x of values) {
+  rows.push([Math.sin(x), Math.cos(x), Math.tan(x), Math.exp(x / 100), Math.cbrt(x), Math.atan2(x, 3)].join(" "));
+  rows.push([Math.log(Math.abs(x)), Math.log10(Math.abs(x)), Math.sinh(x / 1e5), Math.cosh(x / 1e5), Math.expm1(x / 1e6)].join(" "));
+  rows.push([x ** 2, x ** 0.5, Math.abs(x) ** 0.3, x % 1.5].join(" "));
+}
+console.log(rows.join("\n"));
+const edges = JSON.parse('[1, -1, null]');
+console.log(edges[0] ** Infinity, edges[1] ** -Infinity, Number(edges[2]) ** 0.5, (-Infinity) ** 0.5);

tests/harness/corpus-costs/7717-math-fdlibm.ts.json

@@ -0,0 +1 @@
+{ "plain": 3357, "san": 4900 }

tests/harness/corpus-costs/7718-math-fdlibm.cjs.json

@@ -0,0 +1 @@
+{ "plain": 3519, "san": 4613 }