我正致力于从Project Euler解决以下问题，简而言之，该问题涉及在'n'骰子上迭代并更新其值。
一长排骰子-项目Euler问题#641
考虑一排n个骰子，所有骰子都显示1。
每隔一个模具(2，4，6，…(转动一次，使显示的数字增加1。然后每隔三分之一转一圈。第六个骰子现在显示3。然后每隔四个骰子旋转一次，依此类推，直到每n个骰子(仅最后一个骰子(旋转一次。如果要转动的模具显示为6，则将其更改为显示1。
设f(n(为进程结束时显示1的骰子数。给定f(100(=2和f(10^8(=69。
求f(10^36(。
我已经使用numpy在Python中编写了以下代码，但无法准确地弄清楚我对函数输出做了什么错误，以匹配上面的输出。现在f(100)返回1(应该是2(；即使CCD_ 2返回1。
import numpy as np
def f(n):
# establish dice and the value sets for the dice
dice = np.arange(1, n + 1)
dice_values = np.ones(len(dice))
turns = range(2, len(dice) + 1)
print("{a} dice, {b} values, {c} runs to process".format(a=len(dice), b=len(dice_values), c=len(turns)))
# iterate and update the values of each die
# in our array of dice
for turn in turns:
# if die to be processed is 6, update to 1
dice_values[(dice_values == 6) & (dice % turn == 0)] = 1
# update dice_values to if the die's index has no remainder
# from the turn we're processing.
dice_values += dice % turn == 0
# output status
print('Processed every {0} dice'.format(turn))
print('{0}nn'.format(dice_values))
return "f({n}) = {x}".format(n=n, x=len(np.where(dice_values == 1)))
更新11/12/18
@Prune的指导非常有帮助。我的方法如下：
- 找到从1到n的所有正方形
当除以6时，找出所有具有余数为1的因子的平方。
import numpy as np # brute force to find number of factors for each n def factors(n): result = [] i = 1 # This will loop from 1 to int(sqrt(n)) while i * i <= n: # Check if i divides x without leaving a remainder if n % i == 0: result.append(i) if n / i != i: result.append(n / i) i += 1 # Return the list of factors of x return len(result) vect_factors = np.vectorize(factors) # avoid brute forcing all numbers def f(n): # create an array of 1 to n + 1 # find all perfect squares in that range dice = np.arange(1, n + 1)[(np.mod(np.sqrt(np.arange(1, n + 1)), 1) == 0)] # find all squares which have n-factors, which # when divided by 6 have a remainder of 1. dice = dice[np.mod(vect_factors(dice), 6) == 1] return len(dice)
值得注意的是，在我的机器上，我无法运行大于10^10的程序。虽然解决这个问题是理想的，但我觉得我在这个过程中学到的(并决定如何应用(对我来说已经足够了
更新11/13/2018
我会继续花一点时间来优化它，以便更快地处理它。这是更新后的代码库。这将在1分17秒内评估f(10**10(。
import time
from datetime import timedelta
import numpy as np
import math
from itertools import chain, cycle, accumulate
def find_squares(n):
return np.array([n ** 2 for n in np.arange(1, highest = math.sqrt(n) + 1)])
# brute force to find number of factors for each n
def find_factors(n):
def prime_powers(n):
# c goes through 2, 3, 5, then the infinite (6n+1, 6n+5) series
for c in accumulate(chain([2, 1, 2], cycle([2, 4]))):
if c * c > n: break
if n % c: continue
d, p = (), c
while not n % c:
n, p, d = n // c, p * c, d + (p,)
yield (d)
if n > 1: yield ((n,))
r = [1]
for e in prime_powers(n):
r += [a * b for a in r for b in e]
return len(r)
vect_factors = np.vectorize(find_factors)
# avoid brute forcing all numbers
def f(n):
# create an array of 1 to n + 1
# find all perfect squares in that range
start = time.time()
dice = find_squares(n)
# find all squares which have n-factors, which
# when divided by 6 have a remainder of 1.
dice = dice[np.mod(vect_factors(dice), 6) == 1]
diff = (timedelta(seconds=int(time.time() - start))).__str__()
print("{n} has {remain} dice with a value of 1. Computed in {diff}.".format(n=n, remain=len(dice), diff=diff))
我提出了一个x/y问题。修复6=>1翻转将更正您的代码，但它不会在合理的时间内解决出现的问题。要找到f(10^36)，你要处理10^36个骰子，每个10^36次，即使它只是过滤器中的可分割性检查。总共有10^72张支票。我不知道你有什么硬件，但即使是我的多核怪物也不会很快循环10^72次，让人感到舒适。
相反，您需要找出潜在的问题，并尝试为符合描述的整数生成计数。
骰子只是一种在mod 6中计数的设备。我们在计算一个数的除数，包括1和这个数本身。这就是著名的除数函数。
手头的问题并没有要求我们找到所有数字的σ0(n(；它希望我们计数有多少个整数具有σ0(n(=1(mod 6(。这些数字有1、7、13、19。。。除数。
首先，请注意，这是一个奇数。只有除数为奇数的整数是完全平方。看看除数函数；我们如何判断一个数的平方是否具有所需的因子数量1(mod 6(？
这能让你动起来吗？
周末更新
我的代码通过10^18候选人仍然太慢，无法在本日历年完成。它在大约10^7之前表现良好，然后陷入O(N log N(检查步骤。
然而，我在跟踪输出中注意到了更多的限制。主要的一个是描述素数幂的组合会产生一个解。如果我们减少每个功率mod 3，我们有以下内容：
0值不影响结果的有效性1值使该数字无效2值必须成对
此外，这些条件对于声明给定的数字为解既是必要的，也是充分的。因此，可以生成所需的解决方案，而不必费力地遍历所有整数的平方<=10^18.
除其他外，我们只需要10^9以下的素数：一个解的平方根至少需要任何素数因子的2。
我希望这已经足够暗示了。。。您需要构造一个算法来生成某些具有给定乘积上限的受限组合。
正如Thierry在评论中提到的，当你在6点掷骰子时，你会循环回到2。我建议将dice_values[(dice_values == 6) & (dice % turn == 0)] = 1更改为等于0。
你还有一个return "f({n}) = {x}".format(n=n, x=len(np.where(dice_values == 1)))的问题，我会用x=np.count_nonzero(dice_values == 1)代替x=len(np.where(dice_values == 1))来解决
完成这两项更改后，我得到了f(100)=2的输出