Let A[i, j] be zero based on Bezout's identity

           [ii ij]
           [ji jj] is a (k,k) minor of original 3x3 matrix.
766551
def _set_zero(self, i, j, a, b, r, s, t):
        

        L = np.eye(3, dtype='intc')
        L[i, i] = s
        L[i, j] = t
        L[j, i] = -b // r
        L[j, j] = a // r
        self._L.append(L.copy())
        self._A = np.dot(L, self._A)
66810
def apply_M(self, ax, ay):
        
        jac = numpy.array(
            [[self.dx.dot(ax), self.dy.dot(ax)], [self.dx.dot(ay), self.dy.dot(ay)]]
        )

        # jacs and J are of shape (2, 2, k). M must be of the same shape and
        # contain the result of the k 2x2 dot products. Perhaps there's a
        # dot() for this.
        M = numpy.einsum("ijl,jkl->ikl", jac, self.J)
        # M = numpy.array([
        #     [
        #         jac[0][0]*self.J[0][0] + jac[0][1]*self.J[1][0],
        #         jac[0][0]*self.J[0][1] + jac[0][1]*self.J[1][1],
        #     ],
        #     [
        #         jac[1][0]*self.J[0][0] + jac[1][1]*self.J[1][0],
        #         jac[1][0]*self.J[0][1] + jac[1][1]*self.J[1][1],
        #     ],
        #     ])

        # One could use
        #
        #     M = numpy.moveaxis(M, -1, 0)
        #     _, sigma, _ = numpy.linalg.svd(M)
        #
        # but computing the singular values explicitly via
        # <https://scicomp.stackexchange.com/a/14103/3980> is faster and more
        # explicit.
        a = (M[0, 0] + M[1, 1]) / 2
        b = (M[0, 0] - M[1, 1]) / 2
        c = (M[1, 0] + M[0, 1]) / 2
        d = (M[1, 0] - M[0, 1]) / 2

        return a, b, c, d
1061410
public static DMatrixRMaj createFundamental(DMatrixRMaj E, DMatrixRMaj K) {
		DMatrixRMaj K_inv = new DMatrixRMaj(3,3);
		CommonOps_DDRM.invert(K,K_inv);

		DMatrixRMaj F = new DMatrixRMaj(3,3);
		PerspectiveOps.multTranA(K_inv,E,K_inv,F);

		return F;
	}
506915
public static void diag( DMatrix3x3 input , DMatrix3 out ) {
        out.a1 = input.a11;
        out.a2 = input.a22;
        out.a3 = input.a33;
    }
1029866
public Matrix getMatrix(int i0, int i1, int j0, int j1)
    {
        Matrix X = new Matrix(i1 - i0 + 1, j1 - j0 + 1);
        double[][] B = X.getArray();
        try
        {
            for (int i = i0; i <= i1; i++)
            {
                for (int j = j0; j <= j1; j++)
                {
                    B[i - i0][j - j0] = A[i][j];
                }
            }
        }
        catch (ArrayIndexOutOfBoundsException e)
        {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
        return X;
    }
333305
public static double elementMinAbs( DMatrix3x3 a ) {
        double min = Math.abs(a.a11);
        double tmp = Math.abs(a.a12); if( tmp < min ) min = tmp;
        tmp = Math.abs(a.a13); if( tmp < min ) min = tmp;
        tmp = Math.abs(a.a21); if( tmp < min ) min = tmp;
        tmp = Math.abs(a.a22); if( tmp < min ) min = tmp;
        tmp = Math.abs(a.a23); if( tmp < min ) min = tmp;
        tmp = Math.abs(a.a31); if( tmp < min ) min = tmp;
        tmp = Math.abs(a.a32); if( tmp < min ) min = tmp;
        tmp = Math.abs(a.a33); if( tmp < min ) min = tmp;

        return min;
    }
887365
public void set(int i, int j, Matrix<T> B)
    {
        int m = B.rows();
        int n = B.columns();
        for (int ii = 0; ii < m; ii++)
        {
            for (int jj = 0; jj < n; jj++)
            {
                set(i+ii, j+jj, B.get(ii, jj));
            }
        }
    }
830542
public Matrix4f fma4x3(Matrix4fc other, float otherFactor, Matrix4f dest) {
        dest._m00(m00 + other.m00() * otherFactor);
        dest._m01(m01 + other.m01() * otherFactor);
        dest._m02(m02 + other.m02() * otherFactor);
        dest._m03(m03);
        dest._m10(m10 + other.m10() * otherFactor);
        dest._m11(m11 + other.m11() * otherFactor);
        dest._m12(m12 + other.m12() * otherFactor);
        dest._m13(m13);
        dest._m20(m20 + other.m20() * otherFactor);
        dest._m21(m21 + other.m21() * otherFactor);
        dest._m22(m22 + other.m22() * otherFactor);
        dest._m23(m23);
        dest._m30(m30 + other.m30() * otherFactor);
        dest._m31(m31 + other.m31() * otherFactor);
        dest._m32(m32 + other.m32() * otherFactor);
        dest._m33(m33);
        dest._properties(0);
        return dest;
    }
891890
public static void setIdentity( DMatrix4x4 a ) {
        a.a11 = 1; a.a21 = 0; a.a31 = 0; a.a41 = 0;
        a.a12 = 0; a.a22 = 1; a.a32 = 0; a.a42 = 0;
        a.a13 = 0; a.a23 = 0; a.a33 = 1; a.a43 = 0;
        a.a14 = 0; a.a24 = 0; a.a34 = 0; a.a44 = 1;
    }
160651
private int trMedian3(int ISAd, int v1, int v2, int v3) {
        if (SA[ISAd + SA[v1]] > SA[ISAd + SA[v2]]) {
            int t = v1;
            v1 = v2;
            v2 = t;
        }
        if (SA[ISAd + SA[v2]] > SA[ISAd + SA[v3]]) {
            if (SA[ISAd + SA[v1]] > SA[ISAd + SA[v3]]) {
                return v1;
            } else {
                return v3;
            }
        }
        return v2;
    }
222112
public E determinant(E[][] matrix, int r, int s, int i, int j)
    {
        E a, b, c, d;
        a = matrix[r][s];
        b = matrix[r][j];
        c = matrix[i][s];
        d = matrix[i][j];
        E ad  = multiply(a, d); // matrix[r][s], matrix[i][j]);
        E bc = multiply(b, c); // matrix[i][s], matrix[r][j]);
        E ret = subtract(ad, bc);

        return ret;
    }
201410
public static void multAddOuter( double alpha , DMatrix4x4 A , double beta , DMatrix4 u , DMatrix4 v , DMatrix4x4 C ) {
        C.a11 = alpha*A.a11 + beta*u.a1*v.a1;
        C.a12 = alpha*A.a12 + beta*u.a1*v.a2;
        C.a13 = alpha*A.a13 + beta*u.a1*v.a3;
        C.a14 = alpha*A.a14 + beta*u.a1*v.a4;
        C.a21 = alpha*A.a21 + beta*u.a2*v.a1;
        C.a22 = alpha*A.a22 + beta*u.a2*v.a2;
        C.a23 = alpha*A.a23 + beta*u.a2*v.a3;
        C.a24 = alpha*A.a24 + beta*u.a2*v.a4;
        C.a31 = alpha*A.a31 + beta*u.a3*v.a1;
        C.a32 = alpha*A.a32 + beta*u.a3*v.a2;
        C.a33 = alpha*A.a33 + beta*u.a3*v.a3;
        C.a34 = alpha*A.a34 + beta*u.a3*v.a4;
        C.a41 = alpha*A.a41 + beta*u.a4*v.a1;
        C.a42 = alpha*A.a42 + beta*u.a4*v.a2;
        C.a43 = alpha*A.a43 + beta*u.a4*v.a3;
        C.a44 = alpha*A.a44 + beta*u.a4*v.a4;
    }
622077
function vec3lerp(out, a, b, t, oa, ob) {
    var ax = a[oa];
    var ay = a[oa + 1];
    var az = a[oa + 2];
    out[0] = ax + t * (b[ob] - ax);
    out[1] = ay + t * (b[ob + 1] - ay);
    out[2] = az + t * (b[ob + 2] - az);

    return out;
}
231345
public static DMatrixRMaj checkZerosLT(DMatrixRMaj A , int numRows , int numCols) {
        if( A == null ) {
            return new DMatrixRMaj(numRows,numCols);
        }  else if( numRows != A.numRows || numCols != A.numCols ) {
            A.reshape(numRows,numCols);
            A.zero();
        } else {
            for( int i = 0; i < A.numRows; i++ ) {
                int index = i*A.numCols;
                int end = index + Math.min(i,A.numCols);;
                while( index < end ) {
                    A.data[index++] = 0;
                }
            }
        }
        return A;
    }
1002155
def lowerlim(x, y, z, a, b, c):
    
    if x/a + y/b + z/c > 1:
        B = a + b + c - x - y - z
        C = a*b + a*c + b*c - a*y - a*z - b*x - b*z - c*x - c*y
        D = a*b*c - a*b*z - a*c*y - b*c*x
        r = np.roots([1,B,C,D])
        ll = r[~np.iscomplex(r) & (r>0.)]
        return ll[0].real
    else:
        return 0.
285202
def ij_jlk_to_ilk(A, B):
    
    return A.dot(B.reshape(B.shape[0], -1)).reshape(A.shape[0], B.shape[1], B.shape[2])
951934
def s2a(s):
    
    a = np.zeros((3, 3,), 'f')  # make the a matrix
    for i in range(3):
        a[i][i] = s[i]
    a[0][1], a[1][0] = s[3], s[3]
    a[1][2], a[2][1] = s[4], s[4]
    a[0][2], a[2][0] = s[5], s[5]
    return a