Optimize rational arithmetics

[skirpichev]

BPO	43420
Nosy	@tim-one, @rhettinger, @terryjreedy, @mdickinson, @scoder, @skirpichev, @asmeurer
PRs	bpo-43420: Simple optimizations for Fraction's arithmetics #24779
Files	f.patch

^{Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.}

Show more details

GitHub fields:

assignee = None
closed_at = <Date 2021-03-22.02:35:59.048>
created_at = <Date 2021-03-06.14:18:17.372>
labels = ['type-feature', 'library', '3.10']
title = 'Optimize rational arithmetics'
updated_at = <Date 2021-05-17.12:07:39.325>
user = 'https://github.com/skirpichev'

bugs.python.org fields:

activity = <Date 2021-05-17.12:07:39.325>
actor = 'Sergey.Kirpichev'
assignee = 'none'
closed = True
closed_date = <Date 2021-03-22.02:35:59.048>
closer = 'tim.peters'
components = ['Library (Lib)']
creation = <Date 2021-03-06.14:18:17.372>
creator = 'Sergey.Kirpichev'
dependencies = []
files = ['49854']
hgrepos = []
issue_num = 43420
keywords = ['patch']
message_count = 27.0
messages = ['388200', '388210', '388212', '388214', '388215', '388221', '388227', '388232', '388235', '388236', '388249', '388256', '388260', '388328', '388415', '388477', '388563', '388590', '388591', '389270', '389271', '389274', '389276', '389277', '389342', '393807', '393808']
nosy_count = 7.0
nosy_names = ['tim.peters', 'rhettinger', 'terry.reedy', 'mark.dickinson', 'scoder', 'Sergey.Kirpichev', 'asmeurer']
pr_nums = ['24779']
priority = 'normal'
resolution = 'fixed'
stage = 'resolved'
status = 'closed'
superseder = None
type = 'enhancement'
url = 'https://bugs.python.org/issue43420'
versions = ['Python 3.10']

[skirpichev]

fractions.py uses naive algorithms for doing arithmetics.

It may worth implementing less trivial versions for addtion/substraction
and multiplication (e.g. Henrici algorithm and so on), described here:
https://www.eecis.udel.edu/~saunders/courses/822/98f/collins-notes/rnarith.ps
as e.g gmplib does: https://gmplib.org/repo/gmp/file/tip/mpq/aors.c

Some projects (e.g. SymPy here: sympy/sympy#12656) reinvent
the stdlib's Fraction just to add such simple improvements. With big denominators (~10**6)
this really does matter, my local benchmarks suggest the order of magnitude difference for
summation of several such numbers.

[rhettinger]

Here's some code to try out:

from math import gcd
from fractions import Fraction
import operator
import math

class Henrici(Fraction):
    'Reformulate _mul to reduce the size of intermediate products'

# Original has 2 multiplications, 1 gcd calls, and 2 divisions
# This one has 2 multiplications, 2 gcd calls, and 4 divisions  

    def _mul(a, b):
        a_n, a_d = a.numerator, a.denominator
        b_n, b_d = b.numerator, b.denominator
        d1 = math.gcd(a_n, b_d)
        a_n //= d1
        b_d //= d1
        d2 = math.gcd(b_n, a_d)
        b_n //= d2
        a_d //= d2
        result = Fraction(a_n * b_n, a_d * b_d, _normalize=False)
        assert math.gcd(a_n * b_n, a_d * b_d) == 1 and a_d * b_d >= 0
        return result

__mul__, __rmul__ = Fraction._operator_fallbacks(_mul, operator.mul)

assert Henrici(10, 3) * Henrici(6, 5) == Henrici(4, 1)

[rhettinger]

Note that Fraction arithmetic has a huge administrative overhead. The cost of the underlying multiplications and divisions won't dominate the total time until the numerators and denominators are very large.

For the proposed optimization, this implies that cost for the extra Python steps to implement the optimization will be negligible. The benefits of the optimization are similarly attenuated.

-- Update to experimentation code: add guards for the relatively prime case. --

class Henrici(Fraction):
    'Reformulate _mul to reduce the size of intermediate products'
    def _mul(a, b):
        a_n, a_d = a.numerator, a.denominator
        b_n, b_d = b.numerator, b.denominator
        d1 = math.gcd(a_n, b_d)
        if d1 > 1:
            a_n //= d1
            b_d //= d1
        d2 = math.gcd(b_n, a_d)
        if d2 > 1:
            b_n //= d2
            a_d //= d2
        result = Fraction(a_n * b_n, a_d * b_d, _normalize=False)
        assert math.gcd(a_n * b_n, a_d * b_d) == 1 and a_d * b_d >= 0
        return result

[rhettinger]

The cost to the common case for small components is about 20%:

$ python3.10 -m timeit -r11 -s 'from fractions import Fraction as F' -s 'a=F(10,3)' -s 'b=F(6, 5)' 'a * b'
200000 loops, best of 11: 1.8 usec per loop

$ python3.10 -m timeit -r11 -s 'from patched import Fraction as F' -s 'a=F(10,3)' -s 'b=F(6, 5)' 'a * b'
100000 loops, best of 11: 2.14 usec per loop

[rhettinger]

Without the guards the incremental cost drops to 10%.

$ python3.10 -m timeit -r11 -s 'from patched_noguards import Fraction as F' -s 'a=F(10,3)' -s 'b=F(6, 5)' 'a * b'
100000 loops, best of 11: 2.02 usec per loop

[skirpichev]

I have similar timings (for a draft version of PR, see f.patch) as for the last comment, though the small-dens overhead seems to be bigger(~20%):
$ python3.10 -m timeit -r11 -s 'from fractions import Fraction as F' -s 'a=F(10,3)' -s 'b=F(6, 5)' 'a * b'
50000 loops, best of 11: 9.09 usec per loop
$ python3.10 -m timeit -r11 -s 'from patched import Fraction as F' -s 'a=F(10,3)' -s 'b=F(6, 5)' 'a * b'
20000 loops, best of 11: 11.2 usec per loop

On another hand, here are timings for bigger denominators:
$ python3.10 -m timeit -r11 -s 'from fractions import Fraction as F' -s 'import random' -s 'n = [random.randint(1, 1000000) for _ in range(1000)]' -s 'd = [random.randint(1, 1000000) for _ in range(1000)]' -s 'a=list(map(lambda x: F(*x), zip(n, d)))' 'sum(a)'
1 loop, best of 11: 257 msec per loop
$ ... from patched ...
10 loops, best of 11: 33.2 msec per loop

It's not so clear what "are very large" does mean, that could be defined here. BTW, 10**6 denominators are (very!) small for mentioned above use case (CAS package).

[rhettinger]

1 loop, best of 11: 257 msec per loop
$ ... from patched ...
10 loops, best of 11: 33.2 msec per loop

Looks worth it :-)

[mdickinson]

Personally, I'd prefer to keep the code simplicity, and the speed for small inputs here. Python's needs aren't the same as SymPy's needs or SAGE's needs, and not all of the fractions.Fraction use-cases involve summing lots of values with incompatible denominators.

With big denominators (~10**6) this really does matter, my local
benchmarks suggest the order of magnitude difference for
summation of several such numbers.

Could you give some idea of the crossover point for a single addition? At roughly what size numerator/denominator do we start seeing a performance benefit?

[skirpichev]

I'd prefer to keep the code simplicity

It's not going to be complicated so much and algorithms are well known,
but I see your point.

and the speed for small inputs here

Speed loss is not so big, 10-20%.

Python's needs aren't the same as SymPy's needs or SAGE's needs

So, for which needs it will serve?

Sorry, I can't suggest an application, which does use builtin
Fraction's (not sure if even SAGE uses them, as a fallback). SymPy doesn't,
for sure (but it could - it's PythonRational class uses same optimizations,
except for g == 1 branches in _add/_sub, I think).

There is one exception I've found: stdlib's statistics module uses Fraction's
in the _sum() helper, exactly in a paradigm "sum a lot of values".

not all of the fractions.Fraction use-cases involve summing lots of values with incompatible denominators.

No need for a lots of values (i.e. 1000): denominator of the sum will grow very fast, that
why modern CAS use modular GCD algorithms, for example.

Could you give some idea of the crossover point for a single addition?

$ ./python -m timeit -r11 -s 'from fractions import Fraction as F' -s 'a=F(10,31011021112)' -s 'b=F(86,11011021115)' 'a + b'
20000 loops, best of 11: 12.4 usec per loop
$ ./python -m timeit -r11 -s 'from patched import Fraction as F' -s 'a=F(10,31011021112)' -s 'b=F(86,11011021115)' 'a + b'
20000 loops, best of 11: 12.5 usec per loop