Optimize math.lcm(*xs)

[pochmann]

Benchmark for computing the LCM of 1000 random numbers from 1 to 1,000,000:

6.19 ms ± 0.09 ms  current_Python
6.14 ms ± 0.12 ms  current_C
3.46 ms ± 0.16 ms  proposal_Python

Code:

def current_C(xs):
    return lcm(*xs)

def current_Python(xs):
    # re-implementation of the current C implementation
    res = 1
    for x in xs:
       res = lcm(res, x)
    return res

def proposal_Python(xs):
    res = 1
    for x in xs:
        res = lcm(x, res)
    return res

full benchmark code

Attempt This Online!

from math import lcm
from random import randint
from timeit import repeat
from statistics import mean, stdev

def current_C(xs):
    return lcm(*xs)

def current_Python(xs):
    res = 1
    for x in xs:
       res = lcm(res, x)
    return res

def proposal_Python(xs):
    res = 1
    for x in xs:
        res = lcm(x, res)
    return res

funcs = current_C, current_Python, proposal_Python

times = {f: [] for f in funcs}
def stats(f):
    ts = [t * 1e3 for t in times[f][10:]]
    return f'{mean(ts):4.2f} ms ± {stdev(ts):4.2f} ms '

for _ in range(20):
    xs = [randint(1, 10**6) for _ in range(1000)]
    for f in funcs:
        t = min(repeat(lambda: f(xs), number=1))
        times[f].append(t)

for f in sorted(funcs, key=stats, reverse=True):
    print(stats(f), f.__name__)

The difference is using lcm(x, res) instead of lcm(res, x). When computing the LCM of multiple/many values, res tends to be(come) larger than x. And then lcm(x, res) is faster than lcm(res, x).

Why? Because lcm(a, b) is computed as (a // g) * b, where g = gcd(a, b). Let A, B and G be the respective bit lengths, and let's pretend there's no Karatsuba. Then a // g takes AG time and results in a number with A-G bits. Multiplying that with b takes (A-G)B time. So overall AG + (A-G)B = AG + AB - BG = AB + (A-B)G time. So it's better to have B > A.

The simple optimization would be to swap the arguments here:

cpython/Modules/mathmodule.c

Line 786 in 81bf10e

Py_SETREF(res, long_lcm(res, x));

In addition to this "likely we have res > x" rule of thumb, the two-argument lcm could check which one is longer and swap them if b is shorter. And then call gcd(b, a) instead of gcd(a, b) so that gcd doesn't have to swap them right back.

Linked PRs

• gh-102221: speed up math.lcm by swapping numbers #111450
• gh-102221: Optimize math.lcm() for multiple arguments (recursive) #135554
• gh-102221: Optimize math.lcm() for multiple arguments (non-recursive) #135609

[terryjreedy]

@mdickinson @rhettinger any comment?

[OTheDev]

Blindly incorporating your suggestion (i.e. using long_lcm(x, res)) instead of long_lcm(res, x)) and running your script

main:

2.16 ms ± 0.04 ms  current_Python
2.15 ms ± 0.04 ms  current_C
1.21 ms ± 0.02 ms  proposal_Python

pochmann:

2.18 ms ± 0.06 ms  proposal_Python
1.21 ms ± 0.03 ms  current_Python
1.19 ms ± 0.03 ms  current_C

long_lcm(x, res) if res > x, long_lcm(res, x) otherwise:

1.24 ms ± 0.02 ms  current_Python
1.24 ms ± 0.01 ms  proposal_Python
1.20 ms ± 0.02 ms  current_C

[mdickinson]

the [two-argument lcm]could check which one is longer and swap them if b is shorter [...]

Sounds reasonable, provided that it's clearly commented why we're doing the compare and swap in the code.

[mdickinson]

So overall AG + (A-G)B = AG + AB - BG = AB + (A-B)G time.

I think this analysis is a bit too simple in that it neglects that division is likely to have a larger constant factor than multiplication. But I believe the conclusion holds regardless. Suppose AB is the time to multiply an A-bit number by a B-bit number (in whatever units are appropriate), and kAB is the time to divide an A-bit number by a B-bit number. Then we're looking at total time AB + (kA-B)G, so we end up comparing kA - B with kB - A and again conclude that B > A is better.

[pochmann]

Yes, it's both intentionally and unintentionally inaccurate, but I also believe the conclusion holds regardless. The benchmark agrees, and I also tried variations of it (other number sizes and more/fewer numbers), and proposal_Python was always faster than current_Python and almost always faster than current_C (which was faster for very short numbers, where the LCM is short and thus it's only the language speed that matters).

Intuitively, if we do (a // g) * b with a being the shorter one, then we benefit from the shorter one being involved in both operations and the longer one being involved only in the second operation.

[kevin1kevin1k]

Hi everyone, I am trying to contribute to this repo and just stumbled on this issue.
It seems that all of you have reached the same conclusion that the modification is acceptable, given the intuition as well as the benchmark.
If so I can create an PR for this issue; my question is whether new tests or benchmarks should be added (here in the discussion, or into the codebase)
Thanks!

[terryjreedy]

If lcm is properly tested already, then the existing tests should be sufficient to test an altered implementation. That is the purpose of regression tests. My feeling is that the existing speed tests should be sufficient, but note his requirement in #102221 (comment).

[pochmann]

Sorry I didn't get to this yet, still only have a phone, where working on the C code is nearly impossible. Happy if someone else gets it done. I see the PR attempts the second optimization, comparing the numbers and swapping accordingly. But it compares the exact numbers. I'd likely only compare their sizes, I think that matters most and is faster. That's why I spoke of "shorter"/"longer" above. But I'd benchmark both ways.