Counting Divisors of a number with restrictions on the divisors

Question

Given a fraction p/q, I want to count the number of unit fraction solutions

$$\frac{p}{q} = \frac{1}{u} + \frac{1}{v}$$

with some contraints

$$lower \le u \le upper$$

By manipulating the first equation and completing a square you arrive at a formula

$$(pu - q)(pv - q) = q^2$$

So the number of solutions is related to the divisors of q squared. If we only look at the smaller of two divisors (which we will call a), it is a solution if it produces a valid unit fraction (1 / u).

$$a = pu - q \implies u = \frac{a + q}{p}$$

We are interested in a divisor a iff

$$p*lower - q \le a \le p * upper - q, \space a \equiv -q \mod p$$

---

To compute the answer I first get the prime factorisation of q, then I double each power to get the prime factorisation of q squared. Using this I enumerate over all the possible divisors, computing each one by trying all possible combination of powers. When I generate one I check if it matches the above requirements.

def factor_m(p, q, lower, upper):
    num = q
    factors = dict()
    count = 0
    while num % 2 == 0:
        num //= 2
        count += 1
    if count:
        factors[2] = count

    current = 3
    while num > 1:
        count = 0
        while num % current == 0:
            num //= current
            count += 1
        if count:
            factors[current] = count
        current += 2

    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 <= number <= end and number % p == mod:
            total += 1

    return total

I'm calling this function a lot with reasonably big values of q (around 10**14 currently) so anything that saves time here is a massive boost.

Examples:

>>> factor_m(2, 7, 0, 50)  # 3
>>> factor_m(5, 1775025265104, 355005053021, 710010106041)  # 4101
>>> factor_m(737, 1046035200, 1926400, 2838630)  # 1
>>> factor_m(105467, 1231689911361, 11678439, 23356877)  # 0

---

EDIT To guarantee all solutions are distinct I only call the function with values that satisfy the constraint

$$lower \lt upper \le q$$

---

EDIT2 I ran the full program through cProfile along with the program with modifications suggested by Josay. Unfortunately it is slower, I suspect because it almost doubles the number of function calls.

current_code - top 3 functions by time

672574 function calls (614969 primitive calls) in 19.526 seconds

 ncalls  tottime  percall  cumtime  percall filename:lineno(function)
  53273   18.818    0.000   18.893    0.000 fast.py:40(factor_m)
   2957    0.295    0.000    0.298    0.000 fast.py:8(factor_p1)
57606/1    0.149    0.000   19.525   19.525 fast.py:75(f)

With Josay's changes - top 3 functions by time

1163833 function calls (1106228 primitive calls) in 19.835 seconds

ncalls  tottime  percall  cumtime  percall filename:lineno(function) 
53273   18.427    0.000   19.160    0.000 fast.py:61(factor_m)
53273    0.500    0.000    0.650    0.000 fast.py:42(get_factors)
 2957    0.297    0.000    0.301    0.000 fast.py:8(factor_v2)

Answer 1

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 > 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 <= u <= 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 <= number <= 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 <= n. Thus we can stop the search earlier (and the remaining part is known to be a prime).

while n >= d * d:
    while n % d == 0:
        n //= d
        factors[d] += 1
    d += 2
if n >= 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 >= d * d:
        while n % d == 0:
            n //= d
            yield d
        d += 2
    if n >= 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 <= u <= 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 <= a <= 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 <= a <= 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.