totk-pc-recomp/runtime/simd_fp.h

// Adapted from LAS Recomp runtime/simd_fp.h @ da35b83 (MIT, Copyright (c) 2026 Peiikai).
// FP helpers for one precision; included twice by simd.h (F = s/float, d/double). No include guard.
// Scalar entry points take/return T, vector ones VT (all lanes; callers zero the unused half for Q=0).

#define FN(x) SIMD_CAT(x, F)
#define SGN ((U)1 << (NB - 1))
#define EXPM ((((U)1 << EB) - 1) << MB)
#define QB ((U)1 << (MB - 1))
#define DNAN (EXPM | QB)
#define MINN ((U)1 << MB) // smallest normal
#define VK(v) ((VT){} + (T)(v))
#define K(v) ((T)(v))

SI U FN(fb)(T x) { U u; memcpy(&u, &x, sizeof u); return u; }
SI T FN(fv)(U u) { T x; memcpy(&x, &u, sizeof x); return x; }
SI int FN(isnan)(U u) { return (u & ~SGN) > EXPM; }
SI int FN(issnan)(U u) { return FN(isnan)(u) && !(u & QB); }
SI int FN(isqnan)(U u) { return FN(isnan)(u) && (u & QB); }
SI int FN(isinf)(U u) { return (u & ~SGN) == EXPM; }
SI int FN(iszero)(U u) { return !(u & ~SGN); }
SI int FN(isden)(U u) { return !(u & EXPM) && (u & ~SGN); }
// The asm hides the float origin of u: LLVM turns the bit test into an FP class compare otherwise, and with
// MXCSR.DAZ set (FPCR.FZ) the host compare sees a denormal as zero. Same in FN(vflush).
SI T FN(flush)(T x) {
    U u = FN(fb)(x);
    __asm__("" : "+r"(u));
    return FN(fv)(FN(isden)(u) ? u & SGN : u);
}
// FPCR.FZ input flush for helpers that run host FP.
SI T FN(fz)(Cpu *c, T x) { return c->fpcr & FPCR_FZ ? FN(flush)(x) : x; }
SI T FN(neg)(T x) { return FN(fv)(FN(fb)(x) ^ SGN); }

// FPProcessNaN(s): first signalling NaN (quietened), else first quiet NaN, else (invalid op) default NaN
SN T FN(nan1)(Cpu *c, T a) {
    U x = FN(fb)(a);
    if (!(c->fpcr & FPCR_DN) && FN(isnan)(x)) return FN(fv)(x | QB);
    return FN(fv)(DNAN);
}
SN T FN(nan2)(Cpu *c, T a, T b) {
    U x = FN(fb)(a), y = FN(fb)(b);
    if (!(c->fpcr & FPCR_DN)) {
        if (FN(issnan)(x)) return FN(fv)(x | QB);
        if (FN(issnan)(y)) return FN(fv)(y | QB);
        if (FN(isnan)(x)) return a;
        if (FN(isnan)(y)) return b;
    }
    return FN(fv)(DNAN);
}
SN T FN(nan3)(Cpu *c, T a, T b, T d) {
    U x = FN(fb)(a), y = FN(fb)(b), z = FN(fb)(d);
    if (!(c->fpcr & FPCR_DN)) {
        if (FN(issnan)(x)) return FN(fv)(x | QB);
        if (FN(issnan)(y)) return FN(fv)(y | QB);
        if (FN(issnan)(z)) return FN(fv)(z | QB);
        if (FN(isnan)(x)) return a;
        if (FN(isnan)(y)) return b;
        if (FN(isnan)(z)) return d;
    }
    return FN(fv)(DNAN);
}

// ---- FMA flavours. FMA_S / FMA_V: fast paths, host MXCSR is round-to-nearest there (simd_directed() is checked first).
// FMA_X: scalar FMA for the slow paths, any host rounding mode.
// Without hardware FMA (Ivy Bridge, no -mfma) clang turns __builtin_elementwise_fma into a call to the C runtime's
// fmaf/fma for every lane, which is very slow and keeps the thread outside generated code. For float the product of
// two floats is exact in a double (24 x 24 = 48 bits), so the fused op is done in double. A plain
// (float)(x * y + z) rounds twice and can differ by one ulp; rounding the double sum to odd first (Boldo-Melquiond)
// makes the final float rounding exact: TwoSum gives the error e of s = fl(p + z), and when e != 0 and s has an even
// last bit, s is moved one ulp toward the exact value (the odd neighbour). Inf/NaN lanes are left alone (e is NaN
// there, so the ordered compares below are false). Double has no wider type: it keeps the C runtime's fma.
#if defined(__FMA__) || defined(__aarch64__) || NB == 64
#define FMA_S(x, y, z) __builtin_elementwise_fma(x, y, z)
#define FMA_V(x, y, z) __builtin_elementwise_fma(x, y, z)
#define FMA_X(x, y, z) __builtin_elementwise_fma(x, y, z)
#else
typedef double simd_vf64x4 __attribute__((vector_size(32)));
typedef uint64_t simd_vu64x4 __attribute__((vector_size(32)));
SI float fma_rn_s(float x, float y, float z) {
    double p = (double)x * (double)y, zd = (double)z, s = p + zd;
    double bb = s - p, e = (p - (s - bb)) + (zd - bb); // TwoSum: s + e == p + zd exactly
    uint64_t u;
    memcpy(&u, &s, sizeof u);
    if (((e < 0.0) | (e > 0.0)) && !(u & 1)) u += (e > 0.0) == (s > 0.0) ? 1 : (uint64_t)-1;
    memcpy(&s, &u, sizeof s);
    return (float)s;
}
SI vf32 fma_rn_v(vf32 x, vf32 y, vf32 z) {
    simd_vf64x4 p = CVT(x, simd_vf64x4) * CVT(y, simd_vf64x4), zd = CVT(z, simd_vf64x4), s = p + zd;
    simd_vf64x4 bb = s - p, e = (p - (s - bb)) + (zd - bb);
    simd_vu64x4 u = (simd_vu64x4)s;
    simd_vu64x4 fix = (simd_vu64x4)((e < 0.0) | (e > 0.0)) & (simd_vu64x4)((u & 1) == 0);
    simd_vu64x4 up = (simd_vu64x4)((e > 0.0) == (s > 0.0));
    u += fix & ((up & 1) | ~up); // +1 toward the exact value, i.e. +1 if e and s share a sign, else -1
    return CVT((simd_vf64x4)u, vf32);
}
// Any rounding mode: directed modes are exact with a plain double (same-direction double rounding is harmless).
SI float fma_any_s(float x, float y, float z) {
    if (__builtin_expect(_mm_getcsr() & 0x6000u, 0)) return (float)((double)x * (double)y + (double)z);
    return fma_rn_s(x, y, z);
}
#define FMA_S(x, y, z) fma_rn_s(x, y, z)
#define FMA_V(x, y, z) fma_rn_v(x, y, z)
#define FMA_X(x, y, z) fma_any_s(x, y, z)
#endif

// The arithmetic itself, for scalars (KK = K) and vectors (KK = VK). MADD: a + b*d (a = addend).
// RECPS/RSQRTS get their first operand already negated (FPNeg happens before NaN processing).
#define FCALC(op, a, b, d, KK, FMA)                                                                     \
    ((op) == FOP_ADD ? (a) + (b) : (op) == FOP_SUB ? (a) - (b)                                     \
   : (op) == FOP_MUL || (op) == FOP_MULX ? (a) * (b) : (op) == FOP_NMUL ? -((a) * (b))              \
   : (op) == FOP_DIV ? (a) / (b) : (op) == FOP_ABD ? __builtin_elementwise_abs((a) - (b))          \
   : (op) == FOP_SQRT ? __builtin_elementwise_sqrt(a)                                              \
   : (op) == FOP_MADD ? FMA(b, d, a)                                                            \
   : (op) == FOP_RECPS ? FMA(a, b, KK(2))                                                     \
   : (op) == FOP_RSQRTS ? FMA(a, b, KK(3)) * KK(0.5) : (a))

// Lanes whose fast result needs fixing: NaN results (NaN operand or invalid op), RSQRTS overflow.
#define FBAD(op, r) ((op) == FOP_RSQRTS ? ((r) != (r)) | (__builtin_elementwise_abs(r) == K(INFINITY)) : (r) != (r))

// (3 + a*b) / 2 with one rounding (a already negated): halve a normal operand first, so a product beyond
// the format's range (which the fast path rounds to inf, or to the largest normal when rounding down) is exact.
SN T FN(rsqrts)(T a, T b) {
    if (__builtin_elementwise_abs(a) >= FN(fv)(2 * MINN)) return FMA_X(a * K(0.5), b, K(1.5));
    if (__builtin_elementwise_abs(b) >= FN(fv)(2 * MINN)) return FMA_X(a, b * K(0.5), K(1.5));
    return FMA_X(a, b, K(3)) * K(0.5);
}

// ARM result for a lane flagged by FBAD (operands already flushed).
SN T FN(fix)(Cpu *c, int op, T a, T b, T d) {
    U x = FN(fb)(a), y = FN(fb)(b), z = FN(fb)(d);
    int nan = FN(isnan)(x) || FN(isnan)(y);
    switch (op) {
    case FOP_SQRT: return FN(nan1)(c, a); // NaN in, or negative -> default NaN
    case FOP_MADD:
        if (FN(isnan)(x) || FN(isnan)(y) || FN(isnan)(z)) {
            if (FN(isqnan)(x) && ((FN(isinf)(y) && FN(iszero)(z)) || (FN(iszero)(y) && FN(isinf)(z))))
                return FN(fv)(DNAN);
            return FN(nan3)(c, a, b, d);
        }
        return FN(fv)(DNAN);
    case FOP_MULX:
        if (!nan) return FN(fv)(((x ^ y) & SGN) | FN(fb)(K(2)));
        break;
    case FOP_RECPS:
        if (!nan) return K(2);
        break;
    case FOP_RSQRTS:
        if (!nan) {
            if ((FN(isinf)(x) && FN(iszero)(y)) || (FN(iszero)(x) && FN(isinf)(y))) return K(1.5);
            return FN(rsqrts)(a, b);
        }
        break;
    case FOP_NMUL: return FN(neg)(FN(nan2)(c, a, b));
    case FOP_ABD: return FN(fv)(FN(fb)(FN(nan2)(c, a, b)) & ~SGN);
    }
    return FN(nan2)(c, a, b);
}

// FMAX/FMIN/FMAXNM/FMINNM (operands already flushed)
SI T FN(minmax)(Cpu *c, int op, T a, T b) {
    U x = FN(fb)(a), y = FN(fb)(b);
    if (__builtin_expect(FN(isnan)(x) || FN(isnan)(y), 0)) {
        if (op == FOP_MAXNM || op == FOP_MINNM) {
            if (FN(isqnan)(x) && !FN(isnan)(y)) return b;
            if (FN(isqnan)(y) && !FN(isnan)(x)) return a;
        }
        return FN(nan2)(c, a, b);
    }
    if (op == FOP_MAX || op == FOP_MAXNM) return a == b ? FN(fv)(x & y) : a > b ? a : b;
    return a == b ? FN(fv)(x | y) : a < b ? a : b;
}
#define FMINMAX(op) ((op) >= FOP_MAX && (op) <= FOP_MINNM)

// FPCR.FZ / RMode path: flush inputs, round in the FPCR mode, flush tiny results (tininess is detected
// before rounding on ARM: a result rounded up to the smallest normal is recomputed toward zero).
SN T FN(slow)(Cpu *c, int op, T a, T b, T d) {
#pragma STDC FENV_ACCESS ON
    uint64_t fpcr = c->fpcr;
    if (fpcr & FPCR_FZ) a = FN(flush)(a), b = FN(flush)(b), d = FN(flush)(d);
    if (FMINMAX(op)) return FN(minmax)(c, op, a, b);
    unsigned old = simd_rm_enter(fpcr);
    volatile T va = a, vb = b, vd = d;
    volatile T vr = op == FOP_RSQRTS ? FN(rsqrts)(va, vb) : FCALC(op, va, vb, vd, K, FMA_X);
    T r = vr;
    if (FBAD(op, r)) vr = FN(fix)(c, op, a, b, d), r = vr;
    if ((fpcr & FPCR_FZ) && (FN(fb)(r) & ~SGN) <= MINN && !FN(iszero)(FN(fb)(r))) {
        if ((FN(fb)(r) & ~SGN) == MINN) {
            simd_rm_rz();
            volatile T vz = FCALC(op, va, vb, vd, K, FMA_X);
            if ((FN(fb)(vz) & ~SGN) < MINN) r = FN(fv)(FN(fb)(r) & SGN);
        } else {
            r = FN(fv)(FN(fb)(r) & SGN);
        }
    }
    simd_rm_leave(old);
    return r;
}

// Fast-path results that need a second look: FBAD lanes, and +-smallest normal. Under FZ the host runs with
// FTZ|DAZ, which flushes denormal inputs and tiny results like ARM, but x86 detects tininess after rounding:
// an exact result just below the smallest normal that rounds up to it is flushed on ARM only. The slow path
// decides both cases with the operands flushed.
#define FREDO(op, r) (FBAD(op, r) | (__builtin_elementwise_abs(r) == FN(fv)(MINN)))
// The complement of FREDO for every op but RSQRTS, as one compare: |r| ordered and not the smallest normal (LLVM folds
// it to fcmp one, a single ucomiss + jne for scalars and one vcmpps for vectors).
#define FOK(r) ((__builtin_elementwise_abs(r) < FN(fv)(MINN)) | (__builtin_elementwise_abs(r) > FN(fv)(MINN)))
SN T FN(redo)(Cpu *c, int op, T a, T b, T d, T r) {
    if (c->fpcr & (FPCR_FZ | FPCR_RM)) return FN(slow)(c, op, a, b, d);
    return FBAD(op, r) ? FN(fix)(c, op, a, b, d) : r;
}

// Scalar FP op (constant op folds the selection). RMode = nearest runs on the host, with or without FZ. A directed
// RMode is rare, so it shares the result check's branch: the host result, computed in round-to-nearest, is then
// thrown away and redo recomputes in the FPCR mode. Min and max do not round, so RMode does not matter to them.
SI T FN(op)(Cpu *c, int op, T a, T b, T d) {
    if (FMINMAX(op)) {
        if (c->fpcr & FPCR_FZ) a = FN(flush)(a), b = FN(flush)(b);
        return FN(minmax)(c, op, a, b);
    }
    T r = FCALC(op, a, b, d, K, FMA_S);
    if (__builtin_expect((op == FOP_RSQRTS ? FREDO(op, r) : !FOK(r)) | simd_directed(c), 0))
        r = FN(redo)(c, op, a, b, d, r);
    return r;
}

SN VT FN(vslow)(Cpu *c, int op, VT a, VT b, VT d) {
    VT r;
    for (int i = 0; i < NL; i++) r[i] = FN(op)(c, op, a[i], b[i], d[i]);
    return r;
}

SI VTU FN(vsel)(VTS m, VTU x, VTU y) { return (x & (VTU)m) | (y & ~(VTU)m); }

SI VT FN(vflush)(VT a) {
    VTU u = (VTU)a;
    __asm__("" : "+x"(u));
    VTU mag = u & ~((VTU){} + SGN);
    VTS den = (mag != 0) & ((u & EXPM) == 0);
    return (VT)FN(vsel)(den, u & SGN, u);
}

// Vector FP op: native lanes; any lane FREDO flags (or a directed FPCR.RMode, checked in the same branch) reruns all
// lanes through the scalar path. Min and max ignore RMode, as in FN(op).
SI VT FN(vop)(Cpu *c, int op, VT a, VT b, VT d) {
    VT r;
    if (FMINMAX(op)) {
        if (c->fpcr & FPCR_FZ) a = FN(vflush)(a), b = FN(vflush)(b);
        if (__builtin_expect(__builtin_reduce_or((a != a) | (b != b)) != 0, 0)) return FN(vslow)(c, op, a, b, d);
        VTU x = (VTU)a, y = (VTU)b;
        if (op == FOP_MAX || op == FOP_MAXNM)
            r = (VT)FN(vsel)(a > b, x, FN(vsel)(a == b, x & y, y));
        else
            r = (VT)FN(vsel)(a < b, x, FN(vsel)(a == b, x | y, y));
        return r;
    }
    r = FCALC(op, a, b, d, VK, FMA_V);
    if (__builtin_expect((op == FOP_RSQRTS ? __builtin_reduce_or(FREDO(op, r)) != 0 : __builtin_reduce_and(FOK(r)) == 0)
                         | simd_directed(c), 0))
        return FN(vslow)(c, op, a, b, d);
    return r;
}

SI VT FN(vfz)(Cpu *c, VT a) { return c->fpcr & FPCR_FZ ? FN(vflush)(a) : a; }

// ---- compares
// FPCR.FZ needs no flush here: rt_fpcr_sync (runtime/cpu.h) keeps MXCSR.DAZ equal to FPCR.FZ, and x86 compares
// under DAZ read a denormal operand as a zero of its sign, as ARM compares do under FZ.
SI uint64_t FN(cmp)(Cpu *c, T a, T b) { // FCMP NZCV
    (void)c;
    return a == b ? 0x60000000ull : a < b ? 0x80000000ull : a > b ? 0x20000000ull : 0x30000000ull;
}
// FCMEQ/FCMGE/FCMGT/FACGE/FACGT (op 0..4)
SI VTU FN(vcmp)(Cpu *c, int op, VT a, VT b) {
    (void)c;
    if (op >= 3) a = __builtin_elementwise_abs(a), b = __builtin_elementwise_abs(b);
    switch (op) {
    case 0: return (VTU)(a == b);
    case 1: case 3: return (VTU)(a >= b);
    default: return (VTU)(a > b);
    }
}

// ---- rounding to integral (FRINT*): rm 0 N, 1 P, 2 M, 3 Z, 4 A; 5 X / 6 I = FPCR mode
// (rint = ties-to-even: the host stays in round-to-nearest outside the slow paths; no roundeven in MinGW's CRT)
SI T FN(rnd)(T x, int rm) {
    switch (rm) {
    case 0: return __builtin_elementwise_rint(x);
    case 1: return __builtin_elementwise_ceil(x);
    case 2: return __builtin_elementwise_floor(x);
    case 3: return __builtin_elementwise_trunc(x);
    default: return __builtin_elementwise_round(x);
    }
}
SI T FN(frint)(Cpu *c, T x, int rm) {
    x = FN(fz)(c, x);
    if (rm >= 5) rm = (int)((c->fpcr >> 22) & 3);
    T r = FN(rnd)(x, rm);
    if (__builtin_expect(r != r, 0)) r = FN(nan1)(c, x);
    return r;
}
SI VT FN(vfrint)(Cpu *c, VT x, int rm) {
    x = FN(vfz)(c, x);
    if (rm >= 5) rm = (int)((c->fpcr >> 22) & 3);
    VT r;
    switch (rm) {
    case 0: r = __builtin_elementwise_rint(x); break;
    case 1: r = __builtin_elementwise_ceil(x); break;
    case 2: r = __builtin_elementwise_floor(x); break;
    case 3: r = __builtin_elementwise_trunc(x); break;
    default: r = __builtin_elementwise_round(x); break;
    }
    if (__builtin_expect(__builtin_reduce_or(r != r) != 0, 0))
        for (int i = 0; i < NL; i++)
            if (r[i] != r[i]) r[i] = FN(nan1)(c, x[i]);
    return r;
}

// ---- FP -> integer (FCVT[NPMZA][SU], fixed point): saturating, NaN -> 0. ib = 32/64.
SI uint64_t FN(toint)(Cpu *c, T x, int rm, int fbits, int ib, int uns) {
    x = FN(fz)(c, x);
    if (fbits) x = x * (T)(1ull << (fbits - 1)) * K(2);
    T r = FN(rnd)(x, rm);
    if (r != r) return 0;
    if (uns) {
        if (r <= K(0)) return 0;
        if (ib == 32) return r >= K(4294967296.0) ? 0xffffffffull : (uint32_t)r;
        return r >= K(18446744073709551616.0) ? ~0ull : (uint64_t)r;
    }
    if (ib == 32) return r >= K(2147483648.0) ? 0x7fffffffull : r < K(-2147483648.0) ? 0x80000000ull : (uint32_t)(int32_t)r;
    return r >= K(9223372036854775808.0) ? 0x7fffffffffffffffull : r < K(-9223372036854775808.0) ? 0x8000000000000000ull
                                                                                              : (uint64_t)(int64_t)r;
}

// ---- integer -> FP (SCVTF/UCVTF, fixed point), rounded in the FPCR mode
SN T FN(fromint_slow)(Cpu *c, uint64_t v, int ib, int uns, int fbits) {
#pragma STDC FENV_ACCESS ON
    unsigned old = simd_rm_enter(c->fpcr);
    volatile uint64_t vv = v;
    volatile T r;
    if (uns) r = ib == 32 ? (T)(uint32_t)vv : (T)vv;
    else r = ib == 32 ? (T)(int32_t)vv : (T)(int64_t)vv;
    simd_rm_leave(old);
    T x = r;
    if (fbits) x = x / (T)(1ull << (fbits - 1)) / K(2);
    return x;
}
SI T FN(fromint)(Cpu *c, uint64_t v, int ib, int uns, int fbits) {
    if (__builtin_expect(c->fpcr & FPCR_RM, 0)) return FN(fromint_slow)(c, v, ib, uns, fbits);
    T r;
    if (uns) r = ib == 32 ? (T)(uint32_t)v : (T)v;
    else r = ib == 32 ? (T)(int32_t)v : (T)(int64_t)v;
    if (fbits) r = r / (T)(1ull << (fbits - 1)) / K(2);
    return r;
}

// ---- reciprocal estimates (bit-exact ARMv8.0 pseudocode)
SI T FN(recpe)(Cpu *c, T a) {
    U x = FN(fb)(FN(fz)(c, a)), s = x & SGN;
    if (FN(isnan)(x)) return FN(nan1)(c, a);
    if (FN(isinf)(x)) return FN(fv)(s);
    if (FN(iszero)(x)) return FN(fv)(s | EXPM);
    U mag = x & ~SGN;
    if (mag < ((U)1 << (MB - 2))) { // |x| < 2^(emin - 2): result overflows
        int rm = (c->fpcr >> 22) & 3;
        int inf = rm == 0 || (rm == 1 && !s) || (rm == 2 && s);
        return FN(fv)(s | (inf ? EXPM : EXPM - 1));
    }
    int64_t e = (int64_t)(mag >> MB);
    U frac = mag & (((U)1 << MB) - 1);
    if ((c->fpcr & FPCR_FZ) && e >= (int64_t)((1 << EB) - 3)) return FN(fv)(s); // result would be tiny
    if (e == 0) {
        if (!(frac >> (MB - 1))) e = -1, frac = (frac << 2) & (((U)1 << MB) - 1);
        else frac = (frac << 1) & (((U)1 << MB) - 1);
    }
    int scaled = 256 | (int)(frac >> (MB - 8));
    int64_t re = (int64_t)((1 << EB) - 3) - e;
    U est = (U)simd_recip_tab[scaled - 256]; // 8 fraction bits
    U f = est << (MB - 8);
    if (re == 0) f = (f >> 1) | ((U)1 << (MB - 1));
    else if (re == -1) f = (f >> 2) | ((U)1 << (MB - 2)), re = 0;
    return FN(fv)(s | ((U)re << MB) | f);
}
SI T FN(rsqrte)(Cpu *c, T a) {
    U x = FN(fb)(FN(fz)(c, a)), s = x & SGN;
    if (FN(isnan)(x)) return FN(nan1)(c, a);
    if (FN(iszero)(x)) return FN(fv)(s | EXPM);
    if (s) return FN(fv)(DNAN);
    if (FN(isinf)(x)) return K(0);
    int64_t e = (int64_t)(x >> MB);
    U frac = x & (((U)1 << MB) - 1);
    if (e == 0) {
        while (!(frac >> (MB - 1))) frac <<= 1, e--;
        frac = (frac << 1) & (((U)1 << MB) - 1);
    }
    int scaled = (e & 1) ? 128 | (int)(frac >> (MB - 7)) : 256 | (int)(frac >> (MB - 8));
    int64_t re = ((int64_t)(3 * ((1 << (EB - 1)) - 1) - 1) - e) / 2; // (380 - e) / 2, (3068 - e) / 2
    U est = (U)simd_rsqrt_tab[scaled - 128];
    return FN(fv)(((U)re << MB) | (est << (MB - 8)));
}
SI T FN(recpx)(Cpu *c, T a) {
    U x = FN(fb)(FN(fz)(c, a)), s = x & SGN;
    if (FN(isnan)(x)) return FN(nan1)(c, a);
    if (!(x & EXPM)) return FN(fv)(s | (EXPM - MINN));
    return FN(fv)(s | (~x & EXPM));
}

#undef FN
#undef SGN
#undef EXPM
#undef QB
#undef DNAN
#undef MINN
#undef VK
#undef K
#undef FMA_S
#undef FMA_V
#undef FMA_X