First commit

Signed-off-by: Mattias Andrée <maandree@kth.se>

1 files changed, 609 insertions, 0 deletions

diff --git a/equations.c b/equations.c
new file mode 100644
index 0000000..571fc26
--- /dev/null
+++ b/equations.c
@@ -0,0 +1,609 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+CONST_FUNCTION
+static double
+refine_cubic_solution(double t, double a, double b, double c, double d)
+{
+	size_t i;
+	double x;
+	for (i = 0; x = t*t*t*a + t*t*b + t*c + d, !iszeroish(x) && i < 10; i++)
+		t -= x / (t*t*(a*3) + t*(b*2) + c);
+	return t;
+}
+
+
+static size_t
+solve_cubic_with_1_known_root(double *ts, double a, double b, double t1)
+{
+	/*
+	 * (t³ + t²a + tb + c) / (t - t₁) = t² + tp + q
+	 * t³ + t²a + tb + c = (t² + tp + q)(t - t₁)
+	 * t³ + t²a + tb + c = t²(t - t₁) + tp(t - t₁) + q(t - t₁)
+	 * t³ + t²a + tb + c = t³ - t²t₁ + t²p - tpt₁ + tq - qt₁
+	 * t³ + t²a + tb + c = t³ + t²(p - t₁) + t(q - pt₁) - qt₁
+	 *
+	 * ⎧a = p - t₁
+	 * ⎨b = q - pt₁
+	 * ⎩c = -qt₁
+	 *
+	 * ⎧a + t₁ = p
+	 * ⎨b + pt₁ = q
+	 * ⎩c = -qt₁
+	 *
+	 * ⎰a + t₁ = p
+	 * ⎱b + pt₁ = q
+	 *
+	 * ⎰p = a + t₁
+	 * ⎱q = b + pt₁
+	 */
+	double p = a + t1;
+	double q = b + p * t1;
+
+	return solve_quadratic(ts, 1, p, q);
+}
+
+
+static size_t
+solve_cubic_with_2_known_roots(double *ts, double a, double b, double c, double t1, double t2)
+{
+	/*
+	 * (t - t₁)(t - t₂)(t - t₃)
+	 * (t - t₁)t(t - t₃) - (t - t₁)t₂(t - t₃)
+	 * (t - t₁)(t² - t₃t) - (t - t₁)(t₂t - t₂t₃)
+	 * (t - t₁)t² - (t - t₁)t₃t - (t - t₁)t₂t + (t - t₁)t₂t₃
+	 * t³ - t₁t² - t₃t² + t₁t₃t - t₂t² + t₁t₂t + t₂t₃t - t₁t₂t₃
+	 * t³ - t²(t₁ + t₂ + t₃) + t(t₁t₂ + t₁t₃ + t₂t₃) - t₁t₂t₃
+	 * 
+	 * ⎧a = -t₁ - t₂ - t₃
+	 * ⎨b = t₁t₂ + t₁t₃ + t₂t₃
+	 * ⎩c = -t₁t₂t₃
+	 * 
+	 * ⎧t₃ = -t₁ - t₂ - a
+	 * ⎨t₃ = (b - t₁t₂)/(t₁ + t₂)
+	 * ⎩t₃ = -c/t₁t₂
+	 */
+	double t12sum  = t1 + t2;
+	double t12prod = t1 * t2;
+	double t3a = -a - t12sum;
+	double t3b = iszeroish(t12sum) ? t3a : (b - t12prod) / t12sum;
+	double t3c = iszeroish(t12prod) ? t3a : -c / t12prod;
+	if (iszeroish(t3a - t3b) && iszeroish(t3a - t3c)) {
+		ts[0] = (t3a + t3b + t3c) / 3.0;
+		return 0.0 <= ts[0] && ts[0] <= 1.0;
+	}
+	t1 = 0.5 * (t1 + t2);
+	t1 = refine_cubic_solution(t1, 1.0, a, b, c);
+	return solve_cubic_with_1_known_root(ts, a, b, t1);
+}
+
+
+size_t
+solve_cubic(double *ts, double a, double b, double c, double d)
+{
+	double t1, t2, t3;
+	unsigned int solutions = 0;
+	size_t n = 0;
+
+	double u = 1 / (a * -3);
+
+	double D0 = b*b - 3*(a*c);
+	double D1 = b*b*b - (4.5*b)*(a*c) + (13.5*a)*(a*d);
+
+	double kk = D1*D1 - D0*D0*D0;
+
+	double r, i;
+
+	double Kr;
+	if (kk >= 0) {
+		Kr = D1 < 0 ? cbrt(D1 + sqrt(kk)) : cbrt(D1 - sqrt(kk));
+	} else {
+		double ki = sqrt(-kk);
+		double h = pow(D1*D1 + ki*ki, 1.0 / 6.0);
+		double v = (0.1 / 0.3) * atan2(ki, D1);
+
+		Kr = h * cos(v);
+
+		if (1 /* || fabs(cos(v)) > 1.0 - TOLERANCE */) {
+			double Ki = h * sin(v);
+			double t = 1.0 / (Kr * Kr + Ki * Ki);
+			double m = t * D0 + 1;
+			if (iszeroish(Ki - D0 * Ki * t)) {
+				ts[n] = u * (b + Kr + D0 * Kr * t);
+				ts[n] = t1 = refine_cubic_solution(ts[n], a, b, c, d);
+				solutions |= 1;
+				n += 0.0 <= ts[n] && ts[n] <= 1.0;
+			}
+			if (iszeroish(Ki - sqrt(3.0) * Kr)) {
+				ts[n] = u * (b - (0.5 * Kr + sqrt(0.75) * Ki) * m);
+				ts[n] = t2 = refine_cubic_solution(ts[n], a, b, c, d);
+				solutions |= 2;
+				n += 0.0 <= ts[n] && ts[n] <= 1.0;
+			}
+			if (iszeroish(Ki + sqrt(3.0) * Kr)) {
+				ts[n] = u * (b - (0.5 * Kr - sqrt(0.75) * Ki) * m);
+				ts[n] = t3 = refine_cubic_solution(ts[n], a, b, c, d);
+				solutions |= 4;
+				n += 0.0 <= ts[n] && ts[n] <= 1.0;
+			}
+
+			/* Sometimes the maths doesn't work well enough */
+			if (solutions != 7) {
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic push
+# pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
+#elif defined(__clang__)
+# pragma clang diagnostic push
+# pragma clang diagnostic ignored "-Wconditional-uninitialized"
+#endif
+				switch (solutions) {
+				case 1:
+					n += solve_cubic_with_1_known_root(&ts[n], b / a, c / a, t1);
+					break;
+				case 2:
+					n += solve_cubic_with_1_known_root(&ts[n], b / a, c / a, t2);
+					break;
+				case 4:
+					n += solve_cubic_with_1_known_root(&ts[n], b / a, c / a, t3);
+					break;
+				case 7 - 1:
+					n += solve_cubic_with_2_known_roots(&ts[n], b / a, c / a, d / a, t2, t3);
+					break;
+				case 7 - 2:
+					n += solve_cubic_with_2_known_roots(&ts[n], b / a, c / a, d / a, t1, t3);
+					break;
+				case 7 - 4:
+					n += solve_cubic_with_2_known_roots(&ts[n], b / a, c / a, d / a, t1, t2);
+					break;
+				}
+#if defined(__GNUC__) && !defined(__clang__)
+# pragma GCC diagnostic pop
+#elif defined(__clang__)
+# pragma clang diagnostic pop
+#endif
+			}
+			return n;
+		}
+	}
+
+	if (!iszeroish(Kr)) {
+		r = Kr + D0 / Kr;
+		i = u * sqrt(0.75) * (Kr - D0 / Kr);
+
+	real_case:
+		ts[n] = u * (b + r);
+		ts[n] = t1 = refine_cubic_solution(ts[n], a, b, c, d);
+		n += 0.0 <= ts[n] && ts[n] <= 1.0;
+		if (iszeroish(i)) {
+			ts[n] = u * (b - 0.5 * r);
+			ts[n] = refine_cubic_solution(ts[n], a, b, c, d);
+			n += 0.0 <= ts[n] && ts[n] <= 1.0;
+			return n;
+		} else {
+			n += solve_cubic_with_1_known_root(&ts[n], b / a, c / a, t1);
+			return n;
+		}
+	} else {
+		double D0divKr = D0 / Kr;
+		r = Kr + D0divKr;
+		i = u * sqrt(0.75) * (Kr - D0divKr);
+		if (isfinite(D0divKr)) /* NOT -ffast-math compatible */
+			goto real_case;
+
+		ts[n] = u * b;
+		ts[n] = t1 = refine_cubic_solution(ts[n], a, b, c, d);
+		n += 0.0 <= ts[n] && ts[n] <= 1.0;
+		n += solve_cubic_with_1_known_root(&ts[n], b / a, c / a, t1);
+		return n;
+	}
+}
+
+
+size_t
+solve_quadratic(double *ts, double a, double b, double c)
+{
+	size_t n = 0;
+
+	double negahalf_p = -0.5 * (b / a);
+	double q = c / a;
+
+	double D = negahalf_p * negahalf_p - q;
+
+	if (D < 0) {
+		if (D < -TOLERANCE)
+			return 0;
+
+		ts[n] = negahalf_p;
+		n += (0.0 <= ts[n] && ts[n] <= 1.0);
+
+		return n;
+	} else {
+		ts[n] = negahalf_p + sqrt(D);
+		n += (0.0 <= ts[n] && ts[n] <= 1.0);
+
+		ts[n] = negahalf_p - sqrt(D);
+		n += (0.0 <= ts[n] && ts[n] <= 1.0);
+
+		return n;
+	}
+}
+
+
+size_t
+solve_linear(double *ts, double a, double b)
+{
+	ts[0] = -b / a;
+	return (0.0 <= ts[0] && ts[0] <= 1.0);
+}
+
+
+#else
+
+
+#define Re(X) X, 0
+#define Im(X) 0, X
+#define Cx(R, I) R, I
+
+static void
+check_cubic(double r1, double i1, double r2, double i2, double r3, double i3, double k, int strict_real)
+{
+	size_t nexpected = 0, ngot, i, j;
+	double expected[3], got[3];
+	double error, maxerror = 0;
+	double ra, rb, rc;
+	double ia, ib, ic;
+	int found[3], correct;
+
+	if (iszeroish(i1))
+		expected[nexpected++] = r1;
+	if (iszeroish(i2))
+		expected[nexpected++] = r2;
+	if (iszeroish(i3))
+		expected[nexpected++] = r3;
+
+	/*
+	 * k(t₁ - t)(t₂ - t)(t₃ - t) = 0
+	 * k(t₁ - t)t₂(t₃ - t) - k(t₁ - t)t(t₃ - t) = 0
+	 * kt₂(t₁ - t)(t₃ - t) - kt(t₁ - t)(t₃ - t) = 0
+	 * kt₂(t₁(t₃ - t) - t(t₃ - t)) - kt(t₁(t₃ - t) - t(t₃ - t)) = 0
+	 * kt₂(t₁t₃ - t₁t - tt₃ + t²) - kt(t₁t₃ - t₁t - tt₃ + t²) = 0
+	 * kt₂(t₁t₃ - t₁t - tt₃ + t²) + kt(-t₁t₃ + t₁t + tt₃ - t²) = 0
+	 * kt₂t₁t₃ - kt₂t₁t - kt₂tt₃ + kt₂t² - ktt₁t₃ + ktt₁t + kt²t₃ - kt³ = 0
+	 * kt₂t₁t₃ - kt₂t₁t¹ - kt₂t¹t₃ + kt₂t² - kt¹t₁t₃ + kt²t₁ + kt²t₃ - kt³ = 0
+	 * -kt³ + kt₂t² + kt²t₁ + kt²t₃ - kt₂t₁t¹ - kt₂t¹t₃ - kt¹t₁t₃ + kt₂t₁t₃ = 0
+	 * kt³ - kt₂t² - kt²t₁ - kt²t₃ + kt₂t₁t¹ + kt₂t¹t₃ + kt¹t₁t₃ - kt₂t₁t₃ = 0
+	 * t³k - t²k(t₁ + t₂ + t₃) + t¹k(t₂t₁ + t₂t₃ + t₁t₃) - kt₂t₁t₃ = 0
+	 * ⎧a := t₁ + t₂ + t₃      ⎫
+	 * ⎨b := t₁t₂ + t₁t₃ + t₂t₃⎬
+	 * ⎩c := t₁t₂t₃            ⎭
+	 * t³k - t²ka + t¹kb - kc = 0
+	 */
+
+	ra = r1 + r2 + r3;
+	ia = i1 + i2 + i3;
+	rb = r1*r2 - i1*i2 + r1*r3 - i1*i3 + r2*r3 - i2*i3;
+	ib = r1*i2 + i1*r2 + r1*i3 + i1*r3 + r2*i3 + i2*r3;
+	rc = r1*r2*r3 - r1*i2*i3 - i1*r2*i3 - i1*i2*r3;
+	ic = r1*r2*i3 + r1*i2*r3 + i1*r2*r3 - i1*i2*i3;
+
+	if (strict_real) {
+		ASSERT(iszeroish(ia));
+		ASSERT(iszeroish(ib));
+		ASSERT(iszeroish(ic));
+	} else {
+		if (!iszeroish(ia) || !iszeroish(ib) || !iszeroish(ic))
+			return;
+	}
+
+	ngot = solve_cubic(got, k, -k*ra, k*rb, -k*rc);
+
+	for (i = 0; i < sizeof(found) / sizeof(*found); i++)
+		found[i] = (i >= nexpected || expected[i] < 0.0 + TOLERANCE || expected[i] > 1.0 - TOLERANCE);
+	if (!nexpected) {
+		ASSERT(!ngot);
+		return;
+	}
+	for (i = 0; i < ngot; i++) {
+		correct = (got[i] < 0.0 + TOLERANCE || got[i] > 1.0 - TOLERANCE);
+		error = fabs(got[i] - expected[0]);
+		for (j = 0; j < nexpected; j++) {
+			error = fmin(error, fabs(got[i] - expected[j]));
+			if (tolerant_eq(got[i], expected[j]))
+				found[j] = correct = 1;
+		}
+		if (!correct)
+			maxerror = fmax(maxerror, error);
+	}
+	if (!found[0] || !found[1] || !found[2] || !tolerant_eq(maxerror, 0)) {
+		fprintf(stderr, "\n%lg, %lg, %lg, %lg, %lg, %lg, %lg\n", r1, i1, r2, i2, r3, i3, k);
+		fprintf(stderr, "\t(%la, %la, %la, %la, %la, %la, %la)\n", r1, i1, r2, i2, r3, i3, k);
+		fprintf(stderr, "\tmaxerror: %lf (%la)\n", maxerror, maxerror);
+		fprintf(stderr, "\tfound: %i, %i, %i\n", found[0], found[1], found[2]);
+		fprintf(stderr, "\tgot: %zu\n", ngot);
+		for (i = 0; i < ngot; i++)
+			fprintf(stderr, "\t\t%lf (%la)\n", got[i], got[i]);
+		fprintf(stderr, "\texpected: %zu\n", nexpected);
+		for (i = 0; i < nexpected; i++)
+			fprintf(stderr, "\t\t%lf (%la)\n", expected[i], expected[i]);
+	}
+	ASSERT(found[0]);
+	ASSERT(found[1]);
+	ASSERT(found[2]);
+	ASSERT(tolerant_eq(maxerror, 0));
+}
+
+
+static void
+check_quadratic(double r1, double i1, double r2, double i2, double k, int strict_real)
+{
+	size_t nexpected = 0, ngot, i, j;
+	double expected[2], got[2];
+	double error, maxerror = 0;
+	double ra, rb;
+	double ia, ib;
+	int found[2], correct;
+
+	if (iszeroish(i1))
+		expected[nexpected++] = r1;
+	if (iszeroish(i2))
+		expected[nexpected++] = r2;
+
+	/*
+	 * k(t₁ - t)(t₂ - t) = 0
+	 * k(t₁(t₂ - t) - t(t₂ - t)) = 0
+	 * k(t₁t₂ - t¹t₁ - t¹t₂ + t²) = 0
+	 * kt₁t₂ - t¹kt₁ - t¹kt₂ + kt² = 0
+	 * t²k - t¹kt₁ - t¹kt₂ + kt₁t₂ = 0
+	 * t²k - t¹k(t₁ + kt₂) + kt₁t₂ = 0
+	 * ⎰a := t₁ + t₂⎱
+	 * ⎱b := t₁t₂   ⎰
+	 * t²k - t¹ka + kb = 0
+	 */
+
+	ra = r1 + r2;
+	ia = i1 + i2;
+	rb = r1*r2 - i1*i2;
+	ib = r1*i2 + i1*r2;
+
+	if (strict_real) {
+		ASSERT(iszeroish(ia));
+		ASSERT(iszeroish(ib));
+	} else {
+		if (!iszeroish(ia) || !iszeroish(ib))
+			return;
+	}
+
+	ngot = solve_quadratic(got, k, -k*ra, k*rb);
+
+	for (i = 0; i < sizeof(found) / sizeof(*found); i++)
+		found[i] = (i >= nexpected || expected[i] < 0.0 + TOLERANCE || expected[i] > 1.0 - TOLERANCE);
+	if (!nexpected) {
+		ASSERT(!ngot);
+		return;
+	}
+	for (i = 0; i < ngot; i++) {
+		correct = (got[i] < 0.0 + TOLERANCE || got[i] > 1.0 - TOLERANCE);
+		error = fabs(got[i] - expected[0]);
+		for (j = 0; j < nexpected; j++) {
+			error = fmin(error, fabs(got[i] - expected[j]));
+			if (tolerant_eq(got[i], expected[j]))
+				found[j] = correct = 1;
+		}
+		if (!correct)
+			maxerror = fmax(maxerror, error);
+	}
+	if (!found[0] || !found[1] || !tolerant_eq(maxerror, 0)) {
+		fprintf(stderr, "\n%lg, %lg, %lg, %lg, %lg\n", r1, i1, r2, i2, k);
+		fprintf(stderr, "\t(%la, %la, %la, %la, %la)\n", r1, i1, r2, i2, k);
+		fprintf(stderr, "\tmaxerror: %lf (%la)\n", maxerror, maxerror);
+		fprintf(stderr, "\tfound: %i, %i\n", found[0], found[1]);
+		fprintf(stderr, "\tgot: %zu\n", ngot);
+		for (i = 0; i < ngot; i++)
+			fprintf(stderr, "\t\t%lf (%la)\n", got[i], got[i]);
+		fprintf(stderr, "\texpected: %zu\n", nexpected);
+		for (i = 0; i < nexpected; i++)
+			fprintf(stderr, "\t\t%lf (%la)\n", expected[i], expected[i]);
+	}
+	ASSERT(found[0]);
+	ASSERT(found[1]);
+	ASSERT(tolerant_eq(maxerror, 0));
+}
+
+
+static void
+check_linear(double r1, double i1, double k, int strict_real)
+{
+	size_t nexpected = 0, ngot, i;
+	double expected[1] = {999999}, got[1];
+	double ra;
+	double ia;
+
+	if (iszeroish(i1) && 0.0 <= r1 && r1 <= 1.0)
+		expected[nexpected++] = r1;
+
+	/*
+	 * k(t - t₁) = 0
+	 * kt - kt₁ = 0
+	 */
+
+	ra = r1;
+	ia = i1;
+
+	if (strict_real) {
+		ASSERT(iszeroish(ia));
+	} else {
+		if (!iszeroish(ia))
+			return;
+	}
+
+	ngot = solve_linear(got, k, -k*ra);
+
+	ASSERT(ngot == nexpected);
+	for (i = 0; i < ngot; i++)
+		ASSERT(eq(got[i], expected[i]));
+}
+
+
+static void
+check_cubic_randomly(void)
+{
+	double r1 = rand() / (double)RAND_MAX;
+	double i1 = rand() / (double)RAND_MAX;
+	double r2 = rand() / (double)RAND_MAX;
+	double i2 = rand() / (double)RAND_MAX;
+	double r3 = rand() / (double)RAND_MAX;
+	double i3 = rand() / (double)RAND_MAX;
+	double k = rand() / (double)RAND_MAX;
+	r1 = r1 * 3 - 1;
+	i1 = i1 * 3 - 1;
+	r2 = r2 * 3 - 1;
+	i2 = i2 * 3 - 1;
+	r3 = r3 * 3 - 1;
+	i3 = i3 * 3 - 1;
+	k = k * 3 - 1;
+	k += iszeroish(k);
+
+	check_cubic(Cx(r1, i1), Cx(r2, i2), Cx(r3, i3), k, 0);
+	check_cubic(Im(i1), Cx(r2, i2), Cx(r3, i3), k, 0);
+	check_cubic(Re(r1), Cx(r2, i2), Cx(r3, i3), k, 0);
+	check_cubic(Cx(r1, i1), Im(i2), Cx(r3, i3), k, 0);
+	check_cubic(Im(i1), Im(i2), Cx(r3, i3), k, 0);
+	check_cubic(Re(r1), Im(i2), Cx(r3, i3), k, 0);
+	check_cubic(Cx(r1, i1), Re(r2), Cx(r3, i3), k, 0);
+	check_cubic(Im(i1), Re(r2), Cx(r3, i3), k, 0);
+	check_cubic(Re(r1), Re(r2), Cx(r3, i3), k, 0);
+	check_cubic(Cx(r1, i1), Cx(r1, i1), Cx(r3, i3), k, 0);
+	check_cubic(Im(i1), Im(i1), Cx(r3, i3), k, 0);
+	check_cubic(Re(r1), Re(r1), Cx(r3, i3), k, 0);
+	check_cubic(Cx(r1, i1), Cx(r3, i3), Cx(r3, i3), k, 0);
+	check_cubic(Im(i1), Cx(r3, i3), Cx(r3, i3), k, 0);
+	check_cubic(Re(r1), Cx(r3, i3), Cx(r3, i3), k, 0);
+
+	check_cubic(Cx(r1, i1), Cx(r2, i2), Im(i3), k, 0);
+	check_cubic(Im(i1), Cx(r2, i2), Im(i3), k, 0);
+	check_cubic(Re(r1), Cx(r2, i2), Im(i3), k, 0);
+	check_cubic(Cx(r1, i1), Im(i2), Im(i3), k, 0);
+	check_cubic(Im(i1), Im(i2), Im(i3), k, 0);
+	check_cubic(Re(r1), Im(i2), Im(i3), k, 0);
+	check_cubic(Cx(r1, i1), Re(r2), Im(i3), k, 0);
+	check_cubic(Im(i1), Re(r2), Im(i3), k, 0);
+	check_cubic(Re(r1), Re(r2), Im(i3), k, 0);
+	check_cubic(Cx(r1, i1), Cx(r1, i1), Im(i3), k, 0);
+	check_cubic(Im(i1), Im(i1), Im(i3), k, 0);
+	check_cubic(Re(r1), Re(r1), Im(i3), k, 0);
+	check_cubic(Cx(r1, i1), Im(i3), Im(i3), k, 0);
+	check_cubic(Im(i1), Im(i3), Im(i3), k, 0);
+	check_cubic(Re(r1), Im(i3), Im(i3), k, 0);
+
+	check_cubic(Cx(r1, i1), Cx(r2, i2), Re(r3), k, 0);
+	check_cubic(Im(i1), Cx(r2, i2), Re(r3), k, 0);
+	check_cubic(Re(r1), Cx(r2, i2), Re(r3), k, 0);
+	check_cubic(Cx(r1, i1), Im(i2), Re(r3), k, 0);
+	check_cubic(Im(i1), Im(i2), Re(r3), k, 0);
+	check_cubic(Re(r1), Im(i2), Re(r3), k, 0);
+	check_cubic(Cx(r1, i1), Re(r2), Re(r3), k, 0);
+	check_cubic(Im(i1), Re(r2), Re(r3), k, 0);
+	check_cubic(Re(r1), Re(r2), Re(r3), k, 1);
+	check_cubic(Cx(r1, i1), Cx(r1, i1), Re(r3), k, 0);
+	check_cubic(Im(i1), Im(i1), Re(r3), k, 0);
+	check_cubic(Re(r1), Re(r1), Re(r3), k, 1);
+	check_cubic(Cx(r1, i1), Re(r3), Re(r3), k, 0);
+	check_cubic(Im(i1), Re(r3), Re(r3), k, 0);
+	check_cubic(Re(r1), Re(r3), Re(r3), k, 1);
+
+	check_cubic(Cx(r1, i1), Cx(r1, i1), Cx(r1, i1), k, 0);
+	check_cubic(Im(i1), Im(i1), Im(i1), k, 0);
+	check_cubic(Re(r1), Re(r1), Re(r1), k, 1);
+}
+
+
+static void
+check_quadratic_randomly(void)
+{
+	double r1 = rand() / (double)RAND_MAX;
+	double i1 = rand() / (double)RAND_MAX;
+	double r2 = rand() / (double)RAND_MAX;
+	double i2 = rand() / (double)RAND_MAX;
+	double k = rand() / (double)RAND_MAX;
+	r1 = r1 * 3 - 1;
+	i1 = i1 * 3 - 1;
+	r2 = r2 * 3 - 1;
+	i2 = i2 * 3 - 1;
+	k = k * 3 - 1;
+	k += iszeroish(k);
+
+	check_quadratic(Cx(r1, i1), Cx(r2, i2), k, 0);
+	check_quadratic(Im(i1), Cx(r2, i2), k, 0);
+	check_quadratic(Re(r1), Cx(r2, i2), k, 0);
+	check_quadratic(Cx(r1, i1), Im(i2), k, 0);
+	check_quadratic(Im(i1), Im(i2), k, 0);
+	check_quadratic(Re(r1), Im(i2), k, 0);
+	check_quadratic(Cx(r1, i1), Re(r2), k, 0);
+	check_quadratic(Im(i1), Re(r2), k, 0);
+	check_quadratic(Re(r1), Re(r2), k, 1);
+	check_quadratic(Cx(r1, i1), Cx(r1, i1), k, 0);
+	check_quadratic(Im(i1), Im(i1), k, 0);
+	check_quadratic(Re(r1), Re(r1), k, 1);
+}
+
+
+static void
+check_linear_randomly(void)
+{
+	double r1 = rand() / (double)RAND_MAX;
+	double i1 = rand() / (double)RAND_MAX;
+	double k = rand() / (double)RAND_MAX;
+	r1 = r1 * 3 - 1;
+	i1 = i1 * 3 - 1;
+	k = k * 3 - 1;
+	k += iszeroish(k);
+
+	check_linear(Cx(r1, i1), k, 0);
+	check_linear(Im(i1), k, 0);
+	check_linear(Re(r1), k, 0);
+}
+
+
+int
+main(void)
+{
+	size_t i;
+
+	void *r = malloc(1);
+	srand((unsigned)(uintptr_t)r);
+	free(r);
+
+	check_linear(Re(0), 1.0, 1);
+	check_linear(Re(0), 1.5, 1);
+	check_linear(Re(0), 0.5, 1);
+	check_linear(Re(0), -0.5, 1);
+	check_linear(Re(0), -1.0, 1);
+	check_linear(Re(0), -1.5, 1);
+	for (i = 0; i < 100000UL; i++)
+		check_linear_randomly();
+
+	check_quadratic(Re(0), Re(0), 1.0, 1);
+	check_quadratic(Re(0), Re(0), 1.5, 1);
+	check_quadratic(Re(0), Re(0), 0.5, 1);
+	check_quadratic(Re(0), Re(0), -0.5, 1);
+	check_quadratic(Re(0), Re(0), -1.0, 1);
+	check_quadratic(Re(0), Re(0), -1.5, 1);
+	for (i = 0; i < 10000000UL; i++)
+		check_quadratic_randomly();
+
+	check_cubic(Re(0), Re(0), Re(0), 1.0, 1);
+	check_cubic(Re(0), Re(0), Re(0), 1.5, 1);
+	check_cubic(Re(0), Re(0), Re(0), 0.5, 1);
+	check_cubic(Re(0), Re(0), Re(0), -0.5, 1);
+	check_cubic(Re(0), Re(0), Re(0), -1.0, 1);
+	check_cubic(Re(0), Re(0), Re(0), -1.5, 1);
+	for (i = 0; i < 1000000UL; i++)
+		check_cubic_randomly();
+
+	return 0;
+}
+
+
+#endif