From a3df482acd56ca20d49f33e1c63c7bb2e97bbd64 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Sat, 26 Nov 2022 14:14:37 -0500 Subject: [PATCH 1/2] =?UTF-8?q?Use=20non-inlined=20recursion=20for=20?= =?UTF-8?q?=E2=80=A2math=20functions?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/arithm.c | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/builtins/arithm.c b/src/builtins/arithm.c index 0a3c8a4e..581a902f 100644 --- a/src/builtins/arithm.c +++ b/src/builtins/arithm.c @@ -94,15 +94,19 @@ f64 logfact(f64 x) { return lgamma(x+1); } #define P1(N) { if(isArr(x)) { SLOW1("arithm " #N, x); return arith_recm(N##_c1, x); } } B pow_c1(B t, B x) { if (isF64(x)) return m_f64( exp(x.f)); P1( pow); thrM("⋆: Getting exp of non-number"); } B log_c1(B t, B x) { if (isF64(x)) return m_f64( log(x.f)); P1( log); thrM("⋆⁼: Getting log of non-number"); } -#define MATH(n,N) \ - B n##_c1(B t, B x) { if (isF64(x)) return m_f64(n(x.f)); P1(n); thrM("•math." #N ": Argument contained non-number"); } +#undef P1 +static NOINLINE B arith_recm_slow(f64 (*fn)(f64), BB2B rec, B x, char* s) { + if (isF64(x)) return m_f64(fn(x.f)); + if(isArr(x)) return arith_recm(rec, x); + thrF("•math.%S: Argument contained non-number", s); +} +#define MATH(n,N) B n##_c1(B t, B x) { return arith_recm_slow(n, n##_c1, x, #N); } MATH(cbrt,Cbrt) MATH(log2,Log2) MATH(log10,Log10) MATH(log1p,Log1p) MATH(expm1,Expm1) MATH(fact,Fact) MATH(logfact,LogFact) MATH(erf,Erf) MATH(erfc,ErfC) #define TRIG(n,N) MATH(n,N) MATH(a##n,A##n) MATH(n##h,N##h) MATH(a##n##h,A##n##h) TRIG(sin,Sin) TRIG(cos,Cos) TRIG(tan,Tan) #undef TRIG #undef MATH -#undef P1 B lt_c1(B t, B x) { return m_atomUnit(x); } B eq_c1(B t, B x) { if (isAtm(x)) { decA(x); return m_i32(0); } B r = m_i32(RNK(x)); decG(x); return r; } From d911eafd31a52f552ac475f62e39a9b2e6bbe661 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Sat, 26 Nov 2022 14:17:27 -0500 Subject: [PATCH 2/2] =?UTF-8?q?Inverses=20for=20=E2=80=A2math.Fact=20and?= =?UTF-8?q?=20=E2=80=A2math.LogFact?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/arithm.c | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/src/builtins/arithm.c b/src/builtins/arithm.c index 581a902f..41069110 100644 --- a/src/builtins/arithm.c +++ b/src/builtins/arithm.c @@ -90,6 +90,18 @@ GC1f(root, sqrt(xv), "√: Getting square root of non-number") f64 fact(f64 x) { return tgamma(x+1); } f64 logfact(f64 x) { return lgamma(x+1); } +NOINLINE f64 logfact_inv(f64 y) { + if (!(y >= -0.12)) thrM("⁼: required factorial result too small"); + if (y == INFINITY) return y; + f64 x = 4; + PLAINLOOP for (usz i = 0; i < 20; i++) { + f64 x0 = x; + x += (y - logfact(x)) / log(0.52 + x); + if (x == x0) break; + } + return x; +} +f64 fact_inv(f64 y) { return logfact_inv(log(y)); } #define P1(N) { if(isArr(x)) { SLOW1("arithm " #N, x); return arith_recm(N##_c1, x); } } B pow_c1(B t, B x) { if (isF64(x)) return m_f64( exp(x.f)); P1( pow); thrM("⋆: Getting exp of non-number"); } @@ -102,7 +114,7 @@ static NOINLINE B arith_recm_slow(f64 (*fn)(f64), BB2B rec, B x, char* s) { } #define MATH(n,N) B n##_c1(B t, B x) { return arith_recm_slow(n, n##_c1, x, #N); } MATH(cbrt,Cbrt) MATH(log2,Log2) MATH(log10,Log10) MATH(log1p,Log1p) MATH(expm1,Expm1) -MATH(fact,Fact) MATH(logfact,LogFact) MATH(erf,Erf) MATH(erfc,ErfC) +MATH(fact,Fact) MATH(logfact,LogFact) MATH(logfact_inv,LogFact⁼) MATH(fact_inv,Fact⁼) MATH(erf,Erf) MATH(erfc,ErfC) #define TRIG(n,N) MATH(n,N) MATH(a##n,A##n) MATH(n##h,N##h) MATH(a##n##h,A##n##h) TRIG(sin,Sin) TRIG(cos,Cos) TRIG(tan,Tan) #undef TRIG @@ -142,5 +154,7 @@ void arith_init() { INVERSE_PAIR(cosh, acosh) INVERSE_PAIR(tanh, atanh) INVERSE_PAIR(expm1, log1p) + c(BFn,bi_fact)->im = fact_inv_c1; + c(BFn,bi_logfact)->im = logfact_inv_c1; #undef INVERSE_PAIR }