Repository: mpmath/mpmath
Pull request #691 (2023-04-25 08:31:17)

Drop private mpq class, use Rational's, provided by backend
* Missing feature (for pure-python) is an additional caching mechanism.  Not sure it worth.
* Unfortunately, gmpy2 has no released version with support for mpq.as_numer_denom(), hence I use foo.numerator/denominator to unpack components of rationals...  Edit: the numbers.Rational doesn't require this method, so we could not rely on it.

skirpichev: Missing feature (for pure-python) is additional caching mechanism.
oscarbenjamin: This looks good to me.

I don't know how significant the caching mechanism is but it does not look as if there are any places where `MPQ` is used in a tight loop.
skirpichev: @oscarbenjamin, I've some doubts about the ctx.mpq attribute.  It's not documented, but apparently people discover it and use (e.g. #532).  Maybe we should keep this and issue a deprecation warning before removing?

Patch:
diff --git a/mpmath/ctx_fp.py b/mpmath/ctx_fp.py
index f1edd1c..f33a730 100644
--- a/mpmath/ctx_fp.py
+++ b/mpmath/ctx_fp.py
@@ -21,8 +21,6 @@ def __init__(ctx):
 
         ctx._init_aliases()
 
-    _mpq = lambda cls, x: float(x[0])/x[1]
-
     NoConvergence = libmp.NoConvergence
 
     def _get_prec(ctx): return 53
diff --git a/mpmath/ctx_iv.py b/mpmath/ctx_iv.py
index 13e8058..c8a32a4 100644
--- a/mpmath/ctx_iv.py
+++ b/mpmath/ctx_iv.py
@@ -395,12 +395,6 @@ def make_mpc(ctx, v):
         a._mpci_ = v
         return a
 
-    def _mpq(ctx, pq):
-        p, q = pq
-        a = libmp.from_rational(p, q, ctx.prec, round_floor)
-        b = libmp.from_rational(p, q, ctx.prec, round_ceiling)
-        return ctx.make_mpf((a, b))
-
     def convert(ctx, x):
         if isinstance(x, (ctx.mpf, ctx.mpc)):
             return x
diff --git a/mpmath/ctx_mp.py b/mpmath/ctx_mp.py
index 37f9c71..c17017e 100644
--- a/mpmath/ctx_mp.py
+++ b/mpmath/ctx_mp.py
@@ -6,16 +6,16 @@
 import functools
 import re
 
-from . import function_docs, libmp, rational
+from . import function_docs, libmp
 from .ctx_base import StandardBaseContext
-from .libmp import (MPZ_ONE, ComplexResult, bitcount, dps_to_prec, finf, fnan,
-                    fninf, fone, from_rational, fzero, int_types, mpc_add,
-                    mpc_add_mpf, mpc_div, mpc_div_mpf, mpc_mul, mpc_mul_mpf,
-                    mpc_neg, mpc_sub, mpc_sub_mpf, mpc_to_str, mpf_add,
-                    mpf_apery, mpf_catalan, mpf_degree, mpf_div, mpf_e,
-                    mpf_euler, mpf_glaisher, mpf_khinchin, mpf_ln2, mpf_ln10,
-                    mpf_mertens, mpf_mul, mpf_neg, mpf_phi, mpf_pi, mpf_rand,
-                    mpf_sub, mpf_twinprime, repr_dps, to_str)
+from .libmp import (MPQ, MPZ_ONE, ComplexResult, bitcount, dps_to_prec, finf,
+                    fnan, fninf, fone, from_rational, fzero, int_types,
+                    mpc_add, mpc_add_mpf, mpc_div, mpc_div_mpf, mpc_mul,
+                    mpc_mul_mpf, mpc_neg, mpc_sub, mpc_sub_mpf, mpc_to_str,
+                    mpf_add, mpf_apery, mpf_catalan, mpf_degree, mpf_div,
+                    mpf_e, mpf_euler, mpf_glaisher, mpf_khinchin, mpf_ln2,
+                    mpf_ln10, mpf_mertens, mpf_mul, mpf_neg, mpf_phi, mpf_pi,
+                    mpf_rand, mpf_sub, mpf_twinprime, repr_dps, to_str)
 from .libmp.backend import BACKEND
 
 
@@ -41,11 +41,9 @@ def __init__(ctx):
         ctx.trap_complex = False
         ctx.pretty = False
         ctx.types = [ctx.mpf, ctx.mpc, ctx.constant]
-        ctx._mpq = rational.mpq
         ctx.default()
         StandardBaseContext.__init__(ctx)
 
-        ctx.mpq = rational.mpq
         ctx.init_builtins()
 
         ctx.hyp_summators = {}
@@ -300,7 +298,7 @@ def isnan(ctx, x):
             return x._mpf_ == fnan
         if hasattr(x, "_mpc_"):
             return fnan in x._mpc_
-        if isinstance(x, int_types) or isinstance(x, rational.mpq):
+        if isinstance(x, int_types) or isinstance(x, MPQ):
             return False
         x = ctx.convert(x)
         if hasattr(x, '_mpf_') or hasattr(x, '_mpc_'):
@@ -346,8 +344,8 @@ def isnpint(ctx, x):
             return not x.imag and ctx.isnpint(x.real)
         if type(x) in int_types:
             return x <= 0
-        if isinstance(x, ctx.mpq):
-            p, q = x._mpq_
+        if isinstance(x, MPQ):
+            p, q = x.numerator, x.denominator
             if not p:
                 return True
             return q == 1 and p <= 0
@@ -1102,8 +1100,8 @@ def nint_distance(ctx, x):
         typx = type(x)
         if typx in int_types:
             return int(x), ctx.ninf
-        elif typx is rational.mpq:
-            p, q = x._mpq_
+        elif typx is MPQ:
+            p, q = x.numerator, x.denominator
             n, r = divmod(p, q)
             if 2*r >= q:
                 n += 1
diff --git a/mpmath/ctx_mp_python.py b/mpmath/ctx_mp_python.py
index 25fcdb2..af5385f 100644
--- a/mpmath/ctx_mp_python.py
+++ b/mpmath/ctx_mp_python.py
@@ -1,7 +1,7 @@
 import numbers
 
-from . import function_docs, rational
-from .libmp import (MPZ, ComplexResult, bitcount, dps_to_prec, finf, fnan,
+from . import function_docs
+from .libmp import (MPQ, MPZ, ComplexResult, bitcount, dps_to_prec, finf, fnan,
                     fninf, from_Decimal, from_float, from_int, from_man_exp,
                     from_npfloat, from_pickable, from_rational, from_str,
                     fzero, int_types, mpc_abs, mpc_add, mpc_add_mpf,
@@ -94,8 +94,8 @@ def mpf_convert_rhs(cls, x):
         if isinstance(x, int_types): return from_int(x)
         if isinstance(x, float): return from_float(x)
         if isinstance(x, complex_types): return cls.context.mpc(x)
-        if isinstance(x, rational.mpq):
-            p, q = x._mpq_
+        if isinstance(x, MPQ):
+            p, q = x.numerator, x.denominator
             return from_rational(p, q, cls.context.prec)
         if hasattr(x, '_mpf_'): return x._mpf_
         if hasattr(x, '_mpmath_'):
@@ -623,11 +623,11 @@ def convert(ctx, x, strings=True):
             return ctx.make_mpc((from_float(x.real), from_float(x.imag)))
         if type(x).__module__ == 'numpy': return ctx.npconvert(x)
         if isinstance(x, numbers.Rational): # e.g. Fraction
-            try: x = rational.mpq(int(x.numerator), int(x.denominator))
+            try: x = MPQ(int(x.numerator), int(x.denominator))
             except: pass
         prec, rounding = ctx._prec_rounding
-        if isinstance(x, rational.mpq):
-            p, q = x._mpq_
+        if isinstance(x, MPQ):
+            p, q = x.numerator, x.denominator
             return ctx.make_mpf(from_rational(p, q, prec))
         if strings and isinstance(x, str):
             try:
@@ -680,7 +680,7 @@ def isinf(ctx, x):
         if hasattr(x, "_mpc_"):
             re, im = x._mpc_
             return re in (finf, fninf) or im in (finf, fninf)
-        if isinstance(x, int_types) or isinstance(x, rational.mpq):
+        if isinstance(x, int_types) or isinstance(x, MPQ):
             return False
         x = ctx.convert(x)
         return ctx.isinf(x)
@@ -718,7 +718,7 @@ def isnormal(ctx, x):
             if re == fzero: return im_normal
             if im == fzero: return re_normal
             return re_normal and im_normal
-        if isinstance(x, int_types) or isinstance(x, rational.mpq):
+        if isinstance(x, int_types) or isinstance(x, MPQ):
             return bool(x)
         x = ctx.convert(x)
         return ctx.isnormal(x)
@@ -762,8 +762,8 @@ def isint(ctx, x, gaussian=False):
                 im_isint = (iman and iexp >= 0) or im == fzero
                 return re_isint and im_isint
             return re_isint and im == fzero
-        if isinstance(x, rational.mpq):
-            p, q = x._mpq_
+        if isinstance(x, MPQ):
+            p, q = x.numerator, x.denominator
             return p % q == 0
         x = ctx.convert(x)
         return ctx.isint(x, gaussian)
@@ -979,8 +979,6 @@ def _convert_param(ctx, x):
             p = None
             if isinstance(x, tuple):
                 p, q = x
-            elif hasattr(x, '_mpq_'):
-                p, q = x._mpq_
             elif isinstance(x, str) and '/' in x:
                 p, q = x.split('/')
                 p = int(p)
@@ -988,7 +986,7 @@ def _convert_param(ctx, x):
             if p is not None:
                 if not p % q:
                     return p // q, 'Z'
-                return ctx.mpq(p,q), 'Q'
+                return MPQ(p,q), 'Q'
             x = ctx.convert(x)
             if hasattr(x, "_mpc_"):
                 v, im = x._mpc_
@@ -1006,7 +1004,7 @@ def _convert_param(ctx, x):
                 if exp >= 0:
                     return int(man) << exp, 'Z'
                 p, q = int(man), (1<<(-exp))
-                return ctx.mpq(p,q), 'Q'
+                return MPQ(p,q), 'Q'
             x = ctx.make_mpf(v)
             return x, 'R'
         if not exp:
@@ -1057,8 +1055,8 @@ def mag(ctx, x):
             if x:
                 return bitcount(abs(x))
             return ctx.ninf
-        if isinstance(x, rational.mpq):
-            p, q = x._mpq_
+        if isinstance(x, MPQ):
+            p, q = x.numerator, x.denominator
             if p:
                 return 1 + bitcount(abs(p)) - bitcount(q)
             return ctx.ninf
diff --git a/mpmath/functions/bessel.py b/mpmath/functions/bessel.py
index 5bcc7ab..2e6d7fd 100644
--- a/mpmath/functions/bessel.py
+++ b/mpmath/functions/bessel.py
@@ -1,3 +1,4 @@
+from ..libmp.backend import MPQ
 from .functions import defun, defun_wrapped
 
 @defun
@@ -259,7 +260,7 @@ def _anger(ctx,which,v,z,**kwargs):
     v = ctx._convert_param(v)[0]
     z = ctx.convert(z)
     def h(v):
-        b = ctx.mpq_1_2
+        b = MPQ(1,2)
         u = v*b
         m = b*3
         a1,a2,b1,b2 = m-u, m+u, 1-u, 1+u
@@ -288,7 +289,7 @@ def lommels1(ctx, u, v, z, **kwargs):
     v = ctx._convert_param(v)[0]
     z = ctx.convert(z)
     def h(u,v):
-        b = ctx.mpq_1_2
+        b = MPQ(1,2)
         w = ctx.square_exp_arg(z, mult=-0.25)
         return ([u-v+1, u+v+1, z], [-1, -1, u+1], [], [], [1], \
             [b*(u-v+3),b*(u+v+3)], w),
@@ -302,11 +303,11 @@ def lommels2(ctx, u, v, z, **kwargs):
     # Asymptotic expansion (GR p. 947) -- need to be careful
     # not to use for small arguments
     # def h(u,v):
-    #    b = ctx.mpq_1_2
+    #    b = MPQ(1,2)
     #    w = -(z/2)**(-2)
     #    return ([z], [u-1], [], [], [b*(1-u+v)], [b*(1-u-v)], w),
     def h(u,v):
-        b = ctx.mpq_1_2
+        b = MPQ(1,2)
         w = ctx.square_exp_arg(z, mult=-0.25)
         T1 = [u-v+1, u+v+1, z], [-1, -1, u+1], [], [], [1], [b*(u-v+3),b*(u+v+3)], w
         T2 = [2, z], [u+v-1, -v], [v, b*(u+v+1)], [b*(v-u+1)], [], [1-v], w
@@ -425,7 +426,7 @@ def _airyderiv_0(ctx, z, n, ntype, which):
     if ntype == 'Z':
         if n < 0:
             return z
-        r = ctx.mpq_1_3
+        r = MPQ(1,3)
         prec = ctx.prec
         try:
             ctx.prec += 10
@@ -489,8 +490,8 @@ def h():
                     ctx.prec -= extraprec
                     C1 = _airyai_C1(ctx) * 0.5
                     C2 = _airyai_C2(ctx)
-                    T1 = [C1,z],[1,2],[],[],[],[ctx.mpq_5_3],w
-                    T2 = [C2],[1],[],[],[],[ctx.mpq_1_3],w
+                    T1 = [C1,z],[1,2],[],[],[],[MPQ(5,3)],w
+                    T2 = [C2],[1],[],[],[],[MPQ(1,3)],w
                     return T1, T2
             return ctx.hypercomb(h, [], **kwargs)
         else:
@@ -501,7 +502,7 @@ def h(n):
                 ctx.prec += extraprec
                 w = z**3/9
                 ctx.prec -= extraprec
-                q13,q23,q43 = ctx.mpq_1_3, ctx.mpq_2_3, ctx.mpq_4_3
+                q13,q23,q43 = MPQ(1,3), MPQ(2,3), MPQ(4,3)
                 a1=q13; a2=1; b1=(1-n)*q13; b2=(2-n)*q13; b3=1-n*q13
                 T1 = [3, z], [n-q23, -n], [a1], [b1,b2,b3], \
                     [a1,a2], [b1,b2,b3], w
@@ -530,8 +531,8 @@ def h():
                 ctx.prec -= extraprec
                 C1 = _airyai_C1(ctx)
                 C2 = _airyai_C2(ctx)
-                T1 = [C1],[1],[],[],[],[ctx.mpq_2_3],w
-                T2 = [z*C2],[1],[],[],[],[ctx.mpq_4_3],w
+                T1 = [C1],[1],[],[],[],[MPQ(2,3)],w
+                T2 = [z*C2],[1],[],[],[],[MPQ(4,3)],w
                 return T1, T2
         return ctx.hypercomb(h, [], **kwargs)
 
@@ -574,8 +575,8 @@ def h():
                 ctx.prec -= extraprec
                 C1 = _airybi_C1(ctx)*0.5
                 C2 = _airybi_C2(ctx)
-                T1 = [C1,z],[1,2],[],[],[],[ctx.mpq_5_3],w
-                T2 = [C2],[1],[],[],[],[ctx.mpq_1_3],w
+                T1 = [C1,z],[1,2],[],[],[],[MPQ(5,3)],w
+                T2 = [C2],[1],[],[],[],[MPQ(1,3)],w
                 return T1, T2
             return ctx.hypercomb(h, [], **kwargs)
         else:
@@ -585,9 +586,9 @@ def h(n):
                 ctx.prec += extraprec
                 w = z**3/9
                 ctx.prec -= extraprec
-                q13,q23,q43 = ctx.mpq_1_3, ctx.mpq_2_3, ctx.mpq_4_3
-                q16 = ctx.mpq_1_6
-                q56 = ctx.mpq_5_6
+                q13,q23,q43 = MPQ(1,3), MPQ(2,3), MPQ(4,3)
+                q16 = MPQ(1,6)
+                q56 = MPQ(5,6)
                 a1=q13; a2=1; b1=(1-n)*q13; b2=(2-n)*q13; b3=1-n*q13
                 T1 = [3, z], [n-q16, -n], [a1], [b1,b2,b3], \
                     [a1,a2], [b1,b2,b3], w
@@ -606,8 +607,8 @@ def h():
             ctx.prec -= extraprec
             C1 = _airybi_C1(ctx)
             C2 = _airybi_C2(ctx)
-            T1 = [C1],[1],[],[],[],[ctx.mpq_2_3],w
-            T2 = [z*C2],[1],[],[],[],[ctx.mpq_4_3],w
+            T1 = [C1],[1],[],[],[],[MPQ(2,3)],w
+            T2 = [z*C2],[1],[],[],[],[MPQ(4,3)],w
             return T1, T2
         return ctx.hypercomb(h, [], **kwargs)
 
@@ -686,7 +687,7 @@ def h():
         w = z**3/9
         ctx.prec -= extraprec
         T1 = [A], [1], [], [], [], [], 0
-        T2 = [B,z], [-1,2], [], [], [1], [ctx.mpq_4_3,ctx.mpq_5_3], w
+        T2 = [B,z], [-1,2], [], [], [1], [MPQ(4,3),MPQ(5,3)], w
         return T1, T2
     return ctx.hypercomb(h, [], **kwargs)
 
diff --git a/mpmath/functions/functions.py b/mpmath/functions/functions.py
index ba605ac..972f5b9 100644
--- a/mpmath/functions/functions.py
+++ b/mpmath/functions/functions.py
@@ -19,24 +19,6 @@ def __init__(self):
             f, wrap = cls.defined_functions[name]
             cls._wrap_specfun(name, f, wrap)
 
-        self.mpq_1 = self._mpq((1,1))
-        self.mpq_0 = self._mpq((0,1))
-        self.mpq_1_2 = self._mpq((1,2))
-        self.mpq_3_2 = self._mpq((3,2))
-        self.mpq_1_4 = self._mpq((1,4))
-        self.mpq_1_16 = self._mpq((1,16))
-        self.mpq_3_16 = self._mpq((3,16))
-        self.mpq_5_2 = self._mpq((5,2))
-        self.mpq_3_4 = self._mpq((3,4))
-        self.mpq_7_4 = self._mpq((7,4))
-        self.mpq_5_4 = self._mpq((5,4))
-        self.mpq_1_3 = self._mpq((1,3))
-        self.mpq_2_3 = self._mpq((2,3))
-        self.mpq_4_3 = self._mpq((4,3))
-        self.mpq_1_6 = self._mpq((1,6))
-        self.mpq_5_6 = self._mpq((5,6))
-        self.mpq_5_3 = self._mpq((5,3))
-
         self._misc_const_cache = {}
 
         self._aliases.update({
diff --git a/mpmath/functions/hypergeometric.py b/mpmath/functions/hypergeometric.py
index b322146..c83c7b1 100644
--- a/mpmath/functions/hypergeometric.py
+++ b/mpmath/functions/hypergeometric.py
@@ -1,3 +1,4 @@
+from ..libmp.backend import MPQ
 from .functions import defun, defun_wrapped
 
 def _check_need_perturb(ctx, terms, prec, discard_known_zeros):
@@ -291,10 +292,10 @@ def h():
                     w = ctx.sqrt(-z)
                     jw = ctx.j*w
                     u = 1/(4*jw)
-                    c = ctx.mpq_1_2 - b
+                    c = MPQ(1,2) - b
                     E = ctx.exp(2*jw)
-                    T1 = ([-jw,E], [c,-1], [], [], [b-ctx.mpq_1_2, ctx.mpq_3_2-b], [], -u)
-                    T2 = ([jw,E], [c,1], [], [], [b-ctx.mpq_1_2, ctx.mpq_3_2-b], [], u)
+                    T1 = ([-jw,E], [c,-1], [], [], [b-MPQ(1,2), MPQ(3,2)-b], [], -u)
+                    T2 = ([jw,E], [c,1], [], [], [b-MPQ(1,2), MPQ(3,2)-b], [], u)
                     return T1, T2
                 v = ctx.hypercomb(h, [], force_series=True)
                 v = ctx.gamma(b)/(2*ctx.sqrt(ctx.pi))*v
@@ -364,8 +365,8 @@ def _hyp2f1_gosper(ctx,a,b,c,z,**kwargs):
         # things a bit unreadable. The formula is quite messy to begin
         # with, though...
         abz = a*b*z
-        ch = c * ctx.mpq_1_2
-        c1h = (c+1) * ctx.mpq_1_2
+        ch = c * MPQ(1,2)
+        c1h = (c+1) * MPQ(1,2)
         nz = 1-z
         g = z/nz
         abg = a*b*g
@@ -450,9 +451,9 @@ def _hyp2f1(ctx, a_s, b_s, z, **kwargs):
         # Use 1/z transformation
         if absz >= 1.3:
             def h(a,b):
-                t = ctx.mpq_1-c; ab = a-b; rz = 1/z
-                T1 = ([-z],[-a], [c,-ab],[b,c-a], [a,t+a],[ctx.mpq_1+ab],  rz)
-                T2 = ([-z],[-b], [c,ab],[a,c-b], [b,t+b],[ctx.mpq_1-ab],  rz)
+                t = MPQ(1)-c; ab = a-b; rz = 1/z
+                T1 = ([-z],[-a], [c,-ab],[b,c-a], [a,t+a],[MPQ(1)+ab],  rz)
+                T2 = ([-z],[-b], [c,ab],[a,c-b], [b,t+b],[MPQ(1)-ab],  rz)
                 return T1, T2
             v = ctx.hypercomb(h, [a,b], **kwargs)
 
@@ -840,12 +841,12 @@ def _hyp1f2(ctx, a_s, b_s, z, **kwargs):
                 # http://functions.wolfram.com/HypergeometricFunctions/
                 # Hypergeometric1F2/06/02/03/
                 def h(a1,b1,b2):
-                    X = ctx.mpq_1_2*(a1-b1-b2+ctx.mpq_1_2)
+                    X = MPQ(1,2)*(a1-b1-b2+MPQ(1,2))
                     c = {}
                     c[0] = ctx.one
-                    c[1] = 2*(ctx.mpq_1_4*(3*a1+b1+b2-2)*(a1-b1-b2)+b1*b2-ctx.mpq_3_16)
-                    c[2] = 2*(b1*b2+ctx.mpq_1_4*(a1-b1-b2)*(3*a1+b1+b2-2)-ctx.mpq_3_16)**2+\
-                        ctx.mpq_1_16*(-16*(2*a1-3)*b1*b2 + \
+                    c[1] = 2*(MPQ(1,4)*(3*a1+b1+b2-2)*(a1-b1-b2)+b1*b2-MPQ(3,16))
+                    c[2] = 2*(b1*b2+MPQ(1,4)*(a1-b1-b2)*(3*a1+b1+b2-2)-MPQ(3,16))**2+\
+                        MPQ(1,16)*(-16*(2*a1-3)*b1*b2 + \
                         4*(a1-b1-b2)*(-8*a1**2+11*a1+b1+b2-2)-3)
                     s1 = 0
                     s2 = 0
@@ -854,9 +855,9 @@ def h(a1,b1,b2):
                     while 1:
                         if k not in c:
                             uu1 = (3*k**2+(-6*a1+2*b1+2*b2-4)*k + 3*a1**2 - \
-                                (b1-b2)**2 - 2*a1*(b1+b2-2) + ctx.mpq_1_4)
-                            uu2 = (k-a1+b1-b2-ctx.mpq_1_2)*(k-a1-b1+b2-ctx.mpq_1_2)*\
-                                (k-a1+b1+b2-ctx.mpq_5_2)
+                                (b1-b2)**2 - 2*a1*(b1+b2-2) + MPQ(1,4))
+                            uu2 = (k-a1+b1-b2-MPQ(1,2))*(k-a1-b1+b2-MPQ(1,2))*\
+                                (k-a1+b1+b2-MPQ(5,2))
                             c[k] = ctx.one/(2*k)*(uu1*c[k-1]-uu2*c[k-2])
                         w = c[k] * (-z)**(-0.5*k)
                         t1 = (-ctx.j)**k * ctx.mpf(2)**(-k) * w
@@ -919,7 +920,7 @@ def _hyp2f3(ctx, a_s, b_s, z, **kwargs):
                 # http://functions.wolfram.com/HypergeometricFunctions/
                 # Hypergeometric2F3/06/02/03/01/0002/
                 def h(a1,a2,b1,b2,b3):
-                    X = ctx.mpq_1_2*(a1+a2-b1-b2-b3+ctx.mpq_1_2)
+                    X = MPQ(1,2)*(a1+a2-b1-b2-b3+MPQ(1,2))
                     A2 = a1+a2
                     B3 = b1+b2+b3
                     A = a1*a2
@@ -927,8 +928,8 @@ def h(a1,a2,b1,b2,b3):
                     R = b1*b2*b3
                     c = {}
                     c[0] = ctx.one
-                    c[1] = 2*(B - A + ctx.mpq_1_4*(3*A2+B3-2)*(A2-B3) - ctx.mpq_3_16)
-                    c[2] = ctx.mpq_1_2*c[1]**2 + ctx.mpq_1_16*(-16*(2*A2-3)*(B-A) + 32*R +\
+                    c[1] = 2*(B - A + MPQ(1,4)*(3*A2+B3-2)*(A2-B3) - MPQ(3,16))
+                    c[2] = MPQ(1,2)*c[1]**2 + MPQ(1,16)*(-16*(2*A2-3)*(B-A) + 32*R +\
                         4*(-8*A2**2 + 11*A2 + 8*A + B3 - 2)*(A2-B3)-3)
                     s1 = 0
                     s2 = 0
diff --git a/mpmath/functions/orthogonal.py b/mpmath/functions/orthogonal.py
index 74b4fe2..c135864 100644
--- a/mpmath/functions/orthogonal.py
+++ b/mpmath/functions/orthogonal.py
@@ -1,3 +1,4 @@
+from ..libmp.backend import MPQ
 from .functions import defun, defun_wrapped
 
 def _hermite_param(ctx, n, z, parabolic_cylinder):
@@ -8,7 +9,7 @@ def _hermite_param(ctx, n, z, parabolic_cylinder):
     """
     n, ntyp = ctx._convert_param(n)
     z = ctx.convert(z)
-    q = -ctx.mpq_1_2
+    q = -MPQ(1,2)
     # For re(z) > 0, 2F0 -- http://functions.wolfram.com/
     #     HypergeometricFunctions/HermiteHGeneral/06/02/0009/
     # Otherwise, there is a reflection formula
@@ -162,7 +163,7 @@ def pcfu(ctx, a, z, **kwargs):
 
     """
     n, _ = ctx._convert_param(a)
-    return ctx.pcfd(-n-ctx.mpq_1_2, z)
+    return ctx.pcfd(-n-MPQ(1,2), z)
 
 @defun
 def pcfv(ctx, a, z, **kwargs):
@@ -198,8 +199,8 @@ def pcfv(ctx, a, z, **kwargs):
     """
     n, ntype = ctx._convert_param(a)
     z = ctx.convert(z)
-    q = ctx.mpq_1_2
-    r = ctx.mpq_1_4
+    q = MPQ(1,2)
+    r = MPQ(1,4)
     if ntype == 'Q' and ctx.isint(n*2):
         # Faster for half-integers
         def h():
@@ -287,7 +288,7 @@ def h():
         w1 = ctx.fmul(w, -0.25, exact=True)
         w2 = ctx.fmul(w, 0.5, exact=True)
         e = ctx.exp(w1)
-        return [e], [1], [], [], [ctx.mpq_1_2*a+ctx.mpq_1_4], [ctx.mpq_1_2], w2
+        return [e], [1], [], [], [MPQ(1,2)*a+MPQ(1,4)], [MPQ(1,2)], w2
     return ctx.hypercomb(h, [], **kwargs)
 
 @defun
@@ -299,8 +300,8 @@ def h():
         w1 = ctx.fmul(w, -0.25, exact=True)
         w2 = ctx.fmul(w, 0.5, exact=True)
         e = ctx.exp(w1)
-        return [e, z], [1, 1], [], [], [ctx.mpq_1_2*a+ctx.mpq_3_4], \
-            [ctx.mpq_3_2], w2
+        return [e, z], [1, 1], [], [], [MPQ(1,2)*a+MPQ(3,4)], \
+            [MPQ(3,2)], w2
     return ctx.hypercomb(h, [], **kwargs)
 """
 
diff --git a/mpmath/functions/zeta.py b/mpmath/functions/zeta.py
index d977bb1..ddc85cd 100644
--- a/mpmath/functions/zeta.py
+++ b/mpmath/functions/zeta.py
@@ -633,7 +633,7 @@ def _hurwitz_reflection(ctx, s, a, d, atype):
     # We now require a to be standardized
     v = 0
     shift = 0
-    b = a
+    b = ctx.mpf(a)
     while ctx.re(b) > 1:
         b -= 1
         v -= b**negs
@@ -644,7 +644,7 @@ def _hurwitz_reflection(ctx, s, a, d, atype):
         shift += 1
     # Rational reflection formula
     try:
-        p, q = a._mpq_
+        p, q = a.numerator, a.denominator
     except:
         assert a == int(a)
         p = int(a)
diff --git a/mpmath/libmp/__init__.py b/mpmath/libmp/__init__.py
index cc7dafe..bb7abc5 100644
--- a/mpmath/libmp/__init__.py
+++ b/mpmath/libmp/__init__.py
@@ -72,6 +72,6 @@
   isqrt, isqrt_small, isqrt_fast, sqrt_fixed, sqrtrem, ifib, ifac,
   list_primes, isprime, moebius, gcd, eulernum, stirling1, stirling2)
 
-from .backend import (gmpy, sage, BACKEND, STRICT, MPZ,
+from .backend import (gmpy, sage, BACKEND, STRICT, MPZ, MPQ,
   MPZ_ZERO, MPZ_ONE, MPZ_TWO, MPZ_THREE, MPZ_FIVE, int_types,
   HASH_MODULUS, HASH_BITS)
diff --git a/mpmath/libmp/backend.py b/mpmath/libmp/backend.py
index dc6fe04..3f20cb8 100644
--- a/mpmath/libmp/backend.py
+++ b/mpmath/libmp/backend.py
@@ -1,3 +1,4 @@
+from fractions import Fraction
 import os
 import sys
 
@@ -23,6 +24,7 @@
 sage_utils = None
 BACKEND = 'python'
 MPZ = int
+MPQ = Fraction
 HASH_MODULUS = sys.hash_info.modulus
 HASH_BITS = 31 if sys.hash_info.width == 32 else 61
 
@@ -32,6 +34,7 @@
         import gmpy2 as gmpy
         BACKEND = 'gmpy'
         MPZ = gmpy.mpz
+        MPQ = gmpy.mpq
     except ImportError:
         pass
 
@@ -44,6 +47,7 @@
         sage_utils = _sage_utils
         BACKEND = 'sage'
         MPZ = sage.Integer
+        MPQ = lambda x, y=1: sage.Rational((x, y))
     except:
         pass
 
diff --git a/mpmath/libmp/libhyper.py b/mpmath/libmp/libhyper.py
index 6516c87..57f88a9 100644
--- a/mpmath/libmp/libhyper.py
+++ b/mpmath/libmp/libhyper.py
@@ -125,7 +125,7 @@ def make_hyp_summator(key):
             add("%sINT_%i = coeffs[%i]" % (W, i, i))
         elif flag == 'Q':
             ([arat,brat][i >= p]).append(i)
-            add("%sP_%i, %sQ_%i = coeffs[%i]._mpq_" % (W, i, W, i, i))
+            add("%sP_%i, %sQ_%i = coeffs[%i].numerator, coeffs[%i].denominator" % (W, i, W, i, i, i))
         elif flag == 'R':
             ([areal,breal][i >= p]).append(i)
             add("xsign, xm, xe, xbc = coeffs[%i]._mpf_" % i)
diff --git a/mpmath/rational.py b/mpmath/rational.py
deleted file mode 100644
index 118c5bb..0000000
--- a/mpmath/rational.py
+++ /dev/null
@@ -1,229 +0,0 @@
-import numbers
-import operator
-import sys
-
-from .libmp import int_types, mpf_hash, bitcount, from_man_exp, HASH_MODULUS
-
-new = object.__new__
-
-def create_reduced(p, q, _cache={}):
-    key = p, q
-    if key in _cache:
-        return _cache[key]
-    x, y = p, q
-    while y:
-        x, y = y, x % y
-    if x != 1:
-        p //= x
-        q //= x
-    v = new(mpq)
-    v._mpq_ = p, q
-    # Speedup integers, half-integers and other small fractions
-    if q <= 4 and abs(key[0]) < 100:
-        _cache[key] = v
-    return v
-
-class mpq:
-    """
-    Exact rational type, currently only intended for internal use.
-    """
-
-    __slots__ = ["_mpq_"]
-
-    def __new__(cls, p, q=1):
-        if type(p) is tuple:
-            p, q = p
-        elif hasattr(p, '_mpq_'):
-            p, q = p._mpq_
-        return create_reduced(p, q)
-
-    def __repr__(s):
-        return "mpq(%s,%s)" % s._mpq_
-
-    def __str__(s):
-        return "(%s/%s)" % s._mpq_
-
-    def __int__(s):
-        a, b = s._mpq_
-        return a // b
-
-    def __nonzero__(s):
-        return bool(s._mpq_[0])
-
-    __bool__ = __nonzero__
-
-    def __hash__(s):
-        a, b = s._mpq_
-        inverse = pow(b, HASH_MODULUS-2, HASH_MODULUS)
-        if not inverse:
-            h = sys.hash_info.inf
-        else:
-            h = (abs(a) * inverse) % HASH_MODULUS
-        if a < 0: h = -h
-        if h == -1: h = -2
-        return h
-
-    def __eq__(s, t):
-        ttype = type(t)
-        if ttype is mpq:
-            return s._mpq_ == t._mpq_
-        if ttype in int_types:
-            a, b = s._mpq_
-            if b != 1:
-                return False
-            return a == t
-        return NotImplemented
-
-    def __ne__(s, t):
-        ttype = type(t)
-        if ttype is mpq:
-            return s._mpq_ != t._mpq_
-        if ttype in int_types:
-            a, b = s._mpq_
-            if b != 1:
-                return True
-            return a != t
-        return NotImplemented
-
-    def _cmp(s, t, op):
-        ttype = type(t)
-        if ttype in int_types:
-            a, b = s._mpq_
-            return op(a, t*b)
-        if ttype is mpq:
-            a, b = s._mpq_
-            c, d = t._mpq_
-            return op(a*d, b*c)
-        return NotImplementedError
-
-    def __lt__(s, t): return s._cmp(t, operator.lt)
-    def __le__(s, t): return s._cmp(t, operator.le)
-    def __gt__(s, t): return s._cmp(t, operator.gt)
-    def __ge__(s, t): return s._cmp(t, operator.ge)
-
-    def __abs__(s):
-        a, b = s._mpq_
-        if a >= 0:
-            return s
-        v = new(mpq)
-        v._mpq_ = -a, b
-        return v
-
-    def __neg__(s):
-        a, b = s._mpq_
-        v = new(mpq)
-        v._mpq_ = -a, b
-        return v
-
-    def __pos__(s):
-        return s
-
-    def __add__(s, t):
-        ttype = type(t)
-        if ttype is mpq:
-            a, b = s._mpq_
-            c, d = t._mpq_
-            return create_reduced(a*d+b*c, b*d)
-        if ttype in int_types:
-            a, b = s._mpq_
-            v = new(mpq)
-            v._mpq_ = a+b*t, b
-            return v
-        return NotImplemented
-
-    __radd__ = __add__
-
-    def __sub__(s, t):
-        ttype = type(t)
-        if ttype is mpq:
-            a, b = s._mpq_
-            c, d = t._mpq_
-            return create_reduced(a*d-b*c, b*d)
-        if ttype in int_types:
-            a, b = s._mpq_
-            v = new(mpq)
-            v._mpq_ = a-b*t, b
-            return v
-        return NotImplemented
-
-    def __rsub__(s, t):
-        ttype = type(t)
-        if ttype is mpq:
-            a, b = s._mpq_
-            c, d = t._mpq_
-            return create_reduced(b*c-a*d, b*d)
-        if ttype in int_types:
-            a, b = s._mpq_
-            v = new(mpq)
-            v._mpq_ = b*t-a, b
-            return v
-        return NotImplemented
-
-    def __mul__(s, t):
-        ttype = type(t)
-        if ttype is mpq:
-            a, b = s._mpq_
-            c, d = t._mpq_
-            return create_reduced(a*c, b*d)
-        if ttype in int_types:
-            a, b = s._mpq_
-            return create_reduced(a*t, b)
-        return NotImplemented
-
-    __rmul__ = __mul__
-
-    def __div__(s, t):
-        ttype = type(t)
-        if ttype is mpq:
-            a, b = s._mpq_
-            c, d = t._mpq_
-            return create_reduced(a*d, b*c)
-        if ttype in int_types:
-            a, b = s._mpq_
-            return create_reduced(a, b*t)
-        return NotImplemented
-
-    def __rdiv__(s, t):
-        ttype = type(t)
-        if ttype is mpq:
-            a, b = s._mpq_
-            c, d = t._mpq_
-            return create_reduced(b*c, a*d)
-        if ttype in int_types:
-            a, b = s._mpq_
-            return create_reduced(b*t, a)
-        return NotImplemented
-
-    def __pow__(s, t):
-        ttype = type(t)
-        if ttype in int_types:
-            a, b = s._mpq_
-            if t:
-                if t < 0:
-                    a, b, t = b, a, -t
-                v = new(mpq)
-                v._mpq_ = a**t, b**t
-                return v
-            raise ZeroDivisionError
-        return NotImplemented
-
-
-mpq_1 = mpq((1,1))
-mpq_0 = mpq((0,1))
-mpq_1_2 = mpq((1,2))
-mpq_3_2 = mpq((3,2))
-mpq_1_4 = mpq((1,4))
-mpq_1_16 = mpq((1,16))
-mpq_3_16 = mpq((3,16))
-mpq_5_2 = mpq((5,2))
-mpq_3_4 = mpq((3,4))
-mpq_7_4 = mpq((7,4))
-mpq_5_4 = mpq((5,4))
-
-
-# Register with "numbers" ABC
-#   We do not subclass, hence we do not use the @abstractmethod checks. While
-#   this is less invasive it may turn out that we do not actually support
-#   parts of the expected interfaces.  See
-#   https://docs.python.org/3/library/numbers.html for list of abstract methods.
-numbers.Rational.register(mpq)


Test patch:
diff --git a/mpmath/tests/test_basic_ops.py b/mpmath/tests/test_basic_ops.py
index 3b00220..e3daa1a 100644
--- a/mpmath/tests/test_basic_ops.py
+++ b/mpmath/tests/test_basic_ops.py
@@ -6,10 +6,9 @@
 from mpmath import (ceil, fadd, fdiv, floor, fmul, fneg, frac, fsub, inf,
                     isinf, isint, isnan, isnormal, monitor, mp, mpc, mpf, mpi,
                     nan, ninf, nint, nint_distance, pi)
-from mpmath.libmp import (finf, fnan, fninf, fone, from_float, from_int,
+from mpmath.libmp import (MPQ, finf, fnan, fninf, fone, from_float, from_int,
                           from_str, mpf_add, mpf_mul, mpf_sub, round_down,
                           round_nearest, round_up, to_int)
-from mpmath.rational import mpq
 
 
 def test_type_compare():
@@ -180,14 +179,14 @@ def test_hash():
     # Check that overflow doesn't assign equal hashes to large numbers
     assert hash(mpf('1e1000')) != hash('1e10000')
     assert hash(mpc(100,'1e1000')) != hash(mpc(200,'1e1000'))
-    assert hash(mp.mpq(1,3))
-    assert hash(mp.mpq(0,1)) == 0
-    assert hash(mp.mpq(-1,1)) == hash(-1)
-    assert hash(mp.mpq(1,1)) == hash(1)
-    assert hash(mp.mpq(5,1)) == hash(5)
-    assert hash(mp.mpq(1,2)) == hash(0.5)
+    assert hash(MPQ(1,3))
+    assert hash(MPQ(0,1)) == 0
+    assert hash(MPQ(-1,1)) == hash(-1)
+    assert hash(MPQ(1,1)) == hash(1)
+    assert hash(MPQ(5,1)) == hash(5)
+    assert hash(MPQ(1,2)) == hash(0.5)
     assert hash(mpf(1)*2**2000) == hash(2**2000)
-    assert hash(mpf(1)/2**2000) == hash(mpq(1,2**2000))
+    assert hash(mpf(1)/2**2000) == hash(MPQ(1,2**2000))
 
 # Advanced rounding test
 def test_add_rounding():
@@ -382,8 +381,8 @@ def test_isnan_etc():
     assert isnan(mpc(2, 2)) is False
     assert isnan(mpc(nan, inf)) is True
     assert isnan(mpc(inf, inf)) is False
-    assert isnan(mpq((3, 2))) is False
-    assert isnan(mpq((0, 1))) is False
+    assert isnan(MPQ(3, 2)) is False
+    assert isnan(MPQ(0, 1)) is False
     assert isinf(inf) is True
     assert isinf(-inf) is True
     assert isinf(3) is False
@@ -396,8 +395,8 @@ def test_isnan_etc():
     assert isinf(mpc(nan, inf)) is True
     assert isinf(mpc(inf, nan)) is True
     assert isinf(mpc(nan, nan)) is False
-    assert isinf(mpq((3, 2))) is False
-    assert isinf(mpq((0, 1))) is False
+    assert isinf(MPQ(3, 2)) is False
+    assert isinf(MPQ(0, 1)) is False
     pytest.raises(TypeError, lambda: isinf(object()))
     assert isnormal(3) is True
     assert isnormal(3.5) is True
@@ -426,8 +425,8 @@ def test_isnan_etc():
     assert isnormal(mpc(nan, inf)) is False
     assert isnormal(mpc(nan, nan)) is False
     assert isnormal(mpc(inf, inf)) is False
-    assert isnormal(mpq((3, 2))) is True
-    assert isnormal(mpq((0, 1))) is False
+    assert isnormal(MPQ(3, 2)) is True
+    assert isnormal(MPQ(0, 1)) is False
     pytest.raises(TypeError, lambda: isnormal(object()))
     assert isint(3) is True
     assert isint(0) is True
@@ -454,23 +453,23 @@ def test_isnan_etc():
     assert isint(3 + 4j) is False
     assert isint(3 + 4j, gaussian=True) is True
     assert isint(3 + 0j) is True
-    assert isint(mpq((3, 2))) is False
-    assert isint(mpq((3, 9))) is False
-    assert isint(mpq((9, 3))) is True
-    assert isint(mpq((0, 4))) is True
-    assert isint(mpq((1, 1))) is True
-    assert isint(mpq((-1, 1))) is True
+    assert isint(MPQ(3, 2)) is False
+    assert isint(MPQ(3, 9)) is False
+    assert isint(MPQ(9, 3)) is True
+    assert isint(MPQ(0, 4)) is True
+    assert isint(MPQ(1, 1)) is True
+    assert isint(MPQ(-1, 1)) is True
     pytest.raises(TypeError, lambda: isint(object()))
     assert mp.isnpint(0) is True
     assert mp.isnpint(1) is False
     assert mp.isnpint(-1) is True
     assert mp.isnpint(-1.1) is False
     assert mp.isnpint(-1.0) is True
-    assert mp.isnpint(mp.mpq(1, 2)) is False
-    assert mp.isnpint(mp.mpq(-1, 2)) is False
-    assert mp.isnpint(mp.mpq(-3, 1)) is True
-    assert mp.isnpint(mp.mpq(0, 1)) is True
-    assert mp.isnpint(mp.mpq(1, 1)) is False
+    assert mp.isnpint(MPQ(1, 2)) is False
+    assert mp.isnpint(MPQ(-1, 2)) is False
+    assert mp.isnpint(MPQ(-3, 1)) is True
+    assert mp.isnpint(MPQ(0, 1)) is True
+    assert mp.isnpint(MPQ(1, 1)) is False
     assert mp.isnpint(0 + 0j) is True
     assert mp.isnpint(-1 + 0j) is True
     assert mp.isnpint(-1.1 + 0j) is False
@@ -485,6 +484,6 @@ def test_issue_438():
 
 
 def test_ctx_mag():
-    assert mp.mag(mpq(1, 2)) == 0
-    assert mp.mag(mpq(2)) == 2
-    assert mp.mag(mpq(0)) == mpf('-inf')
+    assert mp.mag(MPQ(1, 2)) == 0
+    assert mp.mag(MPQ(2)) == 2
+    assert mp.mag(MPQ(0)) == mpf('-inf')