计算具有特定 属性 的整数元组
Counting integer tuples with specific property
给定两个数 M_max,M_min 计算素数的所有 r 元组 (M_1,...,M_r) 使得：
- M_min < M_i < M_max
- 每个 M_i-1 至少可以被奇素数的一个平方数整除
- lcm(M_1-1,...,M_r-1)=(M_1-1)*(M_2-1)...(M_r-1)/2^(r-1)
改写问题：k 个数 M_1,M_2,...,M_k 中有多少个 r 元组互质？
这是工作代码。在相当未优化的 Python.
# Stupid simple Euclidean primes filter.
def primes ():
found_primes = []
factor_filter = [False, False]
block_size = 2
candidate = 2
while True:
factor_filter.extend([False for _ in factor_filter])
for p in found_primes:
if 2*block_size <= p*p:
break
i = ((block_size + p - 1) // p) * p
while i < 2 * block_size:
factor_filter[i] = 1
i = i + p
block_size += block_size
while (candidate < block_size):
if 0 == factor_filter[candidate]:
yield candidate
found_primes.append(candidate)
candidate = candidate + 1
def factor_range (m, n):
factored = []
remainder = []
for i in range(n - m + 1):
factored.append({})
remainder.append(i+m)
for p in primes():
if n < p*p:
break
else:
count = 1
q = p
while q <= n:
# Find the first thing divisible by q
i = ((m + q - 1) // q) * q
while i <= n:
factored[i-m][p] = count
remainder[i-m] = remainder[i-m] // p
i += q
q = p*q
count += 1
for i in range(len(remainder)):
if 1 < remainder[i]:
factored[i][remainder[i]] = 1
return factored
def int_tuples (m, n, j):
factored = factor_range(m, n)
available = []
for i in range(len(factored)):
if 0 < i and i + m in factored[i]:
for p, count in factored[i-1].items():
if 2 < p and 1 < count:
available.append(i+m)
break
prime_divides = {}
for i in available:
if m < i:
for p, count in factored[i-1-m].items():
if 1 < count or p != 2:
if p not in prime_divides:
prime_divides[p] = set([i])
else:
prime_divides[p].add(i)
options_stack = [(-1, available, False)]
chosen = []
while 0 < len(options_stack):
#print(("...", chosen, options_stack))
i, options, is_pop = options_stack.pop()
i = i+1
if len(chosen) == j:
yield tuple(chosen)
if is_pop:
chosen.pop()
elif len(chosen) + len(options) - i < j:
# Cannot generate an answer.
if is_pop:
chosen.pop()
elif len(options) <= i:
# i reached the end.
if is_pop:
chosen.pop()
else:
options_stack.append((i, options, is_pop))
chosen.append(options[i])
options = options[i:]
for p, count in factored[options[0]-1].items():
if p == 2 and 1 == count:
continue
if p in prime_divides:
options = [option for option in options if option not in prime_divides[p]]
options_stack.append((-1, options, True))
这个想法是找到这样的元组，使得 (M-1)/2 缩减序列彼此成对相对质数。我所做的相当于递归搜索。但是我保留了显式堆栈，以便我可以通过简单地执行 yield.
来制作生成器
我的第一次尝试涉及保留位向量以提高效率。但是较重的数据结构，然后更早地确定“不要进一步搜索”的结论效果更好。
我最终还是决定尝试递归，我很高兴我这么做了。使用缓存是迈向拥有重量级数据结构以避免重复计算的又一步。
# Stupid simple Euclidean primes filter.
def primes ():
found_primes = []
factor_filter = [False, False]
block_size = 2
candidate = 2
while True:
# We get here when we need to double the length of our sieve.
# First add new elements.
factor_filter.extend([False for _ in factor_filter])
# Cross off the ones that divisible by known primes.
for p in found_primes:
if 2*block_size <= p*p:
break
i = ((block_size + p - 1) // p) * p
while i < 2 * block_size:
factor_filter[i] = True
i = i + p
# Record that we did so.
block_size += block_size
# Now that we have sieved, start returning primes.
while (candidate < block_size):
if 0 == factor_filter[candidate]:
yield candidate
found_primes.append(candidate)
candidate = candidate + 1
# Takes a range of numbers, factors all of them at once.
def factor_range (m, n):
# This will hold the factorization.
factored = []
# This will hold the remainder we have not factorized.
remainder = []
# We start with no factors and the number as remainder.
for i in range(n - m + 1):
factored.append({})
remainder.append(i+m)
# We run through the primes up to sqrt(n)
for p in primes():
if n < p*p:
break
else:
# For each power of p, we add it to the factorization and divide it form the remainder.
count = 1
q = p
while q <= n:
# Find the first thing divisible by q
i = ((m + q - 1) // q) * q
while i <= n:
factored[i-m][p] = count
remainder[i-m] = remainder[i-m] // p
i += q
q = p*q
count += 1
# Any left-over remainders must be primes.
for i in range(len(remainder)):
if 1 < remainder[i]:
factored[i][remainder[i]] = 1
return factored
# This function does all of the work.
def _int_tuples (m, n, j):
# BUG FIX HERE:
# If our minimum is odd, extend the boundary down one.
# This does not introduce new solutions, but does mean we
# factorize all of the numbers that we need.
if 1 == m // 2:
m = m-1
# Factorize all numbers.
factored = factor_range(m, n)
# available will be all of our candidate primes.
available = []
for i in range(len(factored)):
if 0 < i and i + m in factored[i]:
for p, count in factored[i-1].items():
if 2 < p and 1 < count:
available.append(i+m)
break
# Construct a data structure to answer, "which other candidate primes
# have p-1 share a non-trivial factor with this one"
prime_divides = {}
for i in available:
if m < i:
for p, count in factored[i-1-m].items():
if 1 < count or p != 2:
if p not in prime_divides:
prime_divides[p] = set([i])
else:
prime_divides[p].add(i)
# PERFORMANCE HACK: We will recurse over lists of remaining candidates.
# Putting candidates that eliminate lots of other candidates first
# makes those lists smaller and that recursion more efficient. We
# also improve the odds that we will be looking at a list we have
# already seen and already have an answer for. And, finally, the
# candidates at the end of our list will largely not filter each
# other out. So we are spending time on things that are more likely to
# give us more solutions.
count_related = {}
for i in available:
total = 0
for p, count in factored[i-1-m].items():
if 1 < count or p != 2:
total += len(prime_divides[p]) - 1
count_related[i] = total
available = sorted(available, key = lambda i: -count_related[i])
# To spend less time filtering, forget about the primes that only
# divide one thing.
for p in list(prime_divides.keys()):
if 1 == len(prime_divides[p]):
i = list(prime_divides[p])[0] # Get the element
prime_divides.pop(p)
factored[i-m-1].pop(p)
# We will keep a cache of already computed results.
# That way we do not have to recalculate them. This trick
# is called "memoizing".
cache = {}
# Recursive search here. options are remaining candidate primes
# and needed are how many more of them we need.
# It returns (count_solutions, data_structure_encoding_solutions)
def search (options, needed):
# This is the key we will cache results under.
key = (tuple(options), needed)
# Return trivial base cases.
if 0 == needed:
return [1, None]
elif len(options) < needed:
return [0, None]
elif key not in cache:
# It is not trivial, and we have never calculated it.
# Actually do work.
answer = []
total = 0
# For each first candidate we can choose.
for i in range(len(options) - needed + 1):
# Remaining candidates are farther along.
next_options = options[i+1:]
# Filter out ones whose p-1 shares a factor with this one.
for p, count in factored[options[i] - m-1].items():
if p == 2 and 1 == count:
# All are even, don't filter on p=2 unless divisible by 4.
continue
# Use list comprehension to do the filtering.
next_options = [option for option in next_options if option not in prime_divides[p]]
# Now the recursive call.
count, result = search(next_options, needed - 1)
# Did we find any?
if 0 < count:
total += count
answer.append([options[i], result])
# cache it
cache[key] = [total, answer]
# and return the cached answer
return cache[key]
# And actually do the recursive search.
return search(available, j)
# The count is the first thing returned.
def count_int_tuples (m, n, j):
return _int_tuples(m, n, j)[0]
# This will yield all of the tuples
def int_tuples (m, n, j):
# the encoded data structure is the second thing returned.
# It has the structure:
# [
# [p1, recursive_encoding of the rest],
# [p2, recursive_encoding of the rest],
# ...
# ]
#
data = _int_tuples(m, n, j)[1]
# chosen is our chosen primes. Since processing always pops first,
# before sticking the next candidate in, it needs to start with
# something, anything. That will be thrown away.
chosen = [None]
# Think of stack as a recursive call stack. Each call has a
# (position, data).
# We start before the beginning of all of our data.
stack = [(-1, data)]
# While there is work to do...
while 0 < len(stack):
# get the next call
i, data = stack.pop()
# clear the previous choice
chosen.pop()
# advance
i = i+1
# If we have another (candidate, more_data) pair to process
if i < len(data):
# Choose our candidate
chosen.append(data[i][0])
# Add walking through current data to the call stack.
stack.append([i, data])
# Did we get a full solution?
if j == len(chosen):
# Return it here, then continue.
# If that is confusing, look up Python iterators.
yield chosen
else:
# Add a recursive call to this data structure.
# Because stacks are Last In First Out (LIFO), this
# will be processed before undoing our choice.
stack.append([-1, data[i][1]])
chosen.append(None)
如果我 运行 print(count_int_tuples(2, 40000, 3)) 在我的笔记本电脑上，只需一秒钟多一点的时间就可以确定有 2,198,808 个长度为 3 的元组满足您的条件。
给定两个数 M_max,M_min 计算素数的所有 r 元组 (M_1,...,M_r) 使得：
- M_min < M_i < M_max
- 每个 M_i-1 至少可以被奇素数的一个平方数整除
- lcm(M_1-1,...,M_r-1)=(M_1-1)*(M_2-1)...(M_r-1)/2^(r-1)
改写问题：k 个数 M_1,M_2,...,M_k 中有多少个 r 元组互质？
这是工作代码。在相当未优化的 Python.
# Stupid simple Euclidean primes filter.
def primes ():
found_primes = []
factor_filter = [False, False]
block_size = 2
candidate = 2
while True:
factor_filter.extend([False for _ in factor_filter])
for p in found_primes:
if 2*block_size <= p*p:
break
i = ((block_size + p - 1) // p) * p
while i < 2 * block_size:
factor_filter[i] = 1
i = i + p
block_size += block_size
while (candidate < block_size):
if 0 == factor_filter[candidate]:
yield candidate
found_primes.append(candidate)
candidate = candidate + 1
def factor_range (m, n):
factored = []
remainder = []
for i in range(n - m + 1):
factored.append({})
remainder.append(i+m)
for p in primes():
if n < p*p:
break
else:
count = 1
q = p
while q <= n:
# Find the first thing divisible by q
i = ((m + q - 1) // q) * q
while i <= n:
factored[i-m][p] = count
remainder[i-m] = remainder[i-m] // p
i += q
q = p*q
count += 1
for i in range(len(remainder)):
if 1 < remainder[i]:
factored[i][remainder[i]] = 1
return factored
def int_tuples (m, n, j):
factored = factor_range(m, n)
available = []
for i in range(len(factored)):
if 0 < i and i + m in factored[i]:
for p, count in factored[i-1].items():
if 2 < p and 1 < count:
available.append(i+m)
break
prime_divides = {}
for i in available:
if m < i:
for p, count in factored[i-1-m].items():
if 1 < count or p != 2:
if p not in prime_divides:
prime_divides[p] = set([i])
else:
prime_divides[p].add(i)
options_stack = [(-1, available, False)]
chosen = []
while 0 < len(options_stack):
#print(("...", chosen, options_stack))
i, options, is_pop = options_stack.pop()
i = i+1
if len(chosen) == j:
yield tuple(chosen)
if is_pop:
chosen.pop()
elif len(chosen) + len(options) - i < j:
# Cannot generate an answer.
if is_pop:
chosen.pop()
elif len(options) <= i:
# i reached the end.
if is_pop:
chosen.pop()
else:
options_stack.append((i, options, is_pop))
chosen.append(options[i])
options = options[i:]
for p, count in factored[options[0]-1].items():
if p == 2 and 1 == count:
continue
if p in prime_divides:
options = [option for option in options if option not in prime_divides[p]]
options_stack.append((-1, options, True))
这个想法是找到这样的元组，使得 (M-1)/2 缩减序列彼此成对相对质数。我所做的相当于递归搜索。但是我保留了显式堆栈，以便我可以通过简单地执行 yield.
我的第一次尝试涉及保留位向量以提高效率。但是较重的数据结构，然后更早地确定“不要进一步搜索”的结论效果更好。
我最终还是决定尝试递归，我很高兴我这么做了。使用缓存是迈向拥有重量级数据结构以避免重复计算的又一步。
# Stupid simple Euclidean primes filter.
def primes ():
found_primes = []
factor_filter = [False, False]
block_size = 2
candidate = 2
while True:
# We get here when we need to double the length of our sieve.
# First add new elements.
factor_filter.extend([False for _ in factor_filter])
# Cross off the ones that divisible by known primes.
for p in found_primes:
if 2*block_size <= p*p:
break
i = ((block_size + p - 1) // p) * p
while i < 2 * block_size:
factor_filter[i] = True
i = i + p
# Record that we did so.
block_size += block_size
# Now that we have sieved, start returning primes.
while (candidate < block_size):
if 0 == factor_filter[candidate]:
yield candidate
found_primes.append(candidate)
candidate = candidate + 1
# Takes a range of numbers, factors all of them at once.
def factor_range (m, n):
# This will hold the factorization.
factored = []
# This will hold the remainder we have not factorized.
remainder = []
# We start with no factors and the number as remainder.
for i in range(n - m + 1):
factored.append({})
remainder.append(i+m)
# We run through the primes up to sqrt(n)
for p in primes():
if n < p*p:
break
else:
# For each power of p, we add it to the factorization and divide it form the remainder.
count = 1
q = p
while q <= n:
# Find the first thing divisible by q
i = ((m + q - 1) // q) * q
while i <= n:
factored[i-m][p] = count
remainder[i-m] = remainder[i-m] // p
i += q
q = p*q
count += 1
# Any left-over remainders must be primes.
for i in range(len(remainder)):
if 1 < remainder[i]:
factored[i][remainder[i]] = 1
return factored
# This function does all of the work.
def _int_tuples (m, n, j):
# BUG FIX HERE:
# If our minimum is odd, extend the boundary down one.
# This does not introduce new solutions, but does mean we
# factorize all of the numbers that we need.
if 1 == m // 2:
m = m-1
# Factorize all numbers.
factored = factor_range(m, n)
# available will be all of our candidate primes.
available = []
for i in range(len(factored)):
if 0 < i and i + m in factored[i]:
for p, count in factored[i-1].items():
if 2 < p and 1 < count:
available.append(i+m)
break
# Construct a data structure to answer, "which other candidate primes
# have p-1 share a non-trivial factor with this one"
prime_divides = {}
for i in available:
if m < i:
for p, count in factored[i-1-m].items():
if 1 < count or p != 2:
if p not in prime_divides:
prime_divides[p] = set([i])
else:
prime_divides[p].add(i)
# PERFORMANCE HACK: We will recurse over lists of remaining candidates.
# Putting candidates that eliminate lots of other candidates first
# makes those lists smaller and that recursion more efficient. We
# also improve the odds that we will be looking at a list we have
# already seen and already have an answer for. And, finally, the
# candidates at the end of our list will largely not filter each
# other out. So we are spending time on things that are more likely to
# give us more solutions.
count_related = {}
for i in available:
total = 0
for p, count in factored[i-1-m].items():
if 1 < count or p != 2:
total += len(prime_divides[p]) - 1
count_related[i] = total
available = sorted(available, key = lambda i: -count_related[i])
# To spend less time filtering, forget about the primes that only
# divide one thing.
for p in list(prime_divides.keys()):
if 1 == len(prime_divides[p]):
i = list(prime_divides[p])[0] # Get the element
prime_divides.pop(p)
factored[i-m-1].pop(p)
# We will keep a cache of already computed results.
# That way we do not have to recalculate them. This trick
# is called "memoizing".
cache = {}
# Recursive search here. options are remaining candidate primes
# and needed are how many more of them we need.
# It returns (count_solutions, data_structure_encoding_solutions)
def search (options, needed):
# This is the key we will cache results under.
key = (tuple(options), needed)
# Return trivial base cases.
if 0 == needed:
return [1, None]
elif len(options) < needed:
return [0, None]
elif key not in cache:
# It is not trivial, and we have never calculated it.
# Actually do work.
answer = []
total = 0
# For each first candidate we can choose.
for i in range(len(options) - needed + 1):
# Remaining candidates are farther along.
next_options = options[i+1:]
# Filter out ones whose p-1 shares a factor with this one.
for p, count in factored[options[i] - m-1].items():
if p == 2 and 1 == count:
# All are even, don't filter on p=2 unless divisible by 4.
continue
# Use list comprehension to do the filtering.
next_options = [option for option in next_options if option not in prime_divides[p]]
# Now the recursive call.
count, result = search(next_options, needed - 1)
# Did we find any?
if 0 < count:
total += count
answer.append([options[i], result])
# cache it
cache[key] = [total, answer]
# and return the cached answer
return cache[key]
# And actually do the recursive search.
return search(available, j)
# The count is the first thing returned.
def count_int_tuples (m, n, j):
return _int_tuples(m, n, j)[0]
# This will yield all of the tuples
def int_tuples (m, n, j):
# the encoded data structure is the second thing returned.
# It has the structure:
# [
# [p1, recursive_encoding of the rest],
# [p2, recursive_encoding of the rest],
# ...
# ]
#
data = _int_tuples(m, n, j)[1]
# chosen is our chosen primes. Since processing always pops first,
# before sticking the next candidate in, it needs to start with
# something, anything. That will be thrown away.
chosen = [None]
# Think of stack as a recursive call stack. Each call has a
# (position, data).
# We start before the beginning of all of our data.
stack = [(-1, data)]
# While there is work to do...
while 0 < len(stack):
# get the next call
i, data = stack.pop()
# clear the previous choice
chosen.pop()
# advance
i = i+1
# If we have another (candidate, more_data) pair to process
if i < len(data):
# Choose our candidate
chosen.append(data[i][0])
# Add walking through current data to the call stack.
stack.append([i, data])
# Did we get a full solution?
if j == len(chosen):
# Return it here, then continue.
# If that is confusing, look up Python iterators.
yield chosen
else:
# Add a recursive call to this data structure.
# Because stacks are Last In First Out (LIFO), this
# will be processed before undoing our choice.
stack.append([-1, data[i][1]])
chosen.append(None)
如果我 运行 print(count_int_tuples(2, 40000, 3)) 在我的笔记本电脑上，只需一秒钟多一点的时间就可以确定有 2,198,808 个长度为 3 的元组满足您的条件。