diff --git a/src/builtins/fold.c b/src/builtins/fold.c index 9a5a3aef..1adb97f4 100644 --- a/src/builtins/fold.c +++ b/src/builtins/fold.c @@ -16,6 +16,51 @@ static i64 bit_diff(u64* x, u64 am) { return r - (i64)(am/2); } +// It's safe to sum a block of integers as long as the current total +// is far enough from +-1ull<<53 (and integer, in dyadic fold). +static const usz sum_small_max = 1<<16; +#define DEF_INT_SUM(T,W,M,A) \ + static i64 sum_small_##T(void* xv, usz ia) { \ + i##A s=0; for (usz i=0; i0 && -lim<=c && c<=lim) { \ + usz e=i0; i0=(i0-1)&~(b-1); \ + c+=sum_small_##T(xp+i0, e-i0); \ + } \ + r = c; \ + } \ + while (i0--) r+=xp[i0]; \ + return r; \ + } +DEF_INT_SUM(i8 ,8 ,32,32) +DEF_INT_SUM(i16,16,32,32) +DEF_INT_SUM(i32,32,52,64) +#undef DEF_SUM +static f64 sum_f64(void* xv, usz i, f64 r) { + while (i--) r += ((f64*)xv)[i]; + return 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 }; + +#define MIN_MAX(T,C) \ + T* xp = xv; T r = xp[0]; \ + for (usz i=1; i) } +DEF_MIN_MAX(i8) DEF_MIN_MAX(i16) DEF_MIN_MAX(i32) DEF_MIN_MAX(f64) +#undef DEF_MIN_MAX +#undef MIN_MAX +static f64 (*const min_fns[])(void*, usz) = { min_i8, min_i16, min_i32, min_f64 }; +static f64 (*const max_fns[])(void*, usz) = { max_i8, max_i16, max_i32, max_f64 }; + B fold_c1(Md1D* d, B x) { B f = d->f; if (isAtm(x) || RNK(x)!=1) thrF("ยด: Argument must be a list (%H โ‰ก โ‰ข๐•ฉ)", x); usz ia = IA(x); @@ -45,23 +90,10 @@ B fold_c1(Md1D* d, B x) { B f = d->f; } if (rtid==n_add) { // + void *xv = tyany_ptr(x); - f64 r; - #define CASE_INT(T,M,A) case el_i##T: { \ - usz b=1<<(M-T); i64 lim = (1ull<<53) - (1ull<0 && -lim<=c && c<=lim) { \ - usz e=i0; i0=(i0-1)&~(b-1); \ - i##A s=0; for (usz i=i0; if; if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=1; for (usz i=ia; i--; ) if (mulOn(c,xp[i]))goto base; decG(x); return m_i32(c); } if (xe==el_f64) { f64* xp = f64any_ptr(x); f64 c=1; for (usz i=ia; i--; ) c*= xp[i]; decG(x); return m_f64(c); } } - #define CASE(T,C) case el_##T: { \ - T* xp = xv; T c = xp[0]; \ - for (usz i=0; i) } // โŒˆ - #undef FC - #undef CASE + if (rtid==n_floor) { f64 r=min_fns[xe-el_i8](tyany_ptr(x), ia); decG(x); return m_f64(r); } // โŒŠ + if (rtid==n_ceil ) { f64 r=max_fns[xe-el_i8](tyany_ptr(x), ia); decG(x); return m_f64(r); } // โŒˆ if (rtid==n_or) { // โˆจ if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool r=0; for (usz i=0; if; goto base; } if (rtid==n_add) { // + - if (xe==el_i8 ) { i8* xp = i8any_ptr (x); i64 c=wi; for (usz i=0; ic) c=xp[i]; decG(x); return m_i32(c); } - if (xe==el_i16) { i16* xp = i16any_ptr(x); i32 c=wi; for (usz i=0; ic) c=xp[i]; decG(x); return m_i32(c); } - if (xe==el_i32) { i32* xp = i32any_ptr(x); i32 c=wi; for (usz i=0; ic) c=xp[i]; decG(x); return m_i32(c); } - } + if (rtid==n_floor) { f64 r=wi; if (ia>0) { f64 m=min_fns[xe-el_i8](tyany_ptr(x), ia); if (m0) { f64 m=max_fns[xe-el_i8](tyany_ptr(x), ia); if (m>r) r=m; } decG(x); return m_f64(r); } // โŒˆ if (rtid==n_or && (wi&1)==wi) { // โˆจ if (xe==el_i8 ) { i8* xp = i8any_ptr (x); bool q=wi; for (usz i=0; i