Counting Divisors of a number with restrictions on the divisors
Code organisation

In order to improve the code performances, if could be a good idea to make the easier to understand and easier to change.

You could:


split your code into smaller functions
add documentation
add tests
reuse existing data structures such as collections.Counter.


We'd get something along the lines of:

from itertools import product
from collections import Counter
import time

def get_factors(n):
    """Get factor for integer n

    Return the factorisation of n as a dictionnary
    mapping prime factors to their exponents."""
    factors = Counter()
    d = 2
    while n % d == 0:
        n //= d
        factors[d] += 1
    d = 3
    while n &gt; 1:
        while n % d == 0:
            n //= d
            factors[d] += 1
        d += 2
    return factors


def factor_m(p, q, lower, upper):
    """Return the number of unit fractions solutions to:
       p   1   1
       - = - + - with lower &lt;= u &lt;= upper
       q   u   v

    The equation is equivalent to ..."""
    factors = get_factors(q)
    start, end = p * lower - q, p * upper - q
    mod = -q % p
    total = 0
    for powers in product(*(range((2 * v) + 1) for v in factors.values())):
        number = 1
        for b, _power in zip(factors, powers):
            number *= b ** _power

        if start &lt;= number &lt;= end and number % p == mod:
            total += 1

    return total

t = time.time()
assert factor_m(2, 7, 0, 50) == 3
assert factor_m(5, 1775025265104, 355005053021, 710010106041) == 4101
assert factor_m(737, 1046035200, 1926400, 2838630) == 1
assert factor_m(105467, 1231689911361, 11678439, 23356877) == 0
print(time.time() - t)


Simple mathematical optimisation

A trivial optimisation can be added to the logic performing the prime decomposition: for all divisors d, d * d &lt;= n. Thus we can stop the search earlier (and the remaining part is known to be a prime).

while n &gt;= d * d:
    while n % d == 0:
        n //= d
        factors[d] += 1
    d += 2
if n &gt;= 1:    # to avoid 1 as a factor
    factors[n] += 1


This makes the code only very slightly faster as the bottleneck is not in that part of the code.

Limit the number of operations performed/limit accesses to factors

When iterating over factors, you get the values to create range of exponents (0, 1, 2, ...) and then, in the loop, you zip this to the keys to generate the prime**exponent value.

You could feed into the product function, the values prime**exponent already computed so that you limit the number of operations (** mostly) and you access factors in a single place.

This would give something like:

def get_divisors_for_square_of(n):
    """"Generate all divisors for n^2.

    Use the factorisation of n and double the exponents."""
    factors = get_factors(q)
    powers_list = [[p**v for v in range(2 * v + 1)] for p, v in factors.items()]
    for powers in product(*powers_list):
       number = 1
       for p in powers:
           number *= p
       yield number


On my setup, this halves the computing time.

Final code

import itertools
import collections

def get_factors(n):
    """Generate factors for integer n."""
    d = 2
    while n % d == 0:
        n //= d
        yield d
    d = 3
    while n &gt;= d * d:
        while n % d == 0:
            n //= d
            yield d
        d += 2
    if n &gt;= 1:    # to avoid 1 as a factor
        yield n


def get_divisors_for_square_of(n):
    """"Generate all divisors for n^2.

    Use the factorisation of n and double the exponents."""
    factors = collections.Counter(get_factors(n))
    powers_list = [[p**v for v in range(2 * v + 1)] for p, v in factors.items()]
    for powers in itertools.product(*powers_list):
       number = 1
       for p in powers:
           number *= p
       yield number


def factor_m(p, q, lower, upper):
    """Return the number of unit fractions solutions to:
       p   1   1
       - = - + - with lower &lt;= u &lt;= upper
       q   u   v

    The equation is equivalent to:
        puv = q(u + v)
    or
        (pu-q) (pv-q) = q^2
    which boils down to finding divisors of q.
    Let a = pu - q, then u = (a + q) / p when it gives an integer
    which is when (a + q ) mod p == 0.
    Also, we want p * lower - q &lt;= a &lt;= p * upper - q.
    ....
    Solution for 2/7 are (u, v) = [(4, 28), (7, 7), (28, 4)].
    """
    start, end = p * lower - q, p * upper - q
    return sum(1
               for a in get_divisors_for_square_of(q)
               if start &lt;= a &lt;= end and (a + q) % p == 0)


if __name__ == '__main__':
    import time
    t = time.time()
    for i in range(5):
        assert factor_m(2, 7, 0, 50) == 3
        assert factor_m(2, 21, 0, 500) == 9
        assert factor_m(5, 1775025265104, 355005053021, 710010106041) == 4101
        assert factor_m(737, 1046035200, 1926400, 2838630) == 1
        assert factor_m(105467, 1231689911361, 11678439, 23356877) == 0
    print(time.time() - t)


I'll update the solution if I find other things.