diff --git a/lib/expr.test.ts b/lib/expr.test.ts index c378db2..707c156 100644 --- a/lib/expr.test.ts +++ b/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 = {}) => 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 = {}) => evaluate(parseExpr(s), env); diff --git a/lib/expr.ts b/lib/expr.ts index c99717b..91a6c81 100644 --- a/lib/expr.ts +++ b/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 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): 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': { diff --git a/lib/glsl.ts b/lib/glsl.ts index 8aa353d..4a82c65 100644 --- a/lib/glsl.ts +++ b/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); }