const float tau = radians(360.0);
const float pi  = radians(180.0);
const float hpi = radians(90.0);
const float phi = sqrt(5.0) * 0.5 + 0.5;

const float goldenAngle = tau / (phi + 1.0);
const mat2 rotateGoldenAngle = mat2(cos(goldenAngle), -sin(goldenAngle), sin(goldenAngle), cos(goldenAngle));

//----------------------------------------------------------------------------//

#define Clamp01(x) clamp(x, 0, 1)
#define Max0(x) max(x, 0)

float MinOf(vec2 x) { return min(x.x, x.y); }
float MinOf(vec3 x) { return min(min(x.x, x.y), x.z); }
float MaxOf(vec2 x) { return max(x.x, x.y); }
float MaxOf(vec3 x) { return max(max(x.x, x.y), x.z); }

vec2 Diagonal(mat2 m) { return vec2(m[0].x, m[1].y); }
vec3 Diagonal(mat3 m) { return vec3(m[0].x, m[1].y, m[2].z); }
vec4 Diagonal(mat4 m) { return vec4(m[0].x, m[1].y, m[2].z, m[3].w); }

// Solid angle and cone angle conversion
#define ConeAngleToSolidAngle(x) (tau * (1.0 - cos(x)))
#define SolidAngleToConeAngle(x) acos(1.0 - (x) / tau)

//----------------------------------------------------------------------------//

// Similar to smoothstep, but using linear interpolation instead of Hermite interpolation.
float LinearStep(float e0, float e1, float x) { return Clamp01((x - e0) / (e1 - e0)); }
vec2  LinearStep(vec2  e0, vec2  e1, vec2  x) { return Clamp01((x - e0) / (e1 - e0)); }
vec3  LinearStep(vec3  e0, vec3  e1, vec3  x) { return Clamp01((x - e0) / (e1 - e0)); }
vec4  LinearStep(vec4  e0, vec4  e1, vec4  x) { return Clamp01((x - e0) / (e1 - e0)); }
vec2  LinearStep(vec2  e0, vec2  e1, float x) { return Clamp01((x - e0) / (e1 - e0)); }
vec3  LinearStep(vec3  e0, vec3  e1, float x) { return Clamp01((x - e0) / (e1 - e0)); }
vec4  LinearStep(vec4  e0, vec4  e1, float x) { return Clamp01((x - e0) / (e1 - e0)); }
vec2  LinearStep(float e0, float e1, vec2  x) { return Clamp01((x - e0) / (e1 - e0)); }
vec3  LinearStep(float e0, float e1, vec3  x) { return Clamp01((x - e0) / (e1 - e0)); }
vec4  LinearStep(float e0, float e1, vec4  x) { return Clamp01((x - e0) / (e1 - e0)); }

//----------------------------------------------------------------------------//

vec3 SphericalToCartesian(vec3 spherical) {
	return vec3(vec2(cos(spherical.y), sin(spherical.y)) * sin(spherical.x), cos(spherical.x)) * spherical.z;
}

//----------------------------------------------------------------------------//

vec3 ViewSpaceFromScreenSpace(vec3 pos, mat4 inverseProjection) {
	pos = pos * 2.0 - 1.0;

	vec3 viewPosition  = vec3(vec2(inverseProjection[0].x, inverseProjection[1].y) * pos.xy + inverseProjection[3].xy, inverseProjection[3].z);
	     viewPosition /= inverseProjection[2].w * pos.z + inverseProjection[3].w;

	return viewPosition;
}
// Return value is NOT normalized!
vec3 CalculateViewDirection(vec2 uv, mat4 inverseProjection) {
	uv = uv * 2.0 - 1.0;

	uv = vec2(inverseProjection[0].x, inverseProjection[1].y) * uv + inverseProjection[3].xy;
	return vec3(uv, inverseProjection[3].z);
}

vec3 ScreenSpaceFromViewSpace(vec3 viewPosition, mat4 projection) {
	vec3 screenPosition = vec3(projection[0].x, projection[1].y, projection[2].z) * viewPosition + projection[3].xyz;
	return screenPosition * (0.5 / -viewPosition.z) + 0.5;
}

//----------------------------------------------------------------------------//

float Bayer2  (vec2 c) { c = 0.5 * floor(c); return fract(1.5 * fract(c.y) + c.x); }
float Bayer4  (vec2 c) { return 0.25 * Bayer2  (0.5 * c) + Bayer2(c); }
float Bayer8  (vec2 c) { return 0.25 * Bayer4  (0.5 * c) + Bayer2(c); }
float Bayer16 (vec2 c) { return 0.25 * Bayer8  (0.5 * c) + Bayer2(c); }
float Bayer32 (vec2 c) { return 0.25 * Bayer16 (0.5 * c) + Bayer2(c); }
float Bayer64 (vec2 c) { return 0.25 * Bayer32 (0.5 * c) + Bayer2(c); }
float Bayer128(vec2 c) { return 0.25 * Bayer64 (0.5 * c) + Bayer2(c); }
float Bayer256(vec2 c) { return 0.25 * Bayer128(0.5 * c) + Bayer2(c); }
float Bayer512(vec2 c) { return 0.25 * Bayer256(0.5 * c) + Bayer2(c); }

//----------------------------------------------------------------------------//

float AddUvMargin(float uv, int resolution) { return uv * (1.0 - 1.0 / resolution) + (0.5 / resolution); }
vec2 AddUvMargin(vec2 uv, ivec2 resolution) { return uv * (1.0 - 1.0 / resolution) + (0.5 / resolution); }
vec3 AddUvMargin(vec3 uv, ivec3 resolution) { return uv * (1.0 - 1.0 / resolution) + (0.5 / resolution); }
vec4 AddUvMargin(vec4 uv, ivec4 resolution) { return uv * (1.0 - 1.0 / resolution) + (0.5 / resolution); }
float RemoveUvMargin(float uv, int resolution) { return (uv - 0.5 / resolution) / (1.0 - 1.0 / resolution); }
vec2 RemoveUvMargin(vec2 uv, ivec2 resolution) { return (uv - 0.5 / resolution) / (1.0 - 1.0 / resolution); }
vec3 RemoveUvMargin(vec3 uv, ivec3 resolution) { return (uv - 0.5 / resolution) / (1.0 - 1.0 / resolution); }
vec4 RemoveUvMargin(vec4 uv, ivec4 resolution) { return (uv - 0.5 / resolution) / (1.0 - 1.0 / resolution); }

//----------------------------------------------------------------------------//

float PackUnorm2x4(vec2 xy) {
	return dot(floor(15.0 * xy + 0.5), vec2(1.0 / 255.0, 16.0 / 255.0));
}
float PackUnorm2x4(float x, float y) { return PackUnorm2x4(vec2(x, y)); }
vec2 UnpackUnorm2x4(float pack) {
	vec2 xy; xy.x = modf(pack * 255.0 / 16.0, xy.y);
	return xy * vec2(16.0 / 15.0, 1.0 / 15.0);
}
float UnpackUnorm2x4X(float pack) { return fract(pack * 255.0 / 16.0) * 16.0 / 15.0; }
float UnpackUnorm2x4Y(float pack) { return floor(pack * 255.0 / 16.0) / 15.0; }

float PackUnorm2x8(vec2 xy) {
	return dot(floor(255.0 * xy + 0.5), vec2(1.0 / 65535.0, 256.0 / 65535.0));
}
float PackUnorm2x8(float x, float y) { return PackUnorm2x8(vec2(x, y)); }
vec2 UnpackUnorm2x8(float pack) {
	vec2 xy; xy.x = modf(pack * 65535.0 / 256.0, xy.y);
	return xy * vec2(256.0 / 255.0, 1.0 / 255.0);
}
float UnpackUnorm2x8X(float pack) { return fract(pack * 65535.0 / 256.0) * 256.0 / 255.0; }
float UnpackUnorm2x8Y(float pack) { return floor(pack * 65535.0 / 256.0) / 255.0; }

float PackSnorm2x8(vec2 xy) {
	vec2 xy2 = floor(mix(128.5 + 127.0 * xy, 127.5 + 127.0 * xy, greaterThan(xy, vec2(0.0))));
	return dot(xy2, vec2(1.0 / 65535.0, 256.0 / 65535.0));
}
float PackSnorm2x8(float x, float y) { return PackSnorm2x8(vec2(x, y)); }
vec2 UnpackSnorm2x8(float pack) {
	vec2 xy; xy.x = modf(pack * 65535.0 / 256.0, xy.y);
	return xy * vec2(256.0 / 127.0, 1.0 / 127.0) - mix(vec2(1.0), vec2(128.0 / 127.0), greaterThan(xy, vec2(127.5 / 256.0, 127.5)));
}
float UnpackSnorm2x8X(float pack) { return UnpackSnorm2x8(pack).x; }
float UnpackSnorm2x8Y(float pack) { return UnpackSnorm2x8(pack).y; }

//----------------------------------------------------------------------------//

// Octahedral Unit Vector encoding
// Intuitive, fast, and has very little error.
vec2 EncodeUnitVector(vec3 vector) {
	// Scale down to octahedron, project onto XY plane
	vector.xy /= abs(vector.x) + abs(vector.y) + abs(vector.z);
	// Reflect -Z hemisphere folds over the diagonals
	return vector.z <= 0.0 ? (1.0 - abs(vector.yx)) * vec2(vector.x >= 0.0 ? 1.0 : -1.0, vector.y >= 0.0 ? 1.0 : -1.0) : vector.xy;
}
vec3 DecodeUnitVector(vec2 encoded) {
	// Exctract Z component
	vec3 vector = vec3(encoded, 1.0 - abs(encoded.x) - abs(encoded.y));
	// Reflect -Z hemisphere folds over the diagonals
	float t = max(-vector.z, 0.0);
	vector.xy += vec2(vector.x >= 0.0 ? -t : t, vector.y >= 0.0 ? -t : t);
	// Normalize and return
	return normalize(vector);
}

// Near-32F range, with only as much precision as is usually needed for color.
vec4 EncodeRGBE8(vec3 rgb) {
	float exponentPart = floor(log2(max(max(rgb.r, rgb.g), max(rgb.b, exp2(-127.0))))); // can remove the clamp to above exp2(-127) if you're sure you're not going to input any values below that
	vec3  mantissaPart = clamp((128.0 / 255.0) * exp2(-exponentPart) * rgb, 0.0, 1.0);
	      exponentPart = clamp(exponentPart / 255.0 + (127.0 / 255.0), 0.0, 1.0);

	return vec4(mantissaPart, exponentPart);
}
vec3 DecodeRGBE8(vec4 rgbe) {
	const float add = log2(255.0 / 128.0) - 127.0;
	return exp2(rgbe.a * 255.0 + add) * rgbe.rgb;
}

//----------------------------------------------------------------------------//

vec4 CalculateSphericalHarmonics(vec3 xyz) {
	const vec2 freqW = vec2(sqrt(1.0 / (4.0 * pi)), sqrt(3.0 / (4.0 * pi)));
	return vec4(freqW.x, freqW.y * xyz.yzx);
}

vec3[4] CalculateSphericalHarmonicCoefficients(vec3 value, vec3 xyz) {
	vec4 harmonics = CalculateSphericalHarmonics(xyz);
	return vec3[4](value * harmonics.x, value * harmonics.y, value * harmonics.z, value * harmonics.w);
}
vec3 ValueFromSphericalHarmonicCoefficients(vec3[4] coefficients, vec3 xyz) {
	vec4 harmonics = CalculateSphericalHarmonics(xyz);
	return coefficients[0] * harmonics.x + coefficients[1] * harmonics.y + coefficients[2] * harmonics.z + coefficients[3] * harmonics.w;
}

vec4[3] CalculateSphericalHarmonicCoefficientsAlt(vec3 value, vec3 xyz) {
	vec4 harmonics = CalculateSphericalHarmonics(xyz);
	return vec4[3](value.r * harmonics, value.g * harmonics, value.b * harmonics);
}
vec3 ValueFromSphericalHarmonicCoefficientsAlt(vec4[3] coefficients, vec3 xyz) {
	vec4 harmonics = CalculateSphericalHarmonics(xyz);
	return vec3(dot(coefficients[0], harmonics), dot(coefficients[1], harmonics), dot(coefficients[2], harmonics));
}

//----------------------------------------------------------------------------//

// fast base 2 halton sequence, hopefully decently fast base 3 halton sequence
vec2 Halton23(int n) {
	float halton2 = float(bitfieldReverse(n) >> 9u) / (exp2(23.0) - 1.0);

	float i = float(n);
	float f = 1.0, r = 0.0;
	while (i > 0.0) {
		f /= 3.0;
		r += f * mod(i, 3);
		i  = floor(i / 3);
	}

	return vec2(halton2, r);
}

/*\
 * returns coords & weights for sampling with bicubic filtering
 * use if you need to do so multiple times in the same location (or different locations an integer number of pixels away)
 * or if you need to do some custom stuff to the individual samples
 *
 * coord needs to be in pixels with integers at pixel centers
 * you can convert to this from normal coords with this:
 * coord = coord * resolution - 0.5
 *
 * use c & m like this:
 * mix(mix(texture(sampler, c.xy), texture(sampler, c.zy).x, m.x),
 *     mix(texture(sampler, c.xw), texture(sampler, c.zw).x, m.x),
 *     m.y);
 * remember that c will also be with integers at pixel centers
 * you can convert this to normal coords like this:
 * c = (c + 0.5) / resolution.xyxy
\*/
void FastBicubicCM(vec2 coord, out vec4 c, out vec2 m) {
	vec2 f = fract(coord);
	coord -= f;

	vec2 ff = f * f;

	vec2 w0 = ff * f;
	vec2 w3 = 1.0 - f; w3 *= w3 * w3;
	vec2 w1 = w3 + 6.0 * f - 2.0 * w0;
	vec2 w2 = 3.0 * w0 + 4.0 - 6.0 * ff;

	vec4 s = vec4(w3, w1) + vec4(w2, w0);
	c = coord.xyxy + vec4(w2, w0) / s;
	c.xy -= 1.0; c.zw += 1.0;

	m = s.zw / (s.xy + s.zw);
}
vec4 TextureBicubic(sampler2D sampler, vec2 coord) {
	ivec2 res = textureSize(sampler, 0);

	coord = coord * res - 0.5;

	vec4 c; vec2 m;
	FastBicubicCM(coord, c, m);

	c = (c + 0.5) / res.xyxy;

	return mix(
		mix(texture(sampler, c.xy), texture(sampler, c.zy), m.x),
		mix(texture(sampler, c.xw), texture(sampler, c.zw), m.x),
		m.y
	);
}
