Finds the roots of a 3 dim polymon of the form x^3 + a1 * x^2 + a2 * x + a3.
    The roots are returned as complex numbers.
424111
def roots_of_cubic_polynom(a1, a2, a3):
    '''
    
    '''
    q = (a1 * a1 - 3.0 * a2) / 9.0
    r = (2 * a1 * a1 * a1 - 9.0 * a1 * a2 + 27.0 * a3) / 54.0
    r2 = r * r
    q3 = q * q * q
    a1d3 = a1 / 3.0

    if r2 - q3 >= 0.0:  # In this case there are 2 complex roots
        # Let a = - sgn(R) * ( |R| + sqrt(R^2 -Q^3) )^(1/3)
        oneThird = 1.0 / 3.0
        a = - sign(r) * (abs(r) + sqrt(r2 - q3)) ** oneThird
        b = q / a if a != 0.0 else 0.0
        apb = a + b
        root1 = complex(apb - a1d3)
        root2 = -0.5 * apb - a1d3 + sqrt(3) / 2.0 * (a1 - a2) * 1j
        root3 = root2.conjugate()
        return root1, root2, root3
    else:  # In this case there are three real roots
        theta = acos(r / sqrt(q3))
        fac = -2.0 * sqrt(q)
        root1 = complex(fac * cos(theta / 3.0) - a1d3)
        root2 = complex(fac * cos((theta + 2.0 * PI) / 3.0) - a1d3)
        root3 = complex(fac * cos((theta - 2.0 * PI) / 3.0) - a1d3)
        return root1, root2, root3
387865
def findXAt(xArr, yArr, yVal, index=0, s=0.0):
    
    if xArr[1] < xArr[0]:
        # numbers must be in ascending order, otherwise method crashes...
        xArr = xArr[::-1]
        yArr = yArr[::-1]

    yArr = yArr - yVal
    if len(yArr) < 5:
        xn = np.linspace(xArr[0], xArr[-1], 5)
        yArr = np.interp(xn, xArr, yArr)
        xArr = xn
    f = interpolate.UnivariateSpline(xArr, yArr, s=s)
    return f.roots()[index]
117216
public static void multAddOuter( double alpha , DMatrix3x3 A , double beta , DMatrix3 u , DMatrix3 v , DMatrix3x3 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.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.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;
    }
201437
function() {
        var derives     = this.derives,
            x           = this.x,
            dimension   = this.dimension,
            h           = this.h

        var i, _y = []

        // Alias: f() <=> this.derives()
        //        Xn  <=> this.x
        //        Yn  <=> this.y
        //        H   <=> this.h
        //        K1  <=> this._k1

        // calculate _k1: K1 = f(Xn, Yn)
        this._k1 = derives(x, this.y)

        // calculate _k2: K2 = f(Xn + 0.5 * H, Yn + 0.5 * H * K1)
        for (i = 0; i < dimension; i++) {
            _y[i] = this.y[i] + h * 0.5 * this._k1[i]
        }
        this._k2 = derives(x + h * 0.5, _y)

        // calculate _k3: K3 = f(Xn + 0.5 * H, Yn + 0.5 * H * K2)
        for (i = 0; i < dimension; i++) {
            _y[i] = this.y[i] + h * 0.5 * this._k2[i]
        }
        this._k3 = derives(x + h * 0.5, _y)

        // calculate _k4: K4 = f(Xn + H, Yn + H * K3)
        for (i = 0; i < dimension; i++) {
            _y[i] = this.y[i] + h * this._k3[i]
        }
        this._k4 = derives(x + h, _y)

        // calculate yNext: Y(n + 1) = Yn + H / 6 * (K1 + 2 * K2 + 2 * K3 + K4)
        for (i = 0; i < dimension; i++) {
            this.y[i] += h / 6 * (this._k1[i] + 2 * this._k2[i] + 2 * this._k3[i] + this._k4[i])
        }

        // calculate xNext: X(n + 1) = Xn + H
        this.x += h

        return this.y
    }
502051
function() {

      // acos(c) = i * log(c - i * sqrt(1 - c^2))

      var a = this['re'];
      var b = this['im'];

      var t1 = new Complex(
              b * b - a * a + 1,
              -2 * a * b)['sqrt']();

      var t2 = new Complex(
              t1['re'] - b,
              t1['im'] + a)['log']();

      return new Complex(Math.PI / 2 - t2['im'], t2['re']);
    }
29153
def gf_poly_mul(p, q):
    ''''''
    # Pre-allocate the result array
    r = bytearray(len(p) + len(q) - 1)
    # Precompute the logarithm of p
    lp = [gf_log[p[i]] for i in xrange(len(p))]
    # Compute the polynomial multiplication (just like the outer product of two vectors, we multiply each coefficients of p with all coefficients of q)
    for j in xrange(len(q)):
        qj = q[j] # optimization: load the coefficient once
        if qj != 0: # log(0) is undefined, we need to check that
            lq = gf_log[qj] # Optimization: precache the logarithm of the current coefficient of q
            for i in xrange(len(p)):
                if p[i] != 0: # log(0) is undefined, need to check that...
                    r[i + j] ^= gf_exp[lp[i] + lq] # equivalent to: r[i + j] = gf_add(r[i+j], gf_mul(p[i], q[j]))
    return r
775333
@Override
    public void solve(ZMatrixRMaj B, ZMatrixRMaj X) {
        if( X.numRows != numCols )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        else if( B.numRows != numRows || B.numCols != X.numCols )
            throw new IllegalArgumentException("Unexpected dimensions for B");

        int BnumCols = B.numCols;
        
        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                int indexB = (i*BnumCols + colB)*2;
                a.data[i*2]   = B.data[indexB];
                a.data[i*2+1] = B.data[indexB+1];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^T)*b = b - u*(gamma*U^T*b)
            for( int n = 0; n < numCols; n++ ) {
                double []u = QR[n];

                double realVV = u[n*2];
                double imagVV = u[n*2+1];

                u[n*2]   = 1;
                u[n*2+1] = 0;

                QrHelperFunctions_ZDRM.rank1UpdateMultR(a, u,0, gammas[n], 0, n, numRows, temp.data);

                u[n*2]   = realVV;
                u[n*2+1] = imagVV;
            }

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_ZDRM.solveU(R.data, a.data, numCols);

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                int indexB = (i*BnumCols + colB)*2;
                X.data[indexB]   = a.data[i*2];
                X.data[indexB+1] = a.data[i*2+1];
            }
        }
    }
443449
def pbis(a):
    """"""
    return(math.cos(3*a - math.pi), (math.sin(3*a - math.pi)))
690707
def Ax(self):
        
        #  0         -self.a3   self.a2
        #  self.a3    0        -self.a1
        # -self.a2    self.a1   0        
        m = np.zeros((len(self.a1), 3, 3))
        m[:, 0, 1] = -self.a3
        m[:, 0, 2] = +self.a2
        m[:, 1, 0] = +self.a3
        m[:, 1, 2] = -self.a1
        m[:, 2, 0] = -self.a2
        m[:, 2, 1] = +self.a1        
        return m
498810
public static void root(ComplexPolar_F64 a , int N , int k , ComplexPolar_F64 result )
    {
        result.r = Math.pow(a.r,1.0/N);
        result.theta = (a.theta + 2.0*k* UtilEjml.PI)/N;
    }
1198789
function() {

      var a = this['re'];
      var b = this['im'];
      var r = this['abs']();

      var re, im;

      if (a >= 0) {

        if (b === 0) {
          return new Complex(Math.sqrt(a), 0);
        }

        re = 0.5 * Math.sqrt(2.0 * (r + a));
      } else {
        re = Math.abs(b) / Math.sqrt(2 * (r - a));
      }

      if (a <= 0) {
        im = 0.5 * Math.sqrt(2.0 * (r - a));
      } else {
        im = Math.abs(b) / Math.sqrt(2 * (r + a));
      }

      return new Complex(re, b < 0 ? -im : im);
    }
13486
def vol(poly1, poly2):
    "
    c_point = central_p(poly1, poly2)
    c_point = (c_point[0], c_point[1], c_point[2])
    vol_therah = 0
    num = len(poly1)
    poly1.append(poly1[0])
    poly2.append(poly2[0])
    for i in range(num - 2):
        # the upper part
        tehrahedron = [c_point, poly1[0], poly1[i+1], poly1[i+2]]
        vol_therah += vol_tehrahedron(tehrahedron)
        # the bottom part
        tehrahedron = [c_point, poly2[0], poly2[i+1], poly2[i+2]]
        vol_therah += vol_tehrahedron(tehrahedron)
    # the middle part
    for i in range(num):
        tehrahedron = [c_point, poly1[i], poly2[i], poly2[i+1]]
        vol_therah += vol_tehrahedron(tehrahedron)
        tehrahedron = [c_point, poly1[i], poly1[i+1], poly2[i]]
        vol_therah += vol_tehrahedron(tehrahedron)
    return vol_therah
644475
def all(A, ax=None):
    
    if isinstance(A, Poly):
        out = numpy.zeros(A.shape, dtype=bool)
        B = A.A
        for key in A.keys:
            out += all(B[key], ax)
        return out

    return numpy.all(A, ax)
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;
    }
152239
def get_support_variables(polynomial):
    
    support = []
    if is_number_type(polynomial):
        return support
    for monomial in polynomial.expand().as_coefficients_dict():
        mon, _ = __separate_scalar_factor(monomial)
        symbolic_support = flatten(split_commutative_parts(mon))
        for s in symbolic_support:
            if isinstance(s, Pow):
                base = s.base
                if is_adjoint(base):
                    base = base.adjoint()
                support.append(base)
            elif is_adjoint(s):
                support.append(s.adjoint())
            elif isinstance(s, Operator):
                support.append(s)
    return support
361452
def _compute_term_3(self, C, mag, R):
        
        return (
            (C['a12'] + C['a13'] * mag + C['a14'] * np.power(mag, 2) +
             C['a15'] * np.power(mag, 3)) * np.power(R, 2)
        )
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;
    }