[seqfan] A071053 (number of odd terms in expansion of (1+x+x^2)^n)

[Joerg Arndt]

Cf.
Shalosh B. Ekhad, N. J. A. Sloane, Doron Zeilberger,
A Meta-Algorithm for Creating Fast Algorithms for Counting
  ON Cells in Odd-Rule Cellular Automata,
http://arxiv.org/abs/1503.01796

The following Pari/GP program computes a(Pi*10^1000) in < 1 second.

-------------------------------------------------------
b(n) = { (2^n - (-1)^n) / 3; }  \\ A001045
a(n)=
{
    if ( n==0, return(1) );

    \\ Use  a( 2^k * t ) = a(t)
    n \= 2^valuation(n,2);
    if ( n==1, return(3) );  \\ Use a(2^k) == 3
    \\ now n is odd

    my ( v1 = valuation(n+1, 2) );
    \\ Use a( 2^k - 1 ) = A001045( 2 + k ):
    if ( n == 2^v1 - 1 ,  return( b( v1 + 2 ) ) );

    my( k2 = 1, k = 0 );
    while ( k2 < n,  k2 <<= 1; k+=1 );
    if ( k2 > n, k2 >>= 1; k-=1 );
    my( t = n - k2 );
    \\ here  n == 2^k + 1 where k maximal

    \\ Use the following:
    \\ a( 2^k + t ) =  3 * a(t)  if  t <= 2^(k-1)
    \\ a( 2^k + 2^(k-1) + t ) =  5 * a(t)  if  t <= 2^(k-2)
    \\ a( 2^k + 2^(k-1) + 2^(k-2) + t ) =  11* a(t)  if  t <= 2^(k-3)
    \\  ... etc. ...
    \\ a( 2^k + ... + 2^(k-s) + t ) = A001045(s+2) * a(t)  if  t <= 2^((k-1)-s)
    my ( s=1 );
    while ( 1 ,
        k2 >>= 1;
        if ( t <= k2 ,  return(  b(s+2) * a(t) ) );
        t -= k2;
        s += 1;
    );
}
-------------------------------------------------------

I'd be happy if anybody could say "this looks correct",
so I could enter the program in the entry.
Tested up to n=10^3 against
\\ A071053: (inefficient)
c(n)= subst( lift( (Mod(1,2)*(1+x+x^2))^n ), x, 1);

Test up to n=10^4 is running.


Best regards,   jj