Kawazoe(@riverplus)さんのCodeIqの問題です。
問題はこちらの御世話になりました。解答例もあります。
->https://blog.goo.ne.jp/r-de-r/e/8af873f3ac9434bb28ea8a73d0d18736
こちらにも情報があります。
->https://togetter.com/li/1094542

縦3個の升を左から並べ右端の調整（後で説明）をしないn列の全体の場合の数をf(n),
k列目までの場合の数を最後の列の状態によって（白黒をoxに、縦を横に変更）
ooo->g(k,1),xoo->g(k,2),oxo->g(k,3),oox->g(k,4),xox->g(k,5)
とすると
f(n)=g(n,1)+g(n,2)+g(n,3)+g(n,4)+g(n,5)
g(n,1)=g(n-1,1)+g(n-1,2)+g(n-1,3)+g(n-1,4)+g(n-1,5)=f(n-1)
g(n,2)=g(n-1,1)+g(n-2,1)+g(n-1,4)、(前が ooo,[ooo,oxo],ooxのどれか)
g(n,3)=途中では g(n-1,1)+g(n-1,2)+g(n-1,4) (前が ooo,xoo,ooxのどれか)
　　　右端では g(n-1,1),([xoo,oxo],[oox,oxo]はだめ)
g(n,4)=g(n,2)
g(n,5)=g(n-1,1) (oooのみ)

上からg(k,_)をg(k-1,_),g(k-2,_)からなる式に変換することを繰り返しつつ、
g(k,1)が現れたらf(k-1)に変換します。

f(n)=f(n-1)+g(n-1,1)*4+g(n-1,2)*4+g(n-2,1)*2
　　....................
　　=f(n-1)+f(n-2)*4+f(n-3)*6+f(n-4)*8+f(n-5)*8+...+f(2)*8+f(1)*4+g(3,2)
f(n-1)=f(n-2)+f(n-3)*4+f(n-4)*6+f(n-5)*8+f(n-6)*8+...+f(2)*8+f(1)*4+g(3,2)
従って
f(n)-f(n-1)=f(n-1)+f(n-2)*3+f(n-3)*2+f(n-4)*2
f(n)=f(n-1)*2+f(n-2)*3+f(n-3)*2+f(n-4)*2
左端から4列目までは別に計算
右端はg(n,3)が違ってきますのでg(n-1,2)=f(n-3)+f(n-4)+g(n-2,2)の2倍を引きます。

function solve(n,a)
    m = 10^6
    g=zeros(Int,5,5)
    f=zeros(Int,5)
    g[1,1],g[1,2],g[1,3],g[1,4],g[1,5] = 1,1,1,1,1
    g[2,1] = g[1,1] + g[1,2] + g[1,3] + g[1,4] + g[1,5]
    g[2,2] = g[1,1] + g[1,4]
    g[2,3] = g[1,1] + g[1,2] + g[1,4]
    g[2,4] = g[1,1] + g[1,2]  # =g[2,2]
    g[2,5] = g[1,1]
    f[1]=g[2,1]
    for i in 3:5
        g[i,1]=g[i-1,1]+g[i-1,2]+g[i-1,3]+g[i-1,4]+g[i-1,5]
        g[i,2]=g[i-1,1]+g[i-2,1]+g[i-1,4]
        g[i,3]=g[i-1,1]+g[i-1,2]+g[i-1,4]
        g[i,4]=g[i-1,1]+g[i-2,1]+g[i-1,2]  # = g[i,2]
        g[i,5]=g[i-1,1]
        f[i-1]=g[i,1]
    end
    if n == 1
      sum = f[1] - g[1,3]
    elseif n < 5
      sum = f[n] - g[n-1,2] * 2
    else
      h = g[3,2] * 2
      for i in 5 : n
        f[5] = f[1] * 2 + f[2] * 2 + f[3] * 3 + f[4] * 2
        if f[5] > 10^11
          f[5] = mod(f[5],m)
        end
        h += (f[2] + f[1]) * 2
        for j in 1 : 4
          f[j] = f[j+1]
        end  
      end
      sum = mod(f[5] - h,m)
    end
    s = a == sum ? "ok " : "no "
    println(s,sum)
end

solve(1,4)
solve(2,11)
solve(3,39)
solve(4,121)
solve(5,387) 
solve(6, 1237)
solve(7,3955)
solve(12,318221)
solve(100,193677)
solve(5190,461261)
solve(9978, 464781)
solve(70013454,184013)
solve(99999991, 885187)

prologでは時間切れとなりますので、g[]のまま解ければ速くなるのではと考え、
掛けられる方の行列の要素に、g[k,_]のほかにg[k-1,1]をくわえれば
行列計算できることに気づき解くことができました。

(..((L1*L3)*L3)*...)=L1*(L3*L3*...)を利用して先にL3のべき乗を計算します。
とはいえ入力値が大きく、普通に計算しますと計算量が多くなりすぎますので、
行列のべき乗の部分はwikipediaにしたがいました。

wikipedia　べき乗
1:指数を n とし、2乗していく値 p := a、結果値 v := 1 とする。
2:n が 0 なら、v を出力して終了する。
3:n の最下位桁が 1 なら、v := v * p とする。
4:n := [n/2] とし（端数切捨て）、 p := p * p として、2. に戻る。

solve/2のL1は[g[1,1],g[1,2],g[1,3],g[1,5],g[[0,1]]、
L2はV:=1に相当する単位行列です。
[[1,2,1,1,0],[1,1,0,0,1],[1,2,0,0,0],[1,0,0,0,0],[1,0,0,0,0]]
の転置行列L3がaに相当します。
その要素はg[k+1,1]=g[k,1]*1+g[k,2]*2+g[k,3]*1+g[k,5]*1+g[k-1,1]*0
に対応する[1,2,1,1,0] 等です。
pow/4で行列のべき乗を処理しています。L2とL3がpに、L4がvに相当します。
行列計算はmuL/3のrot/2で行と列を交換して、
muL0/3、muL1/3、cal/3で行と行の演算として行います。

前に記しましたように最後g[n-1,2]*2が余分です。補正にg[n-1,2]が必要ですので、
solveでL6(最後の列)=L5(ひとつ前の列)*L3としています。
F(n)=g[n,1]+g[n,2]*2+g[n,3]+g[n,5]-g[n-1,2]*2
　　=sum([g[n,1],g[n,2],g[n,3],g[n,5],g[n-1,1]]) +g[n,2]-g[n-1,1]-g[n-1,2]*2
　　=sum(L6)+L6[2]-L6[5]+L5[2]*2

%SWI-Prolog ver 7.42,7.64,9.25
%:-initialization(start).       %Online Prolog Compiler,ideone.    
%start.

cal(L1,L2,R):-
    maplist(mul,L1,L2,L3),sumlist(L3,R1),(R1>10^14->R is R1 mod 10^6;R = R1),!.

mul(N,X,R):-R is N*X.

ht([H|T],H,T).
rot([[]|_],[]).
rot(L,R):-maplist(ht,L,LH,LT),rot(LT,R1),R=[LH|R1].

muL1(_,[],[]).
muL1(L1,[H2|T2],R):-muL1(L1,T2,R1),cal(L1,H2,V),R=[V|R1].

muL0([],_,[]).
muL0([H1|T1],L2,R):-muL0(T1,L2,R1),muL1(H1,L2,V),R=[V|R1],!.

muL(L1,L2,R):-rot(L2,L21),muL0(L1,L21,R).

pow(0,R,_,R).
pow(N,L1,L2,R):-
    divmod(N,2,D,M),muL(L2,L2,L3),(M =:=1->muL(L1,L2,L4);L4=L1),pow(D,L4,L3,R),!.
    
solve(1,4).
solve(N,R):-
    L1=[1,1,1,1,0],L2=[[1,0,0,0,0],[0,1,0,0,0],[0,0,1,0,0],[0,0,0,1,0],[0,0,0,0,1]],
    rot([[1,2,1,1,0],[1,1,0,0,1],[1,2,0,0,0],[1,0,0,0,0],[1,0,0,0,0]],L3),
   (N>2->(N1 is N-2,pow(N1,L2,L3,L4),muL([L1],L4,[L5]));L5=L1),
    muL([L5],L3,[L6]),sumlist(L6,S),nth1(5,L6,X),nth1(2,L6,Y),nth1(2,L5,Z),
    R is (S-X+Y-Z*2) mod 10^6,!.

start:-f(N,A),(solve(N,R)->(R==A->Str=" ok  ";Str=" no  ");
    write(fail)),write(" "),write(Str),writeln(N->R),fail.
start.

f(1,4).
f(2,11).
f(3,39).
f(4,121).
f(5,387). 
f(6, 1237).
f(7,3955).
f(12,318221).
f(100,193677).
f(5190,461261).
f(9978, 464781).
f(70013454,184013).
f(99999991, 885187).