Xtykutl

I'm working on solve the below problem from Project Euler, which in short deals with iterating over 'n' dice and updating their values.

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).

I've written the below code in Python using numpy, but can't exactly figure out what I'm doing wrong to my function output to match the output above. Right now f(100) returns 1 (should be 2); even f(1000) returns 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))) 

UPDATE 11/12/18

@Prune's guidance has been extremely helpful. My methodology is now as follows:

• Find all the squares from 1 to n.


• 
  

  Find all squares with a number of factors which have a remainder of 1, when dividing by 6.

  
  
  

  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)
  
  

  
  

Worth noting - on my machine, I'm unable to run larger than 10^10. While solving this would be ideal, I feel that what I've learned (and determined how to apply) in the process is enough for me.

UPDATE 11/13/2018

I'm continuing to spend a small bit of time trying to optimize this to get it processing more quickly. Here's the updated code base. This evaluates f(10**10) in 1 min and 17 seconds.

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))

I'm raising an x/y issue. Fixing your 6 => 1 flip will correct your code, but it will not solve the presented problem in reasonable time. To find f(10^36), you're processing 10^36 dice 10^36 times each, even if it's merely a divisibility check in the filter. That's a total of 10^72 checks. I don't know what hardware you have, but even my multi-core monster doesn't loop 10^72 times soon enough for comfort.

Instead, you need to figure out the underlying problem and try to generate a count for integers that fit the description.

The dice are merely a device to count something in mod 6. We're counting divisors of a number, including 1 and the number itself. This the (in)famous divisor function.

The problem at hand doesn't ask us to find σ0(n) for all numbers; it wants us to count how many integers have σ0(n) = 1 (mod 6). These are numbers with 1, 7, 13, 19, ... divisors.

First of all, note that this is an odd number. The only integers with an odd number of divisors are perfect squares. Look at the divisor function; how can we tell whether the square of a number will have the desired quantity of factors, 1 (mod 6)?

Does that get you moving?

WEEKEND UPDATE

My code to step through 10^18 candidates is still too slow to finish in this calendar year. It did well up to about 10^7 and then bogged down in the ON(log N) checking steps.

However, there are many more restrictions I've noted in my tracing output.
The main one is in characterizing what combinations of prime powers result in a solution. If we reduce each power mod 3, we have the following:

• 
  0 values do not affect validity of the result.


• 
  1 values make the number invalid.


• 
  2 values must be paired.

Also, these conditions are both necessary and sufficient to declare a given number as a solution. Therefore, it's possible to generate the desired solutions without bothering to step through the squares of all integers <= 10^18.

Among other things, we will need only primes up to 10^9: a solution's square root will need at least 2 of any prime factor.

I hope that's enough hints for now ... you'll need to construct an algorithm to generate certain restricted composite combinations with a given upper limit for the product.

As mentioned by Thierry in the comments, you are looping back to 2 when you flip dice at a 6. I'd suggest just changing dice_values[(dice_values == 6) & (dice % turn == 0)] = 1 to equal 0.

You also have an issue with return "f({n}) = {x}".format(n=n, x=len(np.where(dice_values == 1))) that I'd fix by replacing x=len(np.where(dice_values == 1)) with x=np.count_nonzero(dice_values == 1)

Doing both these changes gave me an output of f(100)=2