#!/usr/bin/env python3
"""kite_exact_parse.py -- the exact parser of the kite bundle's shipped rational connections (no sympy).

The connection JSONs of this directory (kite-connection-A14.json, kite-connection-A14-x3.json; the same strings
are embedded in kite-boundary.py) hold each entry A_ij(d, s) of the 14x14 connection d/ds M = A(d, s) M as a
string in this grammar, and nothing else:
    integer literals, the names s and d, the operators + - * / **, parentheses
(every one of the 392 entries of the two shipped files tokenises in these classes -- the census is in the
bundle's CHANGES.md; a token outside them is refused by name with its position).  Precedence and
associativity are Python's (the strings are Python expressions): ** binds tighter than unary minus, * and /
left to right, + and - left to right; ** takes a non-negative integer exponent.

What the module gives the three scripts (kite-evaluate.py, kite-boundary.py, kite-boundary-general.py):
  parse_entry(text)        -> (P, Q): the entry as an exact bivariate rational function P(d,s)/Q(d,s), the
                              polynomials as dicts {(deg_d, deg_s): Fraction} (fractions.Fraction throughout);
  eps_grade(P, Q, neps)    -> {k: (N_k, D_k)}: the eps-Taylor coefficients of the entry at d = 4 - 2*eps for
                              k = 0..neps as univariate rational functions of s (ascending Fraction lists),
                              by the closed forms c0 = P0/Q0, c1 = (P1 Q0 - P0 Q1)/Q0^2, c2 = (P2 Q0^2 - Q1 N1
                              - Q2 P0 Q0)/Q0^3, ... (exact; nothing is cancelled);
  load_graded(A, graded)   -> {(i, j, k): (pc, qc)}: the eps-graded coefficient lists the numerics consume,
                              READ from a graded JSON (kite-connection-A14-graded.json, -x3-graded.json: the
                              lists the earlier sympy pipeline of these scripts produced, frozen 2026-09-05)
                              and CHECKED entry by entry against the parsed connection -- pc * D_k == qc * N_k
                              as polynomials, exactly, and the two key sets equal -- refused by name otherwise;
  grade_canonical(A, neps) -> the same dict from the parser alone, every entry cancelled to lowest terms with an
                              integer denominator of positive leading coefficient (the form kite-boundary-general.py
                              uses for a connection that ships no graded file; NOT the served lines' form);
  singular_points(A)       -> (points, mixed): the real roots of the s-only irreducible factors of every
                              denominator, exact -- a rational, or a + b*sqrt(n) with n squarefree -- as objects
                              that print in the form the earlier sympy census printed ('7 - 4*sqrt(3)',
                              '4*sqrt(3) + 7', ...), convert to float and to mpf the way that census did, and
                              sort by value; python-flint's fmpz_poly.factor does the factorisation (the only
                              flint use of this module, imported inside this call); s = 0 is always included;
                              'mixed' lists any denominator factor that mixes d and s (none in the shipped files);
  rat(text), mpf_rat(fr)   -> the exact Fraction of a point string ('-7/3', '-4.25') and its mpf at the current
                              working precision by mpmath.libmp.from_rational (the conversion sympy's Rational
                              used), so every endpoint, seed point and rational reference entry is the same
                              binary number as before.

Why the graded lists are read from a file and checked, rather than taken from the parser's own reduced form:
the transport numerics (the mpf coefficient lists and the power-series division of every step) are
representation-sensitive at the last bit -- P/Q and (fP)/(fQ) give the same function but not the same rounding
noise -- and the two-precision and dps-doubling lines of the scripts print that noise as a digit count.  The
earlier pipeline's lists are not the reduced form (in each shipped connection five (i,j,k) entries carry an
uncancelled common factor and thirteen carry a denominator content 2), so byte-identical output at every
documented tier needs those exact lists: they ship once, sha256-pinned by the consuming scripts, and every
start re-derives the entry algebra from the strings and asserts the equality above.  A parser that returned a
wrong coefficient, or a graded list that disagreed with the strings, stops the run by name before any value.

Refusals raise ParseError (a ValueError) with the entry, the position and the token named; the consuming
scripts turn that into their exit code 3 (integrity) before any value is printed.
"""
import json
import re
from fractions import Fraction

__all__ = ["ParseError", "parse_entry", "eps_grade", "load_graded", "grade_canonical", "singular_points",
           "Root", "rat", "mpf_rat", "poly_str"]

NAMES = ("d", "s")          # the two indeterminates of the grammar, in the exponent order of the dict keys
D_SUBST = (Fraction(4), Fraction(-2))   # d = 4 - 2*eps: (constant, eps coefficient)


class ParseError(ValueError):
    """A string outside the shipped grammar, or a graded list that disagrees with the strings."""


# ---------------------------------------------------------------------------
# 1. Tokens and the recursive-descent parser (values are exact rational functions)
# ---------------------------------------------------------------------------
_TOKEN = re.compile(r"\s*(?:(\d+)|([A-Za-z_][A-Za-z_0-9]*)|(\*\*)|([-+*/()]))")


def tokenize(text):
    """The token list of one entry: ('INT', int) | ('NAME', 'd'|'s') | ('OP', '+' '-' '*' '/' '**' '(' ')').
    Anything else is refused by name with its 0-based position."""
    toks, pos, n = [], 0, len(text)
    while pos < n:
        if text[pos:].strip() == "":
            break
        m = _TOKEN.match(text, pos)
        if not m:
            raise ParseError(f"entry {text!r}: character {text[pos]!r} at position {pos} is outside the grammar "
                             f"(integers, the names d and s, + - * / **, parentheses)")
        if m.group(1) is not None:
            toks.append(("INT", int(m.group(1))))
        elif m.group(2) is not None:
            if m.group(2) not in NAMES:
                raise ParseError(f"entry {text!r}: name {m.group(2)!r} at position {m.start(2)} is outside the "
                                 f"grammar (the names are d and s)")
            toks.append(("NAME", m.group(2)))
        elif m.group(3) is not None:
            toks.append(("OP", "**"))
        else:
            toks.append(("OP", m.group(4)))
        pos = m.end()
    return toks


# bivariate polynomials: {(a, b): Fraction} = sum c * d^a * s^b  (zero polynomial = {})
def _padd(p, q):
    r = dict(p)
    for k, v in q.items():
        w = r.get(k, 0) + v
        if w:
            r[k] = w
        elif k in r:
            del r[k]
    return r


def _pneg(p):
    return {k: -v for k, v in p.items()}


def _pmul(p, q):
    r = {}
    for (a1, b1), v1 in p.items():
        for (a2, b2), v2 in q.items():
            k = (a1 + a2, b1 + b2)
            w = r.get(k, 0) + v1 * v2
            if w:
                r[k] = w
            elif k in r:
                del r[k]
    return r


def _ppow(p, n):
    r = {(0, 0): Fraction(1)}
    for _ in range(n):
        r = _pmul(r, p)
    return r


_ONE = {(0, 0): Fraction(1)}


class _Parser:
    def __init__(self, text):
        self.text = text
        self.toks = tokenize(text)
        self.i = 0

    def peek(self):
        return self.toks[self.i] if self.i < len(self.toks) else (None, None)

    def take(self):
        t = self.peek()
        self.i += 1
        return t

    def fail(self, what):
        raise ParseError(f"entry {self.text!r}: {what} at token {self.i} of {len(self.toks)}")

    def parse(self):
        if not self.toks:
            self.fail("empty entry")
        v = self.expr()
        if self.i != len(self.toks):
            self.fail(f"unexpected token {self.peek()!r}")
        return v

    def expr(self):                       # term (('+'|'-') term)*
        p, q = self.term()
        while self.peek() == ("OP", "+") or self.peek() == ("OP", "-"):
            op = self.take()[1]
            p2, q2 = self.term()
            num = _padd(_pmul(p, q2), _pmul(p2, q)) if op == "+" else _padd(_pmul(p, q2), _pneg(_pmul(p2, q)))
            p, q = num, _pmul(q, q2)
        return p, q

    def term(self):                       # factor (('*'|'/') factor)*
        p, q = self.factor()
        while self.peek() == ("OP", "*") or self.peek() == ("OP", "/"):
            op = self.take()[1]
            p2, q2 = self.factor()
            if op == "*":
                p, q = _pmul(p, p2), _pmul(q, q2)
            else:
                if not p2:
                    self.fail("division by the zero polynomial")
                p, q = _pmul(p, q2), _pmul(q, p2)
        return p, q

    def factor(self):                     # ('+'|'-') factor | power
        if self.peek() == ("OP", "-"):
            self.take()
            p, q = self.factor()
            return _pneg(p), q
        if self.peek() == ("OP", "+"):
            self.take()
            return self.factor()
        return self.power()

    def power(self):                      # atom ('**' INT)?
        p, q = self.atom()
        if self.peek() == ("OP", "**"):
            self.take()
            kind, val = self.take()
            if kind != "INT":
                self.fail(f"** needs a non-negative integer exponent, got {(kind, val)!r}")
            p, q = _ppow(p, val), _ppow(q, val)
        return p, q

    def atom(self):                       # INT | NAME | '(' expr ')'
        kind, val = self.take()
        if kind == "INT":
            return ({(0, 0): Fraction(val)} if val else {}), dict(_ONE)
        if kind == "NAME":
            return {((1, 0) if val == "d" else (0, 1)): Fraction(1)}, dict(_ONE)
        if (kind, val) == ("OP", "("):
            v = self.expr()
            if self.take() != ("OP", ")"):
                self.fail("missing ')'")
            return v
        self.fail(f"unexpected token {(kind, val)!r}")


def parse_entry(text):
    """(P, Q): the entry string as the exact rational function P(d, s) / Q(d, s) (bivariate dict polynomials
    over Fraction; Q is the product of the denominators met, never cancelled against P)."""
    return _Parser(text).parse()


# ---------------------------------------------------------------------------
# 2. Univariate helpers (ascending Fraction lists; [] is the zero polynomial)
# ---------------------------------------------------------------------------
def _utrim(p):
    while p and p[-1] == 0:
        p.pop()
    return p


def _uadd(p, q):
    n = max(len(p), len(q))
    return _utrim([(p[i] if i < len(p) else 0) + (q[i] if i < len(q) else 0) for i in range(n)])


def _uscale(p, c):
    return _utrim([c * v for v in p]) if c else []


def _umul(p, q):
    if not p or not q:
        return []
    r = [Fraction(0)] * (len(p) + len(q) - 1)
    for i, a in enumerate(p):
        if a:
            for j, b in enumerate(q):
                r[i + j] += a * b
    return _utrim(r)


def _upow(p, n):
    r = [Fraction(1)]
    for _ in range(n):
        r = _umul(r, p)
    return r


def _udivmod(p, q):
    """Exact polynomial division with remainder over Q."""
    p = list(p)
    if not q:
        raise ZeroDivisionError("polynomial division by zero")
    dq, lq = len(q) - 1, q[-1]
    out = [Fraction(0)] * max(0, len(p) - dq)
    while len(p) - 1 >= dq and p:
        c = p[-1] / lq
        k = len(p) - 1 - dq
        out[k] = c
        for i, b in enumerate(q):
            p[k + i] -= c * b
        _utrim(p)
        if not p:
            break
    return _utrim(out), p


def _ugcd(p, q):
    """Monic gcd over Q (Euclid); gcd(0, 0) = []."""
    p, q = _utrim(list(p)), _utrim(list(q))
    while q:
        _, r = _udivmod(p, q)
        p, q = q, r
    if not p:
        return []
    lc = p[-1]
    return [v / lc for v in p]


def _by_power_of(P, which):
    """P(d, s) as {exponent of d (which=0) or of s (which=1): univariate poly in the other variable}."""
    out = {}
    for (a, b), v in P.items():
        e, o = (a, b) if which == 0 else (b, a)
        lst = out.setdefault(e, [])
        if len(lst) <= o:
            lst.extend([Fraction(0)] * (o + 1 - len(lst)))
        lst[o] += v
    return {e: _utrim(l) for e, l in out.items() if _utrim(list(l))}


def _binom(n, k):
    r = 1
    for i in range(k):
        r = r * (n - i) // (i + 1)
    return r


def _subst_d(P):
    """P(4 - 2*eps, s) as {k: P_k(s)} (univariate Fraction lists in s), exact."""
    c0, c1 = D_SUBST
    out = {}
    for a, ps in _by_power_of(P, 0).items():
        for k in range(a + 1):
            coef = Fraction(_binom(a, k)) * c0 ** (a - k) * c1 ** k
            out[k] = _uadd(out.get(k, []), _uscale(ps, coef))
    return {k: v for k, v in out.items() if v}


def eps_grade(P, Q, neps):
    """{k: (N_k, D_k)} for k = 0..neps: the eps^k coefficient of P/Q at d = 4 - 2*eps as the univariate
    rational function N_k(s) / D_k(s), D_k = Q_0^(k+1), by the recursion N_k = P_k Q_0^k - sum_{m=1..k}
    Q_m N_{k-m} Q_0^(m-1); exact, uncancelled.  Absent k means the coefficient is identically zero."""
    Pk, Qk = _subst_d(P), _subst_d(Q)
    Q0 = Qk.get(0, [])
    if not Q0:
        raise ParseError("the denominator vanishes identically at d = 4 (a pole in eps): outside the grading")
    N = {}
    out = {}
    for k in range(neps + 1):
        nk = _umul(Pk.get(k, []), _upow(Q0, k))
        for m in range(1, k + 1):
            nk = _uadd(nk, _uscale(_umul(_umul(Qk.get(m, []), N[k - m]), _upow(Q0, m - 1)), Fraction(-1)))
        N[k] = nk
        if nk:
            out[k] = (nk, _upow(Q0, k + 1))
    return out


def _reduce(num, den):
    """num/den in lowest terms with an integer, primitive denominator of positive leading coefficient."""
    g = _ugcd(num, den)
    if len(g) > 1:
        num, _ = _udivmod(num, g)
        den, _ = _udivmod(den, g)
    import math
    L = 1
    for v in den:
        L = L * v.denominator // math.gcd(L, v.denominator)
    den_i = [v * L for v in den]
    num = [v * L for v in num]
    G = 0
    for v in den_i:
        G = math.gcd(G, int(v))
    if G > 1:
        den_i = [v / G for v in den_i]
        num = [v / G for v in num]
    if den_i[-1] < 0:
        den_i = [-v for v in den_i]
        num = [-v for v in num]
    return _utrim([Fraction(v) for v in num]), _utrim([Fraction(v) for v in den_i])


def _grade_rows(A, neps):
    """{(i, j, k): (N_k, D_k)} over every nonzero entry of the row-major string matrix A."""
    out = {}
    for i, row in enumerate(A):
        for j, e in enumerate(row):
            P, Q = parse_entry(e)
            if not P:
                continue
            for k, (nk, dk) in eps_grade(P, Q, neps).items():
                out[(i, j, k)] = (nk, dk)
    return out


def grade_canonical(A, neps=2):
    """The parser's own graded lists, cancelled to lowest terms (see the module docstring for why the served
    lines do not use this form)."""
    return {key: _reduce(nk, dk) for key, (nk, dk) in _grade_rows(A, neps).items()}


def load_graded(A, graded, neps=2):
    """The frozen graded coefficient lists of a shipped connection, checked against the strings.
    A: the connection's row-major list of entry strings; graded: the parsed JSON dict (or its path) with
    'grading' == 'd = 4 - 2*eps; eps^0..eps^<neps>' and 'entries' = [{'i','j','k','P': [ints ascending in s],
    'Q': [...]}].  Returns {(i, j, k): (pc, qc)} as Fraction lists in the file's order.  Refuses by name: a
    key set that differs from the parsed one, or any entry with pc * D_k != qc * N_k."""
    if isinstance(graded, str):
        graded = json.load(open(graded))
    want = f"d = 4 - 2*eps; eps^0..eps^{neps}"
    if graded.get("grading") != want:
        raise ParseError(f"graded file declares grading {graded.get('grading')!r}, this run needs {want!r}")
    parsed = _grade_rows(A, neps)
    frozen = {}
    for ent in graded["entries"]:
        key = (int(ent["i"]), int(ent["j"]), int(ent["k"]))
        if key in frozen:
            raise ParseError(f"graded file lists entry {key} twice")
        pc = [Fraction(int(c)) for c in ent["P"]]
        qc = [Fraction(int(c)) for c in ent["Q"]]
        if not _utrim(list(qc)):
            raise ParseError(f"graded entry {key} has a zero denominator")
        frozen[key] = (pc, qc)
    missing = sorted(set(parsed) - set(frozen))
    extra = sorted(set(frozen) - set(parsed))
    if missing or extra:
        raise ParseError(f"graded file disagrees with the parsed connection on the nonzero entries: "
                         f"{len(missing)} parsed but absent from the file {missing[:6]}, "
                         f"{len(extra)} in the file but zero when parsed {extra[:6]}")
    for key, (pc, qc) in frozen.items():
        nk, dk = parsed[key]
        if _umul(_utrim(list(pc)), dk) != _umul(_utrim(list(qc)), nk):
            raise ParseError(f"graded entry {key} is not the parsed rational function: pc*D_k != qc*N_k "
                             f"(pc={[str(c) for c in pc]}, qc={[str(c) for c in qc]})")
    return frozen


# ---------------------------------------------------------------------------
# 3. The singular-point census (exact real roots of the s-only denominator letters)
# ---------------------------------------------------------------------------
def _squarefree_split(n):
    """n = k^2 * m with m squarefree (n > 0): returns (k, m)."""
    k, m, p = 1, n, 2
    while p * p <= m:
        while m % (p * p) == 0:
            m //= p * p
            k *= p
        p += 1
    return k, m


def _rat_str(fr):
    return str(fr.numerator) if fr.denominator == 1 else f"{fr.numerator}/{fr.denominator}"


def _surd_term_str(b, n):
    """b*sqrt(n) printed as sympy prints Mul(Rational, sqrt): 'sqrt(3)', '-sqrt(3)', '4*sqrt(3)', 'sqrt(6)/2',
    '-3*sqrt(2)/2'."""
    sign = "-" if b < 0 else ""
    p, q = abs(b.numerator), b.denominator
    core = f"sqrt({n})" if p == 1 else f"{p}*sqrt({n})"
    return sign + core + (f"/{q}" if q != 1 else "")


class Root:
    """An exact real algebraic number a + b*sqrt(n) (a, b rational; n squarefree > 1, or b = 0)."""
    __slots__ = ("a", "b", "n")

    def __init__(self, a, b=Fraction(0), n=1):
        self.a, self.b, self.n = Fraction(a), Fraction(b), int(n)
        if self.b == 0:
            self.n = 1

    def key(self):
        return (self.a, self.b, self.n)

    def __eq__(self, other):
        return isinstance(other, Root) and self.key() == other.key()

    def __hash__(self):
        return hash(self.key())

    def __repr__(self):
        return f"Root({self.a!s}, {self.b!s}, {self.n})"

    def _mpf_at(self, prec):
        """The value rounded once to prec bits from a 30-guard-bit computation (an mpf at that precision)."""
        import mpmath as mp
        from mpmath.libmp import from_rational, mpf_pos, round_nearest
        with mp.workprec(prec + 30):
            if self.b == 0:
                return mp.mpf(from_rational(self.a.numerator, self.a.denominator, prec, round_nearest))
            x = mp.mpf(from_rational(self.a.numerator, self.a.denominator, prec + 30, round_nearest)) \
                + mp.mpf(from_rational(self.b.numerator, self.b.denominator, prec + 30, round_nearest)) \
                * mp.sqrt(self.n)
            return mp.mpf(mpf_pos(x._mpf_, prec, round_nearest))

    def __float__(self):
        """The double the earlier census produced (float(expr) = the 53-bit evaluation, correctly rounded)."""
        return float(self._mpf_at(53))

    def mpf(self, ndigits):
        """The mpf the earlier census produced from N(expr, ndigits) at the CURRENT working precision:
        evaluated to dps_to_prec(ndigits) bits, then rounded to the working precision."""
        import mpmath as mp
        from mpmath.libmp import dps_to_prec
        return mp.mpf(self._mpf_at(dps_to_prec(ndigits)))

    def __str__(self):
        a, b, n = self.a, self.b, self.n
        if b == 0:
            return _rat_str(a)
        if a == 0:
            return _surd_term_str(b, n)
        ta, tb = _rat_str(a), _surd_term_str(b, n)
        # the order of the two terms: sympy's special case (a positive number first when the surd's
        # coefficient is negative), else ascending by value
        if a > 0 and b < 0:
            first, second = ta, tb
        else:
            surd_val = float(b) * float(n) ** 0.5
            first, second = (ta, tb) if float(a) < surd_val else (tb, ta)
        if second.startswith("-"):
            return f"{first} - {second[1:]}"
        return f"{first} + {second}"


def _roots_of_factor(coeffs):
    """Exact real roots of an irreducible integer polynomial (ascending coefficients) of degree 1 or 2."""
    deg = len(coeffs) - 1
    if deg == 1:
        c0, c1 = coeffs
        return [Root(Fraction(-c0, c1))]
    if deg == 2:
        c, b, a = coeffs
        disc = b * b - 4 * a * c
        if disc <= 0:
            return []
        k, m = _squarefree_split(disc)
        if m == 1:
            r1, r2 = Fraction(-b - k, 2 * a), Fraction(-b + k, 2 * a)
            return [Root(r1), Root(r2)]
        centre, half = Fraction(-b, 2 * a), Fraction(k, 2 * a)
        return [Root(centre, -half, m), Root(centre, half, m)]
    raise ParseError(f"an irreducible s-only denominator factor of degree {deg} ({coeffs}): the shipped "
                     f"connections' letters are linear and quadratic; no exact root form is provided beyond that")


def poly_str(P):
    """A plain printout of a bivariate dict polynomial (used only to name a mixed d,s factor)."""
    terms = []
    for (a, b), v in sorted(P.items(), reverse=True):
        mon = "*".join([f"d**{a}" if a > 1 else "d"] * (a > 0) + [f"s**{b}" if b > 1 else "s"] * (b > 0))
        terms.append(f"{_rat_str(v)}*{mon}" if mon and v != 1 else (mon or _rat_str(v)))
    return " + ".join(terms) if terms else "0"


def singular_points(A):
    """(points, mixed): every real root of every s-only irreducible factor of the denominators of A, with
    s = 0 always included, as Root objects sorted by value (deduplicated); 'mixed' = the residual d,s-mixing
    denominator parts printed by poly_str (the shipped connections have none).  Requires python-flint."""
    from flint import fmpz_poly
    import math
    letters = set()
    mixed = []
    for row in A:
        for e in row:
            P, Q = parse_entry(e)
            if not P:
                continue
            by_d = _by_power_of(Q, 0)
            g = []
            for ps in by_d.values():
                g = _ugcd(g, ps)
            if len(g) > 1:                       # the s-only part of Q, monic
                L = 1
                for v in g:
                    L = L * v.denominator // math.gcd(L, v.denominator)
                gi = [int(v * L) for v in g]
                G = 0
                for v in gi:
                    G = math.gcd(G, v)
                gi = [v // G for v in gi]
                _c, facs = fmpz_poly(gi).factor()
                for f, _mult in facs:
                    letters.add(tuple(int(c) for c in f.coeffs()))
            # the residual R = Q / g: a d-only factor (skipped, as the earlier census skipped d-only letters),
            # else a factor mixing d and s (reported)
            R = {}
            for a, ps in by_d.items():
                q, r = _udivmod(ps, g) if len(g) > 1 else (ps, [])
                if r:
                    raise ParseError("internal: the s-content does not divide a coefficient")
                for b, v in enumerate(q):
                    if v:
                        R[(a, b)] = v
            if any(b > 0 for (_a, b) in R):
                by_s = _by_power_of(R, 1)
                h = []
                for ps in by_s.values():
                    h = _ugcd(h, ps)
                M = {}
                for b, pd in by_s.items():
                    q, r = _udivmod(pd, h) if len(h) > 1 else (pd, [])
                    for a, v in enumerate(q):
                        if v:
                            M[(a, b)] = v
                if any(b > 0 for (_a, b) in M) and any(a > 0 for (a, _b) in M):
                    mixed.append(poly_str(M))
    pts = {Root(0)}
    for coeffs in letters:
        for r in _roots_of_factor(list(coeffs)):
            pts.add(r)
    return sorted(pts, key=float), sorted(set(mixed))


# ---------------------------------------------------------------------------
# 4. Points and rationals
# ---------------------------------------------------------------------------
def rat(text):
    """The exact Fraction of a point / seed / reference string: '-7/3', '-2', '-4.25', '5/2' (an int is
    accepted as is).  Raises ValueError (or ZeroDivisionError for 'p/0') on anything else."""
    if isinstance(text, Fraction):
        return text
    if isinstance(text, int):
        return Fraction(text)
    return Fraction(str(text).strip())


def mpf_rat(fr):
    """The mpf of a Fraction at the CURRENT working precision, correctly rounded (mpmath's from_rational --
    the conversion sympy's Rational used, so every endpoint is the same binary number as before)."""
    import mpmath as mp
    from mpmath.libmp import from_rational, round_nearest
    fr = rat(fr)
    return mp.mpf(from_rational(fr.numerator, fr.denominator, mp.mp.prec, round_nearest))