Fix cube roots: real odd roots of negative bases

lib/expr.test.ts

@@ -162,6 +162,44 @@ describe('toGLSL', () => {
   });
 });
 
+describe('negative base with fractional exponent (real odd roots)', () => {
+  const ev = (s: string, env: Record<string, number> = {}) => evaluate(parseExpr(s), env);
+
+  it('takes the real cube root of a negative base', () => {
+    expect(ev('(-8)^(1/3)')).toBeCloseTo(-2);
+    expect(ev('(-1)^(1/3)')).toBeCloseTo(-1);
+  });
+
+  it('gives a positive result for an even numerator over an odd denominator', () => {
+    expect(ev('(-8)^(2/3)')).toBeCloseTo(4);
+  });
+
+  it('handles negative fractional exponents', () => {
+    expect(ev('(-8)^(-1/3)')).toBeCloseTo(-0.5);
+  });
+
+  it('leaves even roots of a negative base undefined', () => {
+    expect(ev('(-4)^(1/2)')).toBeNaN();
+    expect(ev('(-8)^(1/4)')).toBeNaN();
+  });
+
+  it('does not snap a typed decimal approximation to a nearby rational', () => {
+    // 0.33333 is ~3.3e-6 away from 1/3, outside the 1e-6 tolerance, so this
+    // stays undefined rather than being guessed at as a cube root.
+    expect(ev('(-8)^(0.33333)')).toBeNaN();
+  });
+
+  it('leaves positive bases with fractional exponents unaffected', () => {
+    expect(ev('8^(1/3)')).toBeCloseTo(2);
+    expect(ev('8^(2/3)')).toBeCloseTo(4);
+  });
+
+  it('still handles integer exponents on negative bases', () => {
+    expect(ev('(-2)^3')).toBe(-8);
+    expect(ev('(-2)^2')).toBe(4);
+  });
+});
+
 describe('absolute value bars', () => {
   const ev = (s: string, env: Record<string, number> = {}) => evaluate(parseExpr(s), env);
 

lib/expr.ts

@@ -456,6 +456,55 @@ export const normalcdf = (x: number, mean: number, sd: number): number =>
  */
 export const ISPRIME_MAX = 2048 * 2048 - 1;
 
+/** a^b tolerance for snapping the exponent to a small rational p/q (real odd roots). */
+const POW_RATIONAL_TOL = 1e-6;
+
+/** Largest denominator considered when looking for the exponent's rational form. */
+const POW_RATIONAL_MAX_Q = 12;
+
+/**
+ * Real-valued a^b, matching how graphing calculators (e.g. Desmos) treat a
+ * negative base with a fractional exponent: real odd roots come out real
+ * (e.g. (-8)^(1/3) = -2) instead of NaN, while even roots stay undefined
+ * (e.g. (-4)^(1/2)).
+ *
+ * For a >= 0 this is just Math.pow. For a < 0, Math.pow only agrees with the
+ * "real odd root" convention when b happens to be an exact integer, so
+ * instead we search for the exponent's rational form p/q in lowest terms via
+ * a tolerance search over small denominators (q = 1..POW_RATIONAL_MAX_Q — no
+ * arbitrary-precision rational type needed, just simple fractions like 1/3,
+ * 2/3, 1/5). If q is odd, the root is real: sign * |a|^b, where sign is
+ * negative iff p (the reduced numerator) is odd. If no small-denominator
+ * match is found within tolerance (an irrational-looking exponent) or q is
+ * even (an even root of a negative number), the result is NaN, same as
+ * plain Math.pow.
+ *
+ * The tolerance (1e-6) is deliberately tight: an exponent entered as a
+ * fraction (e.g. "1/3") lands within ~1e-16 of the true rational, so it
+ * always snaps, but a typed decimal approximation like 0.33333 is ~3.3e-6
+ * away from 1/3 — outside tolerance — and is left undefined rather than
+ * silently guessed at.
+ *
+ * Kept in sync with the eq_pow() GLSL twin in glsl.ts (same algorithm, same
+ * tolerance and max denominator, adapted to GLSL's lack of a gcd builtin).
+ */
+export function realPow(a: number, b: number): number {
+  if (a >= 0) return Math.pow(a, b);
+  for (let q = 1; q <= POW_RATIONAL_MAX_Q; q++) {
+    const p = Math.round(b * q);
+    let x = Math.abs(p), y = q;
+    while (y) { const t = x % y; x = y; y = t; } // gcd(|p|, q)
+    const g = x || 1;
+    const pr = p / g, qr = q / g;
+    if (Math.abs(b - pr / qr) < POW_RATIONAL_TOL) {
+      if (qr % 2 === 0) return NaN; // even root of a negative number: undefined
+      const sign = Math.abs(pr) % 2 === 1 ? -1 : 1;
+      return sign * Math.pow(-a, b);
+    }
+  }
+  return NaN; // no small-denominator rational found: irrational-looking exponent
+}
+
 const EVAL_FNS: Record<string, (...xs: number[]) => number> = {
   sin: Math.sin, cos: Math.cos, tan: Math.tan,
   asin: Math.asin, acos: Math.acos, atan: Math.atan, atan2: Math.atan2,
@@ -503,7 +552,7 @@ export function evaluate(e: Expr, env: Record<string, number>): number {
         case '-': return a - b;
         case '*': return a * b;
         case '/': return a / b;
-        case '^': return Math.pow(a, b);
+        case '^': return realPow(a, b);
       }
     }
     case 'call': {

lib/glsl.ts

@@ -64,14 +64,31 @@ float eq_isprime(float x) {
   return 1.0;
 }
 float eq_pow(float a, float b) {
-  // Support negative bases for integer exponents (e.g. (-2)^3).
+  // Support negative bases via the "real odd root" convention, e.g.
+  // (-8)^(1/3) = -2, matching graphing calculators like Desmos. For a < 0,
+  // find the exponent's rational form p/q in lowest terms via a tolerance
+  // search over small denominators (q = 1..12); if q is odd, the root is
+  // real: sign * |a|^b, sign negative iff the reduced numerator p is odd.
+  // No match within tolerance, or an even q (an even root, e.g. (-4)^(1/2)),
+  // leaves the result NaN. Tolerance 1e-6 is tight enough that a typed
+  // decimal like 0.33333 (~3.3e-6 from 1/3) is left undefined rather than
+  // silently snapped. Kept in sync with realPow() in expr.ts (same
+  // algorithm; eq_gcd stands in for the CPU version's inline gcd loop).
   if (a >= 0.0) return pow(a, b);
-  float bi = floor(b + 0.5);
-  if (abs(b - bi) < 1e-6) {
-    float m = pow(-a, b);
-    return mod(bi, 2.0) < 0.5 ? m : -m;
+  for (int q = 1; q <= 12; q++) {
+    float fq = float(q);
+    float p = floor(b * fq + 0.5);
+    float g = eq_gcd(abs(p), fq);
+    if (g < 1.0) g = 1.0;
+    float pr = p / g;
+    float qr = fq / g;
+    if (abs(b - pr / qr) < 1e-6) {
+      if (mod(qr, 2.0) < 0.5) return sqrt(-1.0); // even root: undefined
+      float sign = mod(abs(pr), 2.0) > 0.5 ? -1.0 : 1.0;
+      return sign * pow(-a, b);
+    }
   }
-  return sqrt(-1.0); // NaN: undefined for negative base, fractional exponent
+  return sqrt(-1.0); // no small-denominator rational found
 }
 // Complex arithmetic on vec2(re, im).
 vec2 c_mul(vec2 a, vec2 b) { return vec2(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); }