Compiler Explorer

Source code

#include <math.h>
#include <riscv_vector.h>

#define vec1 vfloat32m1_t
#define vec2 vfloat32m1x2_t
#define vec3 vfloat32m1x3_t
#define vec4 vfloat32m1x4_t
#define vec1b vbool32_t

#define getx(v) __riscv_vget_f32m1(v,0)
#define gety(v) __riscv_vget_f32m1(v,1)
#define getz(v) __riscv_vget_f32m1(v,2)
#define getw(v) __riscv_vget_f32m1(v,3)

#define getxy(v) vec(getx(v), gety(v))
#define getxz(v) vec(getx(v), getz(v))
#define getyz(v) vec(gety(v), getz(v))
#define getyx(v) vec(gety(v), getx(v))

#define getxzy(v) vec(getx(v), getz(v), gety(v))
#define getyzx(v) vec(gety(v), getz(v), getx(v))
#define getzxy(v) vec(getz(v), getx(v), gety(v))

#define VL __riscv_vsetvlmax_e32m1()

inline vec2 vec(vec1 x, vec1 y) { return __riscv_vcreate_v_f32m1x2(x, y); }
inline vec3 vec(vec1 x, vec1 y, vec1 z) { return __riscv_vcreate_v_f32m1x3(x, y, z); }
inline vec3 vec(vec2 xy, vec1 z) { return __riscv_vcreate_v_f32m1x3(getx(xy), gety(xy), z); }
inline vec3 vec(vec1 x, vec2 yz) { return __riscv_vcreate_v_f32m1x3(x, getx(yz), gety(yz)); }
inline vec4 vec(vec1 x, vec1 y, vec1 z, vec1 w) { return __riscv_vcreate_v_f32m1x4(x, y, z, w); }

inline vec2 dup2(vec1 x) { return vec(x, x); }
inline vec3 dup3(vec1 x) { return vec(x, x, x); }
inline vec4 dup4(vec1 x) { return vec(x, x, x, x); }

inline vec1 vec(float x) { return __riscv_vfmv_v_f_f32m1(x, VL); }
inline vec2 vec(float x, float y) { return vec(vec(x), vec(y)); }
inline vec3 vec(float x, float y, float z) { return vec(vec(x), vec(y), vec(z)); }
inline vec4 vec(float x, float y, float z, float w) { return vec(vec(x), vec(y), vec(z), vec(w)); }

inline vec2 vec(vec1 x, float y) { return vec(x, vec(y)); }

struct float2 { const float x, y; constexpr float2(float x, float y) : x(x), y(y) {} };
struct float3 { const float x, y, z; constexpr float3(float x, float y, float z) : x(x), y(y), z(z) {} };
struct float4 { const float x, y, z, w; constexpr float4(float x, float y, float z, float w) : x(x), y(y), z(z), w(w) {} };

#define OP2(op) \
	inline float2 operator op(float a, float2 b) { return float2(a op b.x, a op b.y); } \
	inline float2 operator op(float2 a, float b) { return float2(a.x op b, a.y op b); } \
	inline float2 operator op(float2 a, float2 b) { return float2(a.x op b.x, a.y op b.y); } \
	inline float3 operator op(float a, float3 b) { return float3(a op b.x, a op b.y, a op b.z); } \
	inline float3 operator op(float3 a, float b) { return float3(a.x op b, a.y op b, a.z op b); } \
	inline float3 operator op(float3 a, float3 b) { return float3(a.x op b.x, a.y op b.y, a.z op b.z); }

OP2(+) OP2(-) OP2(/) OP2(*)

#define ALIAS1(n, f) \
	inline auto n(vec1 v) { return __riscv_##f(v, VL); } \
	inline auto n(vec2 v) { return vec(n(getx(v)), n(gety(v))); } \
	inline auto n(vec3 v) { return vec(n(getx(v)), n(gety(v)), n(getz(v))); } \
	inline auto n(vec4 v) { return vec(n(getx(v)), n(gety(v)), n(getz(v)), n(getw(v))); }

#define ALIAS2(n, f) \
	inline auto n(vec1 a, vec1 b) { return __riscv_##f(a, b, VL); } \
	inline auto n(vec2 a, vec2 b) { return vec(n(getx(a),getx(b)), n(gety(a),gety(b))); } \
	inline auto n(vec3 a, vec3 b) { return vec(n(getx(a),getx(b)), n(gety(a),gety(b)), n(getz(a),getz(b))); } \
	inline auto n(vec4 a, vec4 b) { return vec(n(getx(a),getx(b)), n(gety(a),gety(b)), n(getz(a),getz(b)), n(getw(a),getw(b))); }

#define ALIAS2x(n, f) ALIAS2(n, f) \
	inline auto n(vec1 a, float b) { return __riscv_##f(a, b, VL); } \
	inline auto n(vec2 a, float b) { return vec(n(getx(a),b), n(gety(a),b)); } \
	inline auto n(vec3 a, float b) { return vec(n(getx(a),b), n(gety(a),b), n(getz(a),b)); } \
	inline auto n(vec4 a, float b) { return vec(n(getx(a),b), n(gety(a),b), n(getz(a),b), n(getw(a),b)); } \
	inline auto n(vec2 a, float2 b) { return vec(n(getx(a),b.x), n(gety(a),b.y)); } \
	inline auto n(vec3 a, float3 b) { return vec(n(getx(a),b.x), n(gety(a),b.y), n(getz(a),b.z)); } \
	inline auto n(vec4 a, float4 b) { return vec(n(getx(a),b.x), n(gety(a),b.y), n(getz(a),b.z), n(getw(a),b.w)); }

#define ALIAS3(n, f) \
	inline auto n(vec1 a, vec1 b, vec1 c) { return __riscv_##f(a, b, c, VL); } \
	inline auto n(vec2 a, vec2 b, vec2 c) { return vec(n(getx(a),getx(b),getx(c)), n(gety(a),gety(b),gety(c))); } \
	inline auto n(vec3 a, vec3 b, vec3 c) { return vec(n(getx(a),getx(b),getx(c)), n(gety(a),gety(b),gety(c)), n(getz(a),getz(b),getz(c))); } \
	inline auto n(vec4 a, vec4 b, vec4 c) { return vec(n(getx(a),getx(b),getx(c)), n(gety(a),gety(b),gety(c)), n(getz(a),getz(b),getz(c)), n(getw(a),getw(b),getw(c))); }

#define ALIAS3x(n, f) ALIAS3(n, f) \
	inline auto n(vec1 a, float b, vec1 c) { return __riscv_##f(a, b, c, VL); } \
	inline auto n(vec2 a, float b, vec2 c) { return vec(n(getx(a),b,getx(c)), n(gety(a),b,gety(c))); } \
	inline auto n(vec3 a, float b, vec3 c) { return vec(n(getx(a),b,getx(c)), n(gety(a),b,gety(c)), n(getz(a),b,getz(c))); } \
	inline auto n(vec4 a, float b, vec4 c) { return vec(n(getx(a),b,getx(c)), n(gety(a),b,gety(c)), n(getz(a),b,getz(c)), n(getw(a),b,getw(c))); } \
	inline auto n(vec2 a, float2 b, vec2 c) { return vec(n(getx(a),b.x,getx(c)), n(gety(a),b.y,gety(c))); } \
	inline auto n(vec3 a, float3 b, vec3 c) { return vec(n(getx(a),b.x,getx(c)), n(gety(a),b.y,gety(c)), n(getz(a),b.z,getz(c))); } \
	inline auto n(vec4 a, float4 b, vec4 c) { return vec(n(getx(a),b.x,getx(c)), n(gety(a),b.y,gety(c)), n(getz(a),b.z,getz(c)), n(getw(a),b.w,getw(c))); }

ALIAS1(abs, vfabs)
ALIAS1(neg, vfneg)
ALIAS1(sqrt, vfsqrt)
ALIAS2x(div, vfdiv)
ALIAS2x(add, vfadd)
ALIAS2x(sub, vfsub)
ALIAS2x(mul, vfmul)
ALIAS2x(max, vfmax)
ALIAS2x(min, vfmin)
ALIAS2x(sgnj, vfsgnj)
ALIAS2x(sgnjn, vfsgnjn)

ALIAS3x(macc, vfmacc)
ALIAS3x(nmacc, vfnmacc)
ALIAS3x(madd, vfmadd)
ALIAS3x(msub, vfmsub)
ALIAS3x(nmsac, vfnmsac)

inline vec1 rdiv(vec1 a, float b) { return __riscv_vfrdiv(a, b, VL); };
inline vec2 rdiv(vec2 a, float b) { return vec(rdiv(getx(a), b), rdiv(gety(a), b)); };
inline vec3 rdiv(vec3 a, float b) { return vec(rdiv(getx(a), b), rdiv(gety(a), b), rdiv(getz(a), b)); };

inline vec1 rsub(vec1 a, float b) { return __riscv_vfrsub(a, b, VL); };
inline vec2 rsub(vec2 a, float b) { return vec(rsub(getx(a), b), rsub(gety(a), b)); };
inline vec3 rsub(vec3 a, float b) { return vec(rsub(getx(a), b), rsub(gety(a), b), rsub(getz(a), b)); };
inline vec2 rsub(vec2 a, float2 b) { return vec(rsub(getx(a), b.x), rsub(gety(a), b.y)); };
inline vec3 rsub(vec3 a, float3 b) { return vec(rsub(getx(a), b.x), rsub(gety(a), b.y), rsub(getz(a), b.z)); };

#define CMP(f) \
	inline auto f(vec1 a, vec1 b) { return __riscv_vmf##f(a, b, VL); } \
	inline auto f(vec1 a, float b) { return __riscv_vmf##f(a, b, VL); } 
	
CMP(eq) CMP(ne) CMP(lt) CMP(le) CMP(gt) CMP(ge)

inline vec1 vif(vec1b m, vec1 t, vec1 f) { return __riscv_vmerge(f, t, m, VL); }
inline vec2 vif(vec1b m, vec2 t, vec2 f) { return vec(vif(m, getx(t), getx(f)), vif(m, gety(t), gety(f))); }
inline vec3 vif(vec1b m, vec3 t, vec3 f) { return vec(vif(m, getx(t), getx(f)), vif(m, gety(t), gety(f)), vif(m, getz(t), getz(f))); }
inline vec4 vif(vec1b m, vec4 t, vec4 f) { return vec(vif(m, getx(t), getx(f)), vif(m, gety(t), gety(f)), vif(m, getz(t), getz(f)), vif(m, getw(t), getw(f))); }

inline vec1 vif(vec1b m, float t, vec1 f) { return __riscv_vfmerge(f, t, m, VL); }
inline vec2 vif(vec1b m, float2 t, vec2 f) { return vec(vif(m, t.x, getx(f)), vif(m, t.x, gety(f))); }
inline vec3 vif(vec1b m, float3 t, vec3 f) { return vec(vif(m, t.x, getx(f)), vif(m, t.y, gety(f)), vif(m, t.z, getz(f))); }
inline vec4 vif(vec1b m, float4 t, vec4 f) { return vec(vif(m, t.x, getx(f)), vif(m, t.y, gety(f)), vif(m, t.z, getz(f)), vif(m, t.w, getw(f))); }

inline vec1b mor(vec1b a, vec1b b) { return __riscv_vmor(a, b, VL); }
inline vec1b mand(vec1b a, vec1b b) { return __riscv_vmand(a, b, VL); }
inline vec1b mnor(vec1b a, vec1b b) { return __riscv_vmnor(a, b, VL); }
inline vec1b mnand(vec1b a, vec1b b) { return __riscv_vmnand(a, b, VL); }
inline bool all(vec1b m) { return __riscv_vcpop(m, VL) == VL; }
inline bool any(vec1b m) { return __riscv_vfirst(m, VL) >= 0; }
inline bool none(vec1b m) { return __riscv_vfirst(m, VL) < 0; }

inline vec1
dot(vec2 a, vec2 b)
{
	return macc(mul(
		getx(a), getx(b)),
		gety(a), gety(b));
}

inline vec1
dot(vec3 a, vec3 b)
{
	return macc(macc(mul(
		getx(a), getx(b)),
		gety(a), gety(b)),
		getz(a), getz(b));
}

inline vec1
dot(vec2 a, float2 b)
{
	return macc(mul(
		getx(a), b.x),
		b.y, gety(a));
}

inline vec1
dot(vec3 a, float3 b)
{
	return macc(macc(mul(
		getx(a), b.x),
		b.y, gety(a)),
		b.z, getz(a));
}

constexpr float dot(float2 a, float2 b) { return a.x*b.x + a.y*b.y; }
constexpr float dot(float3 a, float3 b) { return a.x*b.x + a.y*b.y + a.z*b.z; }

inline vec1 dot2(vec2 v) { return dot(v,v); }
inline vec1 dot2(vec3 v) { return dot(v,v); }

inline float dot2(float2 v) { return dot(v,v); }
inline float dot2(float3 v) { return dot(v,v); }

inline vec1 length(vec2 v) { return sqrt(dot2(v)); }
inline vec1 length(vec3 v) { return sqrt(dot2(v)); }

inline float length(float2 v) { return sqrt(dot2(v)); }
inline float length(float3 v) { return sqrt(dot2(v)); }

inline vec2 normalize(vec2 v) { return div(v,dup2(length(v))); }
inline vec3 normalize(vec3 v) { return div(v,dup3(length(v))); }

inline float2 normalize(float2 v) { return v/length(v); }
inline float3 normalize(float3 v) { return v/length(v); }

inline vec1 floor(vec1 v) { return __riscv_vfcvt_f(__riscv_vfcvt_x(v,VL),VL); }
inline vec2 floor(vec2 v) { return vec(floor(getx(v)), floor(gety(v))); }
inline vec3 floor(vec3 v) { return vec(floor(getx(v)), floor(gety(v)), floor(getz(v))); }

inline vec1 fract(vec1 v) { return sub(v, floor(v)); }
inline vec2 fract(vec2 v) { return sub(v, floor(v)); }
inline vec3 fract(vec3 v) { return sub(v, floor(v)); }

inline vec3 reflect(vec3 I, vec3 N) { return macc(I, mul(N,-2), dup3(dot(N,I))); }

inline vec1 ndot(vec2 a, vec2 b) { return msub(getx(a), getx(b), mul(gety(a), gety(b))); }
inline vec1 ndot(vec2 a, float2 b) { return msub(getx(a), b.x, mul(gety(a), b.y)); }

#define clamp(x,minv,maxv)  min(max(x, minv), maxv)
inline float sign(float x) { return x == 0 ? 0 : x < 0 ? -1 : 1; }

inline vec1 mix(vec1 x, vec1 y, vec1 a) { return macc(mul(x,rsub(a,1)),y,a); }
inline vec2 mix(vec2 x, vec2 y, vec2 a) { return macc(mul(x,rsub(a,1)),y,a); }
inline vec3 mix(vec3 x, vec3 y, vec3 a) { return macc(mul(x,rsub(a,1)),y,a); }

inline vec1
smoothstep(float a, float b, vec1 v)
{
	v = clamp(div(sub(v, a), b-a), 0, 1);
	return mul(mul(v,v), rsub(add(v,v),3));
}

inline vec1 pow2(vec1 v) { return mul(v, v); }
inline vec1 pow3(vec1 v) { return mul(pow2(v),v); }
inline vec1 pow4(vec1 v) { return pow2(pow2(v)); }
inline vec1 pow5(vec1 v) { vec1 x2 = pow2(v), x3 = mul(x2, v); return mul(x2, x3); }
inline vec1 pow8(vec1 v) { return pow4(pow4(v)); }
inline vec1 pow16(vec1 v) { return pow8(pow8(v)); }

inline vec1 rem(vec1 a, float b) { return nmsac(a, b, floor(div(a, b))); }

inline vec1
sin(vec1 x)
{
	x = sub(rem(add(x,M_PI/2),M_PI),M_PI/2);
	vec1 sign = sub(rem(sub(mul(x,1/M_PI),0.5),2),1);

vec1 x2 = mul(x, x);
	vec1 x3 = mul(x2, x);
	vec1 x5 = mul(x2, x3);
	return sgnj(macc(macc(x, -1./6, x3), 1./120, x5), sign);
}
inline vec2 sin(vec2 v) { return vec(sin(getx(v)), sin(gety(v))); }
inline vec3 sin(vec3 v) { return vec(sin(getx(v)), sin(gety(v)), sin(getz(v))); }

inline vec1 cos(vec1 v) { return sin(add(v,M_PI/2)); }
inline vec2 cos(vec2 v) { return sin(add(v,M_PI/2)); }
inline vec3 cos(vec3 v) { return sin(add(v,M_PI/2)); }

inline vec1
exp(vec1 v)
{
	constexpr float Mln2 = 0.6931471805f;
	constexpr float A = 8388608.0f;
	constexpr float B = 1065353216.0f;
	constexpr float C = 60801.0f;
	return __riscv_vreinterpret_f32m1(__riscv_vfcvt_x(madd(v, A/Mln2, vec(B-C)), VL));
}
inline vec2 exp(vec2 v) { return vec(exp(getx(v)), exp(gety(v))); }
inline vec3 exp(vec3 v) { return vec(exp(getx(v)), exp(gety(v)), exp(getz(v))); }

inline vec2
mat2v(vec2 col1, vec2 col2, vec2 v)
{
	return vec(
		macc(mul(getx(col1),getx(v)), getx(col2), gety(v)),
		macc(mul(gety(col1),getx(v)), gety(col2), gety(v)));
}

inline vec3
mat3v(vec3 col1, vec3 col2, vec3 col3, vec3 v)
{
	return vec(
		macc(macc(mul(getx(col1),getx(v)), getx(col2), gety(v)), getx(col3), getz(v)),
		macc(macc(mul(gety(col1),getx(v)), gety(col2), gety(v)), gety(col3), getz(v)),
		macc(macc(mul(getz(col1),getx(v)), getz(col2), gety(v)), getz(col3), getz(v)));
}

inline vec3
cross(vec3 a, vec3 b)
{
	return vec(
		nmsac(mul(gety(a), getz(b)), getz(a), gety(b)),
		nmsac(mul(getz(a), getx(b)), getx(a), getz(b)),
		nmsac(mul(getx(a), gety(b)), gety(a), getx(b)));
}

inline vec1 sdPlane(vec3 p) { return gety(p); }
inline vec1 sdSphere(vec3 p, float s) { return sub(length(p), s); }

inline vec1
sdBox(vec3 p, float3 b)
{
	vec3 d = sub(abs(p), b);
	return add(min(max(getx(d),max(gety(d),getz(d))),0), length(max(d,0)));
}

inline vec1
sdBoxFrame(vec3 p, float3 b, float e)
{
	p = sub(abs(p),b);
	vec3 q = sub(abs(add(p,e)),e);

return min(min(
	       add(length(max(vec(getx(p),gety(q),getz(q)),0.0)),min(max(getx(p),max(gety(q),getz(q))),0.0)),
	       add(length(max(vec(getx(q),gety(p),getz(q)),0.0)),min(max(getx(q),max(gety(p),getz(q))),0.0))),
	       add(length(max(vec(getx(q),gety(q),getz(p)),0.0)),min(max(getx(q),max(gety(q),getz(p))),0.0)));

}

inline vec1
sdEllipsoid(vec3 p, float3 r)
{
	vec1 k0 = length(mul(p, 1/r)) ;
	vec1 k1 = length(mul(p, 1/(r*r)));
	return div(mul(k0, sub(k0, 1)), k1);
}

inline vec1
sdTorus(vec3 p, float2 t)
{
	return sub(length(vec(sub(length(getxz(p)), t.x),gety(p))), t.y);
}

inline vec1
sdCappedTorus(vec3 p, float2 sc, float ra, float rb)
{
	p = vec(abs(getx(p)), getyz(p));
	vec1 k = vif(gt(mul(getx(p), sc.y), mul(gety(p), sc.x)), dot(getxy(p),sc), length(getxy(p)));
	return sub(sqrt( sub(add(dot2(p), ra*ra), mul(k, 2*ra))), rb);
}

inline vec1
sdHexPrism(vec3 p, float2 h)
{
	constexpr float2 kxy(-0.8660254, 0.5);
	constexpr float kz = 0.57735;
	p = abs(p);
	vec2 pxy = getxy(p);
	p = vec(macc(pxy, kxy*-2, dup2(min(dot(pxy, kxy), 0))), getz(p));
	vec2 d = vec(
		sgnj(length(sub(pxy, vec(clamp(getx(p), -kz*h.x, kz*h.x), h.x))), sub(gety(p), h.x) ),
 		sub(getz(p),h.y));
	return add(min(max(getx(d),gety(d)),0.0), length(max(d,0.0)));
}

inline vec1
sdOctogonPrism(vec3 p, float r, float h)
{
	constexpr float3 k(-0.9238795325,   // sqrt(2+sqrt(2))/2 
	                    0.3826834323,   // sqrt(2-sqrt(2))/2
	                    0.4142135623 ); // sqrt(2)-1 
	// reflections
	p = abs(p);
	vec2 pxy = getxy(p);
	pxy = macc(pxy, -2.0, mul(dup2(min(dot(vec( k.x,k.y),pxy),0.0)), vec( k.x,k.y)));
	pxy = macc(pxy, -2.0, mul(dup2(min(dot(vec(-k.x,k.y),pxy),0.0)), vec(-k.x,k.y)));
	// polygon side
	pxy = sub(pxy, vec(clamp(getx(p), -k.z*r, k.z*r), r));
	vec2 d = vec( sgnj(length(pxy),gety(pxy)), sub(getz(p),h) );
	return add( min(max(getx(d),gety(d)),0.0), length(max(d,0.0)) );
}

inline vec1
sdCapsule(vec3 p, float3 a, float3 b, float r)
{
	vec3 pa = sub(p, a);
	float3 ba = float3(b.x-a.x, b.y-a.y, b.z-a.z);
	vec1 h = clamp( mul(dot(pa,ba),1/dot2(ba)), 0.0, 1.0 );
	return sub(length(nmsac(pa, ba, dup3(h))), r);
}

inline vec1
sdRoundCone(vec3 p, float r1, float r2, float h )
{
	vec2 q = vec(length(getxz(p)), gety(p) );
    
	float b = (r1-r2)/h;
	float a = std::sqrt(1.0-b*b);
	vec1 k = dot(q,float2(-b,a));

vec1 ret = sub( dot(q, float2(a,b) ), r1);
	ret = vif(lt(k, 0),   sub(length(q), r1), ret);
	ret = vif(gt(k, a*h), sub(length(sub(q, float2(0.0,h))), r2), ret);
        
	return ret;
}

static vec1
sdRoundCone(vec3 p, float3 a, float3 b, float r1, float r2)
{
 	// sampling independent computations (only depend on shape)
 	float3  ba = b - a;
	float l2 = dot2(ba);
	float rr = r1 - r2;
	float a2 = l2 - rr*rr;
	float il2 = 1.0/l2;

// sampling dependant computations
	vec3 pa = sub(p, a);
	vec1 y = dot(pa,ba);
	vec1 z = sub(y, l2);
	vec1 x2 = dot2( msub(pa, l2, mul(dup3(y), ba)));
	vec1 y2 = mul(mul(y, y), l2);
	vec1 z2 = mul(mul(z, z), l2);

// single square root!
	vec1 k = sgnj(mul(x2, rr*rr), rr);
	vec1b c1 = gt(sgnj(mul(z2, a2), z), k);
	vec1b c2 = lt(sgnj(mul(y2, a2), y), k);
	vec1b c3 = mnor(c1, c2);
	vec1 sq = vif(c3, mul(mul(x2,il2),a2), add(x2,vif(c1, z2, y2)));
	vec1 ret = sqrt(sq);
	return sub(mul( vif(c3, macc(ret, rr, y), ret), il2), vif(c1, r1, vec(r2)));
}

inline vec1
sdTriPrism(vec3 p, float2 h)
{
	float k = std::sqrt(3.0);
	float hx = h.x*0.5*k;
	vec1 px = mul(getx(p), (1.0f/hx));
	vec1 py = mul(gety(p), (1.0f/hx));
	px = sub(abs(px), 1);
	py = add(py, 1.0/k);
	vec1b m = gt(macc(px,k,py), 0);
	px = vif(m, mul(nmsac(px,k,py), 0.5), px);
	py = vif(m, mul(nmacc(py,k,px), 0.5), px);
	px = sub(px, clamp(px, -2.0, 0.0));
	vec1 d1 = mul(length(vec(px, py)), sgnjn(vec(hx), py));
	vec1 d2 = sub(abs(getz(p)), h.y);
	return add( length(max(vec(d1,d2),0)), min(max(d1,d2),0) );
}

// vertical
inline vec1
sdCylinder(vec3 p, float2 h)
{
	vec2 d = sub(abs(vec(length(getxz(p)),gety(p))),h);
	return add( min(max(getx(d),gety(d)),0.0), length(max(d,0.0)) );
}

// arbitrary orientation
inline vec1
sdCylinder(vec3 p, float3 a, float3 b, float r)
{
	vec3 pa = sub(p, a);
 	float3 ba = b - a;
	float baba = dot2(ba);
	vec1 paba = dot(pa,ba);

vec1 x = sub(length(nmsac(mul(pa, baba), ba, dup3(paba))), r*baba);
	vec1 y = sub(abs(sub(paba, baba*0.5f)), baba*0.5f);
	vec1 x2 = mul(x,x);
	vec1 y2 = mul(mul(y,y),baba);
	vec1 d = vif(lt(max(x,y),0), neg(min(x2,y2)), add(x2, y2));
	return mul(sgnj(sqrt(abs(d)), d), 1/baba);
}

// vertical
inline vec1
sdCone(vec3 p, float2 c, float h)
{
	float2 q = h*float2(c.x,-c.y)/c.y;
	vec2 w = vec( length(getxz(p)), gety(p) );
	vec2 a = nmsac(w, q, dup2(clamp( mul(dot(w,q),1/dot2(q)), 0.0, 1.0 )));
	vec2 b = nmsac(w, q, vec( clamp( mul(getx(w),1/q.x), 0.0, 1.0 ), 1.0 ));
	float k = q.y;
	vec1 d = min(dot2(a), dot2(b));
	vec1 s = max( sgnj(nmsac(mul(getx(w), q.y), q.x, gety(w)), k), sgnj(sub(gety(w), q.y), k)  );
	return sgnj(sqrt(d),s);
}

inline vec1
sdCappedCone(vec3 p, float h, float r1, float r2 )
{
	vec2 q = vec( length(getxz(p)), gety(p) );
	
	float2 k1 = float2(r2,h);
	float2 k2 = float2(r2-r1,2.0*h);
	vec2 ca = vec( sub(getx(q), min(getx(q),vif(lt(gety(q), 0),r1,vec(r2)))), sub(abs(gety(q)),h) );
	vec2 cb = macc(sub(q, k1), k2, dup2(clamp( mul(dot(rsub(q,k1),k2),1/dot2(k2)), 0.0, 1.0 )));
	return sgnj(sqrt( min(dot2(ca),dot2(cb)) ), max(getx(cb), gety(ca)));
}

inline vec1
sdCappedCone(vec3 p, float3 a, float3 b, float ra, float rb)
{
	float rba  = rb-ra;
	float baba = dot(b-a,b-a);
	vec1 papa = dot2(sub(p,a));
	vec1 paba = mul(dot(sub(p,a),b-a),1/baba);

vec1 x = sqrt( nmsac(papa, baba, mul(paba, paba)) );

vec1 cax = max(sub(x, vif(lt(paba,0.5),ra,vec(rb))), 0);
	vec1 cay = sub(abs(sub(paba,0.5)),0.5);

float k = rba*rba + baba;
	vec1 f = clamp( mul(macc(mul(sub(x,ra),rba), baba, paba), 1/k), 0.0, 1.0 );

vec1 cbx = nmsac(sub(x,ra), rba, f);
	vec1 cby = sub(paba, f);

return sgnj(sqrt( min(macc(mul(cax,cax), baba, mul(cay,cay)),
	                      macc(mul(cbx,cbx), baba, mul(cby,cby))) ),max(cbx,cay));
}

// c is the sin/cos of the desired cone angle
inline vec1
sdSolidAngle(vec3 pos, float2 c, float ra)
{
	vec2 p = vec(length(getxz(pos)), gety(pos));
	vec1 l = sub(length(p), ra);
	vec1 m = length(sub(p, mul(dup2(clamp(dot(p,c),0.0,ra)), c)));
	return max(l, sgnj(m, sub(mul(getx(p),c.y), mul(gety(p),c.x))));
}

inline vec1
sdOctahedron(vec3 p, float s)
{
	p = abs(p);
	vec1 m = add(getx(p), add(gety(p), sub(getz(p), s)));
    
	// exact distance
	vec1 mo3 = mul(m,1.0f/3);
	vec1b m1 = lt(getx(p), mo3);
	vec1b m2 = lt(gety(p), mo3);
	vec1b m3 = lt(getz(p), mo3);
	p = vif(m1, p, vif(m2, getyzx(p), getzxy(p)));

vec1 k = clamp(mul(add(sub(getz(p),gety(p)),s), 0.5),0.0,s); 
	vec1 ret = length(vec(getx(p), add(sub(gety(p),s),k), sub(getz(p),k))); 
	return vif(mnor(mor(m1, m2), m3), mul(m,0.57735027f), ret);
}

inline vec1
sdPyramid(vec3 p, float h)
{
	float m2 = h*h + 0.25;

// symmetry
	vec2 pxz = abs(getxz(p));
	pxz = vif(gt(gety(pxz), getx(pxz)), getyx(pxz), getxy(pxz));
	pxz = sub(pxz,0.5);
	p = vec(getx(pxz), gety(p), gety(pxz));

// project into face plane (2D)
	vec3 q = vec(getz(p), msub(gety(p), h, mul(getx(p), 0.5)), madd(getx(p), h, mul(gety(p),0.5)));

vec1 s = max(neg(getx(q)),0.0);
	vec1 t = clamp(mul(nmsac(gety(q), 0.5, getz(p)),1/(m2+0.25)), 0.0, 1.0 );

vec1 a = madd(mul(add(getx(q),s),add(getx(q),s)), m2, mul(gety(q),gety(q)));
	vec1 b = madd(mul(macc(getx(q),0.5,t),macc(getx(q),0.5,t)), m2, mul(nmsac(gety(q),m2, t),nmsac(gety(q),m2,t)));

vec1 d2 = min(gety(q), vif(gt(madd(getx(q),-m2, mul(gety(q),-0.5)),0), 0 , min(a,b)));

// recover 3D and scale, and add sign
	return sqrt(sgnj(mul(macc(d2, getz(q),getz(q)),1/m2), max(getz(q),neg(gety(p)))));
}

// la,lb=semi axis, h=height, ra=corner
inline vec1
sdRhombus(vec3 p, float la, float lb, float h, float ra)
{
	p = abs(p);
	float2 b = float2(la,lb);
	vec1 f = clamp( mul(ndot(sub(mul(getxz(p), 2), b), b),dot(b,b)), -1.0, 1.0 );
	vec2 q = vec( sub(sgnj(length(nmsac(getxz(p),0.5*b,vec(rsub(f,1),add(f,1)))), sub(macc(mul(getx(p),b.y),b.x,getz(p)), b.x*b.y)),ra), sub(gety(p),h));
	return add(min(max(getx(q),gety(q)),0.0), length(max(q,0.0)));
}

inline vec1
sdHorseshoe(vec3 p, float2 c, float r, float le, float2 w)
{
	vec2 pxy = vec(abs(getx(p)), gety(p));
	vec1 pz = getz(p);
	vec1 l = length(pxy);
	pxy = mat2v(vec(-c.x, c.y), 
	            vec(c.y, c.x), pxy);
	pxy = vec(vif(mor(gt(gety(pxy),0.0),gt(getx(pxy),0.0)),getx(pxy),sgnj(l, -c.x)),
	          vif(gt(getx(pxy),0.0),gety(pxy),l) );
	pxy = sub(vec(getx(pxy),abs(sub(gety(pxy),r))),float2(le,0));

vec2 q = vec(add(length(max(pxy,0.0)), min(max(getx(pxy),gety(pxy)),0)),pz);
	vec2 d = sub(abs(q), w);
	return add(min(max(getx(d),gety(d)),0.0), length(max(d,0.0)));
}

inline vec2 opU(vec2 d1, vec2 d2) { return vif(lt(getx(d1), getx(d2)), d1, d2); }

[[riscv::vector_cc]] static
//[[clang::always_inline]]
inline
vec2
map(vec3 pos)
{
	vec1 px = getx(pos);
	vec1 py = gety(pos);
	vec1 pz = getz(pos);
	vec2 res = vec(py, vec(0));

// bounding box
	if (all(mor(mor(lt(px, 0), gt(py, 1.2f)), mor(gt(abs(px), 5), gt(abs(sub(pz,1)), 5))))) return res;

// bounding box
	if (any(mand(lt(abs(add(px,2)), 0.6f), lt(abs(sub(pz,0.5f)), 1.2f)))) {
		res = opU(res, vec(sdSphere(        sub(pos, float3(-2.0,0.25,0.0)), 0.25), 26.9));
		res = opU(res, vec(sdRhombus(getxzy(sub(pos, float3(-2.0,0.25,1.0))), 0.15, 0.25, 0.04, 0.08), 17.0));
	}

// bounding box
	if (any(lt(abs(px), 0.6))) {
		res = opU(res, vec(sdCappedTorus(mul(sub(pos,float3( 0.0,0.30, 1.0)), float3(1,-1,1)), float2(0.866025,-0.5), 0.25, 0.05), 25.0) );
		res = opU(res, vec(sdBoxFrame(       sub(pos,float3( 0.0,0.25, 0.0)), float3(0.3,0.25,0.2), 0.025 ), 16.9 ) );
		res = opU(res, vec(sdCone(           sub(pos,float3( 0.0,0.45,-1.0)), float2(0.6,0.8),0.45 ), 55.0 ) );
		res = opU(res, vec(sdCappedCone(     sub(pos,float3( 0.0,0.25,-2.0)), 0.25, 0.25, 0.1 ), 13.67 ) );
		res = opU(res, vec(sdSolidAngle(     sub(pos,float3( 0.0,0.00,-3.0)), float2(3,4)/5.0, 0.4 ), 49.13 ) );
	}

// bounding box
	if (any(lt(abs(sub(px,1)), 0.6))) {
		res = opU(res, vec(sdTorus(getxzy(sub(pos,float3( 1.0,0.30, 1.0))), float2(0.25,0.05) ), 7.1));
		res = opU(res, vec(sdBox(         sub(pos,float3( 1.0,0.25, 0.0)), float3(0.3,0.25,0.1) ), 3.0));
		res = opU(res, vec(sdCapsule(     sub(pos,float3( 1.0,0.00,-1.0)), float3(-0.1,0.1,-0.1), float3(0.2,0.4,0.2), 0.1 ), 31.9));
		res = opU(res, vec(sdCylinder(    sub(pos,float3( 1.0,0.25,-2.0)), float2(0.15,0.25) ), 8.0));
		res = opU(res, vec(sdHexPrism(    sub(pos,float3( 1.0,0.2, -3.0)), float2(0.2,0.05) ), 18.4));
	}

// bounding box
	if (any(lt(abs(add(px,1)), 0.6))) {
		res = opU(res, vec(sdPyramid(    sub(pos,float3(-1.0,-0.6,-3.0)), 1.0 ), 13.56));
		res = opU(res, vec(sdOctahedron( sub(pos,float3(-1.0,0.15,-2.0)), 0.35 ), 23.56));
		res = opU(res, vec(sdTriPrism(   sub(pos,float3(-1.0,0.15,-1.0)), float2(0.3,0.05) ),43.5));
		res = opU(res, vec(sdEllipsoid(  sub(pos,float3(-1.0,0.25, 0.0)), float3(0.2, 0.25, 0.05) ), 43.17));
		res = opU(res, vec(sdHorseshoe(  sub(pos,float3(-1.0,0.25, 1.0)), float2(cos(1.3),sin(1.3)), 0.2, 0.3, float2(0.03,0.08) ), 11.5));
	}

// bounding box
	if (any(lt(abs(sub(px,2)), 0.6))) {
		res = opU(res, vec(sdOctogonPrism(sub(pos,float3( 2.0,0.2, -3.0)), 0.2, 0.05), 51.8));
		res = opU(res, vec(sdCylinder(    sub(pos,float3( 2.0,0.14,-2.0)), float3(0.1,-0.1,0.0), float3(-0.2,0.35,0.1), 0.08), 31.2));
		res = opU(res, vec(sdCappedCone(  sub(pos,float3( 2.0,0.09,-1.0)), float3(0.1,0.0,0.0), float3(-0.2,0.40,0.1), 0.15, 0.05), 46.1));
		res = opU(res, vec(sdRoundCone(   sub(pos,float3( 2.0,0.15, 0.0)), float3(0.1,0.0,0.0), float3(-0.1,0.35,0.1), 0.15, 0.05), 51.7));
		res = opU(res, vec(sdRoundCone(   sub(pos,float3( 2.0,0.20, 1.0)), 0.2, 0.1, 0.3 ), 37.0));
	}
    
    return res;
}

// https://iquilezles.org/articles/boxfunctions
inline vec2
iBox(vec3 ro, vec3 rd, float3 rad) 
{
	vec3 m = rdiv(rd, 1.0);
	vec3 n = neg(mul(m,ro));
	vec3 k = mul(abs(m), rad);
	vec3 t1 = sub(n, k);
	vec3 t2 = add(n, k);
	return vec( max( max( getx(t1), gety(t1) ), getz(t1) ),
	            min( min( getx(t2), gety(t2) ), getz(t2) ) );
}

inline
vec2 raycast(vec3 ro, vec3 rd )
{
	vec2 res = vec(-1.0,-1.0);

vec1 tmin = vec(1.0);
	vec1 tmax = vec(20.0);

// raytrace floor plane
	vec1 tp1 = div(rsub(gety(ro), 0), gety(rd));
	vec1b m1 = gt(tp1, 0);
	tmax = vif(m1, min(tmax, tp1), tmax);
	res  = vif(m1, vec(tp1, 1), res);

// raymarch primitives   
	vec2 tb = iBox(sub(ro,float3(0.0,0.4,-0.5)), rd, float3(2.5,0.41,3.0) );
	if (none(mand(mand( lt(getx(tb),gety(tb)), gt(gety(tb),0)), lt(getx(tb),tmax))))
		return res;

tmin = max(getx(tb), tmin);
	tmax = min(gety(tb), tmax);

vec1 t = tmin;
	for (size_t i = 0; i < 70 && any(lt(t,tmax)); ++i) {
		vec2 h = map(macc(ro, rd, dup3(t)));
		if (any(lt(abs(getx(h)), mul(t, 0.0001)))) {
			res = vec(t,gety(h));
			break;
		}
		t = add(t, getx(h));
	}
	return res;
}

//------------------------------------------------------------------
  
// https://iquilezles.org/articles/rmshadows
inline vec1
calcSoftshadow(vec3 ro, vec3 rd, float mint, float tmax)
{
	// bounding volume
	vec1 tp = div(rsub(gety(ro), 0.8),gety(rd));
	vec1 vtmax = vif(le(tp,0), tmax, min(tp,tmax));
	vec1 res = vec(1);
	vec1 t = vec(mint);
	vec1b m = eq(vec(0),0);
	for (size_t i = 0; i < 24; ++i) {
		vec1 h = getx(map(macc(ro,rd,dup3(t))));
		vec1 s = clamp(mul(h,mul(t,1/8.0)),0.0,1.0);
		res = vif(m,min(res, s),res);
		t = add(t, clamp(h,0.01,0.2));
		m = mnor(lt(res,0.004), gt(t,vtmax));
		if (none(m)) break;
	}
	res = clamp(res, 0.0, 1.0);
	return mul(mul(res,res), madd(res,-2,vec(3)));
}

// https://iquilezles.org/articles/normalsSDF
inline vec3
calcNormal(vec3 pos)
{
	float2 e = float2(1.0,-1.0)*0.5773*0.0005;
	return normalize(macc(macc(macc(mul(
			dup3(getx(map(add(pos, float3(e.x,e.y,e.y))))), float3(e.x,e.y,e.y)),
			float3(e.y,e.y,e.x), dup3(getx(map(add(pos, float3(e.y,e.y,e.x)))))),
			float3(e.y,e.x,e.y), dup3(getx(map(add(pos, float3(e.y,e.x,e.y)))))),
			float3(e.x,e.x,e.x), dup3(getx(map(add(pos, float3(e.x,e.x,e.x)))))));
}

// https://iquilezles.org/articles/nvscene2008/rwwtt.pdf
inline vec1
calcAO(vec3 pos, vec3 nor)
{
	vec1 occ = vec(0);
	vec1 sca = vec(1);
	vec1b m = eq(vec(0),0);
	for (size_t i = 0; i < 5; ++i) {
		float h = 0.01 + 0.12*i/4.0;
		vec1 d = getx(map(macc(pos,h,nor)));
		occ = vif(m, add(occ,mul(rsub(d,h),sca)), occ);
		sca = mul(sca,0.95);
		m = le(occ,0.35);
		if (none(m)) break;
	}
	return mul(clamp(add(mul(occ,-3.0),1),0,1), add(mul(gety(nor), 0.5), 0.5));
}

// https://iquilezles.org/articles/checkerfiltering
inline vec1
checkersGradBox(vec2 p, vec2 dpdx, vec2 dpdy)
{
    // filter kernel
    vec2 w = add(add(abs(dpdx),abs(dpdy)),0.001);
    // analytical integral (box filter)
    vec2 i = mul(sub(abs(sub(fract(mul(nmsac(p,0.5,w),0.5)),0.5)), abs(sub(fract(mul(macc(p,0.5,w),0.5)),0.5))), mul(w,0.5));
    // xor pattern
    return rsub(mul(mul(getx(i),gety(i)),0.5),0.5);
}

inline vec3
render(vec3 ro, vec3 rd, vec3 rdx, vec3 rdy)
{ 
	// background
	vec3 bg = rsub(dup3(mul(max(gety(rd),0.0),0.3)), float3(0.7, 0.7, 0.9));
	vec3 col;
	
	// raycast scene
	vec2 res = raycast(ro,rd);
	vec1 t = getx(res);
	vec1 m = gety(res);
	vec1b mhit = gt(m,-0.5);
	if (any(mhit)) {
		vec3 pos = macc(ro,dup3(t),rd);
		vec3 nor = vif(lt(m,1.5), float3(0.0,1.0,0.0), calcNormal(pos));
		vec3 ref = reflect(rd, nor);

// material	    
		col = add(mul(sin(add(dup3(mul(m,2.0)), float3(0.0,1.0,2.0))),0.2), 0.2);
		vec1 ks = vec(1);
	    
		vec1b m2 = lt(m,1.5);
		{
			// project pixel footprint into the plane
			vec3 dpdx = mul(dup3(gety(ro)), sub(div(rd,dup3(gety(rd))), div(rdx,dup3(gety(rdx)))));
			vec3 dpdy = mul(dup3(gety(ro)), sub(div(rd,dup3(gety(rd))), div(rdy,dup3(gety(rdy)))));

vec1 f = checkersGradBox( mul(getxz(pos),3), mul(getxz(dpdx),3), mul(getxz(dpdy),3) );
			col = vif(m2, add(mul(dup3(f),float3(0.05,0.05,0.05)),0.15), col);
			ks = vif(m2, 0.4, ks);
		}

// lighting
		vec1 occ = calcAO(pos, nor);
		vec3 lin = vec(0,0,0);

// sun
		{
			vec3 lig = normalize(vec(-0.5, 0.4, -0.6));
			vec3 hal = normalize(sub(lig,rd));
			vec1 dif = clamp(dot(nor,lig), 0.0, 1.0 );
			dif = mul(dif,calcSoftshadow(pos, lig, 0.02, 2.5 ));
			vec1 spe = pow16( clamp(dot(nor,hal), 0.0, 1.0));
			spe = mul(spe,dif);
			spe = mul(spe, add(mul(pow5(clamp(rsub(dot(hal,lig),1),0,1)),0.96), 0.04));
			lin = add(lin, mul(mul(mul(col,2.20),dup3(dif)),float3(1.30,1.00,0.70)));
			lin = add(lin, mul(dup3(mul(spe, 5)), mul(dup3(ks),float3(1.30,1.00,0.70))));
		}
		// sky
		{
			vec1 dif = sqrt(clamp(add(mul(gety(nor),0.5),0.5), 0.0, 1.0 ));
			dif = mul(dif,occ);
			vec1 spe = smoothstep(-0.2, 0.2, gety(ref));
			spe = mul(spe,dif);
			spe = mul(spe, add(mul(pow5(clamp(add(dot(nor,rd),1),0,1)),0.96), 0.04));
			//if( spe>0.001 )
			spe = mul(spe, calcSoftshadow(pos, ref, 0.02, 2.5));
			lin = add(lin, mul(mul(col,0.60),mul(dup3(dif),float3(0.40,0.60,1.15))));
			lin = add(lin, mul(mul(dup3(spe),2), mul(dup3(ks),float3(0.40,0.60,1.30))));
		}
		// back
		{
			vec1 dif = mul(clamp(dot(nor,normalize(float3(0.5,0.0,0.6))),0,1),
			               clamp(rsub(gety(pos),1),0,1));
			dif = mul(dif,occ);
			lin = add(lin, mul(mul(col,0.55),mul(dup3(dif),float3(0.25,0.25,0.25))));
		}
		// sss
		{
			vec1 dif = pow2(clamp(add(dot(nor,rd),1),0,1));
			dif = mul(dif,occ);
			lin = add(lin, mul(mul(col,0.25),mul(dup3(dif),float3(1.00,1.00,1.00))));
		}
	    
		col = lin;

col = mix(col, vec(0.7,0.7,0.9), dup3(rsub(exp(mul(pow3(t),-0.0001)),1)));
	}
	col = vif(mhit, col, bg);

return vec3(clamp(col,0.0,1.0));
}

vec3
mainImage(vec2 fragCoord, vec3 ta, vec3 ro, float2 res)
{
#if 1
	return dup3(getx(map(ro)));
#else
	// camera-to-world transformation
	vec3 cw = normalize(sub(ta,ro));
	vec3 cp = vec(sin(0),cos(0),0.0);
	vec3 cu = normalize(cross(cw,cp));
	vec3 cv =          (cross(cu,cw));

vec3 tot = vec(0,0,0);
	vec2 p = mul(sub(mul(fragCoord,2.0),res),1/res.y);

// focal length
 	vec1 fl = vec(2.5);
        
 	// ray direction
 	vec3 rd = mat3v(cp,cu,cv, normalize(vec(p,fl)));

// ray differentials
 	vec2 px = mul(sub(mul(add(fragCoord,vec(1.0,0.0)),2),res),1/res.y);
 	vec2 py = mul(sub(mul(add(fragCoord,vec(0.0,1.0)),2),res),1/res.y);
 	vec3 rdx = mat3v(cp,cu,cv, normalize(vec(px,fl)));
 	vec3 rdy = mat3v(cp,cu,cv, normalize(vec(py,fl)));
        
 	
 	vec3 col = render(ro, rd, rdx, rdy);
	tot = add(tot,col);
	return tot;

#endif
}

Source code

cbuffer ConstantBuffer : register(b0)
{
    float iTime;
    float2 iResolution;
};

struct PSInput
{
    float4 position : SV_POSITION;
    float2 fragCoord : TEXCOORD0;
};

//////////////////////////////////////////////////////////////////////////
// Utility functions

float dot2(float2 v) { return dot(v, v); }
float dot2(float3 v) { return dot(v, v); }
float ndot(float2 a, float2 b) { return a.x * b.x - a.y * b.y; }

float  fract(float  v) { return v - floor(v); }
float2 fract(float2 v) { return v - floor(v); }
float3 fract(float3 v) { return v - floor(v); }

//////////////////////////////////////////////////////////////////////////
// Signed distance functions

float sdPlane(float3 p)
{
    return p.y;
}

float sdSphere(float3 p, float s)
{
    return length(p) - s;
}

float sdBox(float3 p, float3 b)
{
    float3 d = abs(p) - b;
    return min(max(d.x, max(d.y, d.z)), 0.0) + length(max(d, 0.0));
}

float sdBoxFrame(float3 p, float3 b, float e)
{
    p = abs(p) - b;
    float3 q = abs(p + e) - e;
    return min(min(
        length(max(float3(p.x, q.y, q.z), 0.0)) + min(max(p.x, max(q.y, q.z)), 0.0),
        length(max(float3(q.x, p.y, q.z), 0.0)) + min(max(q.x, max(p.y, q.z)), 0.0)),
        length(max(float3(q.x, q.y, p.z), 0.0)) + min(max(q.x, max(q.y, p.z)), 0.0));
}

float sdEllipsoid(float3 p, float3 r)
{
    float k0 = length(p / r);
    float k1 = length(p / (r * r));
    return k0 * (k0 - 1.0) / k1;
}

float sdTorus(float3 p, float2 t)
{
    float2 q = float2(length(p.xz) - t.x, p.y);
    return length(q) - t.y;
}

float sdCappedTorus(float3 p, float2 sc, float ra, float rb)
{
    p.x = abs(p.x);
    float k = (sc.y * p.x > sc.x * p.y) ? dot(p.xy, sc) : length(p.xy);
    return sqrt(dot(p, p) + ra * ra - 2.0 * ra * k) - rb;
}

float sdHexPrism(float3 p, float2 h)
{
    float3 q = abs(p);
    const float3 k = float3(-0.8660254, 0.5, 0.57735);
    p = abs(p);
    p.xy -= 2.0 * min(dot(k.xy, p.xy), 0.0) * k.xy;
    float2 d = float2(
       length(p.xy - float2(clamp(p.x, -k.z * h.x, k.z * h.x), h.x)) * sign(p.y - h.x),
       p.z - h.y);
    return min(max(d.x, d.y), 0.0) + length(max(d, 0.0));
}

float sdOctogonPrism(float3 p, float r, float h)
{
    const float3 k = float3(-0.9238795325, 0.3826834323, 0.4142135623);
    p = abs(p);
    p.xy -= 2.0 * min(dot(float2(k.x, k.y), p.xy), 0.0) * float2(k.x, k.y);
    p.xy -= 2.0 * min(dot(float2(-k.x, k.y), p.xy), 0.0) * float2(-k.x, k.y);
    p.xy -= float2(clamp(p.x, -k.z * r, k.z * r), r);
    float2 d = float2(length(p.xy) * sign(p.y), p.z - h);
    return min(max(d.x, d.y), 0.0) + length(max(d, 0.0));
}

float sdCapsule(float3 p, float3 a, float3 b, float r)
{
    float3 pa = p - a;
    float3 ba = b - a;
    float h = clamp(dot(pa, ba) / dot(ba, ba), 0.0, 1.0);
    return length(pa - ba * h) - r;
}

float sdRoundCone(float3 p, float r1, float r2, float h)
{
    float2 q = float2(length(p.xz), p.y);
    float b = (r1 - r2) / h;
    float a = sqrt(1.0 - b * b);
    float k = dot(q, float2(-b, a));
    if (k < 0.0) return length(q) - r1;
    if (k > a * h) return length(q - float2(0.0, h)) - r2;
    return dot(q, float2(a, b)) - r1;
}

float sdRoundCone(float3 p, float3 a, float3 b, float r1, float r2)
{
    float3 ba = b - a;
    float l2 = dot(ba, ba);
    float rr = r1 - r2;
    float a2 = l2 - rr * rr;
    float il2 = 1.0 / l2;
    
    float3 pa = p - a;
    float y = dot(pa, ba);
    float z = y - l2;
    float x2 = dot2(pa * l2 - ba * y);
    float y2 = y * y * l2;
    float z2 = z * z * l2;

float k = sign(rr) * rr * rr * x2;
    if (sign(z) * a2 * z2 > k) return sqrt(x2 + z2) * il2 - r2;
    if (sign(y) * a2 * y2 < k) return sqrt(x2 + y2) * il2 - r1;
    return (sqrt(x2 * a2 * il2) + y * rr) * il2 - r1;
}

float sdTriPrism(float3 p, float2 h)
{
    const float k = sqrt(3.0);
    h.x *= 0.5 * k;
    p.xy /= h.x;
    p.x = abs(p.x) - 1.0;
    p.y = p.y + 1.0 / k;
    if (p.x + k * p.y > 0.0)
        p.xy = float2(p.x - k * p.y, -k * p.x - p.y) / 2.0;
    p.x -= clamp(p.x, -2.0, 0.0);
    float d1 = length(p.xy) * sign(-p.y) * h.x;
    float d2 = abs(p.z) - h.y;
    return length(max(float2(d1, d2), 0.0)) + min(max(d1, d2), 0.0);
}

float sdCylinder(float3 p, float2 h)
{
    float2 d = abs(float2(length(p.xz), p.y)) - h;
    return min(max(d.x, d.y), 0.0) + length(max(d, 0.0));
}

float sdCylinder(float3 p, float3 a, float3 b, float r)
{
    float3 pa = p - a;
    float3 ba = b - a;
    float baba = dot(ba, ba);
    float paba = dot(pa, ba);

float x = length(pa * baba - ba * paba) - r * baba;
    float y = abs(paba - baba * 0.5) - baba * 0.5;
    float x2 = x * x;
    float y2 = y * y * baba;
    float d = (max(x, y) < 0.0) ? -min(x2, y2) : (((x > 0.0) ? x2 : 0.0) + ((y > 0.0) ? y2 : 0.0));
    return sign(d) * sqrt(abs(d)) / baba;
}

float sdCone(float3 p, float2 c, float h)
{
    float2 q = h * float2(c.x, -c.y) / c.y;
    float2 w = float2(length(p.xz), p.y);
    float2 a = w - q * clamp(dot(w, q) / dot(q, q), 0.0, 1.0);
    float2 b = w - q * float2(clamp(w.x / q.x, 0.0, 1.0), 1.0);
    float k = sign(q.y);
    float d = min(dot(a, a), dot(b, b));
    float s = max(k * (w.x * q.y - w.y * q.x), k * (w.y - q.y));
    return sqrt(d) * sign(s);
}

float sdCappedCone(float3 p, float h, float r1, float r2)
{
    float2 q = float2(length(p.xz), p.y);
    float2 k1 = float2(r2, h);
    float2 k2 = float2(r2 - r1, 2.0 * h);
    float2 ca = float2(q.x - min(q.x, (q.y < 0.0) ? r1 : r2), abs(q.y) - h);
    float2 cb = q - k1 + k2 * clamp(dot(k1 - q, k2) / dot2(k2), 0.0, 1.0);
    float s = ((cb.x < 0.0 && ca.y < 0.0) ? -1.0 : 1.0);
    return s * sqrt(min(dot2(ca), dot2(cb)));
}

float sdCappedCone(float3 p, float3 a, float3 b, float ra, float rb)
{
    float rba = rb - ra;
    float baba = dot(b - a, b - a);
    float papa = dot(p - a, p - a);
    float paba = dot(p - a, b - a) / baba;

float x = sqrt(papa - paba * paba * baba);
    float cax = max(0.0, x - ((paba < 0.5) ? ra : rb));
    float cay = abs(paba - 0.5) - 0.5;
    float k = rba * rba + baba;
    float f = clamp((rba * (x - ra) + paba * baba) / k, 0.0, 1.0);
    float cbx = x - ra - f * rba;
    float cby = paba - f;
    float s = (cbx < 0.0 && cay < 0.0) ? -1.0 : 1.0;
    return s * sqrt(min(cax * cax + cay * cay * baba, cbx * cbx + cby * cby * baba));
}

float sdSolidAngle(float3 pos, float2 c, float ra)
{
    float2 p = float2(length(pos.xz), pos.y);
    float l = length(p) - ra;
    float m = length(p - c * clamp(dot(p, c), 0.0, ra));
    return max(l, m * sign(c.y * p.x - c.x * p.y));
}

float sdOctahedron(float3 p, float s)
{
    p = abs(p);
    float m = p.x + p.y + p.z - s;
    float3 q;
    if (3.0 * p.x < m) q = p;
    else if (3.0 * p.y < m) q = p.yzx;
    else if (3.0 * p.z < m) q = p.zxy;
    else return m * 0.57735027;
    float k = clamp(0.5 * (q.z - q.y + s), 0.0, s);
    return length(float3(q.x, q.y - s + k, q.z - k));
}

float sdPyramid(float3 p, float h)
{
    float m2 = h * h + 0.25;
    p.xz = abs(p.xz);
    if (p.z > p.x) { float temp = p.x; p.x = p.z; p.z = temp; }
    p.xz -= 0.5;
    float3 q = float3(p.z, h * p.y - 0.5 * p.x, h * p.x + 0.5 * p.y);
    float s = max(-q.x, 0.0);
    float t = clamp((q.y - 0.5 * p.z) / (m2 + 0.25), 0.0, 1.0);
    float a = m2 * (q.x + s) * (q.x + s) + q.y * q.y;
    float b = m2 * (q.x + 0.5 * t) * (q.x + 0.5 * t) + (q.y - m2 * t) * (q.y - m2 * t);
    float d2 = (min(q.y, -q.x * m2 - q.y * 0.5) > 0.0) ? 0.0 : min(a, b);
    return sqrt((d2 + q.z * q.z) / m2) * sign(max(q.z, -p.y));
}

float sdRhombus(float3 p, float la, float lb, float h, float ra)
{
    p = abs(p);
    float2 b = float2(la, lb);
    float f = clamp((ndot(b, b - 2.0 * p.xz)) / dot(b, b), -1.0, 1.0);
    float2 q = float2(length(p.xz - 0.5 * b * float2(1.0 - f, 1.0 + f)) * sign(p.x * b.y + p.z * b.x - b.x * b.y) - ra, p.y - h);
    return min(max(q.x, q.y), 0.0) + length(max(q, 0.0));
}

float sdHorseshoe(float3 p, float2 c, float r, float le, float2 w)
{
    p.x = abs(p.x);
    float l = length(p.xy);
    p.xy = float2(-c.x, c.y) * p.xy; // equivalent to rotation by c (using a 2x2 matrix)
    p.xy = (p.y > 0.0 || p.x > 0.0) ? float2(p.x, p.y) : float2(p.x, l * sign(-c.x));
    p.xy = float2(p.x, abs(p.y - r)) - float2(le, 0.0);
    float2 q = float2(length(max(p.xy, 0.0)) + min(0.0, max(p.x, p.y)), p.z);
    float2 d = abs(q) - w;
    return min(max(d.x, d.y), 0.0) + length(max(d, 0.0));
}

//////////////////////////////////////////////////////////////////////////
// Boolean op

float2 opU(float2 d1, float2 d2)
{
    return (d1.x < d2.x) ? d1 : d2;
}

//////////////////////////////////////////////////////////////////////////
// Scene mapping

float2 map(float3 pos)
{
    float2 res = float2(pos.y, 0.0);

if (pos.y < 0.0 || pos.y > 1.2 || abs(pos.x) > 5.0 || abs(pos.z - 1.0) > 5.0)
        return res;

// bounding box
    if (abs(pos.x + 2.0) < 0.6 && abs(pos.z - 0.5) < 1.2)
    {
        res = opU(res, float2(sdSphere(pos - float3(-2.0, 0.25, 0.0), 0.25), 26.9));
        res = opU(res, float2(sdRhombus((pos - float3(-2.0, 0.25, 1.0)).xzy, 0.15, 0.25, 0.04, 0.08), 17.0));
    }

// bounding box
    if (abs(pos.x) < 0.6)
    {
        res = opU(res, float2(sdCappedTorus((pos - float3(0.0, 0.30, 1.0)) * float3(1, -1, 1), float2(0.866025, -0.5), 0.25, 0.05), 25.0));
        res = opU(res, float2(sdBoxFrame(pos - float3(0.0, 0.25, 0.0), float3(0.3, 0.25, 0.2), 0.025), 16.9));
        res = opU(res, float2(sdCone(pos - float3(0.0, 0.45, -1.0), float2(0.6, 0.8), 0.45), 55.0));
        res = opU(res, float2(sdCappedCone(pos - float3(0.0, 0.25, -2.0), 0.25, 0.25, 0.1), 13.67));
        res = opU(res, float2(sdSolidAngle(pos - float3(0.0, 0.00, -3.0), float2(3, 4) / 5.0, 0.4), 49.13));
    }

// bounding box
    if (abs(pos.x - 1.0) < 0.6)
    {
        res = opU(res, float2(sdTorus((pos - float3(1.0, 0.30, 1.0)).xzy, float2(0.25, 0.05)), 7.1));
        res = opU(res, float2(sdBox(pos - float3(1.0, 0.25, 0.0), float3(0.3, 0.25, 0.1)), 3.0));
        res = opU(res, float2(sdCapsule(pos - float3(1.0, 0.00, -1.0), float3(-0.1, 0.1, -0.1), float3(0.2, 0.4, 0.2), 0.1), 31.9));
        res = opU(res, float2(sdCylinder(pos - float3(1.0, 0.25, -2.0), float2(0.15, 0.25)), 8.0));
        res = opU(res, float2(sdHexPrism(pos - float3(1.0, 0.2, -3.0), float2(0.2, 0.05)), 18.4));
    }

// bounding box
    if (abs(pos.x + 1.0) < 0.6)
    {
        res = opU(res, float2(sdPyramid(pos - float3(-1.0, -0.6, -3.0), 1.0), 13.56));
        res = opU(res, float2(sdOctahedron(pos - float3(-1.0, 0.15, -2.0), 0.35), 23.56));
        res = opU(res, float2(sdTriPrism(pos - float3(-1.0, 0.15, -1.0), float2(0.3, 0.05)), 43.5));
        res = opU(res, float2(sdEllipsoid(pos - float3(-1.0, 0.25, 0.0), float3(0.2, 0.25, 0.05)), 43.17));
        res = opU(res, float2(sdHorseshoe(pos - float3(-1.0, 0.25, 1.0), float2(cos(1.3), sin(1.3)), 0.2, 0.3, float2(0.03, 0.08)), 11.5));
    }

// bounding box
    if (abs(pos.x - 2.0) < 0.6)
    {
        res = opU(res, float2(sdOctogonPrism(pos - float3(2.0, 0.2, -3.0), 0.2, 0.05), 51.8));
        res = opU(res, float2(sdCylinder(pos - float3(2.0, 0.14, -2.0), float3(0.1, -0.1, 0.0), float3(-0.2, 0.35, 0.1), 0.08), 31.2));
        res = opU(res, float2(sdCappedCone(pos - float3(2.0, 0.09, -1.0), float3(0.1, 0.0, 0.0), float3(-0.2, 0.40, 0.1), 0.15, 0.05), 46.1));
        res = opU(res, float2(sdRoundCone(pos - float3(2.0, 0.15, 0.0), float3(0.1, 0.0, 0.0), float3(-0.1, 0.35, 0.1), 0.15, 0.05), 51.7));
        res = opU(res, float2(sdRoundCone(pos - float3(2.0, 0.20, 1.0), 0.2, 0.1, 0.3), 37.0));
    }
    
    return res;
}

//////////////////////////////////////////////////////////////////////////
// Intersection helpers

float2 iBox(float3 ro, float3 rd, float3 rad)
{
    float3 m = 1.0 / rd;
    float3 n = m * ro;
    float3 k = abs(m) * rad;
    float3 t1 = -n - k;
    float3 t2 = -n + k;
    return float2(max(max(t1.x, t1.y), t1.z),
                  min(min(t2.x, t2.y), t2.z));
}

float2 raycast(float3 ro, float3 rd)
{
    float2 res = float2(-1.0, -1.0);
    float tmin = 1.0;
    float tmax = 20.0;
    float tp1 = (0.0 - ro.y) / rd.y;
    if (tp1 > 0.0)
    {
        tmax = min(tmax, tp1);
        res = float2(tp1, 1.0);
    }
    
    float2 tb = iBox(ro - float3(0.0, 0.4, -0.5), rd, float3(2.5, 0.41, 3.0));
    if (tb.x < tb.y && tb.y > 0.0 && tb.x < tmax)
    {
        tmin = max(tb.x, tmin);
        tmax = min(tb.y, tmax);
        float t = tmin;
        [unroll]
        for (int i = 0; i < 70 && t < tmax; i++)
        {
            float2 h = map(ro + rd * t);
            if (abs(h.x) < (0.0001 * t))
            {
                res = float2(t, h.y);
                break;
            }
            t += h.x;
        }
    }
    
    return res;
}

float calcSoftshadow(float3 ro, float3 rd, float mint, float tmax)
{
    float tp = (0.8 - ro.y) / rd.y; 
    if (tp > 0.0) tmax = min(tmax, tp);
    float res = 1.0;
    float t = mint;
    [unroll]
    for (int i = 0; i < 24; i++)
    {
        float h = map(ro + rd * t).x;
        float s = clamp(8.0 * h / t, 0.0, 1.0);
        res = min(res, s);
        t += clamp(h, 0.01, 0.2);
        if (res < 0.004 || t > tmax) break;
    }
    res = clamp(res, 0.0, 1.0);
    return res * res * (3.0 - 2.0 * res);
}

float3 calcNormal(float3 pos)
{
    float2 e = float2(1.0, -1.0) * 0.5773 * 0.0005;
    return normalize(e.xyy * map(pos + float3(e.x, e.x, e.x)).x + 
                     e.yyx * map(pos + float3(e.y, e.y, e.x)).x + 
                     e.yxy * map(pos + float3(e.y, e.x, e.y)).x + 
                     e.xxx * map(pos + float3(e.x, e.x, e.x)).x);
}

float calcAO(float3 pos, float3 nor)
{
    float occ = 0.0;
    float sca = 1.0;
    [unroll]
    for (int i = 0; i < 5; i++)
    {
        float h = 0.01 + 0.12 * float(i) / 4.0;
        float d = map(pos + nor * h).x;
        occ += (h - d) * sca;
        sca *= 0.95;
        if (occ > 0.35) break;
    }
    return clamp(1.0 - 3.0 * occ, 0.0, 1.0) * (0.5 + 0.5 * nor.y);
}

float checkersGradBox(float2 p, float2 dpdx, float2 dpdy)
{
    float2 w = abs(dpdx) + abs(dpdy) + 0.001;
    float2 i = 2.0 * (abs(fract((p - 0.5 * w) * 0.5) - 0.5) - abs(fract((p + 0.5 * w) * 0.5) - 0.5)) / w;
    return 0.5 - 0.5 * i.x * i.y;
}

//////////////////////////////////////////////////////////////////////////
// Rendering

float3 render(float3 ro, float3 rd, float3 rdx, float3 rdy)
{
    float3 col = float3(0.7, 0.7, 0.9) - max(rd.y, 0.0) * 0.3;
    float2 res = raycast(ro, rd);
    float t = res.x;
    float m = res.y;
    if (m > -0.5)
    {
        float3 pos = ro + rd * t;
        float3 nor = (m < 1.5) ? float3(0.0, 1.0, 0.0) : calcNormal(pos);
        float3 ref = reflect(rd, nor);
        col = 0.2 + 0.2 * sin(m * 2.0 + float3(0.0, 1.0, 2.0));
        float ks = 1.0;
        if (m < 1.5)
        {
            float3 dpdx = ro.y * (rd / rd.y - rdx / rdx.y);
            float3 dpdy = ro.y * (rd / rd.y - rdy / rdy.y);
            float f = checkersGradBox(3.0 * pos.xz, 3.0 * dpdx.xz, 3.0 * dpdy.xz);
            col = 0.15 + f * float3(0.05, 0.05, 0.05);
            ks = 0.4;
        }
        float occ = calcAO(pos, nor);
        float3 lin = float3(0.0, 0.0, 0.0);
        {
            float3 lig = normalize(float3(-0.5, 0.4, -0.6));
            float3 hal = normalize(lig - rd);
            float dif = clamp(dot(nor, lig), 0.0, 1.0);
            dif *= calcSoftshadow(pos, lig, 0.02, 2.5);
            float spe = pow(clamp(dot(nor, hal), 0.0, 1.0), 16.0);
            spe *= dif;
            spe *= 0.04 + 0.96 * pow(clamp(1.0 - dot(hal, lig), 0.0, 1.0), 5.0);
            lin += col * 2.20 * dif * float3(1.30, 1.00, 0.70);
            lin += 5.00 * spe * float3(1.30, 1.00, 0.70) * ks;
        }
        {
            float dif = sqrt(clamp(0.5 + 0.5 * nor.y, 0.0, 1.0));
            dif *= occ;
            float spe = smoothstep(-0.2, 0.2, ref.y);
            spe *= dif;
            spe *= 0.04 + 0.96 * pow(clamp(1.0 + dot(nor, rd), 0.0, 1.0), 5.0);
            spe *= calcSoftshadow(pos, ref, 0.02, 2.5);
            lin += col * 0.60 * dif * float3(0.40, 0.60, 1.15);
            lin += 2.00 * spe * float3(0.40, 0.60, 1.30) * ks;
        }
        {
            float dif = clamp(dot(nor, normalize(float3(0.5, 0.0, 0.6))), 0.0, 1.0) * clamp(1.0 - pos.y, 0.0, 1.0);
            dif *= occ;
            lin += col * 0.55 * dif * float3(0.25, 0.25, 0.25);
        }
        {
            float dif = pow(clamp(1.0 + dot(nor, rd), 0.0, 1.0), 2.0);
            dif *= occ;
            lin += col * 0.25 * dif * float3(1.00, 1.00, 1.00);
        }
        col = lin;
        col = lerp(col, float3(0.7, 0.7, 0.9), 1.0 - exp(-0.0001 * t * t * t));
    }
    return saturate(col);
}

//////////////////////////////////////////////////////////////////////////
// Camera setup

float3x3 setCamera(float3 ro, float3 ta, float cr)
{
    float3 cw = normalize(ta - ro);
    float3 cp = float3(sin(cr), cos(cr), 0.0);
    float3 cu = normalize(cross(cw, cp));
    float3 cv = cross(cu, cw);
    return float3x3(cu, cv, cw);
}

//////////////////////////////////////////////////////////////////////////
// Main Image

float4 mainImage(float2 fragCoord)
{
#if 1
return map(fragCoord.xyy).xxxx;
#else
    float time = 32.0 + iTime * 1.5;
    float3 ta = float3(0.25, -0.75, -0.75);
    float3 ro = ta + float3(4.5 * cos(0.1 * time), 2.2, 4.5 * sin(0.1 * time));
    float3x3 ca = setCamera(ro, ta, 0.0);
    float2 p = (2.0 * fragCoord - iResolution) / iResolution.y;
    const float fl = 2.5;
    float3 rd = mul(ca, normalize(float3(p, fl)));
    float2 px = (2.0 * (fragCoord + float2(1.0, 0.0)) - iResolution) / iResolution.y;
    float2 py = (2.0 * (fragCoord + float2(0.0, 1.0)) - iResolution) / iResolution.y;
    float3 rdx = mul(ca, normalize(float3(px, fl)));
    float3 rdy = mul(ca, normalize(float3(py, fl)));
    float3 col = render(ro, rd, rdx, rdy);
    return float4(col, 1.0);
#endif
}

//////////////////////////////////////////////////////////////////////////
// Pixel Shader Entry Point

float4 PSMain(PSInput input) : SV_Target
{
    return mainImage(input.fragCoord);
}

Source code

#include <math.h>
#include <immintrin.h>

typedef __m512i vec1;
struct vec2 { __m512 x,y; vec2(vec1 x, vec1 y):x(x),y(y){} };
struct vec3 { __m512 x,y,z; vec3(vec1 x, vec1 y, vec1 z):x(x),y(y),z(z){} };
struct vec4 { __m512 x,y,z,w; vec4(vec1 x, vec1 y, vec1 z, vec1 w):x(x),y(y),z(z),w(w){} };
typedef __mmask16 vec1b;

#define getx(v) v.x
#define gety(v) v.y
#define getz(v) v.z
#define getw(v) v.w

#define getxy(v) vec(getx(v), gety(v))
#define getxz(v) vec(getx(v), getz(v))
#define getyz(v) vec(gety(v), getz(v))
#define getyx(v) vec(gety(v), getx(v))

#define getxzy(v) vec(getx(v), getz(v), gety(v))
#define getyzx(v) vec(gety(v), getz(v), getx(v))
#define getzxy(v) vec(getz(v), getx(v), gety(v))

inline vec2 vec(vec1 x, vec1 y) { return vec2(x, y); }
inline vec3 vec(vec1 x, vec1 y, vec1 z) { return vec3(x, y, z); }
inline vec3 vec(vec2 xy, vec1 z) { return vec3(getx(xy), gety(xy), z); }
inline vec3 vec(vec1 x, vec2 yz) { return vec3(x, getx(yz), gety(yz)); }
inline vec4 vec(vec1 x, vec1 y, vec1 z, vec1 w) { return vec4(x, y, z, w); }

inline vec2 dup2(vec1 x) { return vec(x, x); }
inline vec3 dup3(vec1 x) { return vec(x, x, x); }
inline vec4 dup4(vec1 x) { return vec(x, x, x, x); }

inline vec1 vec(__m512 x) { return x; }
inline vec1 vec(float x) { return _mm512_set1_ps(x); }
inline vec2 vec(float x, float y) { return vec(vec(x), vec(y)); }
inline vec3 vec(float x, float y, float z) { return vec(vec(x), vec(y), vec(z)); }
inline vec4 vec(float x, float y, float z, float w) { return vec(vec(x), vec(y), vec(z), vec(w)); }

inline vec2 vec(vec1 x, float y) { return vec(x, vec(y)); }

OP2(+) OP2(-) OP2(/) OP2(*)

#define ALIAS1(n, f) \
	inline auto n(vec1 v) { return _mm512_##f##_ps(v); } \
	inline auto n(vec2 v) { return vec(n(getx(v)), n(gety(v))); } \
	inline auto n(vec3 v) { return vec(n(getx(v)), n(gety(v)), n(getz(v))); } \
	inline auto n(vec4 v) { return vec(n(getx(v)), n(gety(v)), n(getz(v)), n(getw(v))); }

#define ALIAS2(n, f) \
	inline auto n(vec1 a, vec1 b) { return _mm512_##f##_ps(a, b); } \
	inline auto n(vec2 a, vec2 b) { return vec(n(getx(a),getx(b)), n(gety(a),gety(b))); } \
	inline auto n(vec3 a, vec3 b) { return vec(n(getx(a),getx(b)), n(gety(a),gety(b)), n(getz(a),getz(b))); } \
	inline auto n(vec4 a, vec4 b) { return vec(n(getx(a),getx(b)), n(gety(a),gety(b)), n(getz(a),getz(b)), n(getw(a),getw(b))); }

#define ALIAS2x(n, f) ALIAS2(n, f) \
	inline auto n(vec1 a, float b) { return _mm512_##f##_ps(a, vec(b)); } \
	inline auto n(vec2 a, float b) { return vec(n(getx(a),b), n(gety(a),b)); } \
	inline auto n(vec3 a, float b) { return vec(n(getx(a),b), n(gety(a),b), n(getz(a),b)); } \
	inline auto n(vec4 a, float b) { return vec(n(getx(a),b), n(gety(a),b), n(getz(a),b), n(getw(a),b)); } \
	inline auto n(vec2 a, float2 b) { return vec(n(getx(a),b.x), n(gety(a),b.y)); } \
	inline auto n(vec3 a, float3 b) { return vec(n(getx(a),b.x), n(gety(a),b.y), n(getz(a),b.z)); } \
	inline auto n(vec4 a, float4 b) { return vec(n(getx(a),b.x), n(gety(a),b.y), n(getz(a),b.z), n(getw(a),b.w)); }

#define ALIAS3(n, f) \
	inline auto n(vec1 a, vec1 b, vec1 c) { return _mm512_##f##_ps(a, b, c); } \
	inline auto n(vec2 a, vec2 b, vec2 c) { return vec(n(getx(a),getx(b),getx(c)), n(gety(a),gety(b),gety(c))); } \
	inline auto n(vec3 a, vec3 b, vec3 c) { return vec(n(getx(a),getx(b),getx(c)), n(gety(a),gety(b),gety(c)), n(getz(a),getz(b),getz(c))); } \
	inline auto n(vec4 a, vec4 b, vec4 c) { return vec(n(getx(a),getx(b),getx(c)), n(gety(a),gety(b),gety(c)), n(getz(a),getz(b),getz(c)), n(getw(a),getw(b),getw(c))); }

#define ALIAS3x(n, f) ALIAS3(n, f) \
	inline auto n(vec1 a, float b, vec1 c) { return _mm512_##f##_ps(a, vec(b), c); } \
	inline auto n(vec2 a, float b, vec2 c) { return vec(n(getx(a),b,getx(c)), n(gety(a),b,gety(c))); } \
	inline auto n(vec3 a, float b, vec3 c) { return vec(n(getx(a),b,getx(c)), n(gety(a),b,gety(c)), n(getz(a),b,getz(c))); } \
	inline auto n(vec4 a, float b, vec4 c) { return vec(n(getx(a),b,getx(c)), n(gety(a),b,gety(c)), n(getz(a),b,getz(c)), n(getw(a),b,getw(c))); } \
	inline auto n(vec2 a, float2 b, vec2 c) { return vec(n(getx(a),b.x,getx(c)), n(gety(a),b.y,gety(c))); } \
	inline auto n(vec3 a, float3 b, vec3 c) { return vec(n(getx(a),b.x,getx(c)), n(gety(a),b.y,gety(c)), n(getz(a),b.z,getz(c))); } \
	inline auto n(vec4 a, float4 b, vec4 c) { return vec(n(getx(a),b.x,getx(c)), n(gety(a),b.y,gety(c)), n(getz(a),b.z,getz(c)), n(getw(a),b.w,getw(c))); }

inline __m512 _mm512_neg_ps(__m512 a) { return -a;} 
inline __m512 _mm512_sgnj_ps(__m512 a, __m512 b) { return b < 0 ? -a : a; } 
inline __m512 _mm512_sgnjn_ps(__m512 a, __m512 b) { return b < 0 ? a : -a; } 
inline __m512 _mm512_fmacc_ps(__m512 a, __m512 b, __m512 c) { return _mm512_fmadd_ps(c, b, a); }

inline __m512 _mm512_fnmacc_ps(__m512 a, __m512 b, __m512 c) { return _mm512_fnmsub_ps(c, b, a); } 
inline __m512 _mm512_fnmsac_ps(__m512 a, __m512 b, __m512 c) { return _mm512_fnmadd_ps(c, b, a); }

ALIAS1(abs, abs)
ALIAS1(neg, neg)
ALIAS1(sqrt, sqrt)
ALIAS2x(div, div)
ALIAS2x(add, add)
ALIAS2x(sub, sub)
ALIAS2x(mul, mul)
ALIAS2x(max, max)
ALIAS2x(min, min)
ALIAS2x(sgnj, sgnj)
ALIAS2x(sgnjn, sgnjn)

ALIAS3x(macc, fmacc)
ALIAS3x(nmacc, fnmacc)
ALIAS3x(madd, fmadd)
ALIAS3x(msub, fmsub)
ALIAS3x(nmsac, fnmsac)

inline auto rsub(vec1 a, float b) { return sub(vec(b),a); }
inline auto rsub(vec2 a, float b) { return sub(vec(b,b),a); }
inline auto rsub(vec3 a, float b) { return sub(vec(b,b,b),a); }
inline auto rsub(vec4 a, float b) { return sub(vec(b,b,b,b),a); }
inline auto rsub(vec2 a, float2 b) { return vec(rsub(a.x,b.x),rsub(a.y,b.y)); }
inline auto rsub(vec3 a, float3 b) { return vec(rsub(a.x,b.x),rsub(a.y,b.y),rsub(a.z,b.z)); }
inline auto rsub(vec4 a, float4 b) { return vec(rsub(a.x,b.x),rsub(a.y,b.y),rsub(a.z,b.z),rsub(a.w,b.w)); }

inline auto rdiv(vec1 a, float b) { return div(vec(b),a); }
inline auto rdiv(vec2 a, float b) { return div(vec(b,b),a); }
inline auto rdiv(vec3 a, float b) { return div(vec(b,b,b),a); }
inline auto rdiv(vec4 a, float b) { return div(vec(b,b,b,b),a); }
inline auto rdiv(vec2 a, float2 b) { return vec(rdiv(a.x,b.x),rdiv(a.y,b.y)); }
inline auto rdiv(vec3 a, float3 b) { return vec(rdiv(a.x,b.x),rdiv(a.y,b.y),rdiv(a.z,b.z)); }
inline auto rdiv(vec4 a, float4 b) { return vec(rdiv(a.x,b.x),rdiv(a.y,b.y),rdiv(a.z,b.z),rdiv(a.w,b.w)); }

#define CMP(n,f) \
	inline auto f(vec1 a, vec1 b) { return _mm512_cmp##f##_ps_mask(a, b); } \
	inline auto f(vec1 a, float b) { return _mm512_cmp##f##_ps_mask(a, vec(b)); }

inline __mmask16 _mm512_cmpgt_ps_mask(__m512 a, __m512 b) { return _mm512_cmplt_ps_mask(b,a); } 
inline __mmask16 _mm512_cmpge_ps_mask(__m512 a, __m512 b) { return _mm512_cmple_ps_mask(b,a); } 
	
CMP(eq,eq) CMP(ne,neq) CMP(lt,lt) CMP(le,le) CMP(gt,gt) CMP(ge,ge)

inline vec1 vif(vec1b m, vec1 t, vec1 f) { return _mm512_mask_blend_ps(m, f, t); }
inline vec2 vif(vec1b m, vec2 t, vec2 f) { return vec(vif(m, getx(t), getx(f)), vif(m, gety(t), gety(f))); }
inline vec3 vif(vec1b m, vec3 t, vec3 f) { return vec(vif(m, getx(t), getx(f)), vif(m, gety(t), gety(f)), vif(m, getz(t), getz(f))); }
inline vec4 vif(vec1b m, vec4 t, vec4 f) { return vec(vif(m, getx(t), getx(f)), vif(m, gety(t), gety(f)), vif(m, getz(t), getz(f)), vif(m, getw(t), getw(f))); }

inline vec1 vif(vec1b m, float t, vec1 f) { return _mm512_mask_blend_ps(m, f, vec(t)); }
inline vec2 vif(vec1b m, float2 t, vec2 f) { return vec(vif(m, t.x, getx(f)), vif(m, t.x, gety(f))); }
inline vec3 vif(vec1b m, float3 t, vec3 f) { return vec(vif(m, t.x, getx(f)), vif(m, t.y, gety(f)), vif(m, t.z, getz(f))); }
inline vec4 vif(vec1b m, float4 t, vec4 f) { return vec(vif(m, t.x, getx(f)), vif(m, t.y, gety(f)), vif(m, t.z, getz(f)), vif(m, t.w, getw(f))); }

inline vec1b mor(vec1b a, vec1b b) { return a|b; }
inline vec1b mand(vec1b a, vec1b b) { return a&b; }
inline vec1b mnor(vec1b a, vec1b b) { return ~(a|b); }
inline bool all(vec1b m) { return m == 16; }
inline bool any(vec1b m) { return m; }
inline bool none(vec1b m) { return !m; }