python - 欧拉计划 #641 Python 3.6 - Numpy

Question

我正在努力解决 Project Euler 中的以下问题，简而言之，它涉及迭代“n”个骰子并更新它们的值。

A Long Row of Dice - project Euler problem #641

Consider a row of n dice all showing 1.

First turn every second die,(2,4,6,…), so that the number showing is increased by 1. Then turn every third die. The sixth die will now show a 3. Then turn every fourth die and so on until every nth die (only the last die) is turned. If the die to be turned is showing a 6 then it is changed to show a 1.

Let f(n) be the number of dice that are showing a 1 when the process finishes. You are given f(100)=2 and f(10^8)=69.

Find f(10^36).

我已经使用 numpy 在 Python 中编写了以下代码，但无法确切地弄清楚我的函数输出做错了什么以匹配上面的输出。现在 f(100) 返回 1(应该是 2)；甚至 f(1000) 也返回 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}\n\n'.format(dice_values))
    return "f({n}) = {x}".format(n=n, x=len(np.where(dice_values == 1))) 

---

18 年 11 月 12 日更新

@Prune 的指导非常有帮助。我现在的方法如下:

• 找出从 1 到 n 的所有方 block 。

• 找出所有除以 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 的数据。虽然解决这个问题是理想的，但我觉得我在这个过程中学到的(并确定如何应用)对我来说已经足够了。

---

2018 年 11 月 13 日更新

我将继续花一点时间尝试优化它以使其处理得更快。这是更新的代码库。这将在 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))

Answer 1

我提出了一个x/y 问题。修复您的 6 => 1 翻转将更正您的 代码，但它不会在合理的时间内解决所提出的问题。要找到 f(10^36)，您需要处理 10^36 个骰子，每个骰子 10^36 次，即使它只是过滤器中的可除性检查。总共有 10^72 张支票。我不知道你有什么硬件，但即使是我的多核怪兽也不会很快循环 10^72 次以保证舒适。

相反，您需要找出根本问题并尝试生成符合描述的整数计数。

在 mod 6 中，骰子只是一种用于计数的设备。我们正在计算一个数字的除数，包括 1 和数字本身。这是著名的 divisor function .

手头的问题并不要求我们找到所有数字的 σ0(n)；它希望我们计算有多少整数具有 σ0(n) = 1 (mod 6)。这些数字的除数为 1、7、13、19 ...。

首先，请注意这是一个奇数。唯一具有奇数个除数的整数是完全平方数。查看除数函数；我们如何判断一个数的平方是否具有所需数量的因子，1 (mod 6)？

这会让你动起来吗？

---

周末更新

我的代码要通过 10^18 个候选人，在这个日历年内完成仍然太慢。它在 10^7 左右表现良好，然后陷入 O(N log N) 检查步骤。

但是，我在跟踪输出中注意到了更多的限制。 最主要的是表征哪些素数幂的组合会导致解决方案。如果我们减少每个功率模 3，我们有以下结果:

• 0 值不影响结果的有效性。
• 1 值使数字无效。
• 2 值必须成对。

此外，这些条件对于将给定数字声明为解决方案既是必要的又是充分的。因此，可以生成所需的解决方案，而不必费心遍历所有 <= 10^18 的整数的平方。

除此之外，我们只需要 10^9 以内的素数:解的平方根至少需要任何素数的 2。

我希望现在这些提示已经足够了……您需要构建一个算法来生成具有给定乘积上限的特定受限复合组合。