def frexp(x):
    """
    Version of frexp that works for numbers with uncertainty, and also
    for regular numbers.
    """

    # 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.).

    aff_func = to_affine_scalar(x)

    if aff_func.derivatives:
        result = math.frexp(aff_func.nominal_value)
        # With frexp(x) = (m, e), dm/dx = 1/(2**e):
        factor = 1/(2**result[1])
        return (
            AffineScalarFunc(
                result[0],
                # Chain rule:
                dict([(var, factor*deriv)
                      for (var, deriv) in aff_func.derivatives.iteritems()])),
            # The exponent is an integer and is supposed to be
            # continuous (small errors):
            result[1])
    else:
        # This function was not called with an AffineScalarFunc
        # argument: there is no need to return numbers with uncertainties:
        return math.frexp(x)

 def rebin(self, x_rebin_fact, y_rebin_fact):
     """ Rebin the data and adjust dims """
     if self.data == None:
         raise Exception('Please read in the file you wish to rebin first')
     (mantis_x, exp_x) = math.frexp(x_rebin_fact)
     (mantis_y, exp_y) = math.frexp(y_rebin_fact)
     # FIXME - this is a floating point comparison, is it always exact?
     if (mantis_x != 0.5 or mantis_y != 0.5):
         raise Exception('Rebin factors not power of 2 not supported (yet)')
     if int(self.dim1 / x_rebin_fact) * x_rebin_fact != self.dim1 or \
        int(self.dim2 / x_rebin_fact) * x_rebin_fact != self.dim2 :
         raise('image size is not divisible by rebin factor - ' + \
               'skipping rebin')
     pass  ## self.data.savespace(1) # avoid the upcasting behaviour
     i = 1
     while i < x_rebin_fact:
         # FIXME - why do you divide by 2? Rebinning should increase counts?
         self.data = ((self.data[:, ::2] + self.data[:, 1::2]) / 2)
         i = i * 2
     i = 1
     while i < y_rebin_fact:
         self.data = ((self.data[::2, :] + self.data[1::2, :]) / 2)
         i = i * 2
     self.resetvals()
     self.dim1 = self.dim1 / x_rebin_fact
     self.dim2 = self.dim2 / y_rebin_fact
     #update header
     self.update_header()

def _product(values):
    """Return product of values as (exponent, mantissa)."""
    errmsg = 'mixed Decimal and float is not supported'
    prod = 1
    for x in values:
        if isinstance(x, float):
            break
        prod *= x
    else:
        return (0, prod)
    if isinstance(prod, Decimal):
        raise TypeError(errmsg)
    # Since floats can overflow easily, we calculate the product as a
    # sort of poor-man's BigFloat. Given that:
    #
    #   x = 2**p * m  # p == power or exponent (scale), m = mantissa
    #
    # we can calculate the product of two (or more) x values as:
    #
    #   x1*x2 = 2**p1*m1 * 2**p2*m2 = 2**(p1+p2)*(m1*m2)
    #
    mant, scale = 1, 0  #math.frexp(prod)  # FIXME
    for y in chain([x], values):
        if isinstance(y, Decimal):
            raise TypeError(errmsg)
        m1, e1 = math.frexp(y)
        m2, e2 = math.frexp(mant)
        scale += (e1 + e2)
        mant = m1*m2
    return (scale, mant)

def frexp(x):
    """
    Version of frexp that works for numbers with uncertainty, and also
    for regular numbers.
    """

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

    aff_func = to_affine_scalar(x)

    if aff_func._linear_part:
        (mantissa, exponent) = math.frexp(aff_func.nominal_value)
        return (
            AffineScalarFunc(
                mantissa,
                # With frexp(x) = (m, e), x = m*2**e, so m = x*2**-e
                # and therefore dm/dx = 2**-e (as e in an integer that
                # does not vary when x changes):
                LinearCombination([2**-exponent, aff_func._linear_part])),
            # The exponent is an integer and is supposed to be
            # continuous (errors must be small):
            exponent)
    else:
        # This function was not called with an AffineScalarFunc
        # argument: there is no need to return numbers with uncertainties:
        return math.frexp(x)

def nearly_equal(float1, float2, places=15):
    """
    Determines whether two floating point numbers are nearly equal (to
    within reasonable rounding errors
    """
    mantissa1, exp1 = math.frexp(float1)
    mantissa2, exp2 = math.frexp(float2)
    return (round(mantissa1, places) == round(mantissa2, places) and
            exp1 == exp2)

 def _not_nearly_equal(self, float1, float2):
     """
     Determines whether two floating point numbers are nearly equal (to
     within reasonable rounding errors
     """
     mantissa1, exp1 = math.frexp(float1)
     mantissa2, exp2 = math.frexp(float2)
     return not ((round(mantissa1, self.nearly_equal_places) ==
                  round(mantissa2, self.nearly_equal_places)) and
                 exp1 == exp2)

    def __init__(self, L):

        self.L = L
        if math.frexp(self.L)[0]!=0.5:
            print "profile size is not a power of 2"
            pdb.set_trace()

        self.J = math.frexp(self.L)[1]-1
        self.data = False
        self.value = dict()

    def SetFog(self,fog):
        projection = (fog.function >> 3) & 1

        if projection:
            if fog.z_far == fog.z_near or fog.z_end == fog.z_start:
                A = 0
                C = 0
            else:
                A = (fog.z_far - fog.z_near)/(fog.z_end - fog.z_start)
                C = (fog.z_start - fog.z_near)/(fog.z_end - fog.z_start)
            b_shift = 0
            b_magnitude = 0

        else:
            if fog.z_far == fog.z_near or fog.z_end == fog.z_start:
                A = 0
                B = 0.5
                C = 0
            else:
                A = fog.z_far*fog.z_near/((fog.z_far - fog.z_near)*(fog.z_end - fog.z_start))
                B = fog.z_far/(fog.z_far - fog.z_near)
                C = fog.z_start/(fog.z_end - fog.z_start)

            if B > 1:
                b_shift = 1 + int(ceil(log(B,2)))
            elif 0 < B < 0.5:
                b_shift = 0
            else:
                b_shift = 1

            A /= 2**b_shift
            b_magnitude = int(2*(B/2**b_shift)*8388638)


        a_mantissa,a_exponent = frexp(A)
        self.fog_param0[0:11] = int(abs(a_mantissa)*2**12) & 0x7FF
        self.fog_param0[11:19] = a_exponent + 126 if A != 0 else 0
        self.fog_param0[19] = a_mantissa < 0

        self.fog_param1[0:24] = b_magnitude
        self.fog_param2[0:5] = b_shift

        c_mantissa,c_exponent = frexp(C)
        self.fog_param3[0:11] = int(abs(c_mantissa)*2**12) & 0x7FF
        self.fog_param3[11:19] = c_exponent + 126 if C != 0 else 0
        self.fog_param3[19] = c_mantissa < 0
        self.fog_param3[20:21] = projection
        self.fog_param3[21:24] = fog.function

        self.fog_color[0:8] = fog.color.b
        self.fog_color[8:16] = fog.color.g
        self.fog_color[16:24] = fog.color.r

def PyNextAfter(x, y):
    """returns the next float after x in the direction of y if possible, else returns x"""
    # if x or y is Nan, we don't do much
    if IsNaN(x) or IsNaN(y):
        return x

    # we can't progress if x == y
    if x == y:
        return x

    # similarly if x is infinity
    if x >= infinity or x <= -infinity:
        return x

    # return small numbers for x very close to 0.0
    if -minFloat < x < minFloat:
        if y > x:
            return x + smallEpsilon
        else:
            return x - smallEpsilon  # we know x != y

    # it looks like we have a normalized number
    # break x down into a mantissa and exponent
    m, e = math.frexp(x)

    # all the special cases have been handled
    if y > x:
        m += epsilon
    else:
        m -= epsilon

    return math.ldexp(m, e)

    def as_integer_ratio(self, a, **args):
        """Convert real number to a (numer, denom) pair. """
        v, n = math.frexp(a)  # XXX: hack, will work only for floats

        for i in xrange(300):
            if v != math.floor(v):
                v, n = 2*v, n - 1
            else:
                break

        numer, denom = int(v), 1

        m = 1 << abs(n)

        if n > 0:
            numer *= m
        else:
            denom = m

        n, d = self.limit_denom(numer, denom, **args)

        if a and not n:
            return numer, denom
        else:
            return n, d

 def _bus_message_received_cb(self, bus, message):
     """
     @param bus: the message bus sending the message
     @param message: the message received
     """
     if message.get_structure().get_name() == 'level':
         s = message.get_structure()
         peak = list(s['peak'])
         decay = list(s['decay'])
         rms = list(s['rms'])
         for l in peak, decay, rms:
             for index, v in enumerate(l):
                 try:
                     v = frexp(v)
                 except (SystemError, OverflowError, ValueError):
                     # It was an invalid value (e.g. -Inf), so clamp to
                     # something appropriate
                     l[index] = -100.0
         if not self.uiState:
             self.warning("effect %s doesn't have a uiState" %
                          self.name)
         else:
             for k, v in ('peak', peak), ('decay', decay), ('rms', rms):
                 self.uiState.set('volume-%s' % k, v)
             if not self.firstVolumeValueReceived:
                 self.uiState.set('volume-volume', self.effect_getVolume())
                 self.firstVolumeValueReceived = True

 def get01(self, x):
     m, e = math.frexp(x - self._base)
     if m >= 0 and e <= _E_MAX:
         v = (e + m) / (2. * _E_MAX)
         return v
     else:
         return 0 if m < 0 else 1

def _hash_float(space, v):
    if not isfinite(v):
        if isinf(v):
            return HASH_INF if v > 0 else -HASH_INF
        return HASH_NAN

    m, e = math.frexp(v)

    sign = 1
    if m < 0:
        sign = -1
        m = -m

    # process 28 bits at a time;  this should work well both for binary
    # and hexadecimal floating point.
    x = r_uint(0)
    while m:
        x = ((x << 28) & HASH_MODULUS) | x >> (HASH_BITS - 28)
        m *= 268435456.0  # 2**28
        e -= 28
        y = r_uint(m)  # pull out integer part
        m -= y
        x += y
        if x >= HASH_MODULUS:
            x -= HASH_MODULUS

    # adjust for the exponent;  first reduce it modulo HASH_BITS
    e = e % HASH_BITS if e >= 0 else HASH_BITS - 1 - ((-1 - e) % HASH_BITS)
    x = ((x << e) & HASH_MODULUS) | x >> (HASH_BITS - e)

    x = intmask(intmask(x) * sign)
    return -2 if x == -1 else x

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 rank_sum_n_sites(measurements, details=False):
    if math.frexp(len(measurements))[0] != 0.5:
        print("rank_sum_n_sites received an input of length %s, which is not equal to the number of genotypes."
              "Quitting." % len(measurements))
        sys.exit()
    output_indices = []
    for genotype in measurements:
        output_indices.append(measurements.keys().index(genotype))
    done = False
    while not done:
        done = True
        for i in range(len(measurements) - 1):
            if ranksums(measurements[measurements.keys()[output_indices[i]]],
                        measurements[measurements.keys()[output_indices[i+1]]])[0] < 0:
                output_indices[i], output_indices[i + 1] = output_indices[i + 1], output_indices[i]
                done = False
    output = []
    output_look_good = []
    number_loci = 0
    for index in output_indices:
        output.append(measurements.keys()[index])
        if len(measurements.keys()[index]) > number_loci:
            number_loci = len(measurements.keys()[index])
    for index in output_indices:
        output_look_good.append(genotype_look_good(measurements.keys()[index], number_loci))
    output_detailed = []
    for genotype in output:
        fitness = measurements[genotype][1:]
        output_detailed.append([genotype, np.mean(fitness)])
    if not details:
        return output
    else:
        return output_detailed

	def __new__(cls, *args):
		if len(args) == 0:
			self = bytes.__new__(cls)
			self._pad = 0
		elif len(args) == 1:
			if isinstance(args[0], IntType):
				self = bytes.__new__(cls, uitob(args[0]))
			else:
				self = bytes.__new__(cls, args[0])
			self._pad = 0
		elif len(args) == 2:
			if isinstance(args[0], str):
				arg = int(args[0], args[1])
				m, e = math.frexp(args[1])
				bpd = e - (1 if m == 0.5 else 0)
				bitlen = len(args[0]) * bpd # num digits * bits per digit
				shift = ((bitlen + 7) & ~7) - bitlen
				self = bytes.__new__(cls, uitob(arg << shift))
				self._pad = shift & 7
			else:
				self = bytes.__new__(cls, args[1])
				self._pad = args[0]
				if pad > 7 or (pad and not self):
					raise ValueError('invalid pad value')
		else:
			raise TypeError('BitString constructor takes 1 or 2 args')
		return self

def encode_double(value):
    assert isinstance(value, float)
    abs_val = _math.fabs(value)
    # Work around repr()'s propensity to use exponential format
    # for doubles with really large or really small absolute values.
    # This extra work is because XML-RPC explicitly DOES NOT support
    # exponential format, just decimal digits separated by a single
    # period.
    if value and (abs_val < 0.0001):
        # Scary magic for small numbers, split the float into it's
        # mantissa and exponent components and use the exponent to
        # guess reasonable format.
        m, e = _math.frexp(abs_val)
        decimal_places = int(e/-3) + 17 # Derived at by pure, old-fashon,
                                        # trail and error.
        if decimal_places > 109:
            # "%.110f" causes an OverflowError.
            min_value = "0." + "0"*108 + "1"
            if abs_val < eval(min_value):
                value = min_value
            else:
                value = "%1.109f" % value
                value = value.rstrip('0')
        else:
            value = "%1.*f" % (decimal_places, value)
            value = value.rstrip('0')
    elif value and (abs_val > 99999999999999984.0):
        value = "%.1f" % value
    else:
        value = repr(value)
    return "<value><double>" + value + "</double></value>\n"

        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 evaluate_func(next_func):

		if next_func[0] == 'x':
			return x
		elif next_func[0] == 'y':
			return y

		elif next_func[0] == 'abs':
			return abs(evaluate_func(next_func[1]))
		elif next_func[0] == 'cos_pi':
			return cos(pi*evaluate_func(next_func[1]))
		elif next_func[0] == 'sin_pi':
			return sin(pi*evaluate_func(next_func[1]))
		elif next_func[0] == 'sqrt':
			return sqrt(abs(pi*evaluate_func(next_func[1])))
		elif next_func[0] == 'diff':
			return evaluate_func(next_func[1])-evaluate_func(next_func[2])
		elif next_func[0] == 'ave':
			return (evaluate_func(next_func[1])+evaluate_func(next_func[2]))/2.0
		elif next_func[0] == 'prod':
			return evaluate_func(next_func[1])*evaluate_func(next_func[2])
		elif next_func[0] == 'frexp':
			return frexp(evaluate_func(next_func[1]))[0]
		elif next_func[0] == 'x':
			return evaluate_func(next_func[1])
		elif next_func[0] == 'y':
			return evaluate_func(next_func[2])

 def __call__(self, x):
     if x == 0:
         return 0, self.mink
     import math
     m, e = math.frexp(x)
     k = max((e - self.bits, self.mink))
     return int(m * 2 ** (e - k)), k