1 module polyplex.math.mathf;
2 import std.traits; // for some introspection (ei isScalar)
3 static import std.math; // for creating public aliases
4 // public imports of math constants
5 public import std.math : E, PI, PI_2, PI_4, M_1_PI, M_2_PI, M_2_SQRTPI, LN10,
6                          LN2, LOG2, LOG2E, LOG2T, LOG10E, SQRT2, SQRT1_2;
7 /// TAU = PI*2
8 public enum real TAU = 6.28318530717958647692;
9 
10 // public aliases of std library math functions to capitalized aliases
11 public alias Abs = std.math.abs;
12 public alias Sin = std.math.sin;
13 public alias Cos = std.math.cos;
14 public alias Tan = std.math.tan;
15 public alias Sqrt = std.math.sqrt;
16 public alias ASin = std.math.asin;
17 public alias ACos = std.math.acos;
18 public alias ATan = std.math.atan;
19 public alias ATan2 = std.math.atan2;
20 public alias SinH = std.math.sinh;
21 public alias CosH = std.math.cosh;
22 public alias TanH = std.math.tanh;
23 public alias ASinH = std.math.asinh;
24 public alias ACosH = std.math.acosh;
25 public alias ATanH = std.math.atanh;
26 public alias Ceil = std.math.ceil;
27 public alias Floor = std.math.floor;
28 public alias Round = std.math.round;
29 public alias Truncate = std.math.trunc;
30 public alias Pow = std.math.pow;
31 public alias Exp = std.math.exp;
32 public alias Exp2 = std.math.exp2;
33 public alias LogNatural = std.math.log;
34 public alias LogBase2 = std.math.log2;
35 public alias LogBase10 = std.math.log10;
36 public alias Fmod = std.math.fmod;
37 public alias Remainder = std.math.remainder;
38 public alias IsFinite = std.math.isFinite;
39 public alias IsIdentical = std.math.isIdentical;
40 public alias IsInfinity = std.math.isInfinity;
41 public alias IsNaN = std.math.isNaN;
42 public alias IsPowerOf2 = std.math.isPowerOf2;
43 public alias Sign = std.math.sgn;
44 
45 
46 /// Checks if two floating-point scalars are approximately
47 /// equal to each other by some epsilon
48 bool ApproxEquals(T)(T a, T b, T eps) pure nothrow if (__traits(isFloating, T)) {
49 	return Abs(a-b) < eps;
50 }
51 unittest {
52 	assert(ApproxEquals(1f, 1.001f, 0.1f));
53 	assert(!ApproxEquals(-1f, 1.001f, 0.1f));
54 	assert(ApproxEquals(10f, 10f, 0.001f));
55 	assert(!ApproxEquals(10f, 10.01f, 0.001f));
56 }
57 
58 /// Converts quantity of degrees to radians
59 T ToRadians(T)(T degrees) pure nothrow {
60 	return (cast(T)(PI)*degrees)/cast(T)180;
61 }
62 /// Converts quantity of radians to degrees
63 T ToDegrees(T)(T radians) pure nothrow {
64 	return (180 * radians)/cast(T)PI;
65 }
66 unittest {
67 	assert(ApproxEquals(ToRadians(90f), ToRadians(ToDegrees(PI_2)), 0.1f));
68 	assert(ApproxEquals(ToRadians(360f+45f)-TAU, PI_4, 0.01f));
69 }
70 
71 /// Minimum of two scalar elements
72 T Min(T)(T scalar_a, T scalar_b) pure nothrow if (__traits(isScalar, T)) {
73 	return ( scalar_a < scalar_b ? scalar_a : scalar_b );
74 }
75 unittest {
76 	assert(Min(5, 10) == 5);
77 	assert(Min(3, -5) == -5);
78 	assert(Min(25.0f, 3) == 3);
79 	assert(!__traits(compiles, Min("a", "b")));
80 	assert(__traits(compiles, Min('a', 'b')));
81 }
82 
83 /// Maximum of two scalar elements
84 T Max(T)(T scalar_a, T scalar_b) pure nothrow if (__traits(isScalar, T)) {
85 	return ( scalar_a > scalar_b ? scalar_a : scalar_b );
86 }
87 unittest {
88 	assert(Max(5, 10) == 10);
89 	assert(Max(3, -5) == 3);
90 	assert(Max(25.0f, 3) == 25.0f);
91 	assert(!__traits(compiles, Max("a", "b")));
92 	assert(__traits(compiles, Max('a', 'b')));
93 }
94 
95 /// Clamps scalar elements
96 T Clamp(T)(T x, T min, T max) pure nothrow if (__traits(isScalar, T)) {
97 	if ( x < min ) return min;
98 	if ( x > max ) return max;
99 	return x;
100 }
101 unittest {
102 	assert(Clamp(3, 0, 10) == 3);
103 	assert(Clamp(3, 5, 10) == 5);
104 	assert(Clamp(3, -2, 2) == 2);
105 }
106 
107 /// Steps scalar `a` on edge `edge`. Returns 0 if a < edge, otherwise 1
108 T Step(T)(T edge, T a) pure nothrow if (__traits(isScalar, T)) {
109 	return (a >= edge);
110 }
111 
112 /// -- All interpolation methods based from
113 ///    http://paulbourke.net/miscellaneous/interpolation/
114 
115 /** Linear interpolation on scalar elements. Interpolates between two scalar
116       values `x` and `y` using a linear gradient `a`
117   Params:
118     x = The minimum element to interpolate by. The result of this function
119           equals `x` when `a = 0`.
120     y = The maximum element to interpolate by. The result of this function
121           equals `x` when `a = 1`.
122     a = The gradient to interpolate between `x` and `y`. Should be a value
123           between `0f` and `1f`, although this function will clamp it to that
124           range. For example, the result of this function will equal (x+y)/2
125           when `a = 0.5`
126 **/
127 T LinearInterpolation(T)(T x, T y, float a) pure nothrow if (__traits(isScalar, T)) {
128 	a = Clamp(a, 0f, 1f);
129 	return cast(T)(x*(1.0f - a) + y*a);
130 }
131 /// Alias for linear interpolation, to follow the same GLSL convention
132 alias Mix = LinearInterpolation;
133 unittest {
134 	assert(LinearInterpolation(0f, 1f, 0.5f) == 0.5f);
135 	assert(LinearInterpolation(0f, 2f, 0.5f) == 1.0f);
136 	assert(LinearInterpolation(0f, 2f, 0.0f) == 0.0f);
137 	assert(LinearInterpolation(0f, 2f, 1.0f) == 2.0f);
138 	assert(LinearInterpolation(0f, 2f, 3.0f) == 2.0f);
139 	assert(LinearInterpolation(0f, 2f, -1.0f) == 0.0f);
140 }
141 
142 /** Cosine interpolation on scalar elements. It interpolates between
143       two scalars `x` and `y` using a cosine-weighted `a` (this function
144       applies the cosine-weight). Check `LinearInterpolation` for more
145       details.
146 **/
147 T CosineInterpolation(T)(T x, T y, float a) pure nothrow if (__traits(isScalar, T)) {
148 	return LinearInterpolation(x, y, (1.0f - Cos(a*PI))/2.0f);
149 }
150 unittest {
151 	assert(__traits(compiles, CosineInterpolation(0, 1, 0.5f)));
152 	assert(__traits(compiles, CosineInterpolation(0f, 1f, 0.5f)));
153 }
154 
155 /** Cube interpolation using Catmull-Rom splines on scalar elements.
156       Interpolates between two scalars `y` and `z` using a gradient `a` that
157       takes `x <-> y` and `z <-> w` into account in order to offer better
158       continuity between segments. Check `LinearInerpolation` for more details.
159 **/
160 T CubicInterpolation(T)(T x, T y, T z, T w, float a) pure nothrow if (__traits(isScalar, T)) {
161 	a = Clamp(a, 0f, 1f);
162 	// Use slope between last and next point as derivative at current point
163 	T cx = cast(T)(-0.5*x + 1.5*y - 1.5*z + 0.5*w);
164 	T cy = cast(T)(x - 2.5*y + 2*z - 0.5*w);
165 	T cz = cast(T)(-0.5*x + 0.5*z);
166 	T cw = y;
167 	return cast(T)(cx*a*a*a + cy*a*a + cz*a + cw);
168 }
169 unittest {
170 	assert(__traits(compiles, CubicInterpolation(0, 1, 2, 3, 0.5f)));
171 	assert(__traits(compiles, CubicInterpolation(0f, 1f, 2f, 3f, 0.5f)));
172 }
173 
174 /** Hermite interpolation on scalar elements as described by GLSL smoothstep.
175       Interpolates between `x` and `y` using gradient `a` that allows
176         smooth transitions as the gradient approaches `1` or `0`. See
177         `LinearInterpolation` for more details.
178 **/
179 T HermiteInterpolation(T)(T x, T y, float a) pure nothrow if (__traits(isScalar, T)) {
180 	T t = cast(T)Clamp((a - x)/(y - x), 0f, 1f);
181 	return t*t*cast(T)(3f - 2f*t);
182 }
183 /// Alias to GLSL name for hermite interpolation
184 alias Smoothstep = HermiteInterpolation;
185 unittest {
186 	assert(__traits(compiles, Smoothstep(0, 1, 0.5f)));
187 	assert(__traits(compiles, Smoothstep(0f, 1f, 0.5f)));
188 }
189 
190 /** Hermite interpolation on scalar elements. Similar to `CubicInterpolation`
191       except that you have control over the tension (tightening of the
192       curvature) and bias (twists the curve about known points).
193   Params:
194     Tension = any value from -1 to 1. -1 is low tension, 1 is high tension
195     Bias = any value from -1 to 1. 0 is no bias, -1 is biased towards first,
196            segment, 1 is biased towards the last segment
197 **/
198 T HermiteInterpolation(T)(T x, T y, T z, T w, float a, float tension, float bias) pure nothrow if (__traits(isScalar, T)) {
199 	a = Clamp(a, 0f, 1f);
200 	float a2 = a*a, a3 = a2*a;
201 	tension = (1f - tension)/2f;
202 
203 	float m_y = (y-z)*(1f+bias)*tension + (z-y)*(1f-bias)*tension;
204 	float m_z = (z-y)*(1f+bias)*tension + (w-z)*(1f-bias)*tension;
205 	float a_x =  2f*a3 - 3f*a2 + 1f;
206 	float a_y =     a3 - 2f*a2 + a;
207 	float a_z =     a3 -    a2;
208 	float a_w = -2f*a3 + 3f*a2;
209 
210 	return cast(T)(a_x*y + a_y*m_y + a_z*m_z + a_w*z);
211 }
212 unittest {
213 	assert(__traits(compiles, HermiteInterpolation(0, 1, 2, 3, 0.5f, 0.2f, 0f)));
214 	assert(__traits(compiles, HermiteInterpolation(0f, 1f, 2f, 3f, 0.5f, 0.2f, 0f)));
215 }
216 
217