Solve a simple maths equation manually
492652
def equation(operation, firstnum, secondnum):
    
    if operation == 'plus':
        return firstnum + secondnum
    elif operation == 'minus':
        return firstnum - secondnum
    elif operation == 'multiply':
        return firstnum * secondnum
    elif operation == 'divide':
        if not secondnum == 0:
            return firstnum / secondnum
        raise ZeroDivisionError("Unable to divide by 0.")
    raise ValueError('Invalid operation provided.')
151450
def equation(self):
        
        mat = np.empty((self.nunknowns, self.model.neq))
        rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
        for icp in range(self.ncp):
            istart = icp * self.nlayers
            ieq = 0
            for e in self.model.elementlist:
                if e.nunknowns > 0:
                    qxin, qyin = e.disvecinflayers(self.xcin[icp], self.ycin[icp], self.layers, aq=self.aqin)
                    qxout, qyout = e.disvecinflayers(self.xcout[icp], self.ycout[icp], self.layers, aq=self.aqout)
                    mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] = \
                        (qxin - qxout) * self.cosnorm[icp] + (qyin - qyout) * self.sinnorm[icp]
                    ieq += e.nunknowns
                else:
                    qxin, qyin = e.disveclayers(self.xcin[icp], self.ycin[icp], self.layers, aq=self.aqin)
                    qxout, qyout = e.disveclayers(self.xcout[icp], self.ycout[icp], self.layers, aq=self.aqout)
                    rhs[istart:istart + self.nlayers] -= (qxin - qxout) * self.cosnorm[icp] + (qyin - qyout) * \
                                                         self.sinnorm[icp]
        return mat, rhs
1039381
def solve(guess_a, guess_b, power, solver='scipy'):
    
    # The problem is 2 dimensional so we need 2 symbols
    x = sp.symbols('x:2', real=True)
    # There is a user specified parameter ``p`` in this problem:
    p = sp.Symbol('p', real=True, negative=False, integer=True)
    # Our system consists of 2-non-linear equations:
    f = [x[0] + (x[0] - x[1])**p/2 - 1,
         (x[1] - x[0])**p/2 + x[1]]
    # We construct our ``SymbolicSys`` instance by passing variables, equations and parameters:
    neqsys = SymbolicSys(x, f, [p])  # (this will derive the Jacobian symbolically)

    # Finally we solve the system using user-specified ``solver`` choice:
    return neqsys.solve([guess_a, guess_b], [power], solver=solver)
657725
def simplify(self):
        
        expr = self.expr
        self.expr = _cancel_mul(expr, self.registry)
        return self
671074
def elasprepFD(self):
    
    self.dx4 = self.dx**4
    self.D = self.E*self.Te**3/(12*(1-self.nu**2))
540522
def solve(self, y, h, t_end):
		
		ts = []
		ys = []
		yi = y
		ti = 0.0
		while ti < t_end:
			ts.append(ti)
			yi = self.step(yi, None, ti, h)
			ys.append(yi)
			ti += h
		return ts, ys
799002
def equation(self):
        
        mat = np.empty((self.nunknowns, self.model.neq))
        # rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
        rhs = self.pc.copy()
        for icp in range(self.ncp):
            istart = icp * self.nlayers
            # rhs[istart:istart+self.nlayers] = self.pc[]
            ieq = 0
            for e in self.model.elementlist:
                if e.nunknowns > 0:
                    mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] = \
                        e.potinflayers(self.xc[icp], self.yc[icp], self.layers)
                    if e == self:
                        mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] -= self.resfac[icp]
                    ieq += e.nunknowns
                else:
                    rhs[istart:istart + self.nlayers] -= \
                        e.potentiallayers(self.xc[icp], self.yc[icp],
                                          self.layers)  # Pretty cool that this works, really
        return mat, rhs
323191
def beta(self):
        '''
        
'''
        t = self.X.transpose()
        XX = dot(t,self.X)
        XY = dot(t,self.y)
        return linalg.solve(XX,XY)
1002257
def expr2dimacssat(ex):
    """"""
    if not ex.simple:
        raise ValueError("expected ex to be simplified")

    litmap, nvars = ex.encode_inputs()

    formula = _expr2sat(ex, litmap)
    if 'xor' in formula:
        if '=' in formula:
            fmt = 'satex'
        else:
            fmt = 'satx'
    elif '=' in formula:
        fmt = 'sate'
    else:
        fmt = 'sat'

    return "p {} {}\n{}".format(fmt, nvars, formula)
544960
def print_equations(self):
        
        print("P_static: ", self.eqn_st)
        print("P_thermal: ", self.eqn_th)
        print("P_anharmonic: ", self.eqn_anh)
        print("P_electronic: ", self.eqn_el)
298723
def equation(self):
        
        mat = np.zeros((self.nunknowns, self.model.neq))
        rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
        rhs[0:self.nlayers - 1] = 0.0
        rhs[self.nlayers - 1] = self.Qc
        ieq = 0
        for e in self.model.elementlist:
            if e.nunknowns > 0:
                head = e.potinflayers(self.xc, self.yc, self.layers) / self.aq.Tcol[self.layers, :]
                mat[0:self.nlayers - 1, ieq:ieq + e.nunknowns] = head[:-1] - head[1:]
                if e == self:
                    for i in range(self.nlayers - 1):
                        mat[i, ieq + i] -= self.resfac[i]
                        mat[i, ieq + i + 1] += self.resfac[i + 1]
                    mat[self.nlayers - 1, ieq:ieq + self.nlayers] = 1.0
                ieq += e.nunknowns
            else:
                head = e.potentiallayers(self.xc, self.yc, self.layers) / self.aq.T[self.layers]
                rhs[0:self.nlayers - 1] -= head[:-1] - head[1:]
        return mat, rhs
609181
def BB(n):
    """"""
    if (n<=1):return Context('0\n1')
    else:
        BB1=BB(n-1)
        AA1=AA(n-1)
        r1 = C1((n-1)*2**(n-2),2**(n-1)) - AA1 - BB1
        r2 = BB1 - C1(2**(n-1),2**(n-1)) - BB1;
        return r1 + r2
120968
def _solve(self, x0, A, l, u, xmin, xmax):
        
        s = pips(self._costfcn, x0, A, l, u, xmin, xmax,
                 self._consfcn, self._hessfcn, self.opt)
        return s
1200752
def _cho_solve(c_and_lower, b, check_finite=True):
    """"""

    L = c_and_lower[0]
    y = cpxl.solve_triangular(L, b, trans=0, lower=True,
                              check_finite=check_finite)
    return cpxl.solve_triangular(L, y, trans=1, lower=True,
                                 check_finite=check_finite)
1318529
function(eqns, solve_for, solutions) { 
        if(eqns instanceof Equation) {
            //if it's zero then we're done
            if(eqns.isZero())
                return [new Symbol(0)];
            //if the lhs = x then we're done
            if(eqns.LHS.equals(solve_for) && !eqns.RHS.contains(solve_for)) {
                return [eqns.RHS];
            }
            //if the rhs = x then we're done
            if(eqns.RHS.equals(solve_for) && !eqns.LHS.contains(solve_for)) {
                return [eqns.LHS];
            }
        }
        //unwrap the vector since what we want are the elements
        if(eqns instanceof core.Vector)
            eqns = eqns.elements;
        solve_for = solve_for || 'x'; //assumes x by default
        //If it's an array then solve it as a system of equations
        if(isArray(eqns)) {
            return sys_solve.apply(undefined, arguments);
        }
        //parse out functions. Fix for issue #300
        //eqns = core.Utils.evaluate(eqns);
        
        solutions = solutions || [];
        //maybe we get lucky
        if(eqns.group === S && eqns.contains(solve_for)) {
            solutions.push(new Symbol(0));
            return solutions;
        }  
        if(eqns.group === CB) {
            var sf = String(solve_for); //everything else belongs to the coeff
            eqns.each(function(x) {
                if(x.contains(sf))
                    solve(x, solve_for, solutions);
            });
 
            return solutions;
        }
        var existing = {}, //mark existing solutions as not to have duplicates
            add_to_result = function(r, has_trig) {
                var r_is_symbol = isSymbol(r);
                if(r === undefined || typeof r === 'number' && isNaN(r))
                    return;
                if(isArray(r)) 
                    solutions = solutions.concat(r);
                else { 
                    if(r.valueOf() !== 'null') {
                        if(!r_is_symbol)
                            r = _.parse(r);
                        //try to convert the number to multiples of pi
                        if(core.Settings.make_pi_conversions && has_trig) {
                            var temp = _.divide(r.clone(), new Symbol(Math.PI)),
                                m = temp.multiplier,
                                a = Math.abs(m.num),
                                b = Math.abs(m.den);
                            if(a < 10 && b < 10)
                                r = _.multiply(temp, new Symbol('pi'));
                        }
                        //convert to a string so we can mark it as a known solution
                        var r_str = r.toString();
                        if(!existing[r_str])
                            solutions.push(r);
                        //mark the answer as seen
                        existing[r_str] = true;
                    }
                }
            };
        //gets points around which to solve. It does that because it builds on the principle that if
        //the sign changes over an interval then there must be a zero on that interval
        
        
        var attempt_Newton = function(symbol) { 
            var has_trig = symbol.hasTrig();
            // we get all the points where a possible zero might exist
            var points1 = get_points(symbol, 0.1);
            var points2 = get_points(symbol, 0.05);
            var points3 = get_points(symbol, 0.01);
            var points = core.Utils.arrayUnique(points1.concat(points2).concat(points3)),
                l = points.length;
            //compile the function and the derivative of the function
            var f = build(symbol.clone()),
                fp = build(_C.diff(symbol.clone()));
            for(var i=0; i<l; i++) {
                var point = points[i];
                add_to_result(Newton(point, f, fp), has_trig);
            }
            solutions.sort();
        };

        var eq = core.Utils.isSymbol(eqns) ? eqns : toLHS(eqns),
            vars = core.Utils.variables(eq),//get a list of all the variables
            numvars = vars.length;//how many variables are we dealing with
        //if we're dealing with a single variable then we first check if it's a 
        //polynomial (including rationals).If it is then we use the Jenkins-Traubb algorithm.     
        //Don't waste time
        if(eq.group === S || eq.group === CB && eq.contains(solve_for))
            return [new Symbol(0)];
        //force to polynomial. We go through each and then we look at what it would 
        //take for its power to be an integer
        //if the power is a fractional we divide by the fractional power
        var fractionals = {},
            cfact;
        var correct_denom = function(symbol) { 
            var original = symbol.clone(); //preserve the original
            if(symbol.symbols) {
                for(var x in symbol.symbols) { 
                    var sym = symbol.symbols[x];
                    var parts = explode(sym, solve_for);
                    var is_sqrt = parts[1].fname === core.Settings.SQRT;
                    var v = Symbol.unwrapSQRT(parts[1]);
                    var p = v.power.clone();
                    //circular logic with sqrt. Since sqrt(x) becomes x^(1/2) which then becomes sqrt(x), this continues forever
                    //this needs to be terminated if p = 1/2
                    if(!isSymbol(p) && !p.equals(1/2)) {
                        if(p.den.gt(1)) { 
                            if(is_sqrt) {
                                symbol = _.subtract(symbol, sym.clone());
                                symbol = _.add(symbol, _.multiply(parts[0].clone(), v));
                                return correct_denom(symbol);
                            }
                            var c = fractionals[p.den];
                            fractionals[p.den] = c ? c++ : 1;
                        }
                        else if(p.sign() === -1){ 
                            var factor = _.parse(solve_for+'^'+Math.abs(p)); //this
                            //unwrap the symbol's denoniator
                            symbol.each(function(y, index) {
                                if(y.contains(solve_for)) {
                                    symbol.symbols[index] = _.multiply(y, factor.clone());
                                } 
                            });
                            fractionals = {};
                            return correct_denom(_.parse(symbol));
                        }
                        else if(sym.group === PL) {
                            var min_p = core.Utils.arrayMin(core.Utils.keys(sym.symbols));
                            if(min_p < 0) {
                                var factor = _.parse(solve_for+'^'+Math.abs(min_p));
                                var corrected = new Symbol(0);
                                original.each(function(x) {
                                    corrected = _.add(corrected, _.multiply(x.clone(), factor.clone()));
                                }, true);
                                return corrected;
                            }
                        }
                    }
                }
            }
            return symbol;
        };
        //rewrites equations/expression in simpler form
        var rewrite = function(rhs, lhs) { 
            lhs = lhs || new Symbol(0);
            rhs = Symbol.unwrapSQRT(_.expand(rhs)); //expand the term expression go get rid of quotients when possible
            var c = 0, //a counter to see if we have all terms with the variable
                l = rhs.length;
            //try to rewrite the whole thing
            if(rhs.group === CP && rhs.contains(solve_for) && rhs.isLinear()) {
                rhs.distributeMultiplier();
                var t = new Symbol(0);
                //first bring all the terms containing the variable to the lhs
                rhs.each(function(x) {
                    if(x.contains(solve_for)) {
                        c++;
                        t = _.add(t, x.clone());
                    }
                    else
                        lhs = _.subtract(lhs, x.clone());
                });
                rhs = t;
                
                //if not all the terms contain the variable so it's in the form
                //a*x^2+x
                if(c !== l)
                    return rewrite(rhs, lhs);
                else { 
                    return [rhs, lhs];
                }
            }
            else if(rhs.group === CB && rhs.contains(solve_for) && rhs.isLinear()) {
                if(rhs.multiplier.lessThan(0)) {
                    rhs.multiplier = rhs.multiplier.multiply(new core.Frac(-1));
                    lhs.multiplier = lhs.multiplier.multiply(new core.Frac(-1));
                }
                if(lhs.equals(0))
                    return new Symbol(0);
                else {
                    var t = new Symbol(1);
                    rhs.each(function(x) { 
                        if(x.contains(solve_for)) 
                            t = _.multiply(t, x.clone());
                        else 
                            lhs = _.divide(lhs, x.clone());
                    });
                    rhs = t;
                    return rewrite(rhs, lhs);
                    
                }
            }   
            else if(!rhs.isLinear() && rhs.contains(solve_for)) { 
                var p = _.parse(rhs.power.clone().invert());
                rhs = _.pow(rhs, p.clone());
                lhs = _.pow(_.expand(lhs), p.clone());
                return rewrite(rhs, lhs);
            }
            else if(rhs.group === FN || rhs.group === S || rhs.group === PL) {
                return [rhs, lhs];
            }
        };

        //separate the equation
        var separate = function(eq) {
            var lhs = new Symbol(0),
                rhs = new Symbol(0);
            eq.each(function(x) {
                if(x.contains(solve_for, true))
                    lhs = _.add(lhs, x.clone());
                else
                    rhs = _.subtract(rhs, x.clone());
            });
            return [lhs, rhs];
        };
        
        //first remove any denominators
        eq = correct_denom(eq);  

        if(eq.equals(0))
            return [eq];
        //correct fractionals. I can only handle one type right now
        var fkeys = core.Utils.keys(fractionals);
        if(fkeys.length === 1) { 
            //make a note of the factor
            cfact = fkeys[0];
            eq.each(function(x, index) {
                if(x.contains(solve_for)) {
                    var parts = explode(x, solve_for);
                    var v = parts[1];
                    var p = v.power;
                    if(p.den.gt(1)) { 
                        v.power = p.multiply(new core.Frac(cfact));
                        eq.symbols[index] = _.multiply(v, parts[0]);
                    }
                } 
            });
            eq = _.parse(eq);
        }
        //polynomial single variable
        if(numvars === 1) { 
            if(eq.isPoly(true)) { 
                var coeffs = core.Utils.getCoeffs(eq, solve_for),
                    deg = coeffs.length - 1,
                    was_calculated = false;
                if(vars[0] === solve_for) {
                    //check to see if all the coefficients are constant
                    if(checkAll(coeffs, function(x) {
                        return x.group !== core.groups.N;
                    })) {
                        var roots = core.Algebra.proots(eq);
                        //if all the roots are integers then return those
                        if(checkAll(roots, function(x) {
                            return !core.Utils.isInt(x);
                        })) {
                            //roots have been calculates
                            was_calculated = true;
                            roots.map(function(x) {
                                add_to_result(new Symbol(x));
                            });
                        }
                    }
                    
                    if(!was_calculated) {
                        //we can solve algebraically for degrees 1, 2, 3. The remainder we switch to Jenkins-
                        if(deg === 1) 
                            add_to_result(_.divide(coeffs[0], coeffs[1].negate()));
                        else if(deg === 2) { 
                            add_to_result(_.expand(quad.apply(undefined, coeffs)));
                            coeffs.push('-');
                            add_to_result(_.expand(quad.apply(undefined, coeffs)));
                        }
                        else if(deg === 3)
                            add_to_result(cubic.apply(undefined, coeffs));
                        else {
                            /*
                            var sym_roots = csolve(eq, solve_for); 
                            if(sym_roots.length === 0)
                                sym_roots = divnconsolve(eq, solve_for);
                            if(sym_roots.length > 0) 
                                add_to_result(sym_roots);
                            else
                            */
                            _A.proots(eq).map(add_to_result);
                        }
                    }
                            
                }
            }
            else { 
                //since it's not a polynomial then we'll try to look for a solution using Newton's method
                //this is not a very broad search but takes the positions that something is better than nothing
                attempt_Newton(eq);
            }
        }
        else { 
            //The idea here is to go through the equation and collect the coefficients
            //place them in an array and call the quad or cubic function to get the results
            if(!eq.hasFunc(solve_for) && eq.isComposite()) { 
                try {
                    var coeffs = core.Utils.getCoeffs(eq, solve_for);
                    
                    var l = coeffs.length,
                        deg = l-1; //the degree of the polynomial
                    //handle the problem based on the degree
                    switch(deg) {
                        case 0:
                            var separated = separate(eq);
                            var lhs = separated[0],
                                rhs = separated[1];
                            if(lhs.group === core.groups.EX) {
                                add_to_result(_.parse(core.Utils.format('log(({0})/({2}))/log({1})', rhs, lhs.value, lhs.multiplier)));
                            }
                            break;
                        case 1:
                            //nothing to do but to return the quotient of the constant and the LT
                            //e.g. 2*x-1
                            add_to_result(_.divide(coeffs[0], coeffs[1].negate()));
                            break;
                        case 2:
                            add_to_result(quad.apply(undefined, coeffs));
                            coeffs.push('-');
                            add_to_result(quad.apply(undefined, coeffs));
                            break;
                        case 3:
                            add_to_result(cubic.apply(undefined, coeffs));
                            break;
                        case 4:
                            add_to_result(quartic.apply(undefined, coeffs));
                            break;
                        default:
                            add_to_result(csolve(eq, solve_for));
                            if(solutions.length === 0)
                                add_to_result(divnconsolve(eq, solve_for));
                    }
                }
                catch(e) { /*something went wrong. EXITING*/; } 
            }
            else {
                try {
                    var rw = rewrite(eq);
                    var lhs = rw[0];
                    var rhs = rw[1];
                    if(lhs.group === FN) {
                        if(lhs.fname === 'abs') {
                            solutions.push(rhs.clone());
                            solutions.push(rhs.negate());
                        }
                        else
                            solutions.push(_.subtract(lhs, rhs));
                    }
                }
                catch(error) {
                    //Let's try this another way
                    try { 
                        //1. if the symbol is in the form a*b*c*... then the solution is zero if 
                        //either a or b or c is zero.
                        if(eq.group === CB)
                            solutions.push(0);
                        else if(eq.group === CP) {
                            var separated = separate(eq);
                            var lhs = separated[0],
                                rhs = separated[1];
                            
                            //reduce the equation
                            if(lhs.group === core.groups.EX && lhs.value === solve_for) {
                                //change the base of both sides
                                var p = lhs.power.clone().invert();
                                solutions.push(_.pow(rhs, p));
                            }
                        }
                    }
                    catch(error){;}
                }
            }
        }
        
        if(cfact) {
            solutions = solutions.map(function(x) {
                return _.pow(x, new Symbol(cfact));
            });
        }
        
        return solutions;
    }
858168
public static int solveQuad (double[] eqn, double[] res) {
        double a = eqn[2];
        double b = eqn[1];
        double c = eqn[0];
        int rc = 0;
        if (a == 0f) {
            if (b == 0f) {
                return -1;
            }
            res[rc++] = -c / b;
        } else {
            double d = b * b - 4f * a * c;
            // d < 0f
            if (d < 0f) {
                return 0;
            }
            d = Math.sqrt(d);
            res[rc++] = (-b + d) / (a * 2f);
            // d != 0f
            if (d != 0f) {
                res[rc++] = (-b - d) / (a * 2f);
            }
        }
        return fixRoots(res, rc);
    }
1026655
def greenhall_BasicSum(J, M, S, F, alpha, d):
    
    first = pow(greenhall_sz(0, F, alpha, d), 2)
    second = (1-float(J)/float(M))*pow(greenhall_sz(float(J)/float(S), F, alpha, d), 2)
    third = 0
    for j in range(1, int(J)):
        third += 2*(1.0-float(j)/float(M))*pow(greenhall_sz(float(j)/float(S), F, alpha, d), 2)
    return first+second+third