Add •math.Sum with Singeli f64 implementation

This commit is contained in:
Marshall Lochbaum 2022-11-23 16:24:15 -05:00
parent 35aef1e5cd
commit 7715cbc158
6 changed files with 77 additions and 22 deletions

View File

@ -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,"˝") \

View File

@ -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);
}

View File

@ -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<ia; i+=b) {
s += sum_small_fns[sel]((u8*)xv + (i<<sel), ia-i<b? ia-i : b);
if (s >= 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<ia; i++) r+=((f64*)xv)[i];
#endif
}
decG(x); return m_f64(r);
}
// Try to keep to i32 product, go to f64 on overflow or non-i32 initial
#define DEF_INT_PROD(T) \
static NOINLINE f64 prod_##T(void* xv, usz i, f64 init) { \

View File

@ -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.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"•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.Sum",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

View File

@ -66,6 +66,10 @@ def __and{a:T, b:T & w256i{T}} = T ~~ emit{[8]f32, '_mm256_and_ps', v2f{a}, v2f{
def __or {a:T, b:T & w256i{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}
@ -132,4 +136,4 @@ def getmask{x:T & w256{T, 64}} = emit{u8, '_mm256_movemask_pd', v2d{x}}
def andIsZero{x:T, y:T & w256i{T}} = emit{u1, '_mm256_testz_si256', x, y}
def any{x:T & w256i{T}} = getmask{x} != 0 # assumes elements of x all have equal bits (avx2 utilizes this for 16 bits)
def all{x:T & w256i{T}} = getmask{x} == (1<<vcount{T})-1 # same assumption
def anyneg{x:T & w256s{T}} = getmask{x}!=0
def anyneg{x:T & w256s{T}} = getmask{x}!=0

View File

@ -4,6 +4,23 @@ include './avx'
include './avx2'
include './mask'
def opsh{op}{v:[4]f64, perm} = op{v, shuf{[4]u64, v, perm}}
def mix{op, v:[4]f64} = { def sh=opsh{op}; sh{sh{v, 4b2301}, 4b1032} }
def reduce_pairwise{op, plog, x:*T, len, init:T} = {
# Pairwise combination to shorten dependency chains
def pairwise{p, i, k} = (if (k==0) { load{p,i} } else {
def l = k-1
op{pairwise{p, i , l},
pairwise{p, i+(1<<l), l}}
})
f:= len >> plog
r:= init
@for (i to f) r = op{r, pairwise{x+(i<<plog), 0, plog}}
@for (x over i from f<<plog to len) r = op{r, x}
r
}
fold_idem{T==f64, op}(x:*T, len:u64) : T = {
def step = 256/width{T}
def V = [step]T
@ -14,28 +31,26 @@ fold_idem{T==f64, op}(x:*T, len:u64) : T = {
assert{len > 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<<l), l}}
}
def pairwise{p, i, k==0} = load{p, i}
def pk = 2 # Combine 1<<pk values in a step
r = load{*V ~~ (x+len-step)}
e:= (len-1)/step
f:= e >> pk
@for (i to f) r = op{r, pairwise{xv+(i<<pk), 0, pk}}
@for (xv over i from f<<pk to e) r = op{r, xv}
r = op{r, shuf{[4]u64, r, 4b2301}}
r = op{r, shuf{[4]u64, r, 4b1032}}
i:= load{*V ~~ (x+len-step)}
r = mix{op, reduce_pairwise{op, 2, xv, (len-1)/step, i}}
}
extract{r, 0}
}
'avx2_fold_min_f64' = fold_idem{f64,min}
'avx2_fold_max_f64' = fold_idem{f64,max}
fold_assoc_0{T==f64, op}(x:*T, len:u64) : T = {
def step = 256/width{T}
def V = [step]T
xv:= *V ~~ x
e:= len / step
i:= load{xv, e} & (V~~maskOf{V, len % step})
r:= reduce_pairwise{op, 2, xv, e, i}
extract{mix{op, r}, 0}
}
'avx2_sum_f64' = fold_assoc_0{f64,+}