pythonmathrecurrenceoeis

Is there a more Pythonic way of coding this recurrance relation re: OEIS A077947


I am working on a thesis regarding Jacobsthal sequences (A001045) and how they can be considered as being composed of some number of distinct sub-sequences. I have made a comment on A077947 indicating this and have included a python program. Unfortunately the program as written leaves a lot to be desired and so of course I wanted to turn to Stack to see if anyone here knows how to improve the code!

Here is the code:

a=1
b=1
c=2
d=5
e=9
f=18
for x in range(0, 100):
  print(a, b, c, d, e, f)
  a = a+(36*64**x)
  b = b+(72*64**x)
  c = c+(144*64**x)
  d = d+(288*64**x)
  e = e+(576*64**x)
  f = f+(1152*64**x) 

I explain the reasoning behind this as follows:

The sequence A077947 is generated by 6 digital root preserving sequences stitched together; per the Python code these sequences initiate at the seed values a-f. The number of iterations required to calculate a given A077947 a(n) is ~n/6. The code when executed returns all the values for A077947 up to range(x), or ~x*6 terms of A077947. I find the repeated digital roots interesting as I look for periodic digital root preservation within sequences as a method to identify patterns within data. For example, digital root preserving sequences enable time series analysis of large datasets when estimating true-or-false status for alarms in large IT ecosystems that undergo maintenance (mod7 environments); such analysis is also related to predicting consumer demand / patterns of behavior. Appropriating those methods of analysis, carving A077947 into 6 digital root preserving sequences was meant to reduce complexity; the Python code reproduces A077947 across 6 "channels" with seed values a-f. This long paragraph boils down to statement, "The digital roots of the terms of the sequence repeat in the pattern (1, 1, 2, 5, 9, 9)." The bigger statement is that all sequences whose digital roots repeat with a pattern can be partitioned/separated into an equal number of distinct sequences and those sequences can be calculated independently. There was a bounty related to this sequence.

This code is ugly but I cannot seem to get the correct answer without coding it this way;

I have not figured out how to write this as a function due to the fact I cannot seem to get the recurrence values to store properly in a function.

So of course if this yields good results we hope to link the discussion to the OEIS references.

Here is a link to the sequence: https://oeis.org/A077947


Solution

  • Here's an alternative way to do it without a second for loop:

    sequences   = [ 1,  1,  2,   5,   9,   18   ]
    multipliers = [ 36, 72, 144, 288, 576, 1152 ]
    
    for x in range(100):
        print(*sequences)
        sequences = [ s + m*64**x for s,m in zip(sequences,multipliers) ]
    

    [EDIT] Looking at the values I noticed that this particular sequence could also be obtained with:

    N[i+1] = 2 * N[i] + (-1,0,1 in rotation)

    or

    N[i+1] = 2 * N[i] + i mod 3 - 1 (assuming a zero based index)

    i  N[i] [-1,0,1]             N[i+1]
    0  1    -1   --> 2*1 - 1 --> 1
    1  1     0   --> 2*1 + 0 --> 2
    2  2     1   --> 2*2 + 1 --> 5
    3  5    -1   --> 2*5 - 1 --> 9
    4  9     0   --> 2*9 + 0 --> 18
    ... 
    

    So a simpler loop to produce the sequence could be:

    n = 1
    for i in range(100):
        print(n)
        n = 2*n + i % 3 - 1
    

    Using the reduce function from functools can make this even more concise:

    from functools import reduce
    
    sequence = reduce(lambda s,i: s + [s[-1]*2 + i%3 - 1],range(20),[1])
    print(sequence)
    
    >>> [1, 1, 2, 5, 9, 18, 37, 73, 146, 293, 585, 1170, 2341, 4681, 9362, 18725, 37449, 74898, 149797, 299593, 599186]
    

    Using your multi-channel approach and my suggested formula this would give:

    sequences   = [ 1,  1,  2,   5,   9,   18   ]
    multipliers = [ 36, 72, 144, 288, 576, 1152 ]
    
    allSequences = reduce(lambda ss,x: ss + [[ s + m*64**x for s,m in zip(ss[-1],multipliers) ]],range(100),[sequences])
    for seq in allSequences: print(*seq) # print 6 by 6
    

    [EDIT2] If all your sequences are going to have a similar pattern (i.e. starting channels, multipliers and calculation formula), you could generalize the printing of such sequences in a function thus only needing one line per sequence:

    def printSeq(calcNext,sequence,multipliers,count):
        for x in range(count):
            print(*sequence)
            sequence = [ calcNext(x,s,m) for s,m in zip(sequence,multipliers) ]
    
    printSeq(lambda x,s,m:s*2+m*64**x,[1,1,2,5,9,18],multipliers=[36,72,144,288,576,1152],count=100)
    

    [EDIT3] Improving on the printSeq function.

    I believe you will not always need an array of multipliers to compute the next value in each channel. An improvement on the function would be to provide a channel index to the lambda function instead of a multiplier. This will allow you to use an an array of multiplier if you need to but will also let you use a more general calculation.

    def printSeq(name,count,calcNext,sequence):
        p = len(sequence)
        for x in range(count):
            print(name, x,":","\t".join(str(s) for s in sequence))
            sequence = [ calcNext(x,s,c,p) for c,s in enumerate(sequence) ]
    

    The lambda function is given 4 parameters and is expected to return the next sequence value for the specified channel:

    s : current sequence value for the channel
    x : iteration number
    c : channel index (zero based)
    p : number of channels
    

    So, using an array inside the formula would express it like this:

    printSeq("A077947",100,lambda x,s,c,p: s + [36,72,144,288,576,1152][c] * 64**x, [1,1,2,5,9,18])
    

    But you could also use a more general formula that is based on the channel index (and number of channels):

    printSeq("A077947",100,lambda x,s,c,p: s + 9 * 2**(p*x+c+2), [1,1,2,5,9,18])
    

    or ( 6 channels based on 2*S + i%3 - 1 ):

    printSeq("A077947",100,lambda x,s,c,p: 64*s + 9*(c%3*2 - (c+2)%3 - 1) ,[1,1,2,5,9,18])
    printSeq("A077947",100,lambda x,s,c,p: 64*s + 9*[-3,1,2][c%3],[1,1,2,5,9,18])
    

    My reasoning here is that if you have a function that can compute the next value based on the current index and value in the sequence, you should be able to define a striding function that will compute the value that is N indexes farther.

    Given F(i,S[i]) --> i+1,S[i+1]

    F2(i,S[i]) --> i+2,S[i+2] = F(F(i,S[i]))
    F3(i,S[i]) --> i+3,S[i+3] = F(F(F(i,S[i])))
    ...
    F6(i,S[i]) --> i+6,S[i+6] = F(F(F(F(F(F(i,S[i]))))))
    ...
    Fn(i,S[i]) --> i+n,S[i+n] = ...
    

    This will always work and should not require an array of multipliers. Most of the time it should be possible to simplify Fn using mere algebra.

    for example A001045 : F(i,S) = i+1, 2*S + (-1)**i

    printSeq("A001045",20,lambda x,s,c,p: 64*s + 21*(-1)**(x*p+c),[0,1,1,3,5,11])
    

    Note that from the 3rd value onward, the next value in that sequence can be computed without knowing the index:

    A001045: F(S) = 2*S + 1 - 2*0**((S+1)%4)