#include "../core.h" #include "../utils/each.h" #include "../builtins.h" #include static f64 pfmod(f64 a, f64 b) { f64 r = fmod(a, b); if (a<0 != b<0 && r!=0) r+= b; return r; } static f64 bqn_or(f64 w, f64 x) { return (w+x)-(w*x); } B add_c2(B, B, B); B ceil_c2(B, B, B); B sub_c2(B, B, B); B div_c2(B, B, B); B not_c2(B, B, B); B pow_c2(B, B, B); B mul_c2(B, B, B); B root_c2(B, B, B); B and_c2(B, B, B); B stile_c2(B, B, B); B or_c2 (B, B, B); B log_c2(B, B, B); B floor_c2(B, B, B); B atan2_c2(B, B, B); typedef void (*AndBytesFn)(u8*, u8*, u64, u64); B leading_axis_arith(FC2 fc2, B w, B x, usz* wsh, usz* xsh, ur mr); #if SINGELI_SIMD #include "../singeli/c/arithdDispatch.c" static const AndBytesFn andBytes_fn = simd_andBytes; #else static void base_andBytes(u8* r, u8* x, u64 repeatedMask, u64 numBytes) { u64* x64 = (u64*)x; usz i; vfor (i = 0; i < numBytes/8; i++) ((u64*)r)[i] = x64[i] & repeatedMask; if (i*8 != numBytes) { u64 v = x64[i]&repeatedMask; for (usz j = 0; j < (numBytes&7); j++) r[i*8 + j] = v>>(j*8); } } static const AndBytesFn andBytes_fn = base_andBytes; #endif B floor_c1(B t, B x); B sub_c1(B t, B x); B fne_c1(B t, B x); B shape_c2(B t, B w, B x); // all divint/floordiv/modint assume integer arguments // floordiv will return float result only on Β―2147483648Γ·Β―1 or nΓ·0, but may not otherwise squeeze integer types; integer argument requirement may be relaxed in the future // divint will return float result if there's a fractional result, or in overflow cases same as floordiv // TODO overflow-checked Singeli code for exact integer divint, and maybe floordiv_AA #define DIVLOOP(RE, WE, EXPR) RE* rp; B r=m_##RE##arrc(&rp, w); usz ia=IA(w); WE* wp=WE##any_ptr(w); vfor(ux i=0; i0 && (wi32&(wi32-1))==0) { if (wi32==1) return i64EachDec(0, x); if (xe==el_bit) return x; // if n>1 (true from the above), 0β€Ώ1 ≑ (2⋆n)|0β€Ώ1 u8 elw = elWidth(xe); u32 mask0 = (u32)wi32; if (mask0 > ((u32)1 << (elw*8-1))) { if (mask0 > 32768) { x=taga(cpyI32Arr(x)); xe=el_i32; elw=4; } else if (mask0 > 128) { x=taga(cpyI16Arr(x)); xe=el_i16; elw=2; } else UD; } u64 mask = (mask0-1)*repeatNum[xe]; usz bytes = IA(x)*elw; u8* rp = m_tyarrc(&r, elw, x, el2t(xe)); andBytes_fn(rp, tyany_ptr(x), mask, bytes); decG(x); if (wi32==2) return taga(cpyBitArr(r)); else if (wi32<256) return taga(cpyI8Arr(r)); // these won't widen, as the code doesn't even get to here if 𝕨 > max possible in 𝕩 else if (wi32<65536) return taga(cpyI16Arr(r)); return r; } else { return modint_SA(wi32, x); } } , /*INT_AS*/ if (q_i32(x)) return modint_AS(w, x); , /*INT_AA*/ return modint_AA(w, x); , /*FLT_SAI*/ if (o2fG(w)==1) { vfor (usz i=0; ixv?xv:wv, 2, {}) AR_I_AA("⌈", ceil , wv>xv?wv:xv, 1, {}) AR_I_AA("+", add, wv+xv, 0, {}) AR_I_AA("-", sub, wv-xv, 0, { if (we==el_c32 && xe==el_i32) { u32* wp = c32any_ptr(w); usz wia = IA(w); u32* rp; r = m_c32arrc(&rp, w); i32* xp = i32any_ptr(x); for (usz i = 0; i < wia; i++) { rp[i] = (u32)((i32)wp[i] - (i32)xp[i]); if (rp[i]>CHR_MAX) thrM("𝕨-𝕩: Invalid character"); // safe - see add } goto dec_ret; } }) #undef AR_I_AA #define AR_I_AS(CHR, NAME, EXPR, DO_AS, EXTRA) NOINLINE B NAME##_AS(B t, B w, B x) { \ B r; u8 we=TI(w,elType); EXTRA \ if (isF64(x)) { usz ia=IA(w); DO_AS(NAME,EXPR) } \ ARITH_SLOW(CHR); return arith_recd(NAME##_c2, w, x); \ dec_ret: decG(w); return r; \ } #define AR_I_SA(CHR, NAME, EXPR, DO_SA, EXTRA) NOINLINE B NAME##_SA(B t, B w, B x) { \ B r; u8 xe=TI(x,elType); EXTRA \ if (isF64(w)) { usz ia=IA(x); DO_SA(NAME,EXPR) } \ ARITH_SLOW(CHR); return arith_recd(NAME##_c2, w, x); \ dec_ret: decG(x); return r; \ } AR_I_SA("-", sub, wv-xv, REG_SA, {}) AR_I_SA("Γ—", mul, wv*xv, REG_SA, {}) AR_I_SA("∧", and, wv*xv, REG_SA, {}) AR_I_SA("∨", or , (wv+xv)-(wv*xv), REG_SA, {}) AR_I_SA("⌊", floor, wv>xv?xv:wv, REG_SA, {}) AR_I_SA("⌈", ceil , wv>xv?wv:xv, REG_SA, {}) AR_I_SA("+", add, wv+xv, REG_SA, { if (isC32(w) && elInt(xe)) { u32 wv = o2cG(w); if (xe!=el_i32) x = taga(cpyI32Arr(x)); i32* xp = i32any_ptr(x); usz xia = IA(x); u32* rp; r = m_c32arrc(&rp, x); for (usz i = 0; i < xia; i++) { rp[i] = (u32)(xp[i]+(i32)wv); if (rp[i]>CHR_MAX) thrM("𝕨+𝕩: Invalid character"); // safe to only check this as wv already must be below CHR_MAX, which is less than U32_MAX/2 } goto dec_ret; } }) #undef AR_I_SA AR_I_AS("-", sub, wv-xv, REG_AS, { if (elChr(we) && isC32(x)) { i32 xv = (i32)o2cG(x); if (we!=el_c32) w = taga(cpyC32Arr(w)); u32* wp = c32any_ptr(w); usz wia = IA(w); i32* rp; r = m_i32arrc(&rp, w); vfor (usz i = 0; i < wia; i++) rp[i] = (i32)wp[i] - xv; goto dec_ret; } }) #undef AR_I_AS #endif // !SINGELI #define add_AS(T, W, X) add_SA(T, X, W) #define mul_AS(T, W, X) mul_SA(T, X, W) #define and_AS(T, W, X) and_SA(T, X, W) #define or_AS(T, W, X) or_SA(T, X, W) #define floor_AS(T, W, X) floor_SA(T, X, W) #define ceil_AS(T, W, X) ceil_SA(T, X, W) #define AR_F_TO_ARR(NAME) return NAME##_c2_arr(t, w, x); #define AR_I_TO_ARR(NAME) \ if (isArr(x)) return isArr(w)? NAME##_AA(t, w, x) : NAME##_SA(t, w, x); \ else if (isArr(w)) return NAME##_AS(t, w, x); #endif // TYPED_ARITH #define AR_I_SCALAR(CHR, NAME, EXPR, MORE) B NAME##_c2(B t, B w, B x) { \ if (isF64(w) & isF64(x)) return m_f64(EXPR); \ MORE; AR_I_TO_ARR(NAME) \ thrM("𝕨"CHR "𝕩: Unexpected argument types"); \ } AR_I_SCALAR("+", add, w.f+x.f, { if (isC32(w) & isF64(x)) { u64 r = (u64)(o2cG(w)+o2i64(x)); if(r>CHR_MAX)thrM("𝕨+𝕩: Invalid character"); return m_c32((u32)r); } if (isF64(w) & isC32(x)) { u64 r = (u64)(o2cG(x)+o2i64(w)); if(r>CHR_MAX)thrM("𝕨+𝕩: Invalid character"); return m_c32((u32)r); } }); AR_I_SCALAR("-", sub, w.f-x.f, { if (isC32(w) & isF64(x)) { u64 r = (u64)((i32)o2cG(w)-o2i64(x)); if(r>CHR_MAX)thrM("𝕨-𝕩: Invalid character"); return m_c32((u32)r); } if (isC32(w) & isC32(x)) return m_f64((i32)(u32)w.u - (i32)(u32)x.u); }) AR_I_SCALAR("Γ—", mul, w.f*x.f, {}) AR_I_SCALAR("∧", and, w.f*x.f, {}) AR_I_SCALAR("∨", or , bqn_or(w.f, x.f), {}) AR_I_SCALAR("⌊", floor, w.f>x.f?x.f:w.f, {}) AR_I_SCALAR("⌈", ceil , w.f>x.f?w.f:x.f, {}) #undef AR_I_SCALAR B not_c2(B t, B w, B x) { return C2(add, m_f64(1), sub_c2(t, w, x)); } #define AR_F_SCALAR(CHR, NAME, EXPR) B NAME##_c2(B t, B w, B x) { \ if (isF64(w) & isF64(x)) return m_f64(EXPR); \ AR_F_TO_ARR(NAME) \ thrM("𝕨"CHR "𝕩: Unexpected argument types"); \ } AR_F_SCALAR("Γ·", div , w.f/(x.f+0)) AR_F_SCALAR("⋆", pow , pow(w.f+0, x.f)) AR_F_SCALAR("√", root , pow(x.f+0, 1.0/(0+w.f))) 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) { assert(k>=0 && n>=2*k); 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 (k= 0) { return j<0? 0 : comb_nat(j, n); } else { if (k<0) return 0; f64 l = -1-n; // l+k == -1-j f64 r = comb_nat(k>= bz; while (1) { a >>= az; u64 d = b - a; b = biw = atan2iw_c2; c(BFn, bi_atan2)->ix = atan2ix_c2; c(BFn, bi_pow)->ix = log_c2; }