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

		final int N = length/2;
		if( N < 4 )
			throw new IllegalArgumentException("Need at least 4 points and not "+N);

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

		double b0=0,b1=0,b2=0,b3=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;
			double x3 = x2*x;
			double x4 = x2*x2;

			sx1 += x; sx2 += x2; sx3 += x3;
			sx4 += x4; sx5 += x4*x;
			sx6 += x4*x2;

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

		A.set(  N,  sx1,sx2,sx3,
				sx1,sx2,sx3,sx4,
				sx2,sx3,sx4,sx5,
				sx3,sx4,sx5,sx6);

		if( !CommonOps_DDF4.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 + A.a14*b3;
		output.b = A.a21*b0 + A.a22*b1 + A.a23*b2 + A.a24*b3;
		output.c = A.a31*b0 + A.a32*b1 + A.a33*b2 + A.a34*b3;
		output.d = A.a41*b0 + A.a42*b1 + A.a43*b2 + A.a44*b3;

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

		final int N = length/2;
		if( N < 4 )
			throw new IllegalArgumentException("Need at least 4 points and not "+N);

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

		double b0=0,b1=0,b2=0,b3=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;
			double x3 = x2*x;
			double x4 = x2*x2;

			sx1 += x; sx2 += x2; sx3 += x3;
			sx4 += x4; sx5 += x4*x;
			sx6 += x4*x2;

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

		A.set(  N,  sx1,sx2,sx3,
				sx1,sx2,sx3,sx4,
				sx2,sx3,sx4,sx5,
				sx3,sx4,sx5,sx6);

		if( !CommonOps_DDF4.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 + A.a14*b3;
		output.b = A.a21*b0 + A.a22*b1 + A.a23*b2 + A.a24*b3;
		output.c = A.a31*b0 + A.a32*b1 + A.a33*b2 + A.a34*b3;
		output.d = A.a41*b0 + A.a42*b1 + A.a43*b2 + A.a44*b3;

		return true;
	}
<p>Fits points to the cubic 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