[ prog / sol / mona ]

prog


LISP Puzzles

18 2020-04-06 10:57

Here is a moderately efficient implementation using O(n^(0.5+epsilon)) space and an iterative computational process. It is based on two simple observations. The first one is that the S sequence can be partitioned into continuous groups with jumps only at the group boundaries. The groups contain no more information than their starting value and length, and can therefore be represented by such a pair, effectively performing run-length encoding on the S sequence. The second one is that, since the gaps in the S sequence are given by the R values, and since R is just a running sum over S, the S-groups can be computed directly from the S sequence, which can thus generate itself without any help from R. Combining these observations yields:

def work_seq (n):
   # heuristic size
   groups = [None] * int (2 * n ** 0.5)
   groups [0] = (2, 1)
   groups [1] = (4, 3)

   group_count  = 2
   index_groups = 1
   used_group   = 1
   lastR        = 3

   start, count = groups [index_groups]
   covered = 2 + count
   S       = start
   length  = count

   while covered < n:
      a = S
      b = S + length
      # R += [a...b)
      lastR += (b * (b - 1) - a * (a - 1)) // 2

      S += length + 1

      if used_group < count:
         length = start + used_group - 1
         used_group += 1
      else:
         index_groups += 1
         start, count = groups [index_groups]
         length = start - 1
         used_group = 1

      groups [group_count] = (S, length)
      group_count += 1
      covered     += length

   # print (groups)
   print ("groups:", group_count, "index:", index_groups)

   a = S
   b = S + length - (covered - n)
   # R += [a...b)
   R = lastR + (b * (b - 1) - a * (a - 1)) // 2

   return R

This is an inner helper function that is only called by its wrapper for n>=3. The function can be easily modified to yield all the R values to its caller. One could also snapshot the state of the computation and resume it later.

Using the convention that R(1)=1 and R(2)=3, the first n of the form 10^k for which R(n) no longer fits into an u64 is ten billion (10^10), so I propose that we use n=10^10 to compare implementations. The value of R(10^10) is
50000940106333921296
While computing it the function stored 141082 group entries, a touch over O(sqrt(n)), out of which 512 were used to compute group lengths, a touch over O(n^0.25).

$ python3 -m flake.hofs seq 10000000000
groups: 141082 index: 512
10000000000: 50000940106333921296 0x2B5E7061C419DA010

The computation takes just over 70 milliseconds. While this can be further improved, keep in mind that this does not use any special libraries like numpy or gmpy2, and that it is run by the default CPython interpreter. A port to C with GMP should easily bring down the time.

$ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_seq (10000000000)'
10 loops, best of 3: 71.3 msec per loop

I'll post some data about the epsilon values in a bit. I hope other Anons will post their timings for R(10^10).

>>16
Thanks for posting timings. Here is confirmation of your R(2*10^9) value:

$ python3 -m flake.hofs seq 2000000000
groups: 63021 index: 340
2000000000: 2000083969706219361 0x1BC1B9C6074AD761

and here is the runtime in my implementation, just over 30 milliseconds:

$ python3 -m timeit -s 'import flake.hofs as mod' 'mod.work_seq (2000000000)'
10 loops, best of 3: 30.5 msec per loop
157


VIP:

do not edit these