Counting Divisors of a number with restrictions on the divisors

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



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)

• Is u <= v? If not, it will "double-count" each not-repeating (commutative) pair of unit fractions, e.g. $$\frac{2}{7} = \frac{1}{4} + \frac{1}{28} = \frac{1}{28} + \frac{1}{4}$$ making u,v have a pair of values [4,28], [28,4] that's counted twice. – Graham Aug 26 '18 at 13:22
• @graham Good question! I'll make an edit to the post with this information. When I call the function I guarantee $$lower \lt upper \le q$$ which should ensure there are no duplicate solutions. – spyr03 Aug 26 '18 at 13:44
• Also, I think you need to add from itertools import product to the top of your code. – Graham Aug 26 '18 at 13:47
• @Graham To be exact, I would set upper = 2 * (7 / 2) = 7, since with 2 fractions the bigger of the two could not be any smaller and still form solutions. So I would call the function like this: factor_m(2, 7, lower, 7) From this, end = 2 * 7 - 7 = 7, so I would check for values of u between start (inclusive) and 7 (exclusive). So I would want to find for the above exactly one solution, the first one. – spyr03 Aug 26 '18 at 14:08

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:

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

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.

• I tried looping until d*d, but I couldn't measure any noticeable difference. Splitting the code out into a function actually worsens the time on my machine. I'll edit in cProfile metrics. – spyr03 Aug 27 '18 at 19:48
• I was tweaking with the code trying to get a better understanding of it before being able to find micro optimisation, better algorithm or mathematical insight but I've thought a lot about this and was unable to find a real way to improve things so far.... – SylvainD Aug 27 '18 at 21:20
• @spyr03 I think things should be faster now. Please let me know... – SylvainD Aug 28 '18 at 16:00