How to find the element of the Digit Sum sequence efficiently? Just out of interest I tried to solve a problem from "Recent" category of Project Euler (  Digit Sum sequence ). But I am unable to think of a way to solve the problem efficiently. The problem is as follows ( in the original question sequence has two ones in beginning , but it does not change the sequence ) :


  The Digit Sum sequence is 1,2,4,8,16,23,28,38,49.... where the $n^{th}$ term of sequence is sum of digits preceding it in the sequence. Find the $10^{15}th$ term of the sequence.


The naive solution can't be implemented because it takes a lot of time. I tried to reduce the problem to a case of matrix exponentiation ( that would takes $O(log ( 10^{15}))$ amount of time ) but could not come up with such a recurrence fitting the linear criteria as recurrence for this sequence is quite peculiar. It can be seen that the sequence is governed by the recurrence :

$$ a_n = a_{n-1} + d( a_{n-1} ) ..... (1 )$$

where $a_n$ is $n^{th}$ term of the sequence and $d$ is a function which when given a natural number as input returns sum of digits of the number ( eg. $\;d(786)=21$ ). My second approach was to try to find some pattern in the sequence. It can be seen that the first few terms of the sequence can be written as

   a_1 = 1  
   a_2 = 1 + d( 1 )
   a_3 = 1 + d( 1 ) + d( 1 + d( 1 ) )
   a_4 = 1 + d( 1 ) + d( 1 + d( 1 ) ) + d( 1 + d( 1 ) + d( 1 + d( 1 ) ) )
   a_5 = 1 + d( 1 ) + d( 1 + d( 1 ) ) + d( 1 + d( 1 ) + d( 1 + d( 1 ) ) ) + d( 1 +  d(  
   1 ) + d( 1 + d( 1 ) ) + d( 1 + d( 1 ) + d( 1 + d( 1 ) ) ) )


From the pattern above it becomes that $n^{th}$ term of the sequence can be  generated by the following method :


Write $2^{n-1}$ $1$'s with addition symbol between them.
Leaving  the first $1$, then apply the function $d$  on the next $2^{0}$ terms then on next $2^{1}$ terms, then on next $2^{2}$ terms and so on.
Then  apply the above method recursively on arguments of each $d$ function applied.


for example if n=3 we perform the following manipulations:

    1 + 1 + 1 + 1 + 1 + 1 + 1 + 1
    1 + d( 1 ) + d( 1 + 1 ) + d( 1 + 1 + 1 + 1 )
    1 + d( 1 ) + d( 1 + d(1) ) + d( 1 + d( 1 ) + d( 1 +d( 1 ) ) )


By dynamic programming I can generate the $n^{th}$ term using the above method in time $O(log ( 2^{10^{15}} ) )$, which again is no better than the naive solution.

EDIT 1
Another thing that can be observed is that $d(a_n) = d(2^{n-1})$. For example $d(a_6)=d(23)=d(32)=5$. But I am unable to make use of this point. I again tried to find a linear recurrence relation ( for matrix exponentiation ) , but I am unable to find it.

EDIT 2 

Following is the graph when the sequence is plotted for smaller range ( first $10^6$ terms of the sequence  are plotted ). 


PS: I know it is not advisable to ask solutions from Project Euler. But I just want a new direction or a hint, as I have been moving in circles for past few days. If that is also unacceptable I can remove the question if suggested.
Since you asked for " a new direction or a hint " and I don't know the answer, I'll leave this here, I hope it's helpful.
some ideas:

It makes sense there would be a pattern mod 9, since

$\forall k &gt; 1,k \in \mathbb{Z}\quad10^k \equiv 1 \mod 9$

Which you can prove by induction.

This means that all numbers are congruent to the sum of their digits mod 9.

Further,
$a_n = d(a_n) \mod 9 \implies $
$$a_{n} = a_{n-1} + d(a_{n-1}) = 2d(a_{n-1}) \mod 9$$

If we keep expanding this recurrence we get

$a_{n} = 2^n \mod 9$

Which explains the pattern mod 9.

It also gets us $a_n = 9k + 2^n$. Every iteration we get a gap that is divisible by 9. How wide are those gaps?

Here's some less than general code:

import matplotlib.pyplot as plt
import numpy as np

#sum digits of n
def sum_digits(n):
    s = 0
    while n:
        s += n % 10
        n //= 10
    return s

#get the sequence to n digits
def calculate(n):
    retval = [1]
    for i in range(n):
        retval.append(retval[-1] + sum_digits(retval[-1]))
    return retval;

#empirically confirm that a_n = 2^n mod 9
def confirmPow2(a):
    count = 0
    for i in a[:10000]:
        if((i%9) != (2**count % 9)):
            print "false"
        count = count + 1

#find gaps divisible by 9 in a subset of a
def find9Gaps(a):
    count = 0
    S = []
    for i in a[:10000]:
         S.append(((2**count ) - i)/9)
         count = count + 1
    return S

#repeatedly sum the digits until they're less than 9...
#gives some interesting patterns
def repeatedDigitSum():
    for i in range(1000, 1100):
         print "=========for ",i
         while i &gt; 9:
                 i = sum_digits(i)
                 print i 


a = calculate(10**6)
b = find9Gaps(a)
plt.plot(range(len(b[:100])), b[:100])
plt.show()


The plot (for the first 100) looks exponential, but I don't think it's perfect.



Here's the output of

&gt;&gt;&gt; plt.plot(range(len(b[5:60])), np.log2(np.array(b[5:60])))
&gt;&gt;&gt; plt.show()




The last thing I have is that it seems if you sum the digits of a number, and then sum the digits of the resulting number, and repeat this, you eventually get that number mod 9.

Makes sense given the above fact about powers of 10 mod 9.

$$n \equiv d(n) \equiv d(d(n)) \equiv \cdots \mod 9$$

It gives an interesting sequence of numbers though.

Edit: Apparently this is called a "digital root".