initial checkin of C++ port (in progress)

git-svn-id: svn://svn.code.sf.net/p/magicseteditor/code/trunk@2 0fc631ac-6414-0410-93d0-97cfa31319b6

src/gfx/polynomial.cpp

@@ -0,0 +1,87 @@
+//+----------------------------------------------------------------------------+
+//| Description:  Magic Set Editor - Program to make Magic (tm) cards          |
+//| Copyright:    (C) 2001 - 2006 Twan van Laarhoven                           |
+//| License:      GNU General Public License 2 or later (see file COPYING)     |
+//+----------------------------------------------------------------------------+
+
+// ----------------------------------------------------------------------------- : Includes
+
+#include <gfx/polynomial.hpp>
+#include <complex>
+
+// ----------------------------------------------------------------------------- : Solving
+
+UInt solveLinear(double a, double b, double* root) {
+	if (a == 0) {
+		if (b == 0) {
+			root[0] = 0;
+			return 1;
+		} else {
+			return 0;
+		}
+	} else {
+		root[0] = -b / a;
+		return 1;
+	}
+}
+
+UInt solveQuadratic(double a, double b, double c, double* roots) {
+	if (a == 0) {
+		return solveLinear(b, c, roots);
+	} else {
+		double d = b*b - 4*a*c;
+		if (d < 0)  return 0;
+		roots[0] = (-b - sqrt(d)) / (2*a);
+		roots[1] = (-b + sqrt(d)) / (2*a);
+		return 2;
+	}
+}
+
+UInt solveCubic(double a, double b, double c, double d, double* roots) {
+	if (a == 0) {
+		return solveQuadratic(b, c, d, roots);
+	} else {
+		return solveCubic(b/a, c/a, d/a, roots);
+	}
+}
+
+// cubic root
+template <typename T>
+inline T curt(T x) { return pow(x, 1.0 / 3); }
+
+UInt solveCubic(double a, double b, double c, double* roots) {
+	double p = b - a*a / 3;
+	double q = c + (2 * a*a*a - 9 * a * b) / 27;
+	if (p == 0 && q == 0) {
+		roots[0] = -a / 3;
+		return 1;
+	}
+	complex<double> u;
+	if (q > 0) {
+		u = curt(q/2 + sqrt(complex<double>(q*q / 4  +  p*p*p / 27)));
+	} else {
+		u = curt(q/2 - sqrt(complex<double>(q*q / 4  +  p*p*p / 27)));
+	}
+	// now for the complex part
+	//              rot1(1,    0)
+	complex<double> rot2(-0.5, sqrt(3.0) / 2);
+	complex<double> rot3(-0.5, -sqrt(3.0) / 2);
+	complex<double> x1 = p / (3.0 * u)         -  u        -  a / 3.0;
+	complex<double> x2 = p / (3.0 * u * rot2)  -  u * rot2 -  a / 3.0;
+	complex<double> x3 = p / (3.0 * u * rot3)  -  u * rot3 -  a / 3.0;
+	// check if the solutions are real
+	UInt count = 0;
+	if (abs(x1.imag()) < 0.00001) {
+		roots[count] = x1.real();
+		count += 1;
+	}
+	if (abs(x2.imag()) < 0.00001) {
+		roots[count] = x2.real();
+		count += 1;
+	}
+	if (abs(x3.imag()) < 0.00001) {
+		roots[count] = x3.real();
+		count += 1;
+	}
+	return count;
+}