"""
Low-level functions for arbitrary-precision floating-point arithmetic.
"""

import math
import random
import re
import sys

from .backend import BACKEND, MPZ, MPZ_FIVE, MPZ_ONE, MPZ_ZERO, gmpy, int_types
from .libintmath import (bctable, bin_to_radix, isqrt, numeral, sqrtrem,
                         stddigits, trailing)


class ComplexResult(ValueError):
    pass

# All supported rounding modes
round_nearest = sys.intern('n')
round_floor = sys.intern('f')
round_ceiling = sys.intern('c')
round_up = sys.intern('u')
round_down = sys.intern('d')

def prec_to_dps(n):
    """Return number of accurate decimals that can be represented
    with a precision of n bits."""
    return max(1, round(int(n)/blog2_10) - 1)

def dps_to_prec(n):
    """Return the number of bits required to represent n decimals
    accurately."""
    return max(1, round((int(n) + 1)*blog2_10))

def repr_dps(n):
    """Return the number of decimal digits required to represent
    a number with n-bit precision so that it can be uniquely
    reconstructed from the representation."""
    return 1 + math.ceil(int(n)/blog2_10)

#----------------------------------------------------------------------------#
#                    Some commonly needed float values                       #
#----------------------------------------------------------------------------#

# Regular number format:
# (-1)**sign * mantissa * 2**exponent, plus mantissa.bit_length()
fzero = (0, MPZ_ZERO, 0, 0)
fone = (0, MPZ_ONE, 0, 1)
fnone = (1, MPZ_ONE, 0, 1)
ftwo = (0, MPZ_ONE, 1, 1)
ften = (0, MPZ_FIVE, 1, 3)
fhalf = (0, MPZ_ONE, -1, 1)

# Arbitrary encoding for special numbers: zero mantissa, nonzero exponent
fnan = (0, MPZ_ZERO, -123, -1)
finf = (0, MPZ_ZERO, -456, -2)
fninf = (1, MPZ_ZERO, -789, -3)

math_float_inf = math.inf
math_float_nan = math.nan
blog2_10 = 3.3219280948873626

float_mant_dig = sys.float_info.mant_dig
float_min_exp = sys.float_info.min_exp
float_max_exp = sys.float_info.max_exp
float_eps = sys.float_info.epsilon
float_max = sys.float_info.max
float_min = sys.float_info.min
float_min_subnormal_exp = float_min_exp - float_mant_dig


#----------------------------------------------------------------------------#
#                                  Rounding                                  #
#----------------------------------------------------------------------------#

# This function can be used to round a mantissa generally. However,
# we will try to do most rounding inline for efficiency.
def round_int(x, n, rnd):
    if rnd == round_nearest:
        if x >= 0:
            t = x >> (n-1)
            if t & 1 and ((t & 2) or (x & h_mask[n<300][n])):
                return (t>>1)+1
            else:
                return t>>1
        else:
            return -round_int(-x, n, rnd)
    if rnd == round_floor:
        return x >> n
    if rnd == round_ceiling:
        return -((-x) >> n)
    if rnd == round_down:
        if x >= 0:
            return x >> n
        return -((-x) >> n)
    if rnd == round_up:
        if x >= 0:
            return -((-x) >> n)
        return x >> n

# These masks are used to pick out segments of numbers to determine
# which direction to round when rounding to nearest.
class h_mask_big:
    def __getitem__(self, n):
        return (MPZ_ONE<<(n-1))-1

h_mask_small = [0]+[((MPZ_ONE<<(_-1))-1) for _ in range(1, 300)]
h_mask = [h_mask_big(), h_mask_small]

# The >> operator rounds to floor. shifts_down[rnd][sign]
# tells whether this is the right direction to use, or if the
# number should be negated before shifting
shifts_down = {round_floor:(1,0), round_ceiling:(0,1),
    round_down:(1,1), round_up:(0,0)}


#----------------------------------------------------------------------------#
#                          Normalization of raw mpfs                         #
#----------------------------------------------------------------------------#

# This function is called almost every time an mpf is created.
# It has been optimized accordingly.

def normalize(sign, man, exp, bc, prec, rnd):
    """
    Create a raw mpf tuple with value (-1)**sign * man * 2**exp and
    normalized mantissa.  The mantissa is rounded according to the specified
    rounding mode if its size exceeds the precision. Trailing zero bits
    are also stripped from the mantissa to ensure that the
    representation is canonical.

    Conditions on the input:
    * The input must represent a regular (finite) number
    * The sign bit must be 0 or 1
    * The mantissa must be nonnegative
    * The exponent must be an integer
    * The bitcount must be exact

    If these conditions are not met, use from_man_exp, mpf_pos, or any
    of the conversion functions to create normalized raw mpf tuples.
    """
    assert type(man) == MPZ
    assert type(bc) in _exp_types
    assert type(exp) in _exp_types
    assert bc == man.bit_length()
    assert man >= 0

    if not man:
        return fzero
    # Cut mantissa down to size if larger than target precision
    n = bc - prec
    if n > 0:
        if rnd == round_nearest:
            t = man >> (n-1)
            if t & 1 and ((t & 2) or (man & h_mask[n<300][n])):
                man = (t>>1)+1
            else:
                man = t>>1
        elif shifts_down[rnd][sign]:
            man >>= n
        else:
            man = -((-man)>>n)
        exp += n
        bc = prec
    # Strip trailing bits
    if not man & 1:
        t = trailing(man)
        man >>= t
        exp += t
        bc -= t
    # Bit count can be wrong if the input mantissa was 1 less than
    # a power of 2 and got rounded up, thereby adding an extra bit.
    # With trailing bits removed, all powers of two have mantissa 1,
    # so this is easy to check for.
    if man == 1:
        bc = 1
    return sign, man, int(exp), int(bc)

_exp_types = (int,)

if gmpy:
    normalize = gmpy._mpmath_normalize

#----------------------------------------------------------------------------#
#                            Conversion functions                            #
#----------------------------------------------------------------------------#

def from_man_exp(man, exp, prec=0, rnd=round_down):
    """Create raw mpf from (man, exp) pair. The mantissa may be signed.
    If no precision is specified, the mantissa is stored exactly."""
    if isinstance(man, int_types):
        man = MPZ(man)
    else:
        raise TypeError("man expected to be an integer")
    sign = 0
    if man < 0:
        sign = 1
        man = -man
    if man < 1024:
        bc = bctable[man]
    else:
        bc = man.bit_length()
    if not prec:
        if not man:
            return fzero
        if not man & 1:
            t = trailing(man)
            return sign, man >> t, int(exp + t), int(bc - t)
        return sign, man, exp, bc
    return normalize(sign, man, exp, bc, prec, rnd)

int_cache = dict((n, from_man_exp(n, 0)) for n in range(-10, 257))

if gmpy:
    from_man_exp = gmpy._mpmath_create

def from_int(n, prec=0, rnd=round_down):
    """Create a raw mpf from an integer. If no precision is specified,
    the mantissa is stored exactly."""
    if not prec:
        if n in int_cache:
            return int_cache[n]
    return from_man_exp(MPZ(n), 0, prec, rnd)

def to_man_exp(s, signed=True):
    """Return (man, exp) of a raw mpf. Raise an error if inf/nan."""
    sign, man, exp, bc = s
    if (not man) and exp:
        raise ValueError("mantissa and exponent are defined "
                         "for finite numbers only")
    if signed and sign:
        man = -man
    return man, exp

def to_int(s, rnd=round_down):
    """Convert a raw mpf to the nearest int. Rounding is done down by
    default (same as int(float) in Python), but can be changed. If the
    input is inf/nan, an exception is raised."""
    sign, man, exp, bc = s
    if (not man) and exp:
        if s == fnan:
            raise ValueError("cannot convert nan to int")
        raise OverflowError("cannot convert infinity to int")
    if exp >= 0:
        if sign:
            return (-man) << exp
        return man << exp
    # Make default rounding fast
    if rnd == round_down:
        if sign:
            return -(man >> (-exp))
        else:
            return man >> (-exp)
    if sign:
        return round_int(-man, -exp, rnd)
    else:
        return round_int(man, -exp, rnd)

def mpf_round_int(s, rnd):
    sign, man, exp, bc = s
    if (not man) and exp:
        return s
    if exp >= 0:
        return s
    mag = exp+bc
    if mag < 1:
        if rnd == round_ceiling:
            if sign: return fzero
            else:    return fone
        elif rnd == round_floor:
            if sign: return fnone
            else:    return fzero
        elif rnd == round_nearest:
            if mag < 0 or man == MPZ_ONE: return fzero
            elif sign: return fnone
            else:      return fone
        else:
            raise NotImplementedError
    return mpf_pos(s, min(bc, mag), rnd)

def mpf_floor(s, prec=0, rnd=round_down):
    v = mpf_round_int(s, round_floor)
    if prec:
        v = mpf_pos(v, prec, rnd)
    return v

def mpf_ceil(s, prec=0, rnd=round_down):
    v = mpf_round_int(s, round_ceiling)
    if prec:
        v = mpf_pos(v, prec, rnd)
    return v

def mpf_nint(s, prec=0, rnd=round_down):
    v = mpf_round_int(s, round_nearest)
    if prec:
        v = mpf_pos(v, prec, rnd)
    return v

def mpf_frac(s, prec=0, rnd=round_down):
    return mpf_sub(s, mpf_floor(s), prec, rnd)

def from_float(x, prec=53, rnd=round_down):
    """Create a raw mpf from a Python float, rounding if necessary.
    If prec >= 53, the result is guaranteed to represent exactly the
    same number as the input. If prec is not specified, use prec=53."""
    # frexp only raises an exception for nan on some platforms
    if x != x: return fnan
    if x == math_float_inf: return finf
    if x == -math_float_inf: return fninf
    m, e = math.frexp(x)
    return from_man_exp(MPZ(m*(1<<53)), e-53, prec, rnd)

def from_npfloat(x, prec=113, rnd=round_down):
    """Create a raw mpf from a numpy float, rounding if necessary.
    If prec >= 113, the result is guaranteed to represent exactly the
    same number as the input. If prec is not specified, use prec=113."""
    y = float(x)
    if x == y: # ldexp overflows for float16
        return from_float(y, prec, rnd)
    import numpy as np
    if np.isfinite(x):
        m, e = np.frexp(x)
        return from_man_exp(MPZ(np.ldexp(m, 113)), int(e)-113, prec, rnd)
    return fnan

def from_Decimal(x, prec=0, rnd=round_down):
    """Create a raw mpf from a Decimal, rounding if necessary.
    If prec is not specified, use the equivalent bit precision
    of the number of significant digits in x."""
    if x.is_nan(): return fnan
    if x.is_infinite(): return fninf if x.is_signed() else finf
    if not prec:
        prec = int(len(x.as_tuple()[1])*blog2_10)
    return from_str(str(x), prec, rnd)

def to_float(s, strict=False, rnd=round_down):
    """
    Convert a raw mpf to a Python float.  The result is exact
    if s.bit_length() <= sys.float_info.mant_dig and no
    underflow/overflow occurs.  Else result is correctly rounded.

    If the magnitude of rounded number is too large to represent as
    a regular float, it will be converted to infinity.  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_nan

    exp2 = exp + bc
    # The smallest normal number is 2^(-1022)=0.1p-1021, and the smallest
    # subnormal is 2^(-1074)=0.1p-1073
    if exp2 <= float_min_subnormal_exp:
        if sign:
            if rnd == round_floor or (rnd == round_nearest
                                      and mpf_cmp(s, (1, MPZ(1), float_min_subnormal_exp
                                                      - 1, 1)) < 0):
                return -float_min * float_eps
            return 0.0
        if rnd == round_ceiling or (rnd == round_nearest
                                    and mpf_cmp(s, (0, MPZ(1), float_min_subnormal_exp
                                                    - 1, 1)) > 0):
            return float_min * float_eps
        return 0.0

    # The largest normal number is 2^1024*(1-2^(-53))=0.111...111p1024
    if exp2 > float_max_exp:
        if sign:
            if rnd == round_down or rnd == round_ceiling:
                return -float_max
            if strict:
                raise OverflowError("math range error")
            return -math_float_inf
        if rnd == round_down or rnd == round_floor:
            return float_max
        if strict:
            raise OverflowError("math range error")
        return math_float_inf

    nbits = float_mant_dig
    if exp2 < float_min_exp:
        # In the subnormal case, compute the exact number of significant bits.
        nbits += exp2 - float_min_exp
        assert 1 <= nbits < float_mant_dig
    if bc > nbits:
        sign, man, exp, bc = normalize(sign, man, exp, bc, nbits, rnd)
    if sign:
        man = -man
    # Should be exact:
    return math.ldexp(man, exp)

def from_rational(p, q, prec, rnd=round_down):
    """Create a raw mpf from a rational number p/q, round if
    necessary."""
    return mpf_div(from_int(p), from_int(q), prec, rnd)

def to_rational(s):
    """Convert a raw mpf to a rational number. Return integers (p, q)
    such that s = p/q exactly."""
    if s == fnan:
        raise ValueError("cannot convert nan to a rational number")
    if s in (finf, fninf):
        raise OverflowError("cannot convert infinity to a rational number")
    sign, man, exp, bc = s
    if sign:
        man = -man
    if exp >= 0:
        return man * (1<<exp), MPZ(1)
    else:
        return man, MPZ(1)<<(-exp)

def to_fixed(s, prec):
    """Convert a raw mpf to a fixed-point big integer"""
    sign, man, exp, bc = s
    offset = exp + prec
    if sign:
        if offset >= 0: return (-man) << offset
        else:           return (-man) >> (-offset)
    else:
        if offset >= 0: return man << offset
        else:           return man >> (-offset)


##############################################################################
##############################################################################

#----------------------------------------------------------------------------#
#                       Arithmetic operations, etc.                          #
#----------------------------------------------------------------------------#

def mpf_rand(prec):
    """Return a raw mpf chosen randomly from [0, 1), with prec bits
    in the mantissa."""
    return from_man_exp(MPZ(random.getrandbits(prec)), -prec, prec, round_floor)

def mpf_eq(s, t):
    """Test equality of two raw mpfs. This is simply tuple comparison
    unless either number is nan, in which case the result is False."""
    if not s[1] or not t[1]:
        if s == fnan or t == fnan:
            return False
    return s == t

def mpf_hash(s):
    # Duplicate the new hash algorithm, introduced in Python 3.2.
    ssign, sman, sexp, sbc = s

    # Handle special numbers
    if not sman:
        if s == fnan: return object.__hash__(s)
        if s == finf: return sys.hash_info.inf
        if s == fninf: return -sys.hash_info.inf

    hash_modulus = sys.hash_info.modulus
    hash_bits = 31 if sys.hash_info.width == 32 else 61
    h = sman % hash_modulus
    if sexp >= 0:
        sexp = sexp % hash_bits
    else:
        sexp = hash_bits - 1 - ((-1 - sexp) % hash_bits)
    h = (h << sexp) % hash_modulus
    if ssign: h = -h
    if h == -1: h = -2
    return int(h)

def mpf_cmp(s, t):
    """Compare the raw mpfs s and t. Return -1 if s < t, 0 if s == t,
    and 1 if s > t. (Same convention as Python's cmp() function.)"""

    # In principle, a comparison amounts to determining the sign of s-t.
    # A full subtraction is relatively slow, however, so we first try to
    # look at the components.
    ssign, sman, sexp, sbc = s
    tsign, tman, texp, tbc = t

    # Handle zeros and special numbers
    if not sman or not tman:
        if s == fzero: return -mpf_sign(t)
        if t == fzero: return mpf_sign(s)
        if s == t: return 0
        # Follow same convention as Python's cmp for float nan
        if t == fnan: return 1
        if s == finf: return 1
        if t == fninf: return 1
        return -1
    # Different sides of zero
    if ssign != tsign:
        if not ssign: return 1
        return -1
    # This reduces to direct integer comparison
    if sexp == texp:
        if sman == tman:
            return 0
        if sman > tman:
            if ssign: return -1
            else:     return 1
        else:
            if ssign: return 1
            else:     return -1
    # Check position of the highest set bit in each number. If
    # different, there is certainly an inequality.
    a = sbc + sexp
    b = tbc + texp
    if ssign:
        if a < b: return 1
        if a > b: return -1
    else:
        if a < b: return -1
        if a > b: return 1

    # Both numbers have the same highest bit. Subtract to find
    # how the lower bits compare.
    delta = mpf_sub(s, t, 5, round_floor)
    if delta[0]:
        return -1
    return 1

def mpf_lt(s, t):
    if s == fnan or t == fnan:
        return False
    return mpf_cmp(s, t) < 0

def mpf_le(s, t):
    if s == fnan or t == fnan:
        return False
    return mpf_cmp(s, t) <= 0

def mpf_gt(s, t):
    if s == fnan or t == fnan:
        return False
    return mpf_cmp(s, t) > 0

def mpf_ge(s, t):
    if s == fnan or t == fnan:
        return False
    return mpf_cmp(s, t) >= 0

def mpf_min_max(seq):
    min = max = seq[0]
    for x in seq[1:]:
        if mpf_lt(x, min): min = x
        if mpf_gt(x, max): max = x
    return min, max

def mpf_pos(s, prec=0, rnd=round_down):
    """Calculate 0+s for a raw mpf (i.e., just round s to the specified
    precision)."""
    if prec:
        sign, man, exp, bc = s
        if (not man) and exp:
            return s
        return normalize(sign, man, exp, bc, prec, rnd)
    return s

def mpf_neg(s, prec=0, rnd=round_down):
    """Negate a raw mpf (return -s), rounding the result to the
    specified precision. The prec argument can be omitted to do the
    operation exactly."""
    sign, man, exp, bc = s
    if not man:
        if exp:
            if s == finf: return fninf
            if s == fninf: return finf
        return s
    if not prec:
        return (1-sign, man, exp, bc)
    return normalize(1-sign, man, exp, bc, prec, rnd)

def mpf_abs(s, prec=0, rnd=round_down):
    """Return abs(s) of the raw mpf s, rounded to the specified
    precision. The prec argument can be omitted to generate an
    exact result."""
    sign, man, exp, bc = s
    if (not man) and exp:
        if s == fninf:
            return finf
        return s
    if not prec:
        if sign:
            return (0, man, exp, bc)
        return s
    return normalize(0, man, exp, bc, prec, rnd)

def mpf_sign(s):
    """Return -1, 0, or 1 (as a Python int, not a raw mpf) depending on
    whether s is negative, zero, or positive. (Nan is taken to give 0.)"""
    sign, man, exp, bc = s
    if not man:
        if s == finf: return 1
        if s == fninf: return -1
        return 0
    return (-1) ** sign

def mpf_add(s, t, prec=0, rnd=round_down, _sub=0):
    """
    Add the two raw mpf values s and t.

    With prec=0, no rounding is performed. Note that this can
    produce a very large mantissa (potentially too large to fit
    in memory) if exponents are far apart.
    """
    ssign, sman, sexp, sbc = s
    tsign, tman, texp, tbc = t
    tsign ^= _sub
    # Standard case: two nonzero, regular numbers
    if sman and tman:
        offset = sexp - texp
        if offset:
            if offset > 0:
                # Outside precision range; only need to perturb
                if offset > 100 and prec:
                    delta = sbc + sexp - tbc - texp
                    if delta > prec + 4:
                        offset = prec + 4
                        sman <<= offset
                        if tsign == ssign: sman += 1
                        else:              sman -= 1
                        return normalize(ssign, sman, sexp-offset,
                                         sman.bit_length(), prec, rnd)
                # Add
                if ssign == tsign:
                    man = tman + (sman << offset)
                # Subtract
                else:
                    if ssign: man = tman - (sman << offset)
                    else:     man = (sman << offset) - tman
                    if man >= 0:
                        ssign = 0
                    else:
                        man = -man
                        ssign = 1
                bc = man.bit_length()
                return normalize(ssign, man, texp, bc, prec or bc, rnd)
            elif offset < 0:
                # Outside precision range; only need to perturb
                if offset < -100 and prec:
                    delta = tbc + texp - sbc - sexp
                    if delta > prec + 4:
                        offset = prec + 4
                        tman <<= offset
                        if ssign == tsign: tman += 1
                        else:              tman -= 1
                        return normalize(tsign, tman, texp-offset,
                            tman.bit_length(), prec, rnd)
                # Add
                if ssign == tsign:
                    man = sman + (tman << -offset)
                # Subtract
                else:
                    if tsign: man = sman - (tman << -offset)
                    else:     man = (tman << -offset) - sman
                    if man >= 0:
                        ssign = 0
                    else:
                        man = -man
                        ssign = 1
                bc = man.bit_length()
                return normalize(ssign, man, sexp, bc, prec or bc, rnd)
        # Equal exponents; no shifting necessary
        if ssign == tsign:
            man = tman + sman
        else:
            if ssign: man = tman - sman
            else:     man = sman - tman
            if man >= 0:
                ssign = 0
            else:
                man = -man
                ssign = 1
        bc = man.bit_length()
        return normalize(ssign, man, texp, bc, prec or bc, rnd)
    # Handle zeros and special numbers
    if _sub:
        t = mpf_neg(t)
    if not sman:
        if sexp:
            if s == t or tman or not texp:
                return s
            return fnan
        if tman:
            return normalize(tsign, tman, texp, tbc, prec or tbc, rnd)
        return t
    if texp:
        return t
    if sman:
        return normalize(ssign, sman, sexp, sbc, prec or sbc, rnd)
    return s

def mpf_sub(s, t, prec=0, rnd=round_down):
    """Return the difference of two raw mpfs, s-t. This function is
    simply a wrapper of mpf_add that changes the sign of t."""
    return mpf_add(s, t, prec, rnd, 1)

def mpf_sum(xs, prec=0, rnd=round_down, absolute=False):
    """
    Sum a list of mpf values efficiently and accurately
    (typically no temporary roundoff occurs). If prec=0,
    the final result will not be rounded either.

    There may be roundoff error or cancellation if extremely
    large exponent differences occur.

    With absolute=True, sums the absolute values.
    """
    man = MPZ(0)
    exp = 0
    max_extra_prec = prec*2 or 1000000  # XXX
    special = None
    for x in xs:
        xsign, xman, xexp, xbc = x
        if xman:
            if xsign and not absolute:
                xman = -xman
            delta = xexp - exp
            if xexp >= exp:
                # x much larger than existing sum?
                # first: quick test
                if (delta > max_extra_prec) and \
                    ((not man) or delta-man.bit_length() > max_extra_prec):
                    man = xman
                    exp = xexp
                else:
                    man += (xman << delta)
            else:
                delta = -delta
                # x much smaller than existing sum?
                if delta-xbc > max_extra_prec:
                    if not man:
                        man, exp = xman, xexp
                else:
                    man = (man << delta) + xman
                    exp = xexp
        elif xexp:
            if absolute:
                x = mpf_abs(x)
            special = mpf_add(special or fzero, x, 1)
    # Will be inf or nan
    if special:
        return special
    return from_man_exp(man, exp, prec, rnd)

def mpf_mul(s, t, prec=0, rnd=round_down):
    """Multiply two raw mpfs"""
    ssign, sman, sexp, sbc = s
    tsign, tman, texp, tbc = t
    sign = ssign ^ tsign
    man = sman*tman
    if man:
        bc = man.bit_length()
        if prec:
            return normalize(sign, man, sexp+texp, bc, prec, rnd)
        else:
            return (sign, man, sexp+texp, bc)
    s_special = (not sman) and sexp
    t_special = (not tman) and texp
    if not s_special and not t_special:
        return fzero
    if fnan in (s, t): return fnan
    if (not tman) and texp: s, t = t, s
    if t == fzero: return fnan
    return {1:finf, -1:fninf}[mpf_sign(s) * mpf_sign(t)]

def gmpy_mpf_mul_int(s, n, prec, rnd=round_down):
    """Multiply by a Python integer."""
    sign, man, exp, bc = s
    if not man:
        return mpf_mul(s, from_int(n), prec, rnd)
    if not n:
        return fzero
    if n < 0:
        sign ^= 1
        n = -n
    man *= n
    return normalize(sign, man, exp, man.bit_length(), prec, rnd)

def python_mpf_mul_int(s, n, prec, rnd=round_down):
    """Multiply by a Python integer."""
    sign, man, exp, bc = s
    if not man:
        return mpf_mul(s, from_int(n), prec, rnd)
    if not n:
        return fzero
    if n < 0:
        sign ^= 1
        n = -n
    man *= n
    # Generally n will be small
    if n < 1024:
        bc += bctable[n] - 1
    else:
        bc += n.bit_length() - 1
    bc += man>>bc
    return normalize(sign, man, exp, bc, prec, rnd)

mpf_mul_int = python_mpf_mul_int

if gmpy:
    mpf_mul_int = gmpy_mpf_mul_int

def mpf_shift(s, n):
    """Quickly multiply the raw mpf s by 2**n without rounding."""
    sign, man, exp, bc = s
    if not man:
        return s
    return sign, man, exp+n, bc

def mpf_frexp(x):
    """Convert x = y*2**n to (y, n) with abs(y) in [0.5, 1) if nonzero"""
    sign, man, exp, bc = x
    if not man:
        return (x, 0)
    return mpf_shift(x, -bc-exp), bc+exp

def mpf_div(s, t, prec, rnd=round_down):
    """Floating-point division"""
    ssign, sman, sexp, sbc = s
    tsign, tman, texp, tbc = t
    if not sman or not tman:
        if s == fzero:
            if t == fzero: raise ZeroDivisionError
            if t == fnan: return fnan
            return fzero
        if t == fzero:
            raise ZeroDivisionError
        s_special = (not sman) and sexp
        t_special = (not tman) and texp
        if s_special and t_special:
            return fnan
        if s == fnan or t == fnan:
            return fnan
        if not t_special:
            if t == fzero:
                return fnan
            return {1:finf, -1:fninf}[mpf_sign(s) * mpf_sign(t)]
        return fzero
    sign = ssign ^ tsign
    if tman == 1:
        return normalize(sign, sman, sexp-texp, sbc, prec, rnd)
    # Same strategy as for addition: if there is a remainder, perturb
    # the result a few bits outside the precision range before rounding
    if not prec:
        extra = max(sbc, tbc) - sbc + tbc + 5
    else:
        extra = prec - sbc + tbc + 5
    if extra < 5:
        extra = 5
    quot, rem = divmod(sman<<extra, tman)
    if rem:
        quot = (quot<<1) + 1
        extra += 1
    bc = quot.bit_length()
    return normalize(sign, quot, sexp-texp-extra, bc, prec or bc, rnd)

def mpf_rdiv_int(n, t, prec, rnd=round_down):
    """Floating-point division n/t with a Python integer as numerator"""
    sign, man, exp, bc = t
    if not n or not man:
        return mpf_div(from_int(n), t, prec, rnd)
    if n < 0:
        sign ^= 1
        n = -n
    extra = prec + bc + 5
    quot, rem = divmod(n<<extra, man)
    if rem:
        quot = (quot<<1) + 1
        extra += 1
        return normalize(sign, quot, -exp-extra, quot.bit_length(), prec, rnd)
    return normalize(sign, quot, -exp-extra, quot.bit_length(), prec, rnd)

def mpf_mod(s, t, prec, rnd=round_down):
    ssign, sman, sexp, sbc = s
    tsign, tman, texp, tbc = t
    if ((not sman) and sexp) or ((not tman) and texp):
        if t == finf or t == fninf:
            return s
        return fnan
    # Important special case: do nothing if t is larger
    if ssign == tsign and texp > sexp+sbc:
        return s
    # Another important special case: this allows us to do e.g. x % 1.0
    # to find the fractional part of x, and it will work when x is huge.
    if tman == 1 and sexp > texp+tbc:
        return fzero
    base = min(sexp, texp)
    sman = (-1)**ssign * sman
    tman = (-1)**tsign * tman
    man = (sman << (sexp-base)) % (tman << (texp-base))
    if man >= 0:
        sign = 0
    else:
        man = -man
        sign = 1
    return normalize(sign, man, base, man.bit_length(), prec, rnd)

reciprocal_rnd = {
  round_down : round_up,
  round_up : round_down,
  round_floor : round_ceiling,
  round_ceiling : round_floor,
  round_nearest : round_nearest
}

negative_rnd = {
  round_down : round_down,
  round_up : round_up,
  round_floor : round_ceiling,
  round_ceiling : round_floor,
  round_nearest : round_nearest
}

def mpf_pow_int(s, n, prec, rnd=round_down):
    """Compute s**n, where s is a raw mpf and n is a Python integer."""
    sign, man, exp, bc = s

    if (not man) and exp:
        if s == finf:
            if n > 0: return s
            if n == 0: return fone
            return fzero
        if s == fninf:
            if n > 0: return [finf, fninf][n & 1]
            if n == 0: return fone
            return fzero
        if n == 0:
            return fone
        return fnan

    n = int(n)
    if n == 0: return fone
    if n == 1: return mpf_pos(s, prec, rnd)
    if n == 2:
        _, man, exp, bc = s
        if not man:
            return fzero
        man = man*man
        if man == 1:
            return (0, MPZ_ONE, exp+exp, 1)
        bc = bc + bc - 2
        bc += bctable[man>>bc]
        return normalize(0, man, exp+exp, bc, prec, rnd)
    if n == -1: return mpf_div(fone, s, prec, rnd)
    if n < 0:
        inverse = mpf_pow_int(s, -n, prec+5, reciprocal_rnd[rnd])
        return mpf_div(fone, inverse, prec, rnd)

    result_sign = sign & n

    # Use exact integer power when the exact mantissa is small
    if man == 1:
        return (result_sign, MPZ_ONE, exp*n, 1)
    if bc*n < 1000:
        man **= n
        return normalize(result_sign, man, exp*n, man.bit_length(), prec, rnd)

    # Use directed rounding all the way through to maintain rigorous
    # bounds for interval arithmetic
    rounds_down = (rnd == round_nearest) or \
        shifts_down[rnd][result_sign]

    # Now we perform binary exponentiation. Need to estimate precision
    # to avoid rounding errors from temporary operations. Roughly log_2(n)
    # operations are performed.
    workprec = prec + 4*n.bit_length() + 4
    _, pm, pe, pbc = fone
    while 1:
        if n & 1:
            pm = pm*man
            pe = pe+exp
            pbc += bc - 2
            pbc = pbc + bctable[pm >> pbc]
            if pbc > workprec:
                if rounds_down:
                    pm = pm >> (pbc-workprec)
                else:
                    pm = -((-pm) >> (pbc-workprec))
                pe += pbc - workprec
                pbc = workprec
            n -= 1
            if not n:
                break
        man = man*man
        exp = exp+exp
        bc = bc + bc - 2
        bc = bc + bctable[man >> bc]
        if bc > workprec:
            if rounds_down:
                man = man >> (bc-workprec)
            else:
                man = -((-man) >> (bc-workprec))
            exp += bc - workprec
            bc = workprec
        n = n // 2

    return normalize(result_sign, pm, pe, pbc, prec, rnd)


def mpf_perturb(x, eps_sign, prec, rnd):
    """
    For nonzero x, calculate x + eps with directed rounding, where
    eps < prec relatively and eps has the given sign (0 for
    positive, 1 for negative).

    With rounding to nearest, this is taken to simply normalize
    x to the given precision.
    """
    if rnd == round_nearest:
        return mpf_pos(x, prec, rnd)
    sign, man, exp, bc = x
    eps = (eps_sign, MPZ_ONE, exp+bc-prec-1, 1)
    if sign:
        away = (rnd in (round_down, round_ceiling)) ^ eps_sign
    else:
        away = (rnd in (round_up, round_ceiling)) ^ eps_sign
    if away:
        return mpf_add(x, eps, prec, rnd)
    else:
        return mpf_pos(x, prec, rnd)


#----------------------------------------------------------------------------#
#                              Radix conversion                              #
#----------------------------------------------------------------------------#

def to_digits_exp(s, dps, base=10):
    """Helper function for representing the floating-point number s as
    a string with dps digits. Returns (sign, string, exponent) where
    sign is '' or '-', string is the digit string in the given base,
    and exponent is the exponent as an int.

    If inexact, the string representation is rounded toward zero."""

    # Extract sign first so it doesn't mess up the string digit count
    if s[0]:
        sign = '-'
        s = mpf_neg(s)
    else:
        sign = ''
    _sign, man, exp, bc = s

    if not man:
        return '', '0'*int(dps), 0

    if base == 10:
        blog2 = blog2_10
    elif pow(2, blog2 := int(math.log2(base))) == base:
        pass
    else:
        raise NotImplementedError

    bitprec = int(dps * blog2) + 10

    # Cut down to size
    # TODO: account for precision when doing this
    exp_from_1 = exp + bc
    if base == 10 and abs(exp_from_1) > 3500:
        from .libelefun import mpf_ln2, mpf_ln10

        # Set b = int(exp * log(2)/log(10))
        # If exp is huge, we must use high-precision arithmetic to
        # find the nearest power of ten
        expprec = exp.bit_length() + 5
        tmp = from_int(exp)
        tmp = mpf_mul(tmp, mpf_ln2(expprec))
        tmp = mpf_div(tmp, mpf_ln10(expprec), expprec)
        b = to_int(tmp)
        s = mpf_div(s, mpf_pow_int(ften, b, bitprec), bitprec)
        _sign, man, exp, bc = s
        exponent = b
    else:
        exponent = 0

    # First, calculate mantissa digits by converting to a binary
    # fixed-point number and then converting that number to
    # a decimal fixed-point number.
    fixprec = max(bitprec - exp - bc, 0)
    fixdps = int(fixprec / blog2 + 0.5)
    sf = to_fixed(s, fixprec)
    sb = bin_to_radix(sf, fixprec, base, fixdps)
    digits = numeral(sb, base=base, size=dps)

    exponent += len(digits) - fixdps - 1
    return sign, digits, exponent

def round_digits(sign, digits, dps, base, rnd=round_down, fixed=False):
    """
    Returns the rounded digits, and the number of places the decimal point was
    shifted.
    """

    assert len(digits) > dps
    assert rnd in (round_nearest, round_up, round_down, round_ceiling,
                   round_floor)

    if rnd == round_ceiling:
        rnd = round_down if sign else round_up
    elif rnd == round_floor:
        rnd = round_up if sign else round_down

    exponent = 0

    if rnd == round_down:
        return digits[:dps], 0
    elif rnd == round_nearest:
        rnd_digs = stddigits[(base//2 + base % 2):base]
    else:
        rnd_digs = stddigits[1:base]

    tie_down = False
    tie_up = False

    if rnd == round_nearest:
        # The first digit after dps is a 5 and we should determine whether we
        # round it up or down.
        if digits[dps] == rnd_digs[0]:
            tie_down = True

            # If the digit we round to is even, we may round down if all the
            # following digits are 0.
            for i in range(dps+1, len(digits)):
                if digits[i] != '0':
                    tie_down = False
                    break
            # If the digit we round to is odd, we round up no matter what.
            if digits[dps-1] in stddigits[1:base:2]:
                tie_down = False

    elif rnd == round_up:
        # If any digit following a 0 is different from zero, we round up.
        if digits[dps] == '0':
            for i in range(dps+1, len(digits)):
                if digits[i] != '0':
                    tie_up = True
                    break

    # Add or subtract a unit to the digit following the one we round to.
    if tie_down:
        digits = digits[:dps] + stddigits[int(digits[dps], base) - 1]
    elif tie_up:
        digits = digits[:dps] + '1'

    # Rounding up kills some instances of "...99999"
    if digits[dps] in rnd_digs:
        digits = digits[:dps]
        i = dps - 1
        dig = stddigits[base-1]
        while i >= 0 and digits[i] == dig:
            i -= 1
        if i >= 0:
            digits = digits[:i] + stddigits[int(digits[i], base) + 1] + \
                '0' * (dps - i - 1)
        else:
            # When rounding up 0.9999... in fixed format, we lose one dps.
            digits = '1' + '0' * (dps - (0 if fixed else 1))
            exponent += 1
    else:
        digits = digits[:dps]

    return digits, exponent


def to_str(s, dps, strip_zeros=True, min_fixed=None, max_fixed=None,
           show_zero_exponent=False, base=10, binary_exp=False,
           rnd=round_nearest):
    """
    Convert a raw mpf to a floating-point literal in the given base
    with at most `dps` digits in the mantissa (not counting extra zeros
    that may be inserted for visual purposes).

    The number will be printed in fixed-point format if the position
    of the leading digit is strictly between min_fixed
    (default = min(-dps/3,-5)) and max_fixed (default = dps).

    To force fixed-point format always, set min_fixed = -inf,
    max_fixed = +inf. To force floating-point format, set
    min_fixed >= max_fixed.

    If binary_exp is True and the base is either 2 or 16, the number will
    be printed in a binary or hexadecimal notation, where the exponent
    separator is the 'p' and the exponent is written in decimal rather than
    hexadecimal or binary.  The number is normalized, i.e. the first
    digit is 1.  This is format of the float.fromhex().

    The literal is formatted so that it can be parsed back to a number
    by from_str, float(), float.fromhex() or Decimal().
    """
    sep = '@' if base > 10 else 'e'

    if binary_exp:
        sep = 'p'
        if base not in (2, 16):
            raise ValueError("binary_exp option could be used for base 2 and 16")

    if rnd not in (round_nearest, round_floor, round_ceiling, round_up,
                   round_down):
        raise ValueError("rnd should be one of " +
                         ", ".join([round_nearest, round_floor, round_ceiling,
                                   round_up, round_down]) + ".")

    if base == 2:
        prefix = "0b"
    elif base == 8:
        prefix = "0o"
    elif base == 16:
        prefix = "0x"
    else:
        prefix = ""

    # Special numbers
    if not s[1]:
        if s == fzero:
            if dps: t = '0.0'
            else:   t = '.0'
            if show_zero_exponent:
                t += sep + '+0'
            return prefix + t
        if s == finf: return 'inf'
        if s == fninf: return '-inf'
        if s == fnan: return 'nan'
        raise ValueError

    if min_fixed is None: min_fixed = min(-(dps//3), -5)
    if max_fixed is None: max_fixed = dps

    # to_digits_exp rounds to floor.
    # This sometimes kills some instances of "...00001"
    sign, digits, exponent = to_digits_exp(s, dps+10, base)

    rnd_digs = stddigits[(base//2 + base%2):base]

    # No digits: show only .0; round exponent to nearest
    if not dps:
        if digits[0] in rnd_digs:
            exponent += 1
        digits = ".0"

    else:
        if binary_exp and base == 16:
            exponent *= 4
            # normalization
            if int(digits[0], 16) > 1:
                shift = math.floor(math.log2(int(digits[0], 16)))
                exponent += shift
                n = int(digits, 16) >> shift
                digits = hex(n)[2:]

        digits, exp_add = round_digits(s[0], digits, dps, base, rnd)
        exponent += exp_add

        # Prettify numbers close to unit magnitude
        if not binary_exp and min_fixed < exponent < max_fixed:
            if exponent < 0:
                digits = ("0"*(-exponent)) + digits
                split = 1
            else:
                split = exponent + 1
                if split > dps:
                    digits += "0"*(split-dps)
            exponent = 0
        else:
            split = 1

        digits = (digits[:split] + "." + digits[split:])

        if strip_zeros:
            # Clean up trailing zeros
            digits = digits.rstrip('0')
            if digits[-1] == ".":
                digits += "0"

    sign += prefix

    if exponent == 0 and dps and not show_zero_exponent: return sign + digits
    return sign + digits + sep + f"{exponent:+}"

def str_to_man_exp(x, base=10):
    """Helper function for from_str."""
    x = x.lower().rstrip('l').replace('_', '')
    # Split into mantissa, exponent
    if base <= 10:
        sep = 'e'
    else:
        sep = '@'
    if pow(2, e2 := int(math.log2(base))) == base and e2 in [1, 4] and x.find('p') >= 0:
        sep = 'p'
    parts = x.split(sep)
    if len(parts) == 1:
        exp = 0
    elif len(parts) == 2:
        x = parts[0]
        exp = int(parts[1])
    else:
        raise ValueError("couldn't convert a str to mpf")
    # Look for radix point in mantissa
    parts = x.split('.')
    if len(parts) == 2:
        a, b = parts[0], parts[1].rstrip('0')
        if sep != 'p':
            exp -= len(b)
        else:
            exp -= len(b)*e2
        if a == '':
            a = '0'
        x = a + b
    int_max_str_digits = 0
    if BACKEND == 'python' and hasattr(sys, 'get_int_max_str_digits'):
        int_max_str_digits = sys.get_int_max_str_digits()
        sys.set_int_max_str_digits(0)
    x = MPZ(x, base)
    if int_max_str_digits:
        sys.set_int_max_str_digits(int_max_str_digits)
    return x, exp

special_str = {'inf':finf, '+inf':finf, '-inf':fninf, 'nan':fnan,
               'oo':finf, '+oo':finf, '-oo':fninf}

def from_str(x, prec=0, rnd=round_down, base=0):
    """Create a raw mpf from a string x in a given base, rounding in the
    specified direction if the input number cannot be represented
    exactly as a binary floating-point number with the given number of
    bits.  The string syntax accepted for float() or float.fromhex()
    is accepted too.

    TODO: the rounding does not work properly for large exponents.
    """
    x = x.lower().strip()
    if x in special_str:
        return special_str[x]

    if not base:
        if x.startswith(('0b', '-0b', '0B', '-0B')):
            base = 2
        elif x.startswith(('0x', '-0x', '0X', '-0X')):
            base = 16
        elif x.startswith(('0o', '-0o')):
            base = 8
        else:
            base = 10

    if '/' in x:
        p, q = x.split('/')
        p, q = p.rstrip('l'), q.rstrip('l')
        return from_rational(int(p, base), int(q, base), prec, rnd)

    man, exp = str_to_man_exp(x, base)

    if base == 10:
        # XXX: appropriate cutoffs & track direction
        # note no factors of 5
        if abs(exp) > 400:
            s = from_int(man, prec+10)
            s = mpf_mul(s, mpf_pow_int(ften, exp, prec+10), prec, rnd)
        else:
            if exp >= 0:
                s = from_int(man * 10**exp, prec, rnd)
            else:
                s = from_rational(man, 10**-exp, prec, rnd)
    elif pow(2, e2 := int(math.log2(base))) == base:
        if x.find('p') < 0:
            s = from_man_exp(man, exp*e2, prec, rnd)
        else:
            s = from_man_exp(man, exp, prec, rnd)
    else:
        raise NotImplementedError
    return s


#----------------------------------------------------------------------------#
#                              String formatting                             #
#----------------------------------------------------------------------------#

_FLOAT_FORMAT_SPECIFICATION_MATCHER = re.compile(r"""
    (?:
        (?P<fill_char>.)?
        (?P<align>[<>=^])
    )?
    (?P<sign>[-+ ]?)
    (?P<no_neg_0>z)?
    (?P<alternate>\#)?
    (?P<zeropad>0(?=0*[1-9]))?
    (?P<width>[0-9]+)?
    (?P<thousands_separators>[,_])?
    (?:\.
        (?=[,_0-9])  # lookahead for digit or separator
        (?P<precision>[0-9]+)?
        (?P<frac_separators>[,_])?
    )?
    (?P<rounding>[UDYZN])?
    (?P<type>[aAbeEfFgG%])?
""", re.DOTALL | re.VERBOSE).fullmatch

_GMPY_ROUND_CHAR_DICT = {
        'U': round_ceiling,
        'D': round_floor,
        'Y': round_up,
        'Z': round_down,
        'N': round_nearest
        }

def calc_padding(nchars, width, align):
    '''
    Computes the left and right padding required to fill the required width,
    according to how the string will be aligned.
    '''
    ntotal = max(nchars, width)

    if align in ('>', '='):
        lpad = ntotal - nchars
        rpad = 0
    elif align == '^':
        lpad = (ntotal - nchars)//2
        rpad = ntotal - nchars - lpad
    else:
        lpad = 0
        rpad = ntotal - nchars

    return (lpad, rpad)


def read_format_spec(format_spec):
    '''
    Reads the format spec into a dictionary.
    This is more or less copied from the CPython implementation for regular
    floats.
    '''

    format_dict = {
        'fill_char': ' ',
        'align': '>',
        'sign': '-',
        'no_neg_0': False,
        'alternate': False,
        'thousands_separators': '',
        'frac_separators': '',
        'width': -1,
        'precision': -1,
        'type': ''
        }

    if match := _FLOAT_FORMAT_SPECIFICATION_MATCHER(format_spec):
        format_dict['fill_char'] = match['fill_char'] or format_dict['fill_char']
        format_dict['align'] = match['align'] or format_dict['align']
        format_dict['sign'] = match['sign'] or format_dict['sign']
        format_dict['no_neg_0'] = bool(match['no_neg_0']) or format_dict['no_neg_0']
        format_dict['alternate'] = bool(match['alternate']) or \
            format_dict['alternate']
        format_dict['thousands_separators'] = match['thousands_separators'] \
            or format_dict['thousands_separators']
        format_dict['width'] = int(match['width'] or format_dict['width'])
        format_dict['precision'] = int(match['precision'] or format_dict['precision'])
        format_dict['frac_separators'] = match['frac_separators'] \
            or format_dict['frac_separators']
        rounding_char = match['rounding']
        format_dict['type'] = match['type'] or format_dict['type']

        if rounding_char is not None:
            format_dict['rounding'] = _GMPY_ROUND_CHAR_DICT[rounding_char]

        if match['zeropad']:
            if not match['align']:
                format_dict['align'] = '='
            if not match['fill_char']:
                format_dict['fill_char'] = '0'

        if format_dict['precision'] < 0 and format_dict['type'].lower() not in ['', 'a', 'b']:
            format_dict['precision'] = 6
    else:
        raise ValueError("Invalid format specifier '{}'".format(format_spec))

    return format_dict


def format_fixed(s, dps, rnd=round_down):
    # First, get the exponent to know how many digits we will need
    base = 10
    _, _, exponent = to_digits_exp(s, 1, base)

    # Now that we have an estimate, compute the correct digits
    # (we do this because the previous computation could yield the wrong
    # exponent by +- 1)
    _, digits, exponent = to_digits_exp(
            s, max(dps+exponent+4, int(s[3]/blog2_10)), base)
    orig_dps = dps
    dps += exponent + 1

    # The number we want to print is lower in magnitude that the requested
    # precision. We should only print 0s.
    if dps < 0:
        int_part = '0'
        frac_part = orig_dps*'0'

    else:
        digits, exp_add = round_digits(s[0], digits, dps, base, rnd, True)
        exponent += exp_add

        # Here we prepend the corresponding 0s to the digits string, according
        # to the value of exponent
        if exponent < 0:
            digits = ("0"*(-exponent)) + digits
            split = 1
        else:
            split = exponent + 1
        int_part = digits[:split]

        # Finally, assemble the digits including the decimal point
        if orig_dps == 0:
            return int_part, ''

        frac_part = digits[split:]

    return int_part, frac_part


def format_scientific(s, dps, rnd=round_down):
    base = 10

    # First, get the exponent to know how many digits we will need
    dps += 1
    _, digits, exponent = to_digits_exp(s, max(dps + 10,
                                               int(s[3]/blog2_10) + 10),
                                        base)
    digits, exp_add = round_digits(s[0], digits, dps, base, rnd)
    exponent += exp_add

    return digits[0], digits[1:], f'e{exponent:+03d}'


def format_hexadecimal(s, dps, rnd=round_down):
    prec = 4*dps + 1 if dps >= 0 else s[1].bit_length()

    if s[1]:
        s = mpf_pos(s, prec, rnd)

        exponent = s[2] + s[3] - 1
        man = s[1] | (1 << s[3] + 2)  # set leading digit (ignored) to 0x9
        man <<= 1 + 4*((s[3] + 3)//4) - s[3]

        frac_digits = hex(man)[3:]
        digits = "1"
    else:
        exponent = 0
        frac_digits = ""
        digits = "0"

    if dps >= 0:
        frac_digits = frac_digits[:dps]
        frac_digits += "0"*(dps - len(frac_digits))
    else:
        # Clean up trailing zeros
        frac_digits = frac_digits.rstrip('0')

    return digits, frac_digits, f'p{exponent:+01d}'


def format_binary(s, dps, rnd=round_down):
    prec = dps + 1 if dps >= 0 else s[1].bit_length()
    s = mpf_pos(s, prec, rnd)

    digits = bin(s[1])[2:]
    digits = digits + '0'*(dps + 1 - len(digits))
    exponent = s[2]
    if s[1]:
        exponent += s[1].bit_length() - 1
    return digits[0], digits[1:], f'p{exponent:+01d}'


_MAP_SPEC_STR = {finf: 'inf', fninf: 'inf', fnan: 'nan'}


def fill_sep(digits, sep, prev, nmod, sep_range):
    return prev + sep.join(digits[pos:pos + sep_range]
                           for pos in range(nmod, len(digits), sep_range))


def format_digits(num, format_dict, prec, rnd, _pretty_repr_dps):
    capitalize = False
    if format_dict['type'] in list('AFGE'):
        capitalize = True

    fmt_type = format_dict['type'].lower()

    percent = False
    if fmt_type == '%':
        percent = True
        fmt_type = 'f'
        num = mpf_mul(num, from_int(100), prec, rnd=round_nearest)

    dps = format_dict['precision']

    int_part = ''
    exponent = ''
    sign = ''

    # Now the general case
    strip_last_zero = False
    strip_zeros = False

    rnd = format_dict.get('rounding', rnd)

    if not fmt_type or fmt_type == 'g':
        if not format_dict['alternate']:
            strip_zeros = True
            if fmt_type == 'g':
                strip_last_zero = True

        if dps < 0:
            dps = repr_dps(prec) if _pretty_repr_dps else prec_to_dps(prec)
        if dps == 0:
            dps = 1

        _, tdigits, exp = to_digits_exp(num, max(53/blog2_10, dps), 10)
        if num[1]:
            _, exp_add = round_digits(num, tdigits, dps, 10, rnd)
            exp += exp_add

        fix0 = 0 if fmt_type else 1
        if -4 <= exp < dps - fix0:
            dps = max(0, dps - exp - 1)
        else:
            fmt_type = 'e'
            dps = max(0, dps - 1)

    if num in _MAP_SPEC_STR:  # special cases
        frac_part = _MAP_SPEC_STR[num]
        if capitalize:
            frac_part = frac_part.upper()

    elif fmt_type == 'e':
        int_part, frac_part, exponent = format_scientific(num, dps, rnd=rnd)
        if strip_zeros:
            frac_part = frac_part.rstrip('0')
        if frac_part or format_dict['alternate']:
            frac_part = '.' + frac_part
        if capitalize:
            exponent = exponent.replace('e', 'E')

    elif fmt_type == 'a':
        int_part, frac_part, exponent = format_hexadecimal(num, dps, rnd=rnd)
        if capitalize:
            int_part = '0X' + int_part
            frac_part = frac_part.upper()
            exponent = exponent.replace('p', 'P')
        else:
            int_part = '0x' + int_part
        if frac_part or format_dict['alternate']:
            frac_part = '.' + frac_part

    elif fmt_type == 'b':
        int_part, frac_part, exponent = format_binary(num, dps, rnd=rnd)
        if frac_part or format_dict['alternate']:
            frac_part = '.' + frac_part

    else:  # fixed-point formats
        int_part, frac_part = format_fixed(num, dps, rnd=rnd)

        if strip_zeros:
            frac_part = frac_part.rstrip('0')
        if not frac_part and not fmt_type:
            frac_part = '0'
        if (frac_part or format_dict['alternate']
                or (dps and not strip_last_zero)):
            frac_part = '.' + frac_part

    sep_range = 3
    sep = format_dict['frac_separators']
    if sep and frac_part:
        frac_part = fill_sep(frac_part, sep, frac_part[0], 1, sep_range)
    digits = frac_part + exponent

    sign = '-' if num[0] else ''
    if sign != '-' and format_dict['sign'] != '-':
        sign = format_dict['sign']
    if fmt_type == 'f' and format_dict['no_neg_0']:
        if int_part == "0" and all(_ in ['0', '.', '_', ',']
                                   for _ in digits):
            if format_dict['sign'] == '-':
                sign = ''
            else:
                sign = format_dict['sign']

    if percent:
        digits += '%'

    sep = format_dict['thousands_separators']
    width = format_dict['width']
    min_leading = width - len(digits) - len(sign)
    if (int_part and fmt_type not in ['a', 'b']
            and format_dict['fill_char'] == '0' and format_dict['align'] == '='
            and min_leading > len(int_part)):
        int_part = int_part.zfill(sep_range*min_leading//(sep_range + 1) + 1
                                  if sep else min_leading)

    # Add the thousands separator every 3 characters.
    split = len(int_part)
    if sep and split > sep_range:
        # the first thousand separator may be located before 3 characters
        nmod = split % sep_range
        if nmod != 0:
            prev = int_part[:nmod] + sep
        else:
            prev = ''
        int_part = fill_sep(int_part, sep, prev, nmod, sep_range)

    return sign, int_part + digits


def format_mpf(num, format_spec, prec, rnd, _pretty_repr_dps):
    format_dict = read_format_spec(format_spec)
    sign, digits = format_digits(num, format_dict, prec, rnd, _pretty_repr_dps)
    nchars = len(digits) + len(sign)
    lpad, rpad = calc_padding(
            nchars, format_dict['width'], format_dict['align'])

    if format_dict['align'] == '=':
        return sign + lpad*format_dict['fill_char'] + digits + \
                rpad*format_dict['fill_char']

    return lpad*format_dict['fill_char'] + sign + digits \
            + rpad*format_dict['fill_char']


def format_mpc(num, format_spec, prec, rnd, _pretty_repr_dps):
    format_dict = read_format_spec(format_spec)

    if format_dict['fill_char'] == '0':
        raise ValueError("Zero padding is not allowed in complex format "
                         "specifier.")
    if format_dict['align'] == '=':
        raise ValueError("'=' alignment flag is not allowed in complex format "
                         "specifier.")
    if format_dict['type'] == '%':
        raise ValueError("'%' formatting type is not allowed in complex "
                         "format specifier.")

    fmt_type = format_dict['type'].lower()
    if not fmt_type:
        format_dict['type'] = 'g'
    sign_re, digits_re = format_digits(num[0], format_dict, prec, rnd, _pretty_repr_dps)
    fmt_sign = format_dict['sign']
    format_dict['sign'] = '+'
    sign_im, digits_im = format_digits(num[1], format_dict, prec, rnd, _pretty_repr_dps)
    digits_im += 'j'

    if not fmt_type:
        if num[0] == fzero:
            sign_re = ''
            digits_re = ''
            if sign_im == '+':
                sign_im = fmt_sign if fmt_sign in [' ', '+'] else ''
        else:
            sign_re = '(' + sign_re
            digits_im += ')'

    nchars = len(sign_re) + len(digits_re) + len(sign_im) + len(digits_im)

    lpad, rpad = calc_padding(nchars, format_dict['width'],
                              format_dict['align'])

    return (lpad*format_dict['fill_char'] + sign_re + digits_re + sign_im
            + digits_im + rpad*format_dict['fill_char'])


#----------------------------------------------------------------------------#
#                                Square roots                                #
#----------------------------------------------------------------------------#


def mpf_sqrt(s, prec, rnd=round_down):
    """
    Compute the square root of a nonnegative mpf value. The
    result is correctly rounded.
    """
    sign, man, exp, bc = s
    if sign:
        raise ComplexResult("square root of a negative number")
    if not man:
        return s
    if exp & 1:
        exp -= 1
        man <<= 1
        bc += 1
    elif man == 1:
        return normalize(sign, man, exp//2, bc, prec, rnd)
    shift = max(4, 2*prec-bc+4)
    shift += shift & 1
    if rnd in 'fd':
        man = isqrt(man<<shift)
    else:
        man, rem = sqrtrem(man<<shift)
        # Perturb up
        if rem:
            man = (man<<1)+1
            shift += 2
    return from_man_exp(man, (exp-shift)//2, prec, rnd)

def mpf_hypot(x, y, prec, rnd=round_down):
    """Compute the Euclidean norm sqrt(x**2 + y**2) of two raw mpfs
    x and y."""
    if y == fzero: return mpf_abs(x, prec, rnd)
    if x == fzero: return mpf_abs(y, prec, rnd)
    hypot2 = mpf_add(mpf_mul(x,x), mpf_mul(y,y), prec+10, rnd)
    return mpf_sqrt(hypot2, prec, rnd)