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 ) ; 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 ) ) return false ; 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