From: Sergey B Kirpichev <skirpichev@gmail.com>
Date: Sun, 14 Mar 2021 03:39:21 +0000 (+0300)
Subject: bpo-43420: Simple optimizations for Fraction's arithmetics
X-Git-Url: http://git.ipfire.org/gitweb.cgi?a=commitdiff_plain;h=430e4764cb554ab4c3fef16990dc55186a18b771;p=thirdparty%2FPython%2Fcpython.git

bpo-43420: Simple optimizations for Fraction's arithmetics

making Fraction class more usable for large arguments (>~ 10**6),
with cost 10-20% for small components case.

Before:
$ ./python -m timeit -s 'from fractions import Fraction as F' \
-s 'a=[F(1, _**3) for _ in range(1, 1000)]' 'sum(a)'
5 loops, best of 5: 81.2 msec per loop

After:
$ ./python -m timeit -s 'from fractions import Fraction as F' \
-s 'a=[F(1, _**3) for _ in range(1, 1000)]' 'sum(a)'
10 loops, best of 5: 23 msec per loop

References:
Knuth, TAOCP, Volume 2, 4.5.1,
https://www.eecis.udel.edu/~saunders/courses/822/98f/collins-notes/rnarith.ps,
https://gmplib.org/ (e.g. https://gmplib.org/repo/gmp/file/tip/mpq/aors.c)
---

diff --git a/Lib/fractions.py b/Lib/fractions.py
index a9fa62dba7b4..9937b960ba84 100644
--- a/Lib/fractions.py
+++ b/Lib/fractions.py
@@ -381,10 +381,79 @@ class Fraction(numbers.Rational):
 
         return forward, reverse
 
+    # Rational arithmetic algorithms: Knuth, TAOCP, Volume 2, 4.5.1.
+    #
+    # Assume input fractions a and b are normalized.
+    #
+    # 1) Consider addition/substraction.
+    #
+    # Let g = gcd(da, db). Then
+    #
+    #              na   nb    na*db ± nb*da
+    #     a ± b == -- ± -- == ------------- ==
+    #              da   db        da*db
+    #
+    #              na*(db//g) ± nb*(da//g)    t
+    #           == ----------------------- == -
+    #                      (da*db)//g         d
+    #
+    # Now, if g > 1, we're working with smaller integers.
+    #
+    # Note, that t, (da//g) and (db//g) are pairwise coprime.
+    #
+    # Indeed, (da//g) and (db//g) share no common factors (they were
+    # removed) and da is coprime with na (since input fractions are
+    # normalized), hence (da//g) and na are coprime.  By symmetry,
+    # (db//g) and nb are coprime too.  Then,
+    #
+    #     gcd(t, da//g) == gcd(na*(db//g), da//g) == 1
+    #     gcd(t, db//g) == gcd(nb*(da//g), db//g) == 1
+    #
+    # Above allows us optimize reduction of the result to lowest
+    # terms.  Indeed,
+    #
+    #     g2 = gcd(t, d) == gcd(t, (da//g)*(db//g)*g) == gcd(t, g)
+    #
+    #                       t//g2                   t//g2
+    #     a ± b == ----------------------- == ----------------
+    #              (da//g)*(db//g)*(g//g2)    (da//g)*(db//g2)
+    #
+    # is a normalized fraction.  This is useful because the unnormalized
+    # denominator d could be much larger than g.
+    #
+    # We should special-case g == 1, since 60.8% of randomly-chosen
+    # integers are coprime:
+    # https://en.wikipedia.org/wiki/Coprime_integers#Probability_of_coprimality
+    #
+    # 2) Consider multiplication
+    #
+    # Let g1 = gcd(na, db) and g2 = gcd(nb, da), then
+    #
+    #            na*nb    na*nb    (na//g1)*(nb//g2)
+    #     a*b == ----- == ----- == -----------------
+    #            da*db    db*da    (db//g1)*(da//g2)
+    #
+    # Note, that after divisions we're multiplying smaller integers.
+    #
+    # Also, the resulting fraction is normalized, because each of
+    # two factors in the numerator is coprime to each of the two factors
+    # in the denominator.
+    #
+    # Indeed, pick (na//g1).  It's coprime with (da//g2), because input
+    # fractions are normalized.  It's also coprime with (db//g1), because
+    # common factors are removed by g1 == gcd(na, db).
+
     def _add_sub_(a, b, pm=int.__add__):
-        da, db = a.denominator, b.denominator
-        return Fraction(pm(a.numerator * db, b.numerator * da),
-                        da * db)
+        na, da = a.numerator, a.denominator
+        nb, db = b.numerator, b.denominator
+        g = math.gcd(da, db)
+        if g == 1:
+            return Fraction(pm(na * db, da * nb), da * db, _normalize=False)
+        else:
+            s = da // g
+            t = pm(na * (db // g), nb * s)
+            g2 = math.gcd(t, g)
+            return Fraction(t // g2, s * (db // g2), _normalize=False)
 
     _add = functools.partial(_add_sub_)
     _add.__doc__ = 'a + b'
@@ -396,14 +465,26 @@ class Fraction(numbers.Rational):
 
     def _mul(a, b):
         """a * b"""
-        return Fraction(a.numerator * b.numerator, a.denominator * b.denominator)
+        na, da = a.numerator, a.denominator
+        nb, db = b.numerator, b.denominator
+        g1 = math.gcd(na, db)
+        g2 = math.gcd(nb, da)
+        return Fraction((na // g1) * (nb // g2),
+                        (db // g1) * (da // g2), _normalize=False)
 
     __mul__, __rmul__ = _operator_fallbacks(_mul, operator.mul)
 
     def _div(a, b):
         """a / b"""
-        return Fraction(a.numerator * b.denominator,
-                        a.denominator * b.numerator)
+        # Same as _mul(), with inversed b.
+        na, da = a.numerator, a.denominator
+        nb, db = b.numerator, b.denominator
+        g1 = math.gcd(na, nb)
+        g2 = math.gcd(db, da)
+        n, d = (na // g1) * (db // g2), (nb // g1) * (da // g2)
+        if nb < 0:
+            n, d = -n, -d
+        return Fraction(n, d, _normalize=False)
 
     __truediv__, __rtruediv__ = _operator_fallbacks(_div, operator.truediv)
 
diff --git a/Misc/NEWS.d/next/Library/2021-03-07-08-03-31.bpo-43420.cee_X5.rst b/Misc/NEWS.d/next/Library/2021-03-07-08-03-31.bpo-43420.cee_X5.rst
new file mode 100644
index 000000000000..02025813a457
--- /dev/null
+++ b/Misc/NEWS.d/next/Library/2021-03-07-08-03-31.bpo-43420.cee_X5.rst
@@ -0,0 +1,2 @@
+Improve performance of class:`fractions.Fraction` arithmetics for large
+components.