Вот простенькая программа подсчёта количества решений уравнения 1/x + 1/y = 1/n в натуральных числах для разных n. Суть в том, что в 11-12 строках проверяется, можно ли выполнить деление нацело, и если да, то мы имеем новое решение.

Я подумал, что для подсчёта количества решений необязательно вычислять значение y, и заменил 11-12 строки на строку:
Не ожидая ничего особенного от такой модификации я обнаружил, что скорость работы возросла почти в два раза, что подтверждает отладочный вывод первой версии:

n = 1000, solutions: 25, time is 0.000
n = 2000, solutions: 32, time is 1.448
n = 3000, solutions: 74, time is 3.853
n = 4000, solutions: 39, time is 7.223
n = 5000, solutions: 32, time is 11.565
n = 6000, solutions: 95, time is 16.873
n = 7000, solutions: 74, time is 23.113
n = 8000, solutions: 46, time is 30.325
n = 9000, solutions: 123, time is 38.495
n = 10000, solutions: 41, time is 47.625
n = 11000, solutions: 74, time is 57.705
n = 12000, solutions: 116, time is 68.748

и второй:

n = 1000, solutions: 25, time is 0.000
n = 2000, solutions: 32, time is 0.783
n = 3000, solutions: 74, time is 2.078
n = 4000, solutions: 39, time is 3.908
n = 5000, solutions: 32, time is 6.263
n = 6000, solutions: 95, time is 9.140
n = 7000, solutions: 74, time is 12.548
n = 8000, solutions: 46, time is 16.473
n = 9000, solutions: 123, time is 20.935
n = 10000, solutions: 41, time is 25.910
n = 11000, solutions: 74, time is 31.405
n = 12000, solutions: 116, time is 37.428
n = 13000, solutions: 74, time is 43.965
n = 14000, solutions: 95, time is 51.035
n = 15000, solutions: 95, time is 58.618
n = 16000, solutions: 53, time is 66.713

Почему так произошло?

Из документации:

divmod(a, b)
Take two (non complex) numbers as arguments and return a pair of numbers consisting of their quotient and remainder when using long division. With mixed operand types, the rules for binary arithmetic operators apply. For plain and long integers, the result is the same as (a // b, a % b). For floating point numbers the result is (q, a % b), where q is usually math.floor(a / b) but may be 1 less than that. In any case q * b + a % b is very close to a, if a % b is non-zero it has the same sign as b, and 0 <= abs(a % b) < abs(b).

Обратили внимание на выделенное жирным шрифтом? divmod делает минимум 2 действия, отсюда и разница по времени около в 2 раза.

Нуу, вообще-то здесь сказано, что результат лишь будет численно равен данному кортежу. Разве не для того существует эта функция, чтобы за одну операцию деления получать частное и остаток? Косвенно ведь это подтверждается определением ассемблерной команды деления div.