Python libmp.normalize函数代码示例

def add_terms(terms, prec, target_prec):
    """
    Helper for evalf_add. Adds a list of (mpfval, accuracy) terms.
    """
    if len(terms) == 1:
        if not terms[0]:
            # XXX: this is supposed to represent a scaled zero
            return mpf_shift(fone, target_prec), -1
        return terms[0]
    max_extra_prec = 2*prec
    sum_man, sum_exp, absolute_error = 0, 0, MINUS_INF
    for x, accuracy in terms:
        if not x:
            continue
        sign, man, exp, bc = x
        if sign:
            man = -man
        absolute_error = max(absolute_error, bc+exp-accuracy)
        delta = exp - sum_exp
        if exp >= sum_exp:
            # x much larger than existing sum?
            # first: quick test
            if (delta > max_extra_prec) and \
                ((not sum_man) or delta-bitcount(abs(sum_man)) > max_extra_prec):
                sum_man = man
                sum_exp = exp
            else:
                sum_man += (man << delta)
        else:
            delta = -delta
            # x much smaller than existing sum?
            if delta-bc > max_extra_prec:
                if not sum_man:
                    sum_man, sum_exp = man, exp
            else:
                sum_man = (sum_man << delta) + man
                sum_exp = exp
    if absolute_error == MINUS_INF:
        return None, None
    if not sum_man:
        # XXX: this is supposed to represent a scaled zero
        return mpf_shift(fone, absolute_error), -1
    if sum_man < 0:
        sum_sign = 1
        sum_man = -sum_man
    else:
        sum_sign = 0
    sum_bc = bitcount(sum_man)
    sum_accuracy = sum_exp + sum_bc - absolute_error
    r = normalize(sum_sign, sum_man, sum_exp, sum_bc, target_prec,
        round_nearest), sum_accuracy
    #print "returning", to_str(r[0],50), r[1]
    return r

def evalf_mul(v, prec, options):
    from sympy.core.core import C

    res = pure_complex(v)
    if res:
        # the only pure complex that is a mul is h*I
        _, h = res
        im, _, im_acc, _ = evalf(h, prec, options)
        return None, im, None, im_acc
    args = list(v.args)

    # see if any argument is NaN or oo and thus warrants a special return
    special = []
    for arg in args:
        arg = evalf(arg, prec, options)
        if arg[0] is None:
            continue
        arg = C.Float._new(arg[0], 1)
        if arg is S.NaN or arg.is_unbounded:
            special.append(arg)
    if special:
        from sympy.core.mul import Mul

        special = Mul(*special)
        return evalf(special, prec + 4, {})

    # With guard digits, multiplication in the real case does not destroy
    # accuracy. This is also true in the complex case when considering the
    # total accuracy; however accuracy for the real or imaginary parts
    # separately may be lower.
    acc = prec

    # XXX: big overestimate
    working_prec = prec + len(args) + 5

    # Empty product is 1
    start = man, exp, bc = MPZ(1), 0, 1

    # First, we multiply all pure real or pure imaginary numbers.
    # direction tells us that the result should be multiplied by
    # I**direction; all other numbers get put into complex_factors
    # to be multiplied out after the first phase.
    last = len(args)
    direction = 0
    args.append(S.One)
    complex_factors = []

    for i, arg in enumerate(args):
        if i != last and pure_complex(arg):
            args[-1] = (args[-1] * arg).expand()
            continue
        elif i == last and arg is S.One:
            continue
        re, im, re_acc, im_acc = evalf(arg, working_prec, options)
        if re and im:
            complex_factors.append((re, im, re_acc, im_acc))
            continue
        elif re:
            (s, m, e, b), w_acc = re, re_acc
        elif im:
            (s, m, e, b), w_acc = im, im_acc
            direction += 1
        else:
            return None, None, None, None
        direction += 2 * s
        man *= m
        exp += e
        bc += b
        if bc > 3 * working_prec:
            man >>= working_prec
            exp += working_prec
        acc = min(acc, w_acc)
    sign = (direction & 2) >> 1
    if not complex_factors:
        v = normalize(sign, man, exp, bitcount(man), prec, rnd)
        # multiply by i
        if direction & 1:
            return None, v, None, acc
        else:
            return v, None, acc, None
    else:
        # initialize with the first term
        if (man, exp, bc) != start:
            # there was a real part; give it an imaginary part
            re, im = (sign, man, exp, bitcount(man)), (0, MPZ(0), 0, 0)
            i0 = 0
        else:
            # there is no real part to start (other than the starting 1)
            wre, wim, wre_acc, wim_acc = complex_factors[0]
            acc = min(acc, complex_accuracy((wre, wim, wre_acc, wim_acc)))
            re = wre
            im = wim
            i0 = 1

        for wre, wim, wre_acc, wim_acc in complex_factors[i0:]:
            # acc is the overall accuracy of the product; we aren't
            # computing exact accuracies of the product.
            acc = min(acc, complex_accuracy((wre, wim, wre_acc, wim_acc)))

            use_prec = working_prec
#.........这里部分代码省略.........

def add_terms(terms, prec, target_prec):
    """
    Helper for evalf_add. Adds a list of (mpfval, accuracy) terms.

    Returns
    -------

    - None, None if there are no non-zero terms;
    - terms[0] if there is only 1 term;
    - scaled_zero if the sum of the terms produces a zero by cancellation
      e.g. mpfs representing 1 and -1 would produce a scaled zero which need
      special handling since they are not actually zero and they are purposely
      malformed to ensure that they can't be used in anything but accuracy
      calculations;
    - a tuple that is scaled to target_prec that corresponds to the
      sum of the terms.

    The returned mpf tuple will be normalized to target_prec; the input
    prec is used to define the working precision.

    XXX explain why this is needed and why one can't just loop using mpf_add
    """
    from sympy.core.core import C

    terms = [t for t in terms if not iszero(t)]
    if not terms:
        return None, None
    elif len(terms) == 1:
        return terms[0]

    # see if any argument is NaN or oo and thus warrants a special return
    special = []
    for t in terms:
        arg = C.Float._new(t[0], 1)
        if arg is S.NaN or arg.is_unbounded:
            special.append(arg)
    if special:
        from sympy.core.add import Add

        rv = evalf(Add(*special), prec + 4, {})
        return rv[0], rv[2]

    working_prec = 2 * prec
    sum_man, sum_exp, absolute_error = 0, 0, MINUS_INF

    for x, accuracy in terms:
        sign, man, exp, bc = x
        if sign:
            man = -man
        absolute_error = max(absolute_error, bc + exp - accuracy)
        delta = exp - sum_exp
        if exp >= sum_exp:
            # x much larger than existing sum?
            # first: quick test
            if (delta > working_prec) and ((not sum_man) or delta - bitcount(abs(sum_man)) > working_prec):
                sum_man = man
                sum_exp = exp
            else:
                sum_man += man << delta
        else:
            delta = -delta
            # x much smaller than existing sum?
            if delta - bc > working_prec:
                if not sum_man:
                    sum_man, sum_exp = man, exp
            else:
                sum_man = (sum_man << delta) + man
                sum_exp = exp
    if not sum_man:
        return scaled_zero(absolute_error)
    if sum_man < 0:
        sum_sign = 1
        sum_man = -sum_man
    else:
        sum_sign = 0
    sum_bc = bitcount(sum_man)
    sum_accuracy = sum_exp + sum_bc - absolute_error
    r = normalize(sum_sign, sum_man, sum_exp, sum_bc, target_prec, rnd), sum_accuracy
    # print "returning", to_str(r[0],50), r[1]
    return r

def evalf_mul(v, prec, options):
    args = v.args
    # With guard digits, multiplication in the real case does not destroy
    # accuracy. This is also true in the complex case when considering the
    # total accuracy; however accuracy for the real or imaginary parts
    # separately may be lower.
    acc = prec
    target_prec = prec
    # XXX: big overestimate
    prec = prec + len(args) + 5
    direction = 0
    # Empty product is 1
    man, exp, bc = MPZ(1), 0, 1
    direction = 0
    complex_factors = []
    # First, we multiply all pure real or pure imaginary numbers.
    # direction tells us that the result should be multiplied by
    # i**direction
    for arg in args:
        re, im, re_acc, im_acc = evalf(arg, prec, options)
        if re and im:
            complex_factors.append((re, im, re_acc, im_acc))
            continue
        elif re:
            (s, m, e, b), w_acc = re, re_acc
        elif im:
            (s, m, e, b), w_acc = im, im_acc
            direction += 1
        else:
            return None, None, None, None
        direction += 2*s
        man *= m
        exp += e
        bc += b
        if bc > 3*prec:
            man >>= prec
            exp += prec
        acc = min(acc, w_acc)
    sign = (direction & 2) >> 1
    v = normalize(sign, man, exp, bitcount(man), prec, round_nearest)
    if complex_factors:
        # make existing real scalar look like an imaginary and
        # multiply by the remaining complex numbers
        re, im = v, (0, MPZ(0), 0, 0)
        for wre, wim, wre_acc, wim_acc in complex_factors:
            # acc is the overall accuracy of the product; we aren't
            # computing exact accuracies of the product.
            acc = min(acc,
                      complex_accuracy((wre, wim, wre_acc, wim_acc)))
            A = mpf_mul(re, wre, prec)
            B = mpf_mul(mpf_neg(im), wim, prec)
            C = mpf_mul(re, wim, prec)
            D = mpf_mul(im, wre, prec)
            re, xre_acc = add_terms([(A, acc), (B, acc)], prec, target_prec)
            im, xim_acc = add_terms([(C, acc), (D, acc)], prec, target_prec)

        if options.get('verbose'):
            print "MUL: wanted", target_prec, "accurate bits, got", acc
        # multiply by i
        if direction & 1:
            return mpf_neg(im), re, acc, acc
        else:
            return re, im, acc, acc
    else:
        # multiply by i
        if direction & 1:
            return None, v, None, acc
        else:
            return v, None, acc, None