main/src/georegression/fitting/curves/FitCurve_F64.java
public static boolean fitMM(double[] data, int offset , int length ,
								PolynomialQuadratic1D_F64 output , @Nullable DMatrix3x3 work ) {
		if( work == null )
			work = new DMatrix3x3();

		final int N = length/2;

		// Unrolled pseudo inverse
		// coef = inv(A^T*A)*A^T*y
		double sx0=N,sx1=0,sx2=0;
		double             sx3=0;
		double             sx4=0;

		double b0=0,b1=0,b2=0;

		int end = offset+length;
		for (int i = offset; i < end; i += 2 ) {
			double x = data[i];
			double y = data[i+1];

			double x2 = x*x;

			sx1 += x;
			sx2 += x2;
			sx3 += x2*x;
			sx4 += x2*x2;

			b0 += y;
			b1 += x*y;
			b2 += x2*y;
		}

		DMatrix3x3 A = work;
		A.set(  sx0,sx1,sx2,
				sx1,sx2,sx3,
				sx2,sx3,sx4);

		if( !CommonOps_DDF3.invert(A,A) ) // TODO use a symmetric inverse. Should be slightly faster
			return false;

		// output = inv(A)*B      Unrolled here for speed
		output.a = A.a11*b0 + A.a12*b1 + A.a13*b2;
		output.b = A.a21*b0 + A.a22*b1 + A.a23*b2;
		output.c = A.a31*b0 + A.a32*b1 + A.a33*b2;

		return true;
	}
public static boolean fitMM(double[] data, int offset , int length ,
								PolynomialQuadratic1D_F64 output , @Nullable DMatrix3x3 work ) {
		if( work == null )
			work = new DMatrix3x3();

		final int N = length/2;

		// Unrolled pseudo inverse
		// coef = inv(A^T*A)*A^T*y
		double sx0=N,sx1=0,sx2=0;
		double             sx3=0;
		double             sx4=0;

		double b0=0,b1=0,b2=0;

		int end = offset+length;
		for (int i = offset; i < end; i += 2 ) {
			double x = data[i];
			double y = data[i+1];

			double x2 = x*x;

			sx1 += x;
			sx2 += x2;
			sx3 += x2*x;
			sx4 += x2*x2;

			b0 += y;
			b1 += x*y;
			b2 += x2*y;
		}

		DMatrix3x3 A = work;
		A.set(  sx0,sx1,sx2,
				sx1,sx2,sx3,
				sx2,sx3,sx4);

		if( !CommonOps_DDF3.invert(A,A) ) // TODO use a symmetric inverse. Should be slightly faster
			return false;

		// output = inv(A)*B      Unrolled here for speed
		output.a = A.a11*b0 + A.a12*b1 + A.a13*b2;
		output.b = A.a21*b0 + A.a22*b1 + A.a23*b2;
		output.c = A.a31*b0 + A.a32*b1 + A.a33*b2;

		return true;
	}
<p>Fits points to the quadratic polynomial. It's recommended that you apply a linear transform to the points
to ensure that they have zero mean and a standard deviation of 1. Then reverse the transform on the result.
A solution is found using the pseudo inverse and inverse by minor matrices.</p>

<p>For a more numerically stable algorithm see {@link #fitQRP(double[], int, int, PolynomialQuadratic1D_F64)}</p>

@param data Interleaved data [input[0], output[0], ....
@param offset first index in data
@param length number of elements in data that are to be read. must be divisible by 2
@param output (Optional) storage for the curve
@return The fitted curve