mirror of https://github.com/Chlumsky/msdfgen.git
79 lines
2.1 KiB
C++
79 lines
2.1 KiB
C++
|
|
#include "equation-solver.h"
|
|
|
|
#define _USE_MATH_DEFINES
|
|
#include <cmath>
|
|
|
|
#define LARGE_RATIO ::msdfgen::real(1e10)
|
|
|
|
#ifndef MSDFGEN_CUBE_ROOT
|
|
#define MSDFGEN_CUBE_ROOT(x) pow((x), ::msdfgen::real(1)/::msdfgen::real(3))
|
|
#endif
|
|
|
|
namespace msdfgen {
|
|
|
|
int solveQuadratic(real x[2], real a, real b, real c) {
|
|
// a == 0 -> linear equation
|
|
if (a == real(0) || fabs(b) > LARGE_RATIO*fabs(a)) {
|
|
// a == 0, b == 0 -> no solution
|
|
if (b == real(0)) {
|
|
if (c == real(0))
|
|
return -1; // 0 == 0
|
|
return 0;
|
|
}
|
|
x[0] = -c/b;
|
|
return 1;
|
|
}
|
|
real dscr = b*b-real(4)*a*c;
|
|
if (dscr > real(0)) {
|
|
dscr = sqrt(dscr);
|
|
x[0] = (-b+dscr)/(a+a);
|
|
x[1] = (-b-dscr)/(a+a);
|
|
return 2;
|
|
} else if (dscr == 0) {
|
|
x[0] = -b/(a+a);
|
|
return 1;
|
|
} else
|
|
return 0;
|
|
}
|
|
|
|
static int solveCubicNormed(real x[3], real a, real b, real c) {
|
|
real a2 = a*a;
|
|
real q = real(1)/real(9)*(a2-real(3)*b);
|
|
real r = real(1)/real(54)*(a*(real(2)*a2-real(9)*b)+real(27)*c);
|
|
real r2 = r*r;
|
|
real q3 = q*q*q;
|
|
a *= real(1)/real(3);
|
|
if (r2 < q3) {
|
|
real t = r/sqrt(q3);
|
|
if (t < real(-1)) t = -1;
|
|
if (t > real(1)) t = 1;
|
|
t = real(1)/real(3)*acos(t);
|
|
q = real(-2)*sqrt(q);
|
|
x[0] = q*cos(t)-a;
|
|
x[1] = q*cos(t+real(2)/real(3)*real(M_PI))-a;
|
|
x[2] = q*cos(t-real(2)/real(3)*real(M_PI))-a;
|
|
return 3;
|
|
} else {
|
|
real u = (r < real(0) ? real(1) : real(-1))*MSDFGEN_CUBE_ROOT(fabs(r)+sqrt(r2-q3));
|
|
real v = u == real(0) ? real(0) : q/u;
|
|
x[0] = (u+v)-a;
|
|
if (u == v || LARGE_RATIO*fabs(u-v) < fabs(u+v)) {
|
|
x[1] = real(-.5)*(u+v)-a;
|
|
return 2;
|
|
}
|
|
return 1;
|
|
}
|
|
}
|
|
|
|
int solveCubic(real x[3], real a, real b, real c, real d) {
|
|
if (a != 0) {
|
|
real bn = b/a;
|
|
if (bn*bn < LARGE_RATIO) // Above this ratio, the numerical error gets larger than if we treated a as zero
|
|
return solveCubicNormed(x, bn, c/a, d/a);
|
|
}
|
|
return solveQuadratic(x, b, c, d);
|
|
}
|
|
|
|
}
|