Python simplify.nsimplify函数代码示例

def get_nullspace(n,q,eigv):
    M = matrix_cache.get_reduced_matrix(n,q, True)
    M.m = M.m._new(M.m.rows, M.m.cols,[nsimplify(v, rational=True) for v in M.m])

    # N = M - lambdaI
    N = M.m - sympy.Matrix.eye(M.get_size())* eigv

    # rationalize entries
    N = N._new(N.rows, N.cols,[nsimplify(v, rational=True) for v in N])

    # find eigenspace - Kernel(M-lambdaI) = Kernel(N)
    EigenSpaceBase = N.nullspace()

    # orthogonal space = Kernel(TransposedEigenVectors)
    TransposedEigenSpaceBase = [list(v) for v in EigenSpaceBase]
    EigenSpaceTransposeMatrix = (sympy.Matrix(numpy.matrix(TransposedEigenSpaceBase)))
    OrthogonalSubspaceBase = EigenSpaceTransposeMatrix.nullspace()

    # need to transpose to put the vectors in columns, and get a matrix for the orthogonal operator
    OrthogonalSubspaceMatrix = (sympy.Matrix(numpy.matrix(OrthogonalSubspaceBase))).transpose()

    # multiply M * Ortohogonal, to get the reduction of M on the orthogonal subspace
    OrthogonalSubspaceMatrix =  OrthogonalSubspaceMatrix._new(OrthogonalSubspaceMatrix.rows, OrthogonalSubspaceMatrix.cols,
                                                              [nsimplify(v, rational=True) for v in OrthogonalSubspaceMatrix])
    EigenValueReducedMatrix = M.m * OrthogonalSubspaceMatrix
    EigenValueReducedMatrix = EigenValueReducedMatrix.rref()

def random_complex_number(a=2, b=-1, c=3, d=1, rational=False):
    """
    Return a random complex number.

    To reduce chance of hitting branch cuts or anything, we guarantee
    b <= Im z <= d, a <= Re z <= c
    """
    A, B = uniform(a, c), uniform(b, d)
    if not rational:
        return A + I*B
    return nsimplify(A, rational=True) + I*nsimplify(B, rational=True)

 def eval(cls, p, q):
     from sympy.simplify.simplify import nsimplify
     if q.is_Number:
         float = not q.is_Rational
         pnew = expand_mul(p)
         if pnew.is_Number:
             float = float or not pnew.is_Rational
             if not float:
                 return pnew % q
             return Float(nsimplify(pnew) % nsimplify(q))
         elif pnew.is_Add and pnew.args[0].is_Number:
             r, p = pnew.as_two_terms()
             p += Mod(r, q)
     return Mod(p, q, evaluate=False)

def random_complex_number(a=2, b=-1, c=3, d=1, rational=False, tolerance=None):
    """
    Return a random complex number.

    To reduce chance of hitting branch cuts or anything, we guarantee
    b <= Im z <= d, a <= Re z <= c

    When rational is True, a rational approximation to a random number
    is obtained within specified tolerance, if any.
    """
    A, B = uniform(a, c), uniform(b, d)
    if not rational:
        return A + I*B
    return (nsimplify(A, rational=True, tolerance=tolerance) +
        I*nsimplify(B, rational=True, tolerance=tolerance))

def _pow_float(inter, power):
    """Evaluates an interval raised to a floating point."""
    power_rational = nsimplify(power)
    num, denom = power_rational.as_numer_denom()
    if num % 2 == 0:
        start = abs(inter.start)**power
        end = abs(inter.end)**power
        if start < 0:
            ret = interval(0, max(start, end))
        else:
            ret = interval(start, end)
        return ret
    elif denom % 2 == 0:
        if inter.end < 0:
            return interval(-float('inf'), float('inf'), is_valid=False)
        elif inter.start < 0:
            return interval(0, inter.end**power, is_valid=None)
        else:
            return interval(inter.start**power, inter.end**power)
    else:
        if inter.start < 0:
            start = -abs(inter.start)**power
        else:
            start = inter.start**power

        if inter.end < 0:
            end = -abs(inter.end)**power
        else:
            end = inter.end**power

        return interval(start, end, is_valid=inter.is_valid)

def test_minpoly_issue_7113():
    # see discussion in https://github.com/sympy/sympy/pull/2234
    from sympy.simplify.simplify import nsimplify
    r = nsimplify(pi, tolerance=0.000000001)
    mp = minimal_polynomial(r, x)
    assert mp == 1768292677839237920489538677417507171630859375*x**109 - \
    2734577732179183863586489182929671773182898498218854181690460140337930774573792597743853652058046464

def _eval_sum_hyper(f, i, a):
    """ Returns (res, cond). Sums from a to oo. """
    from sympy.functions import hyper
    from sympy.simplify import hyperexpand, hypersimp, fraction, simplify
    from sympy.polys.polytools import Poly, factor
    from sympy.core.numbers import Float

    if a != 0:
        return _eval_sum_hyper(f.subs(i, i + a), i, 0)

    if f.subs(i, 0) == 0:
        if simplify(f.subs(i, Dummy('i', integer=True, positive=True))) == 0:
            return S(0), True
        return _eval_sum_hyper(f.subs(i, i + 1), i, 0)

    hs = hypersimp(f, i)
    if hs is None:
        return None

    if isinstance(hs, Float):
        from sympy.simplify.simplify import nsimplify
        hs = nsimplify(hs)

    numer, denom = fraction(factor(hs))
    top, topl = numer.as_coeff_mul(i)
    bot, botl = denom.as_coeff_mul(i)
    ab = [top, bot]
    factors = [topl, botl]
    params = [[], []]
    for k in range(2):
        for fac in factors[k]:
            mul = 1
            if fac.is_Pow:
                mul = fac.exp
                fac = fac.base
                if not mul.is_Integer:
                    return None
            p = Poly(fac, i)
            if p.degree() != 1:
                return None
            m, n = p.all_coeffs()
            ab[k] *= m**mul
            params[k] += [n/m]*mul

    # Add "1" to numerator parameters, to account for implicit n! in
    # hypergeometric series.
    ap = params[0] + [1]
    bq = params[1]
    x = ab[0]/ab[1]
    h = hyper(ap, bq, x)

    return f.subs(i, 0)*hyperexpand(h), h.convergence_statement

    def __new__(cls, *args, **kwargs):
        if iterable(args[0]):
            coords = Tuple(*args[0])
        elif isinstance(args[0], Point):
            coords = args[0].args
        else:
            coords = Tuple(*args)

        if len(coords) != 2:
            raise NotImplementedError("Only two dimensional points currently supported")
        if kwargs.get('evaluate', True):
            coords = [nsimplify(c) for c in coords]

        return GeometryEntity.__new__(cls, *coords)

    def custom_charpoly(self, **flags):
        """
        custom charpoly
        """

        if (self.isSymbolic == True):
            self.m = self.m._new(self.m.rows, self.m.cols,[nsimplify(v, rational=True) for v in self.m])
            max_denom = 0;
            for i in range (0,self.m.rows):
                for j in range (0,self.m.cols):
                    if self.m[i,j] > max_denom:
                        max_denom = self.m[i,j].q
            print max_denom
            self.m *= max_denom
            flags.pop('simplify', None)  # pop unsupported flag
            return self.m.berkowitz_charpoly(Dummy('x'))
        else:
            numpy.rint(self.m)
            return numpy.rint(numpy.poly(self.m))

 def __rpow__(self, other):
     if isinstance(other, (float, int)):
         if not self.is_valid:
             #Don't do anything
             return self
         elif other < 0:
             if self.width > 0:
                 return interval(-float('inf'), float('inf'), is_valid=False)
             else:
                 power_rational = nsimplify(self.start)
                 num, denom = power_rational.as_numer_denom()
                 if denom % 2 == 0:
                     return interval(-float('inf'), float('inf'),
                                     is_valid=False)
                 else:
                     start = -abs(other)**self.start
                     end = start
                     return interval(start, end)
         else:
             return interval(other**self.start, other**self.end)
     elif isinstance(other, interval):
         return other.__pow__(self)
     else:
         return NotImplemented

def polytope_integrate(poly, expr=None, **kwargs):
    """Integrates polynomials over 2/3-Polytopes.

    This function accepts the polytope in `poly` and the function in `expr`
    (uni/bi/trivariate polynomials are implemented) and returns
    the exact integral of `expr` over `poly`.

    Parameters
    ==========

    poly : The input Polygon.

    expr : The input polynomial.

    clockwise : Binary value to sort input points of 2-Polytope clockwise.(Optional)

    max_degree : The maximum degree of any monomial of the input polynomial.(Optional)

    Examples
    ========

    >>> from sympy.abc import x, y
    >>> from sympy.geometry.polygon import Polygon
    >>> from sympy.geometry.point import Point
    >>> from sympy.integrals.intpoly import polytope_integrate
    >>> polygon = Polygon(Point(0,0), Point(0,1), Point(1,1), Point(1,0))
    >>> polys = [1, x, y, x*y, x**2*y, x*y**2]
    >>> expr = x*y
    >>> polytope_integrate(polygon, expr)
    1/4
    >>> polytope_integrate(polygon, polys, max_degree=3)
    {1: 1, x: 1/2, y: 1/2, x*y: 1/4, x*y**2: 1/6, x**2*y: 1/6}
    """
    clockwise = kwargs.get('clockwise', False)
    max_degree = kwargs.get('max_degree', None)

    if clockwise is True:
        if isinstance(poly, Polygon):
            poly = point_sort(poly)
        else:
            raise TypeError("clockwise=True works for only 2-Polytope"
                            "V-representation input")

    if isinstance(poly, Polygon):
        # For Vertex Representation(2D case)
        hp_params = hyperplane_parameters(poly)
        facets = poly.sides
    elif len(poly[0]) == 2:
        # For Hyperplane Representation(2D case)
        plen = len(poly)
        if len(poly[0][0]) == 2:
            intersections = [intersection(poly[(i - 1) % plen], poly[i],
                                          "plane2D")
                             for i in range(0, plen)]
            hp_params = poly
            lints = len(intersections)
            facets = [Segment2D(intersections[i],
                                intersections[(i + 1) % lints])
                      for i in range(0, lints)]
        else:
            raise NotImplementedError("Integration for H-representation 3D"
                                      "case not implemented yet.")
    else:
        # For Vertex Representation(3D case)
        vertices = poly[0]
        facets = poly[1:]
        hp_params = hyperplane_parameters(facets, vertices)

        if max_degree is None:
            if expr is None:
                raise TypeError('Input expression be must'
                                'be a valid SymPy expression')
            return main_integrate3d(expr, facets, vertices, hp_params)

    if max_degree is not None:
        result = {}
        if not isinstance(expr, list) and expr is not None:
            raise TypeError('Input polynomials must be list of expressions')

        if len(hp_params[0][0]) == 3:
            result_dict = main_integrate3d(0, facets, vertices, hp_params,
                                           max_degree)
        else:
            result_dict = main_integrate(0, facets, hp_params, max_degree)

        if expr is None:
            return result_dict

        for poly in expr:
            if poly not in result:
                if poly is S.Zero:
                    result[S.Zero] = S.Zero
                    continue
                integral_value = S.Zero
                monoms = decompose(poly, separate=True)
                for monom in monoms:
                    monom = nsimplify(monom)
                    coeff, m = strip(monom)
                    integral_value += result_dict[m] * coeff
                result[poly] = integral_value
#.........这里部分代码省略.........

    def handle(expr):
        # Handle first reduces to the case
        # expr = 1/d, where d is an add, or d is base**p/2.
        # We do this by recursively calling handle on each piece.
        from sympy.simplify.simplify import nsimplify

        n, d = fraction(expr)

        if expr.is_Atom or (d.is_Atom and n.is_Atom):
            return expr
        elif not n.is_Atom:
            n = n.func(*[handle(a) for a in n.args])
            return _unevaluated_Mul(n, handle(1/d))
        elif n is not S.One:
            return _unevaluated_Mul(n, handle(1/d))
        elif d.is_Mul:
            return _unevaluated_Mul(*[handle(1/d) for d in d.args])

        # By this step, expr is 1/d, and d is not a mul.
        if not symbolic and d.free_symbols:
            return expr

        if ispow2(d):
            d2 = sqrtdenest(sqrt(d.base))**numer(d.exp)
            if d2 != d:
                return handle(1/d2)
        elif d.is_Pow and (d.exp.is_integer or d.base.is_positive):
            # (1/d**i) = (1/d)**i
            return handle(1/d.base)**d.exp

        if not (d.is_Add or ispow2(d)):
            return 1/d.func(*[handle(a) for a in d.args])

        # handle 1/d treating d as an Add (though it may not be)

        keep = True  # keep changes that are made

        # flatten it and collect radicals after checking for special
        # conditions
        d = _mexpand(d)

        # did it change?
        if d.is_Atom:
            return 1/d

        # is it a number that might be handled easily?
        if d.is_number:
            _d = nsimplify(d)
            if _d.is_Number and _d.equals(d):
                return 1/_d

        while True:
            # collect similar terms
            collected = defaultdict(list)
            for m in Add.make_args(d):  # d might have become non-Add
                p2 = []
                other = []
                for i in Mul.make_args(m):
                    if ispow2(i, log2=True):
                        p2.append(i.base if i.exp is S.Half else i.base**(2*i.exp))
                    elif i is S.ImaginaryUnit:
                        p2.append(S.NegativeOne)
                    else:
                        other.append(i)
                collected[tuple(ordered(p2))].append(Mul(*other))
            rterms = list(ordered(list(collected.items())))
            rterms = [(Mul(*i), Add(*j)) for i, j in rterms]
            nrad = len(rterms) - (1 if rterms[0][0] is S.One else 0)
            if nrad < 1:
                break
            elif nrad > max_terms:
                # there may have been invalid operations leading to this point
                # so don't keep changes, e.g. this expression is troublesome
                # in collecting terms so as not to raise the issue of 2834:
                # r = sqrt(sqrt(5) + 5)
                # eq = 1/(sqrt(5)*r + 2*sqrt(5)*sqrt(-sqrt(5) + 5) + 5*r)
                keep = False
                break
            if len(rterms) > 4:
                # in general, only 4 terms can be removed with repeated squaring
                # but other considerations can guide selection of radical terms
                # so that radicals are removed
                if all([x.is_Integer and (y**2).is_Rational for x, y in rterms]):
                    nd, d = rad_rationalize(S.One, Add._from_args(
                        [sqrt(x)*y for x, y in rterms]))
                    n *= nd
                else:
                    # is there anything else that might be attempted?
                    keep = False
                break
            from sympy.simplify.powsimp import powsimp, powdenest

            num = powsimp(_num(rterms))
            n *= num
            d *= num
            d = powdenest(_mexpand(d), force=symbolic)
            if d.is_Atom:
                break

        if not keep:
#.........这里部分代码省略.........