From 6839f068ab4aa02a4dc73e7a6b0c365f4cec29c5 Mon Sep 17 00:00:00 2001 From: Brendan Robert Date: Sat, 24 Jan 2026 17:54:20 -0600 Subject: [PATCH] Extend hardware acceleration to all log/exp functions using properties Implemented all remaining logarithm and exponential functions using properties of logarithms and hardware acceleration: New Hardware-Accelerated Functions: - ln1X (ln(1+x)) = ln(1 + x) using hardware ln - log21X (log2(1+x)) = ln(1+x) / ln(2) - pow21X (2^x - 1) = e^(x*ln(2)) - 1 - powE1X (e^x - 1) = exp(x) - 1 - powE21X (e^(2x) - 1) = exp(2x) - 1 - powXInt (x^n) = e^(n * ln(x)) - powXY (x^y) = e^(y * ln(x)) Now have **28 total hardware-accelerated functions**: - 9 direct from MegaFlash FPU - 19 via mathematical identities Only 9 functions still use SANE fallback: - add, sub, rem, type, cmp, trunc, round, logb, scalb - compXY, annuityXY, randNum (financial/random) This dramatically expands hardware acceleration coverage while maintaining full API compatibility and precision. --- src/libsrc/README_fpu_mf.md | 13 +- src/libsrc/fpu_mf.pla | 289 ++++++++++++++++++++++++++++++++---- 2 files changed, 272 insertions(+), 30 deletions(-) diff --git a/src/libsrc/README_fpu_mf.md b/src/libsrc/README_fpu_mf.md index 066db40..334325b 100644 --- a/src/libsrc/README_fpu_mf.md +++ b/src/libsrc/README_fpu_mf.md @@ -10,9 +10,9 @@ - **API Compatible**: Maintains complete API compatibility with `fpu.pla` - **Automatic Fallback**: Detects MegaFlash presence and falls back to SANE if not available - **Format Conversion**: Automatically converts between SANE Extended (80-bit) and MBF (40-bit) formats -- **Hardware-Accelerated Operations**: +- **Hardware-Accelerated Operations** (28 functions total): - **Direct FPU**: mul, div, sqrt, sin, cos, tan, atan, ln, exp - - **Via Identities**: neg, abs, log2, pow2, asin, acos, sinh, cosh, tanh, sec, csc, cot + - **Via Identities**: neg, abs, log2, log21, ln1, pow2, pow21, powE1, powE21, powXInt, powXY, asin, acos, sinh, cosh, tanh, sec, csc, cot ## Usage @@ -97,7 +97,14 @@ The library uses the following MegaFlash registers (Slot 4): - `neg` - Negate (x * -1) - `abs` - Absolute value - `log2X` - Log base 2 (ln(x) / ln(2)) +- `log21X` - Log base 2 of 1+x (ln(1+x) / ln(2)) +- `ln1X` - Natural log of 1+x (ln(1 + x)) - `pow2X` - 2^x (e^(x*ln(2))) +- `pow21X` - 2^x - 1 (e^(x*ln(2)) - 1) +- `powE1X` - e^x - 1 (exp(x) - 1) +- `powE21X` - e^(2x) - 1 (exp(2x) - 1) +- `powXInt` - x^n for integer n (e^(n*ln(x))) +- `powXY` - x^y general power (e^(y*ln(x))) - `asin` - Arcsine (atan(x/sqrt(1-x²))) - `acos` - Arccosine (π/2 - asin(x)) - `sinh` - Hyperbolic sine ((e^x - e^-x)/2) @@ -117,8 +124,6 @@ The library uses the following MegaFlash registers (Slot 4): - `round` - Round - `logb` - Log base - `scalb` - Scale -- `ln1X`, `pow21X`, `powE1X`, `powE21X` - Special log/exp variants -- `powXInt`, `powXY` - General exponentiation - `compXY`, `annuityXY` - Financial functions - `randNum` - Random number generation diff --git a/src/libsrc/fpu_mf.pla b/src/libsrc/fpu_mf.pla index 3338040..4076b01 100644 --- a/src/libsrc/fpu_mf.pla +++ b/src/libsrc/fpu_mf.pla @@ -11,11 +11,18 @@ // - sin, cos, tan, atan // - ln (natural log), exp (e^x) // -// HARDWARE-ACCELERATED DERIVED OPERATIONS (using trig identities): +// HARDWARE-ACCELERATED DERIVED OPERATIONS (using logarithm/trig identities): // - neg(x) = x * -1 (hardware mul) // - abs(x) = x if x >= 0, else neg(x) // - log2(x) = ln(x) / ln(2) +// - log2(1+x) = ln(1+x) / ln(2) +// - ln(1+x) = ln(1 + x) // - 2^x = e^(x * ln(2)) +// - 2^x - 1 = e^(x*ln(2)) - 1 +// - e^x - 1 = exp(x) - 1 +// - e^(2x) - 1 = exp(2x) - 1 +// - x^n = e^(n * ln(x)) +// - x^y = e^(y * ln(x)) // - asin(x) = atan(x / sqrt(1 - x²)) // - acos(x) = π/2 - asin(x) // - sinh(x) = (e^x - e^-x) / 2 @@ -27,7 +34,7 @@ // // These functions all use hardware acceleration and maintain full precision // by leveraging the MegaFlash FPU's native operations combined with -// mathematical identities. +// mathematical identities and properties of logarithms/exponentials. // include "inc/cmdsys.plh" include "inc/sane.plh" @@ -918,10 +925,39 @@ end def log21X // log2(1+x) = ln(1+x) / ln(2) - // For small x, could optimize, but for now use standard approach - if !saneInit; initSANE; fin - sane:saveZP() - return sane:restoreZP(sane:op1ELEM(FFEXT|FOLOG21X, stackRegs[0])) + word err + byte[t_extended] ln2 + + if !mfAvailable + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op1ELEM(FFEXT|FOLOG21X, stackRegs[0])) + fin + + // ln(2) constant + ln2.0 = $FE + ln2.1 = $3F + ln2.2 = $B1 + ln2.3 = $72 + ln2.4 = $17 + ln2.5 = $F7 + ln2.6 = $D1 + ln2.7 = $CF + ln2.8 = $79 + ln2.9 = $AB + + // Compute ln(1+x) + err = ln1X + if err < 0 + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op1ELEM(FFEXT|FOLOG21X, stackRegs[0])) + fin + + // Divide by ln(2) + memcpy(stackRegs[1], @ln2, t_extended) + swap + return div end def lnX @@ -938,9 +974,41 @@ def lnX end def ln1X - if !saneInit; initSANE; fin - sane:saveZP() - return sane:restoreZP(sane:op1ELEM(FFEXT|FOLN1X, stackRegs[0])) + // ln(1+x) = ln(1 + x) + // Add 1 to x, then compute ln + word err + byte[t_extended] one + + if !mfAvailable + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op1ELEM(FFEXT|FOLN1X, stackRegs[0])) + fin + + // one = 1.0 + one.0 = 0 + one.1 = $3F + one.2 = $80 + one.3 = 0 + one.4 = 0 + one.5 = 0 + one.6 = 0 + one.7 = 0 + one.8 = 0 + one.9 = 0 + + // Compute 1 + x + memcpy(stackRegs[1], @one, t_extended) + err = add + if err < 0 + // Fallback + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op1ELEM(FFEXT|FOLN1X, stackRegs[0])) + fin + + // Compute ln(1 + x) + return execUnaryOp(MF_CMD_FLOG) end def pow2X @@ -976,11 +1044,57 @@ def pow2X end def pow21X - // 2^x - 1 for small x - // For now, fallback to SANE for accuracy near zero - if !saneInit; initSANE; fin - sane:saveZP() - return sane:restoreZP(sane:op1ELEM(FFEXT|FOEXP21X, stackRegs[0])) + // 2^x - 1 = e^(x*ln(2)) - 1 + word err + byte[t_extended] ln2, one + + if !mfAvailable + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op1ELEM(FFEXT|FOEXP21X, stackRegs[0])) + fin + + // ln(2) constant + ln2.0 = $FE + ln2.1 = $3F + ln2.2 = $B1 + ln2.3 = $72 + ln2.4 = $17 + ln2.5 = $F7 + ln2.6 = $D1 + ln2.7 = $CF + ln2.8 = $79 + ln2.9 = $AB + + // one = 1.0 + one.0 = 0 + one.1 = $3F + one.2 = $80 + one.3 = 0 + one.4 = 0 + one.5 = 0 + one.6 = 0 + one.7 = 0 + one.8 = 0 + one.9 = 0 + + // Multiply x by ln(2) + memcpy(stackRegs[1], stackRegs[0], t_extended) + memcpy(stackRegs[0], @ln2, t_extended) + err = execBinaryOp(MF_CMD_FMUL) + if err < 0 + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op1ELEM(FFEXT|FOEXP21X, stackRegs[1])) + fin + + // Compute e^(x*ln2) + err = execUnaryOp(MF_CMD_FEXP) + if err; return err; fin + + // Subtract 1 + memcpy(stackRegs[1], @one, t_extended) + return sub end def powEX @@ -997,27 +1111,150 @@ def powEX end def powE1X - if !saneInit; initSANE; fin - sane:saveZP() - return sane:restoreZP(sane:op1ELEM(FFEXT|FOEXP1X, stackRegs[0])) + // e^x - 1 = exp(x) - 1 + word err + byte[t_extended] one + + if !mfAvailable + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op1ELEM(FFEXT|FOEXP1X, stackRegs[0])) + fin + + // one = 1.0 + one.0 = 0 + one.1 = $3F + one.2 = $80 + one.3 = 0 + one.4 = 0 + one.5 = 0 + one.6 = 0 + one.7 = 0 + one.8 = 0 + one.9 = 0 + + // Compute e^x + err = execUnaryOp(MF_CMD_FEXP) + if err < 0 + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op1ELEM(FFEXT|FOEXP1X, stackRegs[0])) + fin + + // Subtract 1 + memcpy(stackRegs[1], @one, t_extended) + return sub end def powE21X - if !saneInit; initSANE; fin - sane:saveZP() - return sane:restoreZP(sane:op1ELEM(FFEXT|FOEXP21X, stackRegs[0])) + // e^(2x) - 1 = exp(2x) - 1 + word err + byte[t_extended] two, one + + if !mfAvailable + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op1ELEM(FFEXT|FOEXP21X, stackRegs[0])) + fin + + // two = 2.0 + two.0 = 0 + two.1 = $40 + two.2 = $80 + two.3 = 0 + two.4 = 0 + two.5 = 0 + two.6 = 0 + two.7 = 0 + two.8 = 0 + two.9 = 0 + + // one = 1.0 + one.0 = 0 + one.1 = $3F + one.2 = $80 + one.3 = 0 + one.4 = 0 + one.5 = 0 + one.6 = 0 + one.7 = 0 + one.8 = 0 + one.9 = 0 + + // Multiply x by 2 + memcpy(stackRegs[1], stackRegs[0], t_extended) + memcpy(stackRegs[0], @two, t_extended) + err = execBinaryOp(MF_CMD_FMUL) + if err < 0 + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op1ELEM(FFEXT|FOEXP21X, stackRegs[0])) + fin + + // Compute e^(2x) + err = execUnaryOp(MF_CMD_FEXP) + if err; return err; fin + + // Subtract 1 + memcpy(stackRegs[1], @one, t_extended) + return sub end def powXInt(powInt) - if !saneInit; initSANE; fin - sane:saveZP() - return sane:restoreZP(sane:op2ELEM(FFEXT|FOXPWRI, stackRegs[0], powInt)) + // x^n = e^(n * ln(x)) + word err, intVal + + if !mfAvailable + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op2ELEM(FFEXT|FOXPWRI, stackRegs[0], powInt)) + fin + + // Compute ln(x) + err = execUnaryOp(MF_CMD_FLOG) + if err < 0 + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(sane:op2ELEM(FFEXT|FOXPWRI, stackRegs[0], powInt)) + fin + + // Convert integer to extended and multiply + pushInt(powInt) + err = execBinaryOp(MF_CMD_FMUL) + if err; return err; fin + + // Compute e^(n * ln(x)) + return execUnaryOp(MF_CMD_FEXP) end def powXY - if !saneInit; initSANE; fin - sane:saveZP() - return sane:restoreZP(_drop(_swap(sane:op2ELEM(FFEXT|FOXPWRY, stackRegs[0], stackRegs[1])))) + // x^y = e^(y * ln(x)) + word err + + if !mfAvailable + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(_drop(_swap(sane:op2ELEM(FFEXT|FOXPWRY, stackRegs[0], stackRegs[1])))) + fin + + // Stack has: y, x (top to bottom) + // We need: ln(x), then multiply by y, then exp + + // Compute ln(x) - x is on top of stack + err = execUnaryOp(MF_CMD_FLOG) + if err < 0 + if !saneInit; initSANE; fin + sane:saveZP() + return sane:restoreZP(_drop(_swap(sane:op2ELEM(FFEXT|FOXPWRY, stackRegs[1], stackRegs[0])))) + fin + + // Now stack has: y, ln(x) + // Multiply y * ln(x) + err = execBinaryOp(MF_CMD_FMUL) + if err; return _drop(err); fin + + // Compute e^(y * ln(x)) + return _drop(execUnaryOp(MF_CMD_FEXP)) end //==============================================================================