difference of squares factorization
542361
function(symbol, factors) { 
                try {
                    var remove_square = function(x) {
                        return core.Utils.block('POSITIVE_MULTIPLIERS', function() {
                            return Symbol.unwrapPARENS(math.sqrt(math.abs(x)));
                        }, true);
                    }; 
                    var separated = core.Utils.separate(symbol.clone());
                    var obj_array = [];
                    //get the unique variables
                    for(var x in separated) {
                        if(x !== 'constants') {
                            obj_array.push(separated[x]);
                        }
                    }
                    obj_array.sort(function(a, b) {
                        return b.power - a.power;
                    });

                    //if we have the same number of variables as unique variables then we can apply the difference of squares
                    if(obj_array.length === 2) { 
                        var a, b;
                        a = obj_array.pop();
                        b = obj_array.pop();
                        if(a.isComposite() && b.power.equals(2)) {
                            //remove the square from b
                            b = remove_square(b);
                            var f = __.Factor.factor(_.add(a, separated.constants));
                            if(f.power.equals(2)) {
                                f.toLinear();
                                factors.add(_.subtract(f.clone(), b.clone()));
                                factors.add(_.add(f, b));
                                symbol = new Symbol(1);
                            }
                        }
                        else {
                            a = a.powSimp();
                            b = b.powSimp();
                            
                            if((a.group === S || a.fname === '') && a.power.equals(2) && (b.group === S || b.fname === '') && b.power.equals(2)) {
                                if(a.multiplier.lessThan(0)) {
                                    var t = b; b = a; a = t;
                                }
                                if(a.multiplier.greaterThan(0)) {
                                    a = remove_square(a);
                                    b = remove_square(b);
                                }
                                factors.add(_.subtract(a.clone(), b.clone()));
                                factors.add(_.add(a, b));
                                symbol = new Symbol(1);
                            }
                        }
                    }
                }
                catch(e){;}
                
                return symbol;
            }
882706
def get_factors_iterative1(n):
    

    todo, res = [(n, 2, [])], []
    while todo:
        n, i, combi = todo.pop()
        while i * i <= n:
            if n % i == 0:
                res += combi + [i, n//i],
                todo.append((n//i, i, combi+[i])),
            i += 1
    return res
827941
public static LinearSolverDense<DMatrixRMaj> general(int numRows , int numCols ) {
        if( numRows == numCols )
            return linear(numRows);
        else
            return leastSquares(numRows,numCols);
    }
434552
@SuppressWarnings("unchecked")
    public static double cosineSimilarity(DoubleVector a, DoubleVector b) {
        double dotProduct = 0.0;
        double aMagnitude = a.magnitude();
        double bMagnitude = b.magnitude();

        // Check whether both vectors support fast iteration over their non-zero
        // values.  If so, use only the non-zero indices to speed up the
        // computation by avoiding zero multiplications
        if (a instanceof Iterable && b instanceof Iterable) {
            // Check whether we can easily determine how many non-zero values
            // are in each vector.  This value is used to select the iteration
            // order, which affects the number of get(value) calls.
            boolean useA =
                (a.length() < b.length() ||
                 (a instanceof SparseVector && b instanceof SparseVector) &&
                 ((SparseVector)a).getNonZeroIndices().length <
                 ((SparseVector)b).getNonZeroIndices().length);

            // Choose the smaller of the two to use in computing the dot
            // product.  Because it would be more expensive to compute the
            // intersection of the two sets, we assume that any potential
            // misses would be less of a performance hit.
            if (useA) {
                DoubleVector t = a;
                a = b;
                b = t;
            }

            for (DoubleEntry e : ((Iterable<DoubleEntry>)b)) {
                int index = e.index();                    
                double aValue = a.get(index);
                double bValue = e.value();
                dotProduct += aValue * bValue;
            }
        }

        // Check whether both vectors are sparse.  If so, use only the non-zero
        // indices to speed up the computation by avoiding zero multiplications
        else if (a instanceof SparseVector && b instanceof SparseVector) {
            SparseVector svA = (SparseVector)a;
            SparseVector svB = (SparseVector)b;
            int[] nzA = svA.getNonZeroIndices();
            int[] nzB = svB.getNonZeroIndices();

            // Choose the smaller of the two to use in computing the dot
            // product.  Because it would be more expensive to compute the
            // intersection of the two sets, we assume that any potential
            // misses would be less of a performance hit.
            if (a.length() < b.length() ||
                nzA.length < nzB.length) {
                DoubleVector t = a;
                a = b;
                b = t;
            }

            for (int nz : nzB) {
                double aValue = a.get(nz);
                double bValue = b.get(nz);
                dotProduct += aValue * bValue;
            }
        }

        // Check if the second vector is sparse.  If so, use only the non-zero
        // indices of b to speed up the computation by avoiding zero
        // multiplications.
        else if (b instanceof SparseVector) {
            SparseVector svB = (SparseVector)b;
            for (int nz : svB.getNonZeroIndices())
                dotProduct += b.get(nz) * a.get(nz);
        }

        // Check if the first vector is sparse.  If so, use only the non-zero
        // indices of a to speed up the computation by avoiding zero
        // multiplications.
        else if (a instanceof SparseVector) {
            SparseVector svA = (SparseVector)a;
            for (int nz : svA.getNonZeroIndices())
                dotProduct += b.get(nz) * a.get(nz);
        }

        // Otherwise, just assume both are dense and compute the full amount
        else {
            // Swap the vectors such that the b is the shorter vector and a is
            // the longer vector, or of equal length.   In the case that the two
            // vectors of unequal length, this will prevent any calls to out of
            // bounds values in the smaller vector.
            if (a.length() < b.length()) {
                DoubleVector t = a;
                a = b;
                b = t;
            }

            for (int i = 0; i < b.length(); i++) {
                double aValue = a.get(i);
                double bValue = b.get(i);
                dotProduct += aValue * bValue;
            }
        }
        return (aMagnitude == 0 || bMagnitude == 0)
            ? 0 : dotProduct / (aMagnitude * bMagnitude);
    }

    /**
     * 
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    @SuppressWarnings("unchecked")
    public static double cosineSimilarity(IntegerVector a, IntegerVector b) {
        check(a,b);

        int dotProduct = 0;
        double aMagnitude = a.magnitude();
        double bMagnitude = b.magnitude();

        // Check whether both vectors support fast iteration over their non-zero
        // values.  If so, use only the non-zero indices to speed up the
        // computation by avoiding zero multiplications
        if (a instanceof Iterable && b instanceof Iterable) {
            // Check whether we can easily determine how many non-zero values
            // are in each vector.  This value is used to select the iteration
            // order, which affects the number of get(value) calls.
            boolean useA =
                (a instanceof SparseVector && b instanceof SparseVector)
                && ((SparseVector)a).getNonZeroIndices().length <
                   ((SparseVector)b).getNonZeroIndices().length;

            // Choose the smaller of the two to use in computing the dot
            // product.  Because it would be more expensive to compute the
            // intersection of the two sets, we assume that any potential
            // misses would be less of a performance hit.
            if (useA) {
                IntegerVector t = a;
                a = b;
                b = t;
            }

            for (IntegerEntry e : ((Iterable<IntegerEntry>)b)) {
                int index = e.index();                    
                int aValue = a.get(index);
                int bValue = e.value();
                dotProduct += aValue * bValue;
            }            
        }

        // Check whether both vectors are sparse.  If so, use only the non-zero
        // indices to speed up the computation by avoiding zero multiplications
        else if (a instanceof SparseVector && b instanceof SparseVector) {
            SparseVector svA = (SparseVector)a;
            SparseVector svB = (SparseVector)b;
            int[] nzA = svA.getNonZeroIndices();
            int[] nzB = svB.getNonZeroIndices();
            // Choose the smaller of the two to use in computing the dot
            // product.  Because it would be more expensive to compute the
            // intersection of the two sets, we assume that any potential
            // misses would be less of a performance hit.
            if (nzA.length < nzB.length) {
                for (int nz : nzA) {
                    int aValue = a.get(nz);
                    int bValue = b.get(nz);
                    dotProduct += aValue * bValue;
                }
            }
            else {
                for (int nz : nzB) {
                    int aValue = a.get(nz);
                    int bValue = b.get(nz);
                    dotProduct += aValue * bValue;
                }
            }
        }

        // Otherwise, just assume both are dense and compute the full amount
        else {
            for (int i = 0; i < b.length(); i++) {
                int aValue = a.get(i);
                int bValue = b.get(i);
                dotProduct += aValue * bValue;
            }
        }
        return (aMagnitude == 0 || bMagnitude == 0)
            ? 0 : dotProduct / (aMagnitude * bMagnitude);
    }

    /**
     * 
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double cosineSimilarity(Vector a, Vector b) {
        return 
            (a instanceof IntegerVector && b instanceof IntegerVector)
            ? cosineSimilarity((IntegerVector)a, (IntegerVector)b)
            : cosineSimilarity(Vectors.asDouble(a), Vectors.asDouble(b));
    }

    /**
     * Returns the Pearson product-moment correlation coefficient of the two
     * arrays.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double correlation(double[] arr1, double[] arr2) {
        check(arr1, arr2);

        // REMINDER: this could be made more effecient by not looping
        double xSum = 0;
        double ySum = 0;
        for (int i = 0; i < arr1.length; ++i) {
            xSum += arr1[i];
            ySum += arr2[i];
        }
        
        double xMean = xSum / arr1.length;
        double yMean = ySum / arr1.length;
    
        double numerator = 0, xSqSum = 0, ySqSum = 0;
        for (int i = 0; i < arr1.length; ++i) {
            double x = arr1[i] - xMean;
            double y = arr2[i] - yMean;
            numerator += x * y;
            xSqSum += (x * x);
            ySqSum += (y * y);
        }
        if (xSqSum == 0 || ySqSum == 0)
            return 0;

        return numerator / Math.sqrt(xSqSum * ySqSum);
    }

    /**
     * Returns the Pearson product-moment correlation coefficient of the two
     * arrays.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double correlation(int[] arr1, int[] arr2) {
        check(arr1, arr2);

        // REMINDER: this could be made more effecient by not looping
        long xSum = 0;
        long ySum = 0;
        for (int i = 0; i < arr1.length; ++i) {
            xSum += arr1[i];
            ySum += arr2[i];
        }
        
        double xMean = xSum / (double)(arr1.length);
        double yMean = ySum / (double)(arr1.length);
    
        double numerator = 0, xSqSum = 0, ySqSum = 0;
        for (int i = 0; i < arr1.length; ++i) {
            double x = arr1[i] - xMean;
            double y = arr2[i] - yMean;
            numerator += x * y;
            xSqSum += (x * x);
            ySqSum += (y * y);
        }
        return numerator / Math.sqrt(xSqSum * ySqSum);
    }

    /**
     * Returns the Pearson product-moment correlation coefficient of the two
     * {@code Vector}s.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double correlation(DoubleVector arr1, DoubleVector arr2) {
        check(arr1, arr2);

        check(arr1, arr2);

        // REMINDER: this could be made more effecient by not looping
        double xSum = 0;
        double ySum = 0;
        for (int i = 0; i < arr1.length(); ++i) {
            xSum += arr1.get(i);
            ySum += arr2.get(i);
        }
        
        double xMean = xSum / arr1.length();
        double yMean = ySum / arr1.length();
    
        double numerator = 0, xSqSum = 0, ySqSum = 0;
        for (int i = 0; i < arr1.length(); ++i) {
            double x = arr1.get(i) - xMean;
            double y = arr2.get(i) - yMean;
            numerator += x * y;
            xSqSum += (x * x);
            ySqSum += (y * y);
        }
        return numerator / Math.sqrt(xSqSum * ySqSum);
    }

    /**
     * Returns the Pearson product-moment correlation coefficient of the two
     * {@code Vector}s.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double correlation(IntegerVector arr1, DoubleVector arr2) {
        check(arr1, arr2);

        // REMINDER: this could be made more effecient by not looping
        double xSum = 0;
        double ySum = 0;
        for (int i = 0; i < arr1.length(); ++i) {
            xSum += arr1.get(i);
            ySum += arr2.get(i);
        }
        
        double xMean = xSum / arr1.length();
        double yMean = ySum / arr1.length();
    
        double numerator = 0, xSqSum = 0, ySqSum = 0;
        for (int i = 0; i < arr1.length(); ++i) {
            double x = arr1.get(i) - xMean;
            double y = arr2.get(i) - yMean;
            numerator += x * y;
            xSqSum += (x * x);
            ySqSum += (y * y);
        }
        return numerator / Math.sqrt(xSqSum * ySqSum);
    }

    /**
     * Returns the Pearson product-moment correlation coefficient of the two
     * {@code Vector}s.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double correlation(Vector a, Vector b) {
        return correlation(Vectors.asDouble(a), Vectors.asDouble(b));
    }

    /**
     * Returns the euclidian distance between two arrays of {code double}s.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double euclideanDistance(double[] a, double[] b) {
        check(a, b);        
        double sum = 0;
        for (int i = 0; i < a.length; ++i)
            sum += Math.pow((a[i] - b[i]), 2);
        return Math.sqrt(sum);
    }

    /**
     * Returns the euclidian distance between two arrays of {code double}s.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double euclideanDistance(int[] a, int[] b) {
        check(a, b);
        
        long sum = 0;
        for (int i = 0; i < a.length; ++i)
            sum += Math.pow(a[i] - b[i], 2);
        return Math.sqrt(sum);
    }

    /**
     * Returns the euclidian distance between two {@code DoubleVector}s.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double euclideanDistance(DoubleVector a, DoubleVector b) {
        check(a, b);
        
        if (a instanceof SparseVector && b instanceof SparseVector) {
            SparseVector svA = (SparseVector)a;
            SparseVector svB = (SparseVector)b;

            int[] aNonZero = svA.getNonZeroIndices();
            int[] bNonZero = svB.getNonZeroIndices();
            TIntSet union = new TIntHashSet(aNonZero);
            union.addAll(bNonZero);
            
            double sum = 0;
            int[] nzIndices = union.toArray();
            for (int nz : nzIndices) {
                double x = a.get(nz);
                double y = b.get(nz);
                double diff = x - y;
                sum += diff * diff;
            }
            return Math.sqrt(sum);

        } else if (b instanceof SparseVector) {
            // If b is sparse, use a special case where we use the cached
            // magnitude of a and the sparsity of b to avoid most of the
            // computations.
            SparseVector sb = (SparseVector) b;
            int[] bNonZero = sb.getNonZeroIndices();
            double sum = 0;

            // Get the magnitude for a.  This value will often only be computed
            // once for the first vector once since the DenseVector caches the
            // magnitude, thus saving a large amount of computation.
            double aMagnitude = Math.pow(a.magnitude(), 2);

            // Compute the difference between the nonzero values of b and the
            // corresponding values for a.
            for (int index : bNonZero) {
                double value = a.get(index);
                // Decrement a's value at this index from it's magnitude.
                aMagnitude -= Math.pow(value, 2);
                sum += Math.pow(value - b.get(index), 2);
            }

            // Since the rest of b's values are 0, the difference between a and
            // b for these values is simply the magnitude of indices which have
            // not yet been traversed in a.  This corresponds to the modified
            // magnitude that was computed.
            sum += aMagnitude;

            return (sum < 0d) ? 0 : Math.sqrt(sum);
        }

        double sum = 0;
        for (int i = 0; i < a.length(); ++i)
            sum += Math.pow((a.get(i) - b.get(i)), 2);
        return Math.sqrt(sum);
    }

    /**
     * Returns the euclidian distance between two {@code DoubleVector}s.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double euclideanDistance(IntegerVector a, IntegerVector b) {
        check(a, b);
         
        if (a instanceof SparseVector && b instanceof SparseVector) {
            SparseVector svA = (SparseVector)a;
            SparseVector svB = (SparseVector)b;

            int[] aNonZero = svA.getNonZeroIndices();
            int[] bNonZero = svB.getNonZeroIndices();
            HashSet<Integer> sparseIndicesA = new HashSet<Integer>(
                    aNonZero.length);
            double sum = 0;
            for (int nonZero : aNonZero) {
                sum += Math.pow((a.get(nonZero) - b.get(nonZero)), 2);
                sparseIndicesA.add(nonZero);
            }

            for (int nonZero : bNonZero)
                if (!sparseIndicesA.contains(bNonZero))
                    sum += Math.pow(b.get(nonZero), 2);
            return sum;
        }

        double sum = 0;
        for (int i = 0; i < a.length(); ++i)
            sum += Math.pow((a.get(i) - b.get(i)), 2);
        return Math.sqrt(sum);
    }

    /**
     * Returns the euclidian distance between two {@code Vector}s.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double euclideanDistance(Vector a, Vector b) {
        return euclideanDistance(Vectors.asDouble(a), Vectors.asDouble(b));
    }

    /**
     * Returns the euclidian similiarty between two arrays of values.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double euclideanSimilarity(int[] a, int[] b) {
        return 1 / (1 + euclideanDistance(a,b));
    }

    /**
     * Returns the euclidian similiarty between two arrays of values.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double euclideanSimilarity(double[] a, double[] b) {
        return 1 / (1 + euclideanDistance(a,b));
    }

    /**
     * Returns the euclidian similiarty between two arrays of values.
     *
     * @throws IllegalArgumentException when the length of the two vectors are
     *                                   not the same.
     */
    public static double euclideanSimilarity(Vector a, Vector b) {
        return 1 / (1 + euclideanDistance(a,b));
    }

    /**
     * Computes the <a href="http://en.wikipedia.org/wiki/Jaccard_index">Jaccard
     * index</a> of the two sets of elements.
     */
    public static double jaccardIndex(Set<?> a, Set<?> b) {        
        int intersection = 0;
        for (Object o : a) {
            if (b.contains(o))
                intersection++;
        }

        double union = a.size() + b.size() - intersection;
        return intersection / union;
    }

    /**
     * Computes the <a href="http://en.wikipedia.org/wiki/Jaccard_index">Jaccard
     * index</a> comparing the similarity both arrays when viewed as sets of
     * samples.
     */
    public static double jaccardIndex(double[] a, double[] b) {        
        Set<Double> intersection = new HashSet<Double>();
        Set<Double> union = new HashSet<Double>();
        for (double d : a) {
            intersection.add(d);
            union.add(d);
        }
        Set<Double> tmp = new HashSet<Double>();
        for (double d : b) {
            tmp.add(d);
            union.add(d);
        }

        intersection.retainAll(tmp);
        return ((double)(intersection.size())) / union.size();
    }

    /**
     * Computes the <a href="http://en.wikipedia.org/wiki/Jaccard_index">Jaccard
     * index</a> comparing the similarity both arrays when viewed as sets of
     * samples.
     */
    public static double jaccardIndex(int[] a, int[] b) {
        // The BitSets should be faster than a HashMap since it's back by an
        // array and operations are just logical bit operations and require no
        // auto-boxing.  However, if a or b contains large values, then the cost
        // of creating the necessary size for the BitSet may outweigh its
        // performance.  At some point, it would be useful to profile the two
        // methods and their associated worst cases. -jurgens
        BitSet c = new BitSet();
        BitSet d = new BitSet();
        BitSet union = new BitSet();
        for (int i : a) {
            c.set(i);
            union.set(i);
        }
        for (int i : b) {
            d.set(i);
            union.set(i);
        }
        
        // get the intersection
        c.and(d); 
        return ((double)(c.cardinality())) / union.cardinality();
    }
207420
def pairwise_euclid_distance(A, B):
    
    batchA = tf.shape(A)[0]
    batchB = tf.shape(B)[0]

    sqr_norm_A = tf.reshape(tf.reduce_sum(tf.pow(A, 2), 1), [1, batchA])
    sqr_norm_B = tf.reshape(tf.reduce_sum(tf.pow(B, 2), 1), [batchB, 1])
    inner_prod = tf.matmul(B, A, transpose_b=True)

    tile_1 = tf.tile(sqr_norm_A, [batchB, 1])
    tile_2 = tf.tile(sqr_norm_B, [1, batchA])
    return (tile_1 + tile_2 - 2 * inner_prod)
243161
def solve_factorized_aug(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """"""
    m = M    

    """Primal variable"""
    x = z[0:N]

    """Multiplier for equality constraints"""
    nu = z[N:N+P]

    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """RHS for condensed system"""
    b1 = -rd - mydot(G.T, mydot(SIG, rp2)) + mydot(G.T, rc/s)
    b2 = -rp1
    b = np.hstack((b1, b2))
    dxnu = mysolve(LU, b)
    dx = dxnu[0:N]
    dnu = dxnu[N:]

    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
242601
def _compute_projection_filters(G, sf, estimated_source):
    
    # epsilon
    eps = np.finfo(np.float).eps

    # shapes
    (nsampl, nchan) = estimated_source.shape
    # handles the case where we are calling this with only one source
    # G should be nsrc X nsrc X nchan X nchan X filters_len X filters_len
    # and sf should be nsrc X nchan X filters_len
    if len(G.shape) == 4:
        G = G[None, None, ...]
        sf = sf[None, ...]
    nsrc = G.shape[0]
    filters_len = G.shape[-1]

    # zero pad estimates and put chan in first dimension
    estimated_source = _zeropad(estimated_source.T, filters_len - 1, axis=1)

    # compute its FFT
    n_fft = int(2**np.ceil(np.log2(nsampl + filters_len - 1.)))
    sef = scipy.fftpack.fft(estimated_source, n=n_fft)

    # compute the cross-correlations between sources and estimates
    D = np.zeros((nsrc, nchan, filters_len, nchan))
    for (j, cj, c) in itertools.product(
        list(range(nsrc)), list(range(nchan)), list(range(nchan))
    ):
        ssef = sf[j, cj] * np.conj(sef[c])
        ssef = np.real(scipy.fftpack.ifft(ssef))
        D[j, cj, :, c] = np.hstack((ssef[0], ssef[-1:-filters_len:-1]))

    # reshape matrices to build the filters
    D = D.reshape(nsrc * nchan * filters_len, nchan)
    G = _reshape_G(G)

    # Distortion filters
    try:
        C = np.linalg.solve(G + eps*np.eye(G.shape[0]), D).reshape(
            nsrc, nchan, filters_len, nchan
        )
    except np.linalg.linalg.LinAlgError:
        C = np.linalg.lstsq(G, D)[0].reshape(
            nsrc, nchan, filters_len, nchan
        )

    # if we asked for one single reference source,
    # return just a nchan X filters_len matrix
    if nsrc == 1:
        C = C[0]
    return C
750785
def diff_non_uni(self, y, coords, dim, coefs):
        """"""

        yd = np.zeros_like(y)

        ndims = len(y.shape)
        multi_slice = [slice(None, None)] * ndims
        ref_multi_slice = [slice(None, None)] * ndims

        for i, x in enumerate(coords):
            weights = coefs[i]["coefficients"]
            offsets = coefs[i]["offsets"]
            ref_multi_slice[dim] = i

            for off, w in zip(offsets, weights):
                multi_slice[dim] = i + off
                yd[ref_multi_slice] += w * y[multi_slice]

        return yd
650218
def solve_factorized_schur(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """"""
    m = M    

    """Primal variable"""
    x = z[0:N]

    """Multiplier for equality constraints"""
    nu = z[N:N+P]

    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """Assemble right hand side of augmented system"""
    r1 = rd + mydot(G.T, mydot(SIG, rp2)) - mydot(G.T, rc/s)
    r2 = rp1

    """Unpack LU-factors"""
    LU_S, LU_H = LU

    """Assemble right hand side for normal equation"""
    b = r2 - mydot(A, mysolve(LU_H, r1))  

    """Solve for dnu"""
    dnu = mysolve(LU_S, b)
       
    """Solve for dx"""
    dx = mysolve(LU_H, -(r1 + mydot(A.T, dnu)))    
    
    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
1049602
def factorize_modp(an = 0)
      f = self
      f = f.sqfree_zp if an == 0
      an += 1 while f.deg_at(an) <= 0
      f = f.pp(an)
      lc_of_f = f.lc_at(an)
      f = f.monic_int(an)
      facts = f._factorize(an)
      facts = facts.map! do |f2|
        f3 = f2.project(self.class) { |c, j| c * lc_of_f**j[an] }
        f3.pp(an)
      end
      facts = facts.fact_all_u(self)
      facts.mcorrect_lc!(self, an) do |fac|
        fac.factorize_modp(an + 1)
      end
      facts
    end
1327131
private void _symmetricTred2 ()
  {
    // This is derived from the Algol procedures tred2 by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.
    for (int j = 0; j < m_nDim; j++)
    {
      m_aEVd[j] = m_aEigenVector[m_nDim - 1][j];
    }

    // Householder reduction to tridiagonal form.

    for (int i = m_nDim - 1; i > 0; i--)
    {
      // Scale to avoid under/overflow.
      double scale = 0.0;
      double h = 0.0;
      for (int k = 0; k < i; k++)
        scale += MathHelper.abs (m_aEVd[k]);
      if (scale == 0.0)
      {
        m_aEVe[i] = m_aEVd[i - 1];
        for (int j = 0; j < i; j++)
        {
          m_aEVd[j] = m_aEigenVector[i - 1][j];
          m_aEigenVector[i][j] = 0.0;
          m_aEigenVector[j][i] = 0.0;
        }
      }
      else
      {
        // Generate Householder vector.

        for (int k = 0; k < i; k++)
        {
          m_aEVd[k] /= scale;
          h += m_aEVd[k] * m_aEVd[k];
        }
        double f = m_aEVd[i - 1];
        double g = Math.sqrt (h);
        if (f > 0)
        {
          g = -g;
        }
        m_aEVe[i] = scale * g;
        h -= f * g;
        m_aEVd[i - 1] = f - g;
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] = 0.0;
        }

        // Apply similarity transformation to remaining columns.
        for (int j = 0; j < i; j++)
        {
          f = m_aEVd[j];
          m_aEigenVector[j][i] = f;
          g = m_aEVe[j] + m_aEigenVector[j][j] * f;
          for (int k = j + 1; k <= i - 1; k++)
          {
            g += m_aEigenVector[k][j] * m_aEVd[k];
            m_aEVe[k] += m_aEigenVector[k][j] * f;
          }
          m_aEVe[j] = g;
        }
        f = 0.0;
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] /= h;
          f += m_aEVe[j] * m_aEVd[j];
        }
        final double hh = f / (h + h);
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] -= hh * m_aEVd[j];
        }
        for (int j = 0; j < i; j++)
        {
          f = m_aEVd[j];
          g = m_aEVe[j];
          for (int k = j; k <= i - 1; k++)
          {
            m_aEigenVector[k][j] -= (f * m_aEVe[k] + g * m_aEVd[k]);
          }
          m_aEVd[j] = m_aEigenVector[i - 1][j];
          m_aEigenVector[i][j] = 0.0;
        }
      }
      m_aEVd[i] = h;
    }

    // Accumulate transformations.
    for (int i = 0; i < m_nDim - 1; i++)
    {
      m_aEigenVector[m_nDim - 1][i] = m_aEigenVector[i][i];
      m_aEigenVector[i][i] = 1.0;
      final double h = m_aEVd[i + 1];
      if (h != 0.0)
      {
        for (int k = 0; k <= i; k++)
        {
          m_aEVd[k] = m_aEigenVector[k][i + 1] / h;
        }
        for (int j = 0; j <= i; j++)
        {
          double g = 0.0;
          for (int k = 0; k <= i; k++)
          {
            g += m_aEigenVector[k][i + 1] * m_aEigenVector[k][j];
          }
          for (int k = 0; k <= i; k++)
          {
            m_aEigenVector[k][j] -= g * m_aEVd[k];
          }
        }
      }
      for (int k = 0; k <= i; k++)
      {
        m_aEigenVector[k][i + 1] = 0.0;
      }
    }
    for (int j = 0; j < m_nDim; j++)
    {
      m_aEVd[j] = m_aEigenVector[m_nDim - 1][j];
      m_aEigenVector[m_nDim - 1][j] = 0.0;
    }
    m_aEigenVector[m_nDim - 1][m_nDim - 1] = 1.0;
    m_aEVe[0] = 0.0;
  }
687027
function(poly1, poly2)
{
    var mini1 = findMSA(poly2, poly1.tPlanes);
    if(mini1 == -1) return NONE;
    var min1 = last_MSA_min;
    
    var mini2 = findMSA(poly1, poly2.tPlanes);
    if(mini2 == -1) return NONE;
    var min2 = last_MSA_min;
    
    // There is overlap, find the penetrating verts
    if(min1 > min2)
        return findVerts(poly1, poly2, poly1.tPlanes[mini1].n, min1);
    else
        return findVerts(poly1, poly2, vneg(poly2.tPlanes[mini2].n), min2);
}
1136279
def sq_dist(X1,X2=None):
    
    if X2==None:
        X2 = X1
    assert X1.shape[1]==X2.shape[1], 'dimensions do not match'

    n = X1.shape[0]
    m = X2.shape[0]
    d = X1.shape[1]
    # (X1 - X2)**2 = X1**2 + X2**2 - 2X1X2
    X1sq = sp.reshape((X1**2).sum(1),n,1)
    X2sq = sp.reshape((X2**2).sum(1),m,1)

    K = sp.tile((X1*X1).sum(1),(m,1)).T + sp.tile((X2*X2).sum(1),(n,1)) - 2*sp.dot(X1,X2.T)
    return K
614211
def _D_constraint(self, neg_pairs, w):
    
    diff = neg_pairs[:, 0, :] - neg_pairs[:, 1, :]
    diff_sq = diff * diff
    dist = np.sqrt(diff_sq.dot(w))
    sum_deri1 = np.einsum('ij,i', diff_sq, 0.5 / np.maximum(dist, 1e-6))
    sum_deri2 = np.einsum(
        'ij,ik->jk',
        diff_sq,
        diff_sq / (-4 * np.maximum(1e-6, dist**3))[:,None]
    )
    sum_dist = dist.sum()
    return (
      np.log(sum_dist),
      sum_deri1 / sum_dist,
      sum_deri2 / sum_dist - np.outer(sum_deri1, sum_deri1) / (sum_dist * sum_dist)
    )
23777
function() { 
            var a = this.clone(),
                i = 1,
                b = a.clone().diff(),
                c = a.clone().gcd(b),
                w = a.divide(c)[0];
            var output = Polynomial.fromArray([new Frac(1)], a.variable);
            while(!c.equalsNumber(1)) { 
                var y = w.gcd(c); 
                var z = w.divide(y)[0];
                //one of the factors may have shown up since it's square but smaller than the 
                //one where finding
                if(!z.equalsNumber(1) && i>1) {
                    var t = z.clone();
                    for(var j=1; j<i; j++)
                        t.multiply(z.clone());
                    z = t;
                }
                output = output.multiply(z); 
                i++;
                w = y;
                c = c.divide(y)[0];
            }
            return [output, w, i];
        }
569946
function factorize(n, addFactor) {
    let remainder = 0;
    let i = 0;
    let multiplier = FOUND_PRIMES[i];

    //go through prime multipliers
    while (multiplier <= n) {
        do {
            remainder = n % multiplier;
            if (!remainder) {
                n /= multiplier;
                addFactor(multiplier);
            }
        } while (!remainder);
        multiplier = FOUND_PRIMES[++i];
    }
}
776373
@Override
	public void computeHessian(DMatrixRMaj jacLeft , DMatrixRMaj jacRight) {
		A.reshape(jacLeft.numCols,jacLeft.numCols);
		B.reshape(jacLeft.numCols,jacRight.numCols);
		D.reshape(jacRight.numCols,jacRight.numCols);

		// take advantage of the inner product's symmetry when possible to reduce
		// the number of calculations
		MatrixMultProduct_DDRM.inner_reorder_lower(jacLeft,A);
		CommonOps_DDRM.symmLowerToFull(A);
		CommonOps_DDRM.multTransA(jacLeft,jacRight,B);
		MatrixMultProduct_DDRM.inner_reorder_lower(jacRight,D);
		CommonOps_DDRM.symmLowerToFull(D);
	}