From 51840bc853510c3d39a0db797f39cf6fde8ff091 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Tue, 22 Nov 2022 16:47:53 -0500 Subject: [PATCH 01/11] =?UTF-8?q?=E2=80=A2math.GCD=20and=20=E2=80=A2math.L?= =?UTF-8?q?CM=20for=20u64=20arguments?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins.h | 2 +- src/builtins/arithd.c | 38 ++++++++++++++++++++++++++++++++++++++ src/builtins/arithm.c | 4 ++-- src/builtins/sysfn.c | 2 +- 4 files changed, 42 insertions(+), 4 deletions(-) diff --git a/src/builtins.h b/src/builtins.h index 5d79619f..6f881265 100644 --- a/src/builtins.h +++ b/src/builtins.h @@ -17,7 +17,7 @@ /*internal.c*/M(itype,"•internal.Type") M(elType,"•internal.ElType") M(refc,"•internal.Refc") M(isPure,"•internal.IsPure") A(info,"•internal.Info") M(heapDump,"•internal.HeapDump") \ /*internal.c*/M(squeeze,"•internal.Squeeze") M(deepSqueeze,"•internal.DeepSqueeze") D(eequal,"•internal.EEqual") A(internalTemp,"•internal.Temp") \ /*internal.c*/D(variation,"•internal.Variation") A(listVariations,"•internal.ListVariations") M(clearRefs,"•internal.ClearRefs") M(unshare,"•internal.Unshare") \ -/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") +/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(gcd,"•math.GCD") D(lcm,"•math.LCM") #define FOR_PM1(A,M,D) \ /*md1.c*/A(tbl,"⌜") A(each,"¨") A(fold,"´") A(scan,"`") A(const,"˙") A(swap,"˜") A(cell,"˘") A(insert,"˝") \ diff --git a/src/builtins/arithd.c b/src/builtins/arithd.c index 211d9955..82258fb2 100644 --- a/src/builtins/arithd.c +++ b/src/builtins/arithd.c @@ -366,4 +366,42 @@ B atan2_c2(B t, B w, B x) { thrM("•math.Atan2: Unexpected argument types"); } +static u64 gcd_u64(u64 a, u64 b) { + if (a == 0) return b; + if (b == 0) return a; + u8 az = CTZ(a); + u8 bz = CTZ(b); + u8 sh = az>= bz; + while (a > 0) { + a >>= az; + u64 d = b - a; + az = CTZ(d); + b = b Date: Tue, 22 Nov 2022 17:27:01 -0500 Subject: [PATCH 02/11] Hyperbolic trig --- src/builtins.h | 2 +- src/builtins/arithm.c | 22 ++++++++++++++-------- src/builtins/sysfn.c | 2 +- 3 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/builtins.h b/src/builtins.h index 6f881265..e27d4550 100644 --- a/src/builtins.h +++ b/src/builtins.h @@ -17,7 +17,7 @@ /*internal.c*/M(itype,"•internal.Type") M(elType,"•internal.ElType") M(refc,"•internal.Refc") M(isPure,"•internal.IsPure") A(info,"•internal.Info") M(heapDump,"•internal.HeapDump") \ /*internal.c*/M(squeeze,"•internal.Squeeze") M(deepSqueeze,"•internal.DeepSqueeze") D(eequal,"•internal.EEqual") A(internalTemp,"•internal.Temp") \ /*internal.c*/D(variation,"•internal.Variation") A(listVariations,"•internal.ListVariations") M(clearRefs,"•internal.ClearRefs") M(unshare,"•internal.Unshare") \ -/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(gcd,"•math.GCD") D(lcm,"•math.LCM") +/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") D(gcd,"•math.GCD") D(lcm,"•math.LCM") #define FOR_PM1(A,M,D) \ /*md1.c*/A(tbl,"⌜") A(each,"¨") A(fold,"´") A(scan,"`") A(const,"˙") A(swap,"˜") A(cell,"˘") A(insert,"˝") \ diff --git a/src/builtins/arithm.c b/src/builtins/arithm.c index 61c4f0a5..1fd45f74 100644 --- a/src/builtins/arithm.c +++ b/src/builtins/arithm.c @@ -92,12 +92,12 @@ GC1f(root, sqrt(xv), "√: Getting square root of non-number") #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"); } -B sin_c1(B t, B x) { if (isF64(x)) return m_f64( sin(x.f)); P1( sin); thrM("•math.Sin: Argument contained non-number"); } -B cos_c1(B t, B x) { if (isF64(x)) return m_f64( cos(x.f)); P1( cos); thrM("•math.Cos: Argument contained non-number"); } -B tan_c1(B t, B x) { if (isF64(x)) return m_f64( tan(x.f)); P1( tan); thrM("•math.Tan: Argument contained non-number"); } -B asin_c1(B t, B x) { if (isF64(x)) return m_f64( asin(x.f)); P1( asin); thrM("•math.Asin: Argument contained non-number"); } -B acos_c1(B t, B x) { if (isF64(x)) return m_f64( acos(x.f)); P1( acos); thrM("•math.Acos: Argument contained non-number"); } -B atan_c1(B t, B x) { if (isF64(x)) return m_f64( atan(x.f)); P1( atan); thrM("•math.Atan: Argument contained 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"); } +#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 MATH +#undef TRIG #undef P1 B lt_c1(B t, B x) { return m_atomUnit(x); } @@ -109,8 +109,8 @@ static B mathNS; B getMathNS() { if (mathNS.u == 0) { #define F(X) inc(bi_##X), - Body* d = m_nnsDesc("sin","cos","tan","asin","acos","atan","atan2","gcd","lcm"); - mathNS = m_nns(d, F(sin)F(cos)F(tan)F(asin)F(acos)F(atan)F(atan2)F(gcd)F(lcm)); + Body* d = m_nnsDesc("sin","cos","tan","asin","acos","atan","atan2","sinh","cosh","tanh","asinh","acosh","atanh","gcd","lcm"); + mathNS = m_nns(d, F(sin)F(cos)F(tan)F(asin)F(acos)F(atan)F(atan2)F(sinh)F(cosh)F(tanh)F(asinh)F(acosh)F(atanh)F(gcd)F(lcm)); #undef F gc_add(mathNS); } @@ -130,4 +130,10 @@ void arith_init() { c(BFn,bi_asin)->im = sin_c1; c(BFn,bi_acos)->im = cos_c1; c(BFn,bi_atan)->im = tan_c1; + c(BFn,bi_sinh)->im = asinh_c1; + c(BFn,bi_cosh)->im = acosh_c1; + c(BFn,bi_tanh)->im = atanh_c1; + c(BFn,bi_asinh)->im = sinh_c1; + c(BFn,bi_acosh)->im = cosh_c1; + c(BFn,bi_atanh)->im = tanh_c1; } diff --git a/src/builtins/sysfn.c b/src/builtins/sysfn.c index 1c2a0359..50ee562b 100644 --- a/src/builtins/sysfn.c +++ b/src/builtins/sysfn.c @@ -1551,7 +1551,7 @@ u32* dsv_text[] = { U"•file.MapBytes",U"•file.Modified",U"•file.Name",U"•file.Parent",U"•file.Remove",U"•file.Rename",U"•file.Size",U"•file.Type", U"•internal.ClearRefs",U"•internal.DeepSqueeze",U"•internal.EEqual",U"•internal.ElType",U"•internal.HeapDump",U"•internal.Info",U"•internal.IsPure",U"•internal.ListVariations",U"•internal.Refc",U"•internal.Squeeze",U"•internal.Temp",U"•internal.Type",U"•internal.Unshare",U"•internal.Variation", - U"•math.Acos",U"•math.Asin",U"•math.Atan",U"•math.Atan2",U"•math.Cos",U"•math.GCD",U"•math.LCM",U"•math.Sin",U"•math.Tan", + U"•math.Acos",U"•math.Acosh",U"•math.Asin",U"•math.Asinh",U"•math.Atan",U"•math.Atan2",U"•math.Atanh",U"•math.Cos",U"•math.Cosh",U"•math.GCD",U"•math.LCM",U"•math.Sin",U"•math.Sinh",U"•math.Tan",U"•math.Tanh", U"•rand.Deal",U"•rand.Range",U"•rand.Subset", U"•term.CharB",U"•term.CharN",U"•term.ErrRaw",U"•term.Flush",U"•term.OutRaw",U"•term.RawMode", NULL From 58cd6f0f4e54ff27fe74efda4a81574dbafb2ccc Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Tue, 22 Nov 2022 18:03:59 -0500 Subject: [PATCH 03/11] Other useful math: log2 log10 log1p expm1 hypot --- src/builtins.h | 2 +- src/builtins/arithd.c | 13 ++++++++----- src/builtins/arithm.c | 9 ++++++--- src/builtins/sysfn.c | 2 +- 4 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/builtins.h b/src/builtins.h index e27d4550..a12d8ebd 100644 --- a/src/builtins.h +++ b/src/builtins.h @@ -17,7 +17,7 @@ /*internal.c*/M(itype,"•internal.Type") M(elType,"•internal.ElType") M(refc,"•internal.Refc") M(isPure,"•internal.IsPure") A(info,"•internal.Info") M(heapDump,"•internal.HeapDump") \ /*internal.c*/M(squeeze,"•internal.Squeeze") M(deepSqueeze,"•internal.DeepSqueeze") D(eequal,"•internal.EEqual") A(internalTemp,"•internal.Temp") \ /*internal.c*/D(variation,"•internal.Variation") A(listVariations,"•internal.ListVariations") M(clearRefs,"•internal.ClearRefs") M(unshare,"•internal.Unshare") \ -/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") D(gcd,"•math.GCD") D(lcm,"•math.LCM") +/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(hypot,"•math.Hypot") M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") M(cbrt,"•math.Cbrt") M(log2,"•math.Log2") M(log10,"•math.Log10") M(log1p,"•math.Log1p") M(expm1,"•math.Expm1") D(gcd,"•math.GCD") D(lcm,"•math.LCM") #define FOR_PM1(A,M,D) \ /*md1.c*/A(tbl,"⌜") A(each,"¨") A(fold,"´") A(scan,"`") A(const,"˙") A(swap,"˜") A(cell,"˘") A(insert,"˝") \ diff --git a/src/builtins/arithd.c b/src/builtins/arithd.c index 82258fb2..56b2d9ad 100644 --- a/src/builtins/arithd.c +++ b/src/builtins/arithd.c @@ -360,11 +360,14 @@ AR_F_SCALAR("|", stile, pfmod(x.f, w.f)) AR_F_SCALAR("⋆⁼",log , log(x.f)/log(w.f)) #undef AR_F_SCALAR -B atan2_c2(B t, B w, B x) { - if (isNum(w) && isNum(x)) return m_f64(atan2(x.f, w.f)); - P2(atan2) - thrM("•math.Atan2: Unexpected argument types"); -} +#define MATH(n,N) \ + B n##_c2(B t, B w, B x) { \ + if (isNum(w) && isNum(x)) return m_f64(n(x.f, w.f)); \ + P2(n) \ + thrM("•math." #N ": Unexpected argument types"); \ + } +MATH(atan2,Atan2) MATH(hypot,Hypot) +#undef MATH static u64 gcd_u64(u64 a, u64 b) { if (a == 0) return b; diff --git a/src/builtins/arithm.c b/src/builtins/arithm.c index 1fd45f74..881e9c7b 100644 --- a/src/builtins/arithm.c +++ b/src/builtins/arithm.c @@ -94,10 +94,11 @@ B pow_c1(B t, B x) { if (isF64(x)) return m_f64( exp(x.f)); P1( pow); thrM(" 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"); } +MATH(cbrt,Cbrt) MATH(log2,Log2) MATH(log10,Log10) MATH(log1p,Log1p) MATH(expm1,Expm1) #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 MATH #undef TRIG +#undef MATH #undef P1 B lt_c1(B t, B x) { return m_atomUnit(x); } @@ -109,8 +110,8 @@ static B mathNS; B getMathNS() { if (mathNS.u == 0) { #define F(X) inc(bi_##X), - Body* d = m_nnsDesc("sin","cos","tan","asin","acos","atan","atan2","sinh","cosh","tanh","asinh","acosh","atanh","gcd","lcm"); - mathNS = m_nns(d, F(sin)F(cos)F(tan)F(asin)F(acos)F(atan)F(atan2)F(sinh)F(cosh)F(tanh)F(asinh)F(acosh)F(atanh)F(gcd)F(lcm)); + Body* d = m_nnsDesc("sin","cos","tan","asin","acos","atan","atan2","sinh","cosh","tanh","asinh","acosh","atanh","cbrt","log2","log10","log1p","expm1","hypot","gcd","lcm"); + mathNS = m_nns(d, F(sin)F(cos)F(tan)F(asin)F(acos)F(atan)F(atan2)F(sinh)F(cosh)F(tanh)F(asinh)F(acosh)F(atanh)F(cbrt)F(log2)F(log10)F(log1p)F(expm1)F(hypot)F(gcd)F(lcm)); #undef F gc_add(mathNS); } @@ -136,4 +137,6 @@ void arith_init() { c(BFn,bi_asinh)->im = sinh_c1; c(BFn,bi_acosh)->im = cosh_c1; c(BFn,bi_atanh)->im = tanh_c1; + c(BFn,bi_expm1)->im = log1p_c1; + c(BFn,bi_log1p)->im = expm1_c1; } diff --git a/src/builtins/sysfn.c b/src/builtins/sysfn.c index 50ee562b..db507622 100644 --- a/src/builtins/sysfn.c +++ b/src/builtins/sysfn.c @@ -1551,7 +1551,7 @@ u32* dsv_text[] = { U"•file.MapBytes",U"•file.Modified",U"•file.Name",U"•file.Parent",U"•file.Remove",U"•file.Rename",U"•file.Size",U"•file.Type", U"•internal.ClearRefs",U"•internal.DeepSqueeze",U"•internal.EEqual",U"•internal.ElType",U"•internal.HeapDump",U"•internal.Info",U"•internal.IsPure",U"•internal.ListVariations",U"•internal.Refc",U"•internal.Squeeze",U"•internal.Temp",U"•internal.Type",U"•internal.Unshare",U"•internal.Variation", - U"•math.Acos",U"•math.Acosh",U"•math.Asin",U"•math.Asinh",U"•math.Atan",U"•math.Atan2",U"•math.Atanh",U"•math.Cos",U"•math.Cosh",U"•math.GCD",U"•math.LCM",U"•math.Sin",U"•math.Sinh",U"•math.Tan",U"•math.Tanh", + U"•math.Acos",U"•math.Acosh",U"•math.Asin",U"•math.Asinh",U"•math.Atan",U"•math.Atan2",U"•math.Atanh",U"•math.Cbrt",U"•math.Cos",U"•math.Cosh",U"•math.Expm1",U"•math.GCD",U"•math.Hypot",U"•math.LCM",U"•math.Log10",U"•math.Log1p",U"•math.Log2",U"•math.Sin",U"•math.Sinh",U"•math.Tan",U"•math.Tanh", U"•rand.Deal",U"•rand.Range",U"•rand.Subset", U"•term.CharB",U"•term.CharN",U"•term.ErrRaw",U"•term.Flush",U"•term.OutRaw",U"•term.RawMode", NULL From f6109b39858a46e1741ee42cf167f2d1db5ea43e Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Tue, 22 Nov 2022 19:20:13 -0500 Subject: [PATCH 04/11] Handle no-op and boolean cases of monadic arithmetic quickly --- src/builtins/arithm.c | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/builtins/arithm.c b/src/builtins/arithm.c index 881e9c7b..f4751ea8 100644 --- a/src/builtins/arithm.c +++ b/src/builtins/arithm.c @@ -22,11 +22,12 @@ B bit_negate(B x) { // consumes return r; } -#define GC1i(SYMB,NAME,FEXPR,IBAD,IEXPR,BX,SQF) B NAME##_c1(B t, B x) { \ +#define GC1i(SYMB,NAME,FEXPR,IBAD,IEXPR,SQF,TMIN,RMIN) B NAME##_c1(B t, B x) { \ if (isF64(x)) { f64 v = x.f; return m_f64(FEXPR); } \ if (RARE(!isArr(x))) thrM(SYMB ": Expected argument to be a number"); \ u8 xe = TI(x,elType); \ - i64 sz = IA(x); BX \ + if (xe<=TMIN) return RMIN; \ + i64 sz = IA(x); \ if (xe==el_i8) { i8 MAX=I8_MAX; i8 MIN=I8_MIN; i8* xp=i8any_ptr(x); i8* rp; B r=m_i8arrc(&rp,x); \ for (i64 i = 0; i < sz; i++) { i8 v = xp[i]; if (RARE(IBAD)) { decG(r); goto base; } rp[i] = IEXPR; } \ decG(x); (void)MIN;(void)MAX; return r; \ @@ -56,14 +57,12 @@ B add_c1(B t, B x) { return x; } -GC1i("-", sub, -v, v== MIN, -v, {}, 0) // change icond to v==-v to support ¯0 (TODO that won't work for i8/i16) -GC1i("|", stile, fabs(v), v== MIN, v<0?-v:v,{}, 0) -GC1i("⌊", floor, floor(v), 0, v, {}, 1) -GC1i("⌈", ceil, ceil(v), 0, v, {}, 1) -GC1i("×", mul, v==0?0:v>0?1:-1, 0, v==0?0:v>0?1:-1,{}, 1) -GC1i("¬", not, 1-v, v<=-MAX, 1-v, { - if(xe==el_bit) return bit_negate(x); -}, 0) +GC1i("-", sub, -v, v== MIN, -v, 0, el_bit, bit_sel(x,m_f64(0),m_f64(-1))) // change icond to v==-v to support ¯0 (TODO that won't work for i8/i16) +GC1i("|", stile, fabs(v), v== MIN, v<0?-v:v,0, el_bit, x) +GC1i("⌊", floor, floor(v), 0, v, 1, el_i32, x) +GC1i("⌈", ceil, ceil(v), 0, v, 1, el_i32, x) +GC1i("×", mul, v==0?0:v>0?1:-1, 0,v==0?0:v>0?1:-1,1, el_bit, x) +GC1i("¬", not, 1-v, v<=-MAX, 1-v, 0, el_bit, bit_negate(x)) #define GC1f(N, F, MSG) B N##_c1(B t, B x) { \ if (isF64(x)) { f64 xv=o2fG(x); return m_f64(F); } \ From 6eb21bd68cf69c62c560fc95c9e333131d0f65c8 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Tue, 22 Nov 2022 20:51:11 -0500 Subject: [PATCH 05/11] Factorial and combinations functions --- src/builtins.h | 2 +- src/builtins/arithd.c | 36 +++++++++++++++++++++++++++++++++++- src/builtins/arithm.c | 8 +++++--- src/builtins/sysfn.c | 2 +- 4 files changed, 42 insertions(+), 6 deletions(-) diff --git a/src/builtins.h b/src/builtins.h index a12d8ebd..c7b02f7e 100644 --- a/src/builtins.h +++ b/src/builtins.h @@ -17,7 +17,7 @@ /*internal.c*/M(itype,"•internal.Type") M(elType,"•internal.ElType") M(refc,"•internal.Refc") M(isPure,"•internal.IsPure") A(info,"•internal.Info") M(heapDump,"•internal.HeapDump") \ /*internal.c*/M(squeeze,"•internal.Squeeze") M(deepSqueeze,"•internal.DeepSqueeze") D(eequal,"•internal.EEqual") A(internalTemp,"•internal.Temp") \ /*internal.c*/D(variation,"•internal.Variation") A(listVariations,"•internal.ListVariations") M(clearRefs,"•internal.ClearRefs") M(unshare,"•internal.Unshare") \ -/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(hypot,"•math.Hypot") M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") M(cbrt,"•math.Cbrt") M(log2,"•math.Log2") M(log10,"•math.Log10") M(log1p,"•math.Log1p") M(expm1,"•math.Expm1") D(gcd,"•math.GCD") D(lcm,"•math.LCM") +/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(hypot,"•math.Hypot") M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") M(cbrt,"•math.Cbrt") M(log2,"•math.Log2") M(log10,"•math.Log10") M(log1p,"•math.Log1p") M(expm1,"•math.Expm1") M(fact,"•math.Fact") D(comb,"•math.Comb") D(gcd,"•math.GCD") D(lcm,"•math.LCM") #define FOR_PM1(A,M,D) \ /*md1.c*/A(tbl,"⌜") A(each,"¨") A(fold,"´") A(scan,"`") A(const,"˙") A(swap,"˜") A(cell,"˘") A(insert,"˝") \ diff --git a/src/builtins/arithd.c b/src/builtins/arithd.c index 56b2d9ad..ae81f553 100644 --- a/src/builtins/arithd.c +++ b/src/builtins/arithd.c @@ -360,13 +360,47 @@ AR_F_SCALAR("|", stile, pfmod(x.f, w.f)) AR_F_SCALAR("⋆⁼",log , log(x.f)/log(w.f)) #undef AR_F_SCALAR +static f64 comb_nat(f64 k, f64 n, f64 j) { + if (j < k) k = j; + if (k > 514) return INFINITY; + f64 p = 1; + for (usz i=0; i<(usz)k; i++) { + p*= (n-i) / (k-i); + if (p == INFINITY) return p; + } + return round(p); +} +static f64 comb(f64 k, f64 n) { // n choose k + f64 j = n - k; // j+k == n + bool jint = j == round(j); + if (k == round(k)) { + if (jint) { + if (n >= 0) { + if (!(k>=0 && j>=0)) return 0; // Negative phrasing to catch NaN + return comb_nat(k, n, j); + } else { + if (k<0) { + if (j<0) return 0; + f64 t=k; k=j; j=t; // Swap so k is non-negative + } + f64 r = comb_nat(k, -1-j, -1-n); + return k<(1ull<<53) && ((i64)k&1)? -r : r; + } + } + if (k < 0) return 0; + } else if (jint) { + if (j < 0) return 0; + } + return exp(lgamma(n+1) - lgamma(k+1) - lgamma(j+1)); +} + #define MATH(n,N) \ B n##_c2(B t, B w, B x) { \ if (isNum(w) && isNum(x)) return m_f64(n(x.f, w.f)); \ P2(n) \ thrM("•math." #N ": Unexpected argument types"); \ } -MATH(atan2,Atan2) MATH(hypot,Hypot) +MATH(atan2,Atan2) MATH(hypot,Hypot) MATH(comb,Comb) #undef MATH static u64 gcd_u64(u64 a, u64 b) { diff --git a/src/builtins/arithm.c b/src/builtins/arithm.c index f4751ea8..f3d90315 100644 --- a/src/builtins/arithm.c +++ b/src/builtins/arithm.c @@ -88,12 +88,14 @@ GC1f( div, 1/xv, "÷: Getting reciprocal of non-number") GC1f(root, sqrt(xv), "√: Getting square root of non-number") #undef GC1f +f64 fact(f64 x) { return tgamma(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"); } -MATH(cbrt,Cbrt) MATH(log2,Log2) MATH(log10,Log10) MATH(log1p,Log1p) MATH(expm1,Expm1) +MATH(cbrt,Cbrt) MATH(log2,Log2) MATH(log10,Log10) MATH(log1p,Log1p) MATH(expm1,Expm1) MATH(fact,Fact) #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 @@ -109,8 +111,8 @@ static B mathNS; B getMathNS() { if (mathNS.u == 0) { #define F(X) inc(bi_##X), - Body* d = m_nnsDesc("sin","cos","tan","asin","acos","atan","atan2","sinh","cosh","tanh","asinh","acosh","atanh","cbrt","log2","log10","log1p","expm1","hypot","gcd","lcm"); - mathNS = m_nns(d, F(sin)F(cos)F(tan)F(asin)F(acos)F(atan)F(atan2)F(sinh)F(cosh)F(tanh)F(asinh)F(acosh)F(atanh)F(cbrt)F(log2)F(log10)F(log1p)F(expm1)F(hypot)F(gcd)F(lcm)); + Body* d = m_nnsDesc("sin","cos","tan","asin","acos","atan","atan2","sinh","cosh","tanh","asinh","acosh","atanh","cbrt","log2","log10","log1p","expm1","hypot","fact","comb","gcd","lcm"); + mathNS = m_nns(d, F(sin)F(cos)F(tan)F(asin)F(acos)F(atan)F(atan2)F(sinh)F(cosh)F(tanh)F(asinh)F(acosh)F(atanh)F(cbrt)F(log2)F(log10)F(log1p)F(expm1)F(hypot)F(fact)F(comb)F(gcd)F(lcm)); #undef F gc_add(mathNS); } diff --git a/src/builtins/sysfn.c b/src/builtins/sysfn.c index db507622..17330017 100644 --- a/src/builtins/sysfn.c +++ b/src/builtins/sysfn.c @@ -1551,7 +1551,7 @@ u32* dsv_text[] = { U"•file.MapBytes",U"•file.Modified",U"•file.Name",U"•file.Parent",U"•file.Remove",U"•file.Rename",U"•file.Size",U"•file.Type", U"•internal.ClearRefs",U"•internal.DeepSqueeze",U"•internal.EEqual",U"•internal.ElType",U"•internal.HeapDump",U"•internal.Info",U"•internal.IsPure",U"•internal.ListVariations",U"•internal.Refc",U"•internal.Squeeze",U"•internal.Temp",U"•internal.Type",U"•internal.Unshare",U"•internal.Variation", - U"•math.Acos",U"•math.Acosh",U"•math.Asin",U"•math.Asinh",U"•math.Atan",U"•math.Atan2",U"•math.Atanh",U"•math.Cbrt",U"•math.Cos",U"•math.Cosh",U"•math.Expm1",U"•math.GCD",U"•math.Hypot",U"•math.LCM",U"•math.Log10",U"•math.Log1p",U"•math.Log2",U"•math.Sin",U"•math.Sinh",U"•math.Tan",U"•math.Tanh", + U"•math.Acos",U"•math.Acosh",U"•math.Asin",U"•math.Asinh",U"•math.Atan",U"•math.Atan2",U"•math.Atanh",U"•math.Cbrt",U"•math.Comb",U"•math.Cos",U"•math.Cosh",U"•math.Expm1",U"•math.Fact",U"•math.GCD",U"•math.Hypot",U"•math.LCM",U"•math.Log10",U"•math.Log1p",U"•math.Log2",U"•math.Sin",U"•math.Sinh",U"•math.Tan",U"•math.Tanh", U"•rand.Deal",U"•rand.Range",U"•rand.Subset", U"•term.CharB",U"•term.CharN",U"•term.ErrRaw",U"•term.Flush",U"•term.OutRaw",U"•term.RawMode", NULL From 9a1268ea3a1fa28e084594ca1ba782a609e1c1bd Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Wed, 23 Nov 2022 09:49:17 -0500 Subject: [PATCH 06/11] =?UTF-8?q?Simplify=20=E2=80=A2math.Comb=20implement?= =?UTF-8?q?ation?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins/arithd.c | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/builtins/arithd.c b/src/builtins/arithd.c index ae81f553..8f0b6ad1 100644 --- a/src/builtins/arithd.c +++ b/src/builtins/arithd.c @@ -360,8 +360,8 @@ AR_F_SCALAR("|", stile, pfmod(x.f, w.f)) AR_F_SCALAR("⋆⁼",log , log(x.f)/log(w.f)) #undef AR_F_SCALAR -static f64 comb_nat(f64 k, f64 n, f64 j) { - if (j < k) k = j; +static f64 comb_nat(f64 k, f64 n) { + assert(k>=0 && n>=2*k); if (k > 514) return INFINITY; f64 p = 1; for (usz i=0; i<(usz)k; i++) { @@ -375,15 +375,13 @@ static f64 comb(f64 k, f64 n) { // n choose k bool jint = j == round(j); if (k == round(k)) { if (jint) { + if (k= 0) { - if (!(k>=0 && j>=0)) return 0; // Negative phrasing to catch NaN - return comb_nat(k, n, j); + return j<0? 0 : comb_nat(j, n); } else { - if (k<0) { - if (j<0) return 0; - f64 t=k; k=j; j=t; // Swap so k is non-negative - } - f64 r = comb_nat(k, -1-j, -1-n); + if (k<0) return 0; + f64 l = -1-n; // l+k == -1-j + f64 r = comb_nat(k Date: Wed, 23 Nov 2022 11:34:33 -0500 Subject: [PATCH 07/11] =?UTF-8?q?=E2=80=A2math.Erf,=20=E2=80=A2math.ErfC,?= =?UTF-8?q?=20and=20=E2=80=A2math.LogFact?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins.h | 2 +- src/builtins/arithm.c | 37 ++++++++++++++++++------------------- src/builtins/sysfn.c | 2 +- 3 files changed, 20 insertions(+), 21 deletions(-) diff --git a/src/builtins.h b/src/builtins.h index c7b02f7e..77e8a1bb 100644 --- a/src/builtins.h +++ b/src/builtins.h @@ -17,7 +17,7 @@ /*internal.c*/M(itype,"•internal.Type") M(elType,"•internal.ElType") M(refc,"•internal.Refc") M(isPure,"•internal.IsPure") A(info,"•internal.Info") M(heapDump,"•internal.HeapDump") \ /*internal.c*/M(squeeze,"•internal.Squeeze") M(deepSqueeze,"•internal.DeepSqueeze") D(eequal,"•internal.EEqual") A(internalTemp,"•internal.Temp") \ /*internal.c*/D(variation,"•internal.Variation") A(listVariations,"•internal.ListVariations") M(clearRefs,"•internal.ClearRefs") M(unshare,"•internal.Unshare") \ -/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(hypot,"•math.Hypot") M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") M(cbrt,"•math.Cbrt") M(log2,"•math.Log2") M(log10,"•math.Log10") M(log1p,"•math.Log1p") M(expm1,"•math.Expm1") M(fact,"•math.Fact") D(comb,"•math.Comb") D(gcd,"•math.GCD") D(lcm,"•math.LCM") +/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(hypot,"•math.Hypot") M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") M(cbrt,"•math.Cbrt") M(log2,"•math.Log2") M(log10,"•math.Log10") M(log1p,"•math.Log1p") M(expm1,"•math.Expm1") M(fact,"•math.Fact") D(comb,"•math.Comb") M(logfact,"•math.LogFact") M(erf,"•math.Erf") M(erfc,"•math.ErfC") D(gcd,"•math.GCD") D(lcm,"•math.LCM") #define FOR_PM1(A,M,D) \ /*md1.c*/A(tbl,"⌜") A(each,"¨") A(fold,"´") A(scan,"`") A(const,"˙") A(swap,"˜") A(cell,"˘") A(insert,"˝") \ diff --git a/src/builtins/arithm.c b/src/builtins/arithm.c index f3d90315..7e6b40f4 100644 --- a/src/builtins/arithm.c +++ b/src/builtins/arithm.c @@ -88,14 +88,16 @@ GC1f( div, 1/xv, "÷: Getting reciprocal of non-number") GC1f(root, sqrt(xv), "√: Getting square root of non-number") #undef GC1f -f64 fact(f64 x) { return tgamma(x+1); } +f64 fact(f64 x) { return tgamma(x+1); } +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"); } -MATH(cbrt,Cbrt) MATH(log2,Log2) MATH(log10,Log10) MATH(log1p,Log1p) MATH(expm1,Expm1) MATH(fact,Fact) +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 @@ -111,8 +113,8 @@ static B mathNS; B getMathNS() { if (mathNS.u == 0) { #define F(X) inc(bi_##X), - Body* d = m_nnsDesc("sin","cos","tan","asin","acos","atan","atan2","sinh","cosh","tanh","asinh","acosh","atanh","cbrt","log2","log10","log1p","expm1","hypot","fact","comb","gcd","lcm"); - mathNS = m_nns(d, F(sin)F(cos)F(tan)F(asin)F(acos)F(atan)F(atan2)F(sinh)F(cosh)F(tanh)F(asinh)F(acosh)F(atanh)F(cbrt)F(log2)F(log10)F(log1p)F(expm1)F(hypot)F(fact)F(comb)F(gcd)F(lcm)); + Body* d = m_nnsDesc("sin","cos","tan","asin","acos","atan","atan2","sinh","cosh","tanh","asinh","acosh","atanh","cbrt","log2","log10","log1p","expm1","hypot","fact","logfact","erf","erfc","comb","gcd","lcm"); + mathNS = m_nns(d, F(sin)F(cos)F(tan)F(asin)F(acos)F(atan)F(atan2)F(sinh)F(cosh)F(tanh)F(asinh)F(acosh)F(atanh)F(cbrt)F(log2)F(log10)F(log1p)F(expm1)F(hypot)F(fact)F(logfact)F(erf)F(erfc)F(comb)F(gcd)F(lcm)); #undef F gc_add(mathNS); } @@ -124,20 +126,17 @@ void arith_init() { c(BFn,bi_mul)->ident = c(BFn,bi_div)->ident = c(BFn,bi_and)->ident = c(BFn,bi_eq)->ident = c(BFn,bi_ge)->ident = c(BFn,bi_pow)->ident = c(BFn,bi_not)->ident = m_i32(1); c(BFn,bi_floor)->ident = m_f64(1.0/0.0); c(BFn,bi_ceil )->ident = m_f64(-1.0/0.0); - + + #define INVERSE_PAIR(F,G) \ + c(BFn,bi_##F)->im = G##_c1; \ + c(BFn,bi_##G)->im = F##_c1; c(BFn,bi_sub)->im = sub_c1; - c(BFn,bi_sin)->im = asin_c1; - c(BFn,bi_cos)->im = acos_c1; - c(BFn,bi_tan)->im = atan_c1; - c(BFn,bi_asin)->im = sin_c1; - c(BFn,bi_acos)->im = cos_c1; - c(BFn,bi_atan)->im = tan_c1; - c(BFn,bi_sinh)->im = asinh_c1; - c(BFn,bi_cosh)->im = acosh_c1; - c(BFn,bi_tanh)->im = atanh_c1; - c(BFn,bi_asinh)->im = sinh_c1; - c(BFn,bi_acosh)->im = cosh_c1; - c(BFn,bi_atanh)->im = tanh_c1; - c(BFn,bi_expm1)->im = log1p_c1; - c(BFn,bi_log1p)->im = expm1_c1; + INVERSE_PAIR(sin, asin) + INVERSE_PAIR(cos, acos) + INVERSE_PAIR(tan, atan) + INVERSE_PAIR(sinh, asinh) + INVERSE_PAIR(cosh, acosh) + INVERSE_PAIR(tanh, atanh) + INVERSE_PAIR(expm1, log1p) + #undef INVERSE_PAIR } diff --git a/src/builtins/sysfn.c b/src/builtins/sysfn.c index 17330017..8adbe7a0 100644 --- a/src/builtins/sysfn.c +++ b/src/builtins/sysfn.c @@ -1551,7 +1551,7 @@ u32* dsv_text[] = { U"•file.MapBytes",U"•file.Modified",U"•file.Name",U"•file.Parent",U"•file.Remove",U"•file.Rename",U"•file.Size",U"•file.Type", U"•internal.ClearRefs",U"•internal.DeepSqueeze",U"•internal.EEqual",U"•internal.ElType",U"•internal.HeapDump",U"•internal.Info",U"•internal.IsPure",U"•internal.ListVariations",U"•internal.Refc",U"•internal.Squeeze",U"•internal.Temp",U"•internal.Type",U"•internal.Unshare",U"•internal.Variation", - U"•math.Acos",U"•math.Acosh",U"•math.Asin",U"•math.Asinh",U"•math.Atan",U"•math.Atan2",U"•math.Atanh",U"•math.Cbrt",U"•math.Comb",U"•math.Cos",U"•math.Cosh",U"•math.Expm1",U"•math.Fact",U"•math.GCD",U"•math.Hypot",U"•math.LCM",U"•math.Log10",U"•math.Log1p",U"•math.Log2",U"•math.Sin",U"•math.Sinh",U"•math.Tan",U"•math.Tanh", + U"•math.Acos",U"•math.Acosh",U"•math.Asin",U"•math.Asinh",U"•math.Atan",U"•math.Atan2",U"•math.Atanh",U"•math.Cbrt",U"•math.Comb",U"•math.Cos",U"•math.Cosh",U"•math.Erf",U"•math.ErfC",U"•math.Expm1",U"•math.Fact",U"•math.GCD",U"•math.Hypot",U"•math.LCM",U"•math.Log10",U"•math.Log1p",U"•math.Log2",U"•math.LogFact",U"•math.Sin",U"•math.Sinh",U"•math.Tan",U"•math.Tanh", U"•rand.Deal",U"•rand.Range",U"•rand.Subset", U"•term.CharB",U"•term.CharN",U"•term.ErrRaw",U"•term.Flush",U"•term.OutRaw",U"•term.RawMode", NULL From 7715cbc1580bae58b8a298192bb093bf8838e8a1 Mon Sep 17 00:00:00 2001 From: Marshall Lochbaum Date: Wed, 23 Nov 2022 16:24:15 -0500 Subject: [PATCH 08/11] =?UTF-8?q?Add=20=E2=80=A2math.Sum=20with=20Singeli?= =?UTF-8?q?=20f64=20implementation?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/builtins.h | 2 +- src/builtins/arithm.c | 4 +-- src/builtins/fold.c | 36 ++++++++++++++++++++++++++ src/builtins/sysfn.c | 2 +- src/singeli/src/avx.singeli | 6 ++++- src/singeli/src/fold.singeli | 49 +++++++++++++++++++++++------------- 6 files changed, 77 insertions(+), 22 deletions(-) diff --git a/src/builtins.h b/src/builtins.h index 77e8a1bb..b7b21ebf 100644 --- a/src/builtins.h +++ b/src/builtins.h @@ -17,7 +17,7 @@ /*internal.c*/M(itype,"•internal.Type") M(elType,"•internal.ElType") M(refc,"•internal.Refc") M(isPure,"•internal.IsPure") A(info,"•internal.Info") M(heapDump,"•internal.HeapDump") \ /*internal.c*/M(squeeze,"•internal.Squeeze") M(deepSqueeze,"•internal.DeepSqueeze") D(eequal,"•internal.EEqual") A(internalTemp,"•internal.Temp") \ /*internal.c*/D(variation,"•internal.Variation") A(listVariations,"•internal.ListVariations") M(clearRefs,"•internal.ClearRefs") M(unshare,"•internal.Unshare") \ -/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(hypot,"•math.Hypot") M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") M(cbrt,"•math.Cbrt") M(log2,"•math.Log2") M(log10,"•math.Log10") M(log1p,"•math.Log1p") M(expm1,"•math.Expm1") M(fact,"•math.Fact") D(comb,"•math.Comb") M(logfact,"•math.LogFact") M(erf,"•math.Erf") M(erfc,"•math.ErfC") D(gcd,"•math.GCD") D(lcm,"•math.LCM") +/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(hypot,"•math.Hypot") M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") M(cbrt,"•math.Cbrt") M(log2,"•math.Log2") M(log10,"•math.Log10") M(log1p,"•math.Log1p") M(expm1,"•math.Expm1") M(fact,"•math.Fact") D(comb,"•math.Comb") M(logfact,"•math.LogFact") M(erf,"•math.Erf") M(erfc,"•math.ErfC") D(gcd,"•math.GCD") D(lcm,"•math.LCM") M(sum,"•math.Sum") #define FOR_PM1(A,M,D) \ /*md1.c*/A(tbl,"⌜") A(each,"¨") A(fold,"´") A(scan,"`") A(const,"˙") A(swap,"˜") A(cell,"˘") A(insert,"˝") \ diff --git a/src/builtins/arithm.c b/src/builtins/arithm.c index 7e6b40f4..daebb007 100644 --- a/src/builtins/arithm.c +++ b/src/builtins/arithm.c @@ -113,8 +113,8 @@ static B mathNS; B getMathNS() { if (mathNS.u == 0) { #define F(X) inc(bi_##X), - Body* d = m_nnsDesc("sin","cos","tan","asin","acos","atan","atan2","sinh","cosh","tanh","asinh","acosh","atanh","cbrt","log2","log10","log1p","expm1","hypot","fact","logfact","erf","erfc","comb","gcd","lcm"); - mathNS = m_nns(d, F(sin)F(cos)F(tan)F(asin)F(acos)F(atan)F(atan2)F(sinh)F(cosh)F(tanh)F(asinh)F(acosh)F(atanh)F(cbrt)F(log2)F(log10)F(log1p)F(expm1)F(hypot)F(fact)F(logfact)F(erf)F(erfc)F(comb)F(gcd)F(lcm)); + Body* d = m_nnsDesc("sin","cos","tan","asin","acos","atan","atan2","sinh","cosh","tanh","asinh","acosh","atanh","cbrt","log2","log10","log1p","expm1","hypot","fact","logfact","erf","erfc","comb","gcd","lcm","sum"); + mathNS = m_nns(d, F(sin)F(cos)F(tan)F(asin)F(acos)F(atan)F(atan2)F(sinh)F(cosh)F(tanh)F(asinh)F(acosh)F(atanh)F(cbrt)F(log2)F(log10)F(log1p)F(expm1)F(hypot)F(fact)F(logfact)F(erf)F(erfc)F(comb)F(gcd)F(lcm)F(sum)); #undef F gc_add(mathNS); } diff --git a/src/builtins/fold.c b/src/builtins/fold.c index 8937050a..378a76d4 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -9,6 +9,8 @@ // COULD implement fast numeric -´ // ∨ on boolean-valued integers, stopping at 1 +// •math.Sum: +´ with faster and more precise SIMD code for i32, f64 + #include "../core.h" #include "../builtins.h" @@ -66,6 +68,40 @@ static f64 sum_f64(void* xv, usz i, f64 r) { static i64 (*const sum_small_fns[])(void*, usz) = { sum_small_i8, sum_small_i16, sum_small_i32 }; static f64 (*const sum_fns[])(void*, usz, f64) = { sum_i8, sum_i16, sum_i32, sum_f64 }; +B sum_c1(B t, B x) { + if (isAtm(x) || RNK(x)!=1) thrF("•math.Sum: Argument must be a list (%H ≡ ≢𝕩)", x); + usz ia = IA(x); + if (ia==0) return m_f64(0); + u8 xe = TI(x,elType); + if (!elNum(xe)) { + x = any_squeeze(x); xe = TI(x,elType); + if (!elNum(xe)) thrF("•math.Sum: Argument elements must be numbers", x); + } + f64 r; + void* xv = tyany_ptr(x); + if (xe == el_bit) { + r = bit_sum(xv, ia); + } else if (xe <= el_i32) { + u8 sel = xe - el_i8; + i64 s = 0; r = 0; + i64 m = 1ull<<48; + usz b = sum_small_max; + for (usz i=0; i= m) { r+=m; s-=m; } + if (s <= -m) { r-=m; s+=m; } + } + r += s; + } else { + #if SINGELI + r = avx2_sum_f64(xv, ia); + #else + r=0; for (usz i=0; i> plog + r:= init + @for (i to f) r = op{r, pairwise{x+(i< 0} r = load{xv} if (len > 1) { - if (len > 2) r = op{r, shuf{[4]u64, r, 4b2222}} - r = op{r, shuf{[4]u64, r, 4b1111}} + if (len > 2) r = opsh{op}{r, 4b2222} + r = opsh{op}{r, 4b1111} } } else { - # Pairwise combination to shorten dependency chains - def pairwise{p, i, k} = { - def l = k-1 - op{pairwise{p, i , l}, - pairwise{p, i+(1<> pk - @for (i to f) r = op{r, pairwise{xv+(i< Date: Wed, 23 Nov 2022 16:43:30 -0500 Subject: [PATCH 09/11] Formatting --- src/builtins.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/builtins.h b/src/builtins.h index b7b21ebf..c1691095 100644 --- a/src/builtins.h +++ b/src/builtins.h @@ -17,7 +17,10 @@ /*internal.c*/M(itype,"•internal.Type") M(elType,"•internal.ElType") M(refc,"•internal.Refc") M(isPure,"•internal.IsPure") A(info,"•internal.Info") M(heapDump,"•internal.HeapDump") \ /*internal.c*/M(squeeze,"•internal.Squeeze") M(deepSqueeze,"•internal.DeepSqueeze") D(eequal,"•internal.EEqual") A(internalTemp,"•internal.Temp") \ /*internal.c*/D(variation,"•internal.Variation") A(listVariations,"•internal.ListVariations") M(clearRefs,"•internal.ClearRefs") M(unshare,"•internal.Unshare") \ -/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(hypot,"•math.Hypot") M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") M(cbrt,"•math.Cbrt") M(log2,"•math.Log2") M(log10,"•math.Log10") M(log1p,"•math.Log1p") M(expm1,"•math.Expm1") M(fact,"•math.Fact") D(comb,"•math.Comb") M(logfact,"•math.LogFact") M(erf,"•math.Erf") M(erfc,"•math.ErfC") D(gcd,"•math.GCD") D(lcm,"•math.LCM") M(sum,"•math.Sum") +/* arithm.c*/M(sin,"•math.Sin") M(cos,"•math.Cos") M(tan,"•math.Tan") M(asin,"•math.Asin") M(acos,"•math.Acos") M(atan,"•math.Atan") D(atan2,"•math.Atan2") D(hypot,"•math.Hypot") \ +/* arithm.c*/M(sinh,"•math.Sinh") M(cosh,"•math.Cosh") M(tanh,"•math.Tanh") M(asinh,"•math.Asinh") M(acosh,"•math.Acosh") M(atanh,"•math.Atanh") \ +/* arithm.c*/M(cbrt,"•math.Cbrt") M(log2,"•math.Log2") M(log10,"•math.Log10") M(log1p,"•math.Log1p") M(expm1,"•math.Expm1") M(fact,"•math.Fact") D(comb,"•math.Comb") M(logfact,"•math.LogFact") \ +/* arithm.c*/M(erf,"•math.Erf") M(erfc,"•math.ErfC") D(gcd,"•math.GCD") D(lcm,"•math.LCM") M(sum,"•math.Sum") #define FOR_PM1(A,M,D) \ /*md1.c*/A(tbl,"⌜") A(each,"¨") A(fold,"´") A(scan,"`") A(const,"˙") A(swap,"˜") A(cell,"˘") A(insert,"˝") \ From cf5ce07f51ed96f9401244eb1fb377b9c2a54c45 Mon Sep 17 00:00:00 2001 From: dzaima Date: Thu, 24 Nov 2022 00:21:19 +0200 Subject: [PATCH 10/11] simplify singeli xor/and/or definitions --- src/singeli/src/avx.singeli | 10 +++------- src/singeli/src/sse3.singeli | 6 +++--- 2 files changed, 6 insertions(+), 10 deletions(-) diff --git a/src/singeli/src/avx.singeli b/src/singeli/src/avx.singeli index 76434442..f27c7878 100644 --- a/src/singeli/src/avx.singeli +++ b/src/singeli/src/avx.singeli @@ -61,15 +61,11 @@ def iota{T & w256{T,16}} = make{T,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15} def iota{T & w256{T,8}} = make{T,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31} # bit arith -def __xor{a:T, b:T & w256i{T}} = T ~~ emit{[8]f32, '_mm256_xor_ps', v2f{a}, v2f{b}} -def __and{a:T, b:T & w256i{T}} = T ~~ emit{[8]f32, '_mm256_and_ps', v2f{a}, v2f{b}} -def __or {a:T, b:T & w256i{T}} = T ~~ emit{[8]f32, '_mm256_or_ps', v2f{a}, v2f{b}} +def __xor{a:T, b:T & w256{T}} = T ~~ emit{[8]f32, '_mm256_xor_ps', v2f{a}, v2f{b}} +def __and{a:T, b:T & w256{T}} = T ~~ emit{[8]f32, '_mm256_and_ps', v2f{a}, v2f{b}} +def __or {a:T, b:T & w256{T}} = T ~~ emit{[8]f32, '_mm256_or_ps', v2f{a}, v2f{b}} def __not{a:T & w256u{T}} = a ^ broadcast{T, ~cast{eltype{T},0}} -def __xor{a:T, b:T & T==[4]f64} = emit{[4]f64, '_mm256_xor_pd', a, b} -def __and{a:T, b:T & T==[4]f64} = emit{[4]f64, '_mm256_and_pd', a, b} -def __or {a:T, b:T & T==[4]f64} = emit{[4]f64, '_mm256_or_pd', a, b} - # float comparison local def f32cmpAVX{a,b,n} = [8]u32 ~~ emit{[8]f32, '_mm256_cmp_ps', a, b, n} local def f64cmpAVX{a,b,n} = [4]u64 ~~ emit{[4]f64, '_mm256_cmp_pd', a, b, n} diff --git a/src/singeli/src/sse3.singeli b/src/singeli/src/sse3.singeli index 39d95049..8d0298ab 100644 --- a/src/singeli/src/sse3.singeli +++ b/src/singeli/src/sse3.singeli @@ -59,9 +59,9 @@ def iota{T & w128{T,16}} = make{T,0,1,2,3,4,5,6,7} def iota{T & w128{T,8}} = make{T,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15} # bit arith -def __xor{a:T, b:T & w128i{T}} = T ~~ emit{[4]f32, '_mm_xor_ps', v2f{a}, v2f{b}} -def __and{a:T, b:T & w128i{T}} = T ~~ emit{[4]f32, '_mm_and_ps', v2f{a}, v2f{b}} -def __or {a:T, b:T & w128i{T}} = T ~~ emit{[4]f32, '_mm_or_ps', v2f{a}, v2f{b}} +def __xor{a:T, b:T & w128{T}} = T ~~ emit{[4]f32, '_mm_xor_ps', v2f{a}, v2f{b}} +def __and{a:T, b:T & w128{T}} = T ~~ emit{[4]f32, '_mm_and_ps', v2f{a}, v2f{b}} +def __or {a:T, b:T & w128{T}} = T ~~ emit{[4]f32, '_mm_or_ps', v2f{a}, v2f{b}} def __not{a:T & w128u{T}} = a ^ broadcast{T, ~cast{eltype{T},0}} # signed comparison From 34e37d1b1a4cfc2481de5e79ef233fea59730cf4 Mon Sep 17 00:00:00 2001 From: dzaima Date: Thu, 24 Nov 2022 01:09:50 +0200 Subject: [PATCH 11/11] increase gc root buffer size --- src/opt/gc.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/opt/gc.c b/src/opt/gc.c index 0b9c2cda..6dda6768 100644 --- a/src/opt/gc.c +++ b/src/opt/gc.c @@ -14,11 +14,11 @@ void gc_addFn(vfn f) { gc_roots[gc_rootSz++] = f; } -Value* gc_rootObjs[256]; +Value* gc_rootObjs[512]; u32 gc_rootObjSz; void gc_add(B x) { assert(isVal(x)); - if (gc_rootObjSz>=256) err("Too many GC root objects"); + if (gc_rootObjSz>=512) err("Too many GC root objects"); gc_rootObjs[gc_rootObjSz++] = v(x); }