aboutsummaryrefslogtreecommitdiffstats
path: root/equations.c
diff options
context:
space:
mode:
Diffstat (limited to '')
-rw-r--r--equations.c609
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
generated by cgit v1.2.3-70-g09d2 (git 2.47.0) at 2024-11-28 04:49:45 +0000