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)
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. \$\endgroup\$from itertools import product
to the top of your code. \$\endgroup\$