What number can come next only depends on the previous two numbers. So here's an iterative solution that only keeps O(1) numbers around, it's about twice as fast as Tomer's for N=2000 (even without any optimizations, I even call gcd all the time):
from collections import Counter
from math import gcd
def count_seq(N):
ctr = Counter([(7, 7)])
for _ in range(N):
temp = Counter()
for (a, b), count in ctr.items():
for c in range(1, 7):
if c != a and c != b and gcd(b, c) == 1:
temp[b, c] += count
ctr = temp
return sum(ctr.values())
My ctr is a hash map whose keys are pairs representing the last two numbers, and the values tell how many valid sequences end in those two numbers. Initialized with (7, 7) because that allows all extensions.
For fun and max performance, hardcoding all states and transitions, this is about 14 times faster than my above solution for N=2000 and solves N=100,000 in about 10 seconds (the result is a number with 45,077 digits):
def count_seq(N):
if N < 2:
return 6 if N else 1
x12 = x13 = x14 = x15 = x16 = x21 = x23 = x25 = x31 = x32 = x34 = x35 = x41 = x43 = x45 = x51 = x52 = x53 = x54 = x56 = x61 = x65 = 1
for _ in range(N - 2):
x12, x13, x14, x15, x16, x21, x23, x25, x31, x32, x34, x35, x41, x43, x45, x51, x52, x53, x54, x56, x61, x65, = x31+x41+x51+x61, x21+x41+x51+x61, x21+x31+x51+x61, x21+x31+x41+x61, x21+x31+x41+x51, x32+x52, x12+x52, x12+x32, x23+x43+x53, x13+x43+x53, x13+x23+x53, x13+x23+x43, x34+x54, x14+x54, x14+x34, x25+x35+x45+x65, x15+x35+x45+x65, x15+x25+x45+x65, x15+x25+x35+x65, x15+x25+x35+x45, x56, x16
return x12 + x13 + x14 + x15 + x16 + x21 + x23 + x25 + x31 + x32 + x34 + x35 + x41 + x43 + x45 + x51 + x52 + x53 + x54 + x56 + x61 + x65
Here, for example x23 is the number of valid sequences ending in the numbers 2 and 3. There's room for further microoptimizations (using additional y variables to alternate between x and y instead of building+unpacking tuples, or exploiting common subsums like x21+x41, which is used three times), but I'll leave it at this.
Oh, actually... as one might've seen with Fibonacci numbers, we can represent the transitions as a 22×22 matrix and then quickly compute the Nth (or N-2th) power of that matrix with exponentiation by squaring. That should be even much faster.
Meh... implemented the matrix power method now, and sadly it's slower. I guess it really only helps if you only need some kind of truncated values, like only the first or last 100 digits and the number of digits, or the number of sequences modulo some number. Otherwise, while the matrix power method does need fewer operations, they include multiplications of large numbers instead of just additions, which is slower. Anyway, here's my implementation (Try it online!):
from math import gcd
import numpy as np
def count_seq(N):
if N < 2:
return 6 if N else 1
states = [(a, b) for a in range(1, 7) for b in range(1, 7) if a != b and gcd(a, b) == 1]
A = [[1 if b2 == b and a != c else 0
for b2, c in states]
for a, b in states]
return (np.matrix(A, dtype=object) ** (N-2)).sum()
As demonstration, here's a modification that computes the last three digits. For N=100,000 it takes about 0.26 seconds and for N=1,000,000,000,000,000 it takes about one second.
def count_seq(N):
if N < 2:
return 6 if N else 1
modulo = 1000
class IntModulo:
def __init__(self, value):
self.value = value
def __add__(self, other):
return IntModulo((self.value + other.value) % modulo)
def __mul__(self, other):
return IntModulo((self.value * other.value) % modulo)
def __int__(self):
return self.value
states = [(a, b) for a in range(1, 7) for b in range(1, 7) if a != b and gcd(a, b) == 1]
A = [[IntModulo(1 if b2 == b and a != c else 0)
for b2, c in states]
for a, b in states]
return int((np.matrix(A, dtype=object) ** (N-2)).sum())
