Compute Math.pow(2, n) for integer exponents from the result's bits

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

@@ -969,6 +969,64 @@ export function emitPathUrlLibCall(host: LlvmEmitterContext, e: LibCallExpr): Ll
   return host.emitGenericLibCall(e);
 }
 
+/** Math.pow(2, n): for an integer n from -1074 to 1023 the result is the
+ * power of two whose bits are its exponent field (or, below -1022, one
+ * mantissa bit of a subnormal), the value pow returns. Other exponents call
+ * scr_math_pow. */
+function emitPowerOfTwo(host: LlvmEmitterContext, n: string): string {
+  const B = host.B;
+  const check = B.newLabel("pow2.check");
+  const exact = B.newLabel("pow2.exact");
+  const call = B.newLabel("pow2.call");
+  const join = B.newLabel("pow2.join");
+  const low = B.tmp();
+  const high = B.tmp();
+  const inRange = B.tmp();
+  B.line(`${low} = fcmp oge double ${n}, ${f64Lit(-1074)}`);
+  B.line(`${high} = fcmp ole double ${n}, ${f64Lit(1023)}`);
+  B.line(`${inRange} = and i1 ${low}, ${high}`);
+  B.condBr(inRange, check, call);
+  B.startBlock(check);
+  const int = B.tmp();
+  const back = B.tmp();
+  const integral = B.tmp();
+  B.line(`${int} = fptosi double ${n} to i32`);
+  B.line(`${back} = sitofp i32 ${int} to double`);
+  B.line(`${integral} = fcmp oeq double ${back}, ${n}`);
+  B.condBr(integral, exact, call);
+  B.startBlock(exact);
+  const normal = B.tmp();
+  const biased = B.tmp();
+  const biasedWide = B.tmp();
+  const normalBits = B.tmp();
+  const shift = B.tmp();
+  const shiftWide = B.tmp();
+  const subnormalBits = B.tmp();
+  const bits = B.tmp();
+  const power = B.tmp();
+  B.line(`${normal} = icmp sge i32 ${int}, -1022`);
+  B.line(`${biased} = add i32 ${int}, 1023`);
+  B.line(`${biasedWide} = zext i32 ${biased} to i64`);
+  B.line(`${normalBits} = shl i64 ${biasedWide}, 52`);
+  B.line(`${shift} = add i32 ${int}, 1074`);
+  B.line(`${shiftWide} = zext i32 ${shift} to i64`);
+  B.line(`${subnormalBits} = shl i64 1, ${shiftWide}`);
+  B.line(`${bits} = select i1 ${normal}, i64 ${normalBits}, i64 ${subnormalBits}`);
+  B.line(`${power} = bitcast i64 ${bits} to double`);
+  const exactEnd = B.currentLabel();
+  B.br(join);
+  B.startBlock(call);
+  host.declare(`declare double @scr_math_pow(double, double)`);
+  const called = B.tmp();
+  B.line(`${called} = call double @scr_math_pow(double ${f64Lit(2)}, double ${n})`);
+  const callEnd = B.currentLabel();
+  B.br(join);
+  B.startBlock(join);
+  const out = B.tmp();
+  B.line(`${out} = phi double [ ${power}, %${exactEnd} ], [ ${called}, %${callEnd} ]`);
+  return out;
+}
+
 export function emitPrimitiveLibCall(host: LlvmEmitterContext, e: LibCallExpr): LlValue {
   const B = host.B;
   if (e.fn === "math.min" || e.fn === "math.max") {
@@ -986,6 +1044,10 @@ export function emitPrimitiveLibCall(host: LlvmEmitterContext, e: LibCallExpr):
     const v = host.emitExpr(e.args[0]!);
     return { name: emitTruncF64(host, v.name), type: e.type };
   }
+  if (e.fn === "math.pow" && e.args[0]!.kind === "numLit" && e.args[0]!.value === 2) {
+    const v = host.emitExpr(e.args[1]!);
+    return { name: emitPowerOfTwo(host, v.name), 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/test/ts7/baselines/order-parity/tests/corpus/6863-math-pow-base-two.ts.json

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

packages/runtime/src/scr_lib.c

@@ -6821,6 +6821,17 @@ double scr_math_round(double x) {
  * powers, follow libm. */
 double scr_math_pow(double base, double exponent) {
   if (isnan(exponent) || (isinf(exponent) && fabs(base) == 1.0)) return NAN;
+  /* 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) {
+    int32_t n = (int32_t)exponent;
+    if ((double)n == exponent) {
+      uint64_t bits = n >= -1022 ? (uint64_t)(n + 1023) << 52 : UINT64_C(1) << (n + 1074);
+      double result;
+      memcpy(&result, &bits, sizeof result);
+      return result;
+    }
+  }
   return pow(base, exponent);
 }
 

tests/corpus/6863-math-pow-base-two.ts

@@ -0,0 +1,51 @@
+// Math.pow(2, n) with a literal base of 2 and with a base of 2 in a variable,
+// over every integer exponent from -1080 to 1030: subnormal results down to
+// 2 ** -1074, zero below it, Infinity above 1023. Non-integer, signed-zero,
+// infinite and NaN exponents, and bases next to 2, take the general path.
+function bits(x: number): string {
+  const view = new DataView(new ArrayBuffer(8));
+  view.setFloat64(0, x);
+  return view.getBigUint64(0).toString(16).padStart(16, "0");
+}
+
+function show(label: string, x: number): void {
+  console.log(label, bits(x), Object.is(x, -0) ? "-0" : String(x));
+}
+
+let base = 2;
+if (process.argv.length > 99) base = 3;
+
+let literalDigest = 0n;
+let variableDigest = 0n;
+let mismatches = 0;
+for (let n = -1080; n <= 1030; n++) {
+  const literal = Math.pow(2, n);
+  const variable = Math.pow(base, n);
+  literalDigest = (literalDigest * 31n + BigInt(`0x${bits(literal)}`)) & 0xffffffffffffffffn;
+  variableDigest = (variableDigest * 31n + BigInt(`0x${bits(variable)}`)) & 0xffffffffffffffffn;
+  if (!Object.is(literal, variable)) mismatches++;
+  if (n <= -1070 || (n >= -1025 && n <= -1019) || (n >= -2 && n <= 2) || (n >= 51 && n <= 54) || n >= 1021) {
+    show(`2 ** ${n}`, literal);
+  }
+}
+console.log("digests", literalDigest.toString(16), variableDigest.toString(16), "mismatches", mismatches);
+
+const exponents = [0.5, -0.5, 1.5, -0, 0, 1e-300, -1e-300, 1023.0000001, 1022.9999999, -1074.5, -1073.5, Infinity, -Infinity, NaN, 52.999999999999996, 2 ** 31, -(2 ** 31)];
+for (const e of exponents) {
+  show(`literal ${e}`, Math.pow(2, e));
+  show(`variable ${e}`, Math.pow(base, e));
+}
+for (const b of [-2, 2.0000000000000004, 1.9999999999999998, 4, 0.5]) {
+  show(`${b} ** 10`, Math.pow(b, 10));
+  show(`${b} ** -3`, Math.pow(b, -3));
+}
+
+// In a loop that masks values, as a k-mer counter does.
+let key = 0;
+let total = 0;
+for (let i = 0; i < 100000; i++) {
+  const length = (i % 18) + 1;
+  key = (key * 4 + (i % 4)) % Math.pow(2, 2 * length);
+  total += key - Math.floor(key / Math.pow(2, length)) * Math.pow(2, length);
+}
+console.log("masked total", total);

tests/harness/corpus-costs/6863-math-pow-base-two.ts.json

@@ -0,0 +1 @@
+{ "plain": 2110, "san": 4860 }