Python math.ldexp函数代码示例

def get_heavy_correction(beta, heavytype):
    m = get_heavyq_mass(beta, heavytype)

    def prop(p):
        W = 1-np.cos(p)
        Weplus = W + 1./2 + np.sqrt(W + 1./4)
        Wemins = W + 1./2 - np.sqrt(W + 1./4)
        num = np.sin(p) * np.sin(p*t) + m*(1-Wemins) * np.cos(p*t)
        #num = 2*np.sin(p) * np.sin(p*t) #+ m*(1-Wemins) * np.cos(p*t)
        den = - (1-Weplus) + m**2*(1-Wemins)
        return num/den /(2*np.pi)

    Q = ((1+m**2)/(1-m**2))**2
    T  = (1-Q)/2 + np.sqrt(3*Q+Q**2)/2
    W0 = (1+Q)/2 - np.sqrt(3*Q+Q**2)/2
    m1 = np.log(T+np.sqrt(T**2-1))
    K = 2.0/((1-m**2)*(1.0+np.sqrt(Q/(1.0+4*W0))))

    Nt = 65
    it = np.arange(Nt)
    Ct = np.empty(Nt)
    R  = {}
    minDouble  = math.ldexp(1.0, -1022)
    smallEpsilon  = math.ldexp(1.0, -1074)
    epsilon  = math.ldexp(1.0, -53)
    for t in it:
        Ct[t], err = integrate.quad(prop,-np.pi,np.pi)
        if Ct[t] < err or (K*np.exp(-m1*t)) < err:
            R[t] = 1.
        else:
            R[t] = Ct[t] / (K*np.exp(-m1*t))

    return R

def _write_float(f, x):
	import math
	if x < 0:
		sign = 0x8000
		x = x * -1
	else:
		sign = 0
	if x == 0:
		expon = 0
		himant = 0
		lomant = 0
	else:
		fmant, expon = math.frexp(x)
		if expon > 16384 or fmant >= 1:		# Infinity or NaN
			expon = sign|0x7FFF
			himant = 0
			lomant = 0
		else:					# Finite
			expon = expon + 16382
			if expon < 0:			# denormalized
				fmant = math.ldexp(fmant, expon)
				expon = 0
			expon = expon | sign
			fmant = math.ldexp(fmant, 32)
			fsmant = math.floor(fmant)
			himant = long(fsmant)
			fmant = math.ldexp(fmant - fsmant, 32)
			fsmant = math.floor(fmant)
			lomant = long(fsmant)
	_write_short(f, expon)
	_write_long(f, himant)
	_write_long(f, lomant)

	def set_ref(self, refval):
		x = self.teach()
		m,e = math.frexp(refval / 5e-9)
		if m < 0:
			s = 0x80
			m = -m
		else:
			s = 0
		m = math.ldexp(m, -1)
		for i in range(14,20):
			m = math.ldexp(m, 8)
			h = int(m)
			m -= h
			x[i] = s | h
			s = 0
		
		e -= 31
		if e < 0:
			e += 256
		x[20] = e

		y=bytearray("LN")
		for i in x:
			if i == 10 or i == 27 or i == 13 or i == 43:
				y.append(27)
			y.append(i)
		d.wr(y)

        def msum(iterable):
            """Full precision summation.  Compute sum(iterable) without any
            intermediate accumulation of error.  Based on the 'lsum' function
            at http://code.activestate.com/recipes/393090/

            """
            tmant, texp = 0, 0
            for x in iterable:
                mant, exp = math.frexp(x)
                mant, exp = int(math.ldexp(mant, mant_dig)), exp - mant_dig
                if texp > exp:
                    tmant <<= texp-exp
                    texp = exp
                else:
                    mant <<= exp-texp
                tmant += mant
            # Round tmant * 2**texp to a float.  The original recipe
            # used float(str(tmant)) * 2.0**texp for this, but that's
            # a little unsafe because str -> float conversion can't be
            # relied upon to do correct rounding on all platforms.
            tail = max(len(bin(abs(tmant)))-2 - mant_dig, etiny - texp)
            if tail > 0:
                h = 1 << (tail-1)
                tmant = tmant // (2*h) + bool(tmant & h and tmant & 3*h-1)
                texp += tail
            return math.ldexp(tmant, texp)

def sqrt(x):
    sqrt_special = [
        [inf-infj, 0-infj, 0-infj, infj, infj, inf+infj, nan+infj],
        [inf-infj, None, None, None, None, inf+infj, nan+nanj],
        [inf-infj, None, 0-0j, 0+0j, None, inf+infj, nan+nanj],
        [inf-infj, None, 0-0j, 0+0j, None, inf+infj, nan+nanj],
        [inf-infj, None, None, None, None, inf+infj, nan+nanj],
        [inf-infj, complex(float("inf"), -0.0), complex(float("inf"), -0.0), inf, inf, inf+infj, inf+nanj],
        [inf-infj, nan+nanj, nan+nanj, nan+nanj, nan+nanj, inf+infj, nan+nanj]
    ]

    z = _make_complex(x)

    if math.isinf(z.real) or math.isinf(z.imag):
        return sqrt_special[_special_type(z.real)][_special_type(z.imag)]

    abs_x, abs_y = abs(z.real), abs(z.imag)
    if abs_x < _DBL_MIN and abs_y < _DBL_MIN:
        if abs_x > 0 or abs_y > 0:
            abs_x = math.ldexp(abs_x, _CM_SCALE_UP)
            s = math.ldexp(math.sqrt(abs_x +
                                     math.hypot(abs_x,
                                                math.ldexp(abs_y,
                                                           _CM_SCALE_UP))),
                           _CM_SCALE_DOWN)
        else:
            return complex(0, z.imag)
    else:
        abs_x /= 8
        s = 2 * math.sqrt(abs_x + math.hypot(abs_x, abs_y/8))

    if z.real >= 0:
        return complex(s, math.copysign(abs_y/(2*s), z.imag))
    return complex(abs_y/(2*s), math.copysign(s, z.imag))

def _write_float(f, x):
    import math
    if x < 0:
        sign = 32768
        x = x * -1
    else:
        sign = 0
    if x == 0:
        expon = 0
        himant = 0
        lomant = 0
    else:
        fmant, expon = math.frexp(x)
        if expon > 16384 or fmant >= 1 or fmant != fmant:
            expon = sign | 32767
            himant = 0
            lomant = 0
        else:
            expon = expon + 16382
            if expon < 0:
                fmant = math.ldexp(fmant, expon)
                expon = 0
            expon = expon | sign
            fmant = math.ldexp(fmant, 32)
            fsmant = math.floor(fmant)
            himant = long(fsmant)
            fmant = math.ldexp(fmant - fsmant, 32)
            fsmant = math.floor(fmant)
            lomant = long(fsmant)
    _write_ushort(f, expon)
    _write_ulong(f, himant)
    _write_ulong(f, lomant)

def to_float(s, strict=False):
    """
    Convert a raw mpf to a Python float. The result is exact if the
    bitcount of s is <= 53 and no underflow/overflow occurs.

    If the number is too large or too small to represent as a regular
    float, it will be converted to inf or 0.0. Setting strict=True
    forces an OverflowError to be raised instead.
    """
    sign, man, exp, bc = s
    if not man:
        if s == fzero: return 0.0
        if s == finf: return math_float_inf
        if s == fninf: return -math_float_inf
        return math_float_inf/math_float_inf
    if sign:
        man = -man
    try:
        if bc < 100:
            return math.ldexp(man, exp)
        # Try resizing the mantissa. Overflow may still happen here.
        n = bc - 53
        m = man >> n
        return math.ldexp(m, exp + n)
    except OverflowError:
        if strict:
            raise
        # Overflow to infinity
        if exp + bc > 0:
            if sign:
                return -math_float_inf
            else:
                return math_float_inf
        # Underflow to zero
        return 0.0

def ldexp(x, y):
    # The code below is inspired by uncertainties.wrap().  It is
    # simpler because only 1 argument is given, and there is no
    # delegation to other functions involved (as for __mul__, etc.).

    # Another approach would be to add an additional argument to
    # uncertainties.wrap() so that some arguments are automatically
    # considered as constants.

    aff_func = to_affine_scalar(x)  # y must be an integer, for math.ldexp

    if aff_func.derivatives:
        factor = 2**y
        return AffineScalarFunc(
            math.ldexp(aff_func.nominal_value, y),
            # Chain rule:
            dict([(var, factor*deriv)
                  for (var, deriv) in aff_func.derivatives.iteritems()]))
    else:
        # This function was not called with an AffineScalarFunc
        # argument: there is no need to return numbers with uncertainties:

        # aff_func.nominal_value is not passed instead of x, because
        # we do not have to care about the type of the return value of
        # math.ldexp, this way (aff_func.nominal_value might be the
        # value of x coerced to a difference type [int->float, for
        # instance]):
        return math.ldexp(x, y)

	def __init__(self, pj, lo):
		super(float70, self).__init__(pj, lo, lo + 7)
		x = pj.m.rd(self.lo)
		if x & 0x80:
			s = -1
			x ^= 0x80
		else:
			s = 1
		o = -7
		m = 0
		for i in range(6):
			m +=  math.ldexp(x, o)
			o -= 8
			x = pj.m.rd(self.lo + 1 + i)
		e =  pj.m.s8(self.lo + 6)
		self.lcmt = "m %g e %g" % (m, e)
		v = s * math.ldexp(m, e)
		self.val = v

		x = "%.9e" % v
		if x.find(".") == -1 and x.find("e") == -1:
			x += "."
		self.fmt = x
		self.typ = ".FLOAT"
		self.compact = False
		pj.set_label(lo, "c_" + x)
		print(self, self.fmt)

    def test_705836(self):
        # SF bug 705836.  "<f" and ">f" had a severe rounding bug, where a carry
        # from the low-order discarded bits could propagate into the exponent
        # field, causing the result to be wrong by a factor of 2.
        for base in range(1, 33):
            # smaller <- largest representable float less than base.
            delta = 0.5
            while base - delta / 2.0 != base:
                delta /= 2.0
            smaller = base - delta
            # Packing this rounds away a solid string of trailing 1 bits.
            packed = struct.pack("<f", smaller)
            unpacked = struct.unpack("<f", packed)[0]
            # This failed at base = 2, 4, and 32, with unpacked = 1, 2, and
            # 16, respectively.
            self.assertEqual(base, unpacked)
            bigpacked = struct.pack(">f", smaller)
            self.assertEqual(bigpacked, string_reverse(packed))
            unpacked = struct.unpack(">f", bigpacked)[0]
            self.assertEqual(base, unpacked)

        # Largest finite IEEE single.
        big = (1 << 24) - 1
        big = math.ldexp(big, 127 - 23)
        packed = struct.pack(">f", big)
        unpacked = struct.unpack(">f", packed)[0]
        self.assertEqual(big, unpacked)

        # The same, but tack on a 1 bit so it rounds up to infinity.
        big = (1 << 25) - 1
        big = math.ldexp(big, 127 - 24)
        self.assertRaises(OverflowError, struct.pack, ">f", big)

def float_unpack(Q, size):
    """Convert a 16-bit, 32-bit, or 64-bit integer created
    by float_pack into a Python float."""
    if size == 8:
        MIN_EXP = -1021  # = sys.float_info.min_exp
        MAX_EXP = 1024  # = sys.float_info.max_exp
        MANT_DIG = 53  # = sys.float_info.mant_dig
        BITS = 64
    elif size == 4:
        MIN_EXP = -125  # C's FLT_MIN_EXP
        MAX_EXP = 128  # FLT_MAX_EXP
        MANT_DIG = 24  # FLT_MANT_DIG
        BITS = 32
    elif size == 2:
        MIN_EXP = -13
        MAX_EXP = 16
        MANT_DIG = 11
        BITS = 16
    else:
        raise ValueError("invalid size value")

    if not objectmodel.we_are_translated():
        # This tests generates wrong code when translated:
        # with gcc, shifting a 64bit int by 64 bits does
        # not change the value.
        if Q >> BITS:
            raise ValueError("input '%r' out of range '%r'" % (Q, Q >> BITS))

    # extract pieces with assumed 1.mant values
    one = r_ulonglong(1)
    sign = rarithmetic.intmask(Q >> BITS - 1)
    exp = rarithmetic.intmask((Q & ((one << BITS - 1) - (one << MANT_DIG - 1))) >> MANT_DIG - 1)
    mant = Q & ((one << MANT_DIG - 1) - 1)

    if exp == MAX_EXP - MIN_EXP + 2:
        # nan or infinity
        if mant == 0:
            result = rfloat.INFINITY
        else:
            # preserve at most 52 bits of mant value, but pad w/zeros
            exp = r_ulonglong(0x7FF) << 52
            sign = r_ulonglong(sign) << 63
            if MANT_DIG < 53:
                mant = r_ulonglong(mant) << (53 - MANT_DIG)
            if mant == 0:
                result = rfloat.NAN
            else:
                uint = exp | mant | sign
                result = longlong2float(cast(LONGLONG, uint))
            return result
    elif exp == 0:
        # subnormal or zero
        result = math.ldexp(mant, MIN_EXP - MANT_DIG)
    else:
        # normal: add implicit one value
        mant += one << MANT_DIG - 1
        result = math.ldexp(mant, exp + MIN_EXP - MANT_DIG - 1)
    return -result if sign else result

def c_log(x, y):
    # The usual formula for the real part is log(hypot(z.real, z.imag)).
    # There are four situations where this formula is potentially
    # problematic:
    #
    # (1) the absolute value of z is subnormal.  Then hypot is subnormal,
    # so has fewer than the usual number of bits of accuracy, hence may
    # have large relative error.  This then gives a large absolute error
    # in the log.  This can be solved by rescaling z by a suitable power
    # of 2.
    #
    # (2) the absolute value of z is greater than DBL_MAX (e.g. when both
    # z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX)
    # Again, rescaling solves this.
    #
    # (3) the absolute value of z is close to 1.  In this case it's
    # difficult to achieve good accuracy, at least in part because a
    # change of 1ulp in the real or imaginary part of z can result in a
    # change of billions of ulps in the correctly rounded answer.
    #
    # (4) z = 0.  The simplest thing to do here is to call the
    # floating-point log with an argument of 0, and let its behaviour
    # (returning -infinity, signaling a floating-point exception, setting
    # errno, or whatever) determine that of c_log.  So the usual formula
    # is fine here.

    # XXX the following two lines seem unnecessary at least on Linux;
    # the tests pass fine without them
    if not isfinite(x) or not isfinite(y):
        return log_special_values[special_type(x)][special_type(y)]

    ax = fabs(x)
    ay = fabs(y)

    if ax > CM_LARGE_DOUBLE or ay > CM_LARGE_DOUBLE:
        real = math.log(math.hypot(ax/2., ay/2.)) + M_LN2
    elif ax < DBL_MIN and ay < DBL_MIN:
        if ax > 0. or ay > 0.:
            # catch cases where hypot(ax, ay) is subnormal
            real = math.log(math.hypot(math.ldexp(ax, DBL_MANT_DIG),
                                       math.ldexp(ay, DBL_MANT_DIG)))
            real -= DBL_MANT_DIG*M_LN2
        else:
            # log(+/-0. +/- 0i)
            raise ValueError("math domain error")
            #real = -INF
            #imag = atan2(y, x)
    else:
        h = math.hypot(ax, ay)
        if 0.71 <= h and h <= 1.73:
            am = max(ax, ay)
            an = min(ax, ay)
            real = log1p((am-1)*(am+1) + an*an) / 2.
        else:
            real = math.log(h)
    imag = math.atan2(y, x)
    return (real, imag)

def sqrt(x):
    """
       Return the square root of x.

       This has the same branch cut as log().
    """
    #   Method: use symmetries to reduce to the case when x = z.real and y
    #   = z.imag are nonnegative.  Then the real part of the result is
    #   given by
    #
    #       s = sqrt((x + hypot(x, y))/2)
    #
    #   and the imaginary part is
    #
    #        d = (y/2)/s
    #
    #   If either x or y is very large then there's a risk of overflow in
    #   computation of the expression x + hypot(x, y).  We can avoid this
    #   by rewriting the formula for s as:
    #
    #       s = 2*sqrt(x/8 + hypot(x/8, y/8))
    #
    #   This costs us two extra multiplications/divisions, but avoids the
    #   overhead of checking for x and y large.
    #   If both x and y are subnormal then hypot(x, y) may also be
    #   subnormal, so will lack full precision.  We solve this by rescaling
    #   x and y by a sufficiently large power of 2 to ensure that x and y
    #   are normal.
    s, d, ax, ay = .0, .0, math.fabs(x.real), math.fabs(x.imag)

    ret = _SPECIAL_VALUE(x, _sqrt_special_values)
    if ret is not None:
        return ret

    if x.real == .0 and x.imag == .0:
        _real = .0
        _imag = x.imag
        return complex(_real,_imag)

    if ax < sys.float_info.min and ay < sys.float_info.min and (ax > 0. or ay > 0.):
        #here we catch cases where hypot(ax, ay) is subnormal
        ax = math.ldexp(ax, _CM_SCALE_UP);
        s = math.ldexp(math.sqrt(ax + math.hypot(ax, math.ldexp(ay, _CM_SCALE_UP))),_CM_SCALE_DOWN)
    else:
        ax /= 8.0;
        s = 2.0*math.sqrt(ax + math.hypot(ax, ay/8.0));

    d = ay/(2.0*s)

    if x.real >= .0:
        _real = s;
        _imag = math.copysign(d, x.imag)
    else:
        _real = d;
        _imag = math.copysign(s, x.imag)

    return complex(_real,_imag)

def TakeData():
    t0 = time.clock();
    result = "Test started at " + str(datetime.now()) + "\n";

    LIA1.write("STRT");
    #LIA1.write("STRD"); #start scan after 0.5 sec
    t = time.clock() - t0;
    result += "Scan started at t = " + str(t) + "\n";

    time.sleep(0.5); #32s max
    LIA1.write("PAUS"); #pause
    t = time.clock() - t0;
    result += "Scan ended at t = " + str(t) + "\n";
    LIA1.ask("SPTS ?");
    num = LIA1.ask("SPTS ?"); #number of points
    result += "There are " + num + " points\n";
    t = time.clock() - t0;
    result += "Transfer started at t = " + str(t) + "\n";

    #LIA1.term_chars = "\0";
    CH1 = LIA1.ask_for_values("TRCA ? 1,0," + num);
    CH2 = LIA1.ask_for_values("TRCA ? 2,0," + num);
    LIA1.values_format = visa.single;
    BIN1 = LIA1.ask("TRCL ? 1,0," + num);
    BIN2 = LIA1.ask("TRCL ? 2,0," + num);
    t = time.clock() - t0;
    result += "Transfer ended at t = " + str(t) + "\n";

    BIN1 = list(BIN1);
    BIN2 = list(BIN2);
    CONV1 = list(CH1);
    CONV2 = list(CH2);
    if len(BIN1) != 4*int(num) or len(BIN2) != 4*int(num):
        print "ERROR in num!" + num + "\t" + str(len(BIN1));
        
    for i in range(0,int(num)):
        mantissa = ord(BIN1[4*i+1])*256+ord(BIN1[4*i]);
        if mantissa >= 32768:
            mantissa -= 65536;
        CONV1[i] = math.ldexp(mantissa, ord(BIN1[4*i+2])-124);

        mantissa = ord(BIN2[4*i+1])*256+ord(BIN2[4*i]);
        if mantissa >= 32768:
            mantissa -= 65536;
        CONV2[i] = math.ldexp(mantissa, ord(BIN2[4*i+2])-124);
            
        show = str(CH1[i]) + "\t" + str(CONV1[i]) + \
               "\t" + str(CH2[i]) + "\t" + str(CONV2[i]);
        print show;
        result += show + "\n";

    t = time.clock() - t0;
    result += "Evrything ended at t = " + str(t) + "\n";

    output.write(result);

 def __float__(s):
     """Convert s to a Python float. OverflowError will be raised
     if the magnitude of s is too large."""
     try:
         return math.ldexp(s.man, s.exp)
     # Handle case when mantissa has too many bits (will still
     # overflow if exp is too large)
     except OverflowError:
         n = s.bc - 64
         m = s.man >> n
         return math.ldexp(m, s.exp + n)