オイラープロジェクト321～330 - c21coterie @ ウィキ

問い324

以下の問題があります。

3x3xn の塔を 2x1x1 のブロックで埋める場合の数を f(n) とする。
ブロックは好きな方向に回転することができる；しかし、塔自身を回転・反転等させたものは区別して数える。
たとえば（q=100000007 として）：
f(2) = 229,
f(4) = 117805,
f(10) mod q = 96149360,
f(10^3) mod q = 24806056,
f(10^6) mod q = 30808124.
f(10^10000) mod 100000007 を求めよ。

問題の原文サイトはこちら。
http://projecteuler.net/problem=324
プロジェクトオイラーというサイトで数学の問題をプログラムで解こうという練習問題を扱ったサイトです。

自分で考え出した解法

自前アルゴリズムでこの問題を解きました。
今のところ手元のパソコンに6時間30分ほど計算させて正しい答えを出せるコードまで高速化できましたがそこから先に進めていません。
1分ルール（プロジェクトオイラーに掲載されている問題は正しいコードを書けば一分以内で問題の答えが出るというルールがある）等夢のまた夢ですね。
ｎＣｒとパスカルの三角形より高度な授業を受けたことがない私ではここら辺が限界なのかもしれません。

この問題ですが何か非常に高速な漸化式があるのでしょうか？
数式を思いつかないのでとりあえず動的計画法で解こうとしてみました。

まず１段のブロックタワーを作り、１段を２つ繋げて２段を作ります。
２段を２つ結合で4段、４段を２つ結合して８段、８段を２つ結合して１６段、、として10＾10000段に近い２＾ｎの段数まで到達し後は10＾10000を２進数に変換２進数で1のついたブロックの塊だけを結合して組み合わせ数を掛け算していけば規則的に計算できます。

例えば１０段のブロックタワーで計算するなら
10を２進数表示すると1010で２段と８段の結合となります。
まず最初に２段が求まった時２段で底が凸凹してないものの組み合わせを答えの計算用に設定します。
次に８段のブロックタワーが出来た時、答えとして設定した２段の上に８段でのせれるものをのせその組み合わせ数を計算します。
後は１０段が出来たので上面が凸凹してない物の組み合わせ数を答えとして出力します。

速度は遅いが正しい答えを出すコードも掲載しておきます。
コード記述者情報
小学校の算数までしかできない。ガイキチ。頭おかしい。ひたすらキモイやつ。大勢からの笑い物と創価学会員の間で評判の堀江伸一こと私によるコードです。

今回の問題について解法とコードは自前で考えたものであることを記述しておきます。
他の評判についてはともかく噂どおり小学校の算数までしかできない人間であるかどうかはコードを読んだ方の判断にお任せします。

正答した後に知ったのですが正しいコードは漸化式から一般項を導き出して一瞬で解くそうです。

少しだけ独り言
知恵袋でＢＡを獲得しました
http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q14110242764?fr=chie_my_notice_ba
勉強が不得手な私がこの質問でＢＡなんかもらってもよいのだろうか？

以下自前コード

#include<stdio.h>
#include<string.h>
#include<bitset>
#include<iostream>
 
const int bitSize=34000;
unsigned long long int memo[2][512][512];
std::bitset<bitSize> bits10no10000zyo;
unsigned long long int perm[2][512][512],ans[512],nextAns[512];
const unsigned long long int mod=100000007;
 
 
void my_add(std::bitset<bitSize> add){
   int c=0;
       for(int i=0;i<bitSize;i++){
               int a=bits10no10000zyo.test(i);
               int b=add.test(i);
               bits10no10000zyo.set(i,(a+b+c)%2);
               c=(a+b+c)>1;
       }
}

void mult10(){
        std::bitset<bitSize> temp=bits10no10000zyo<<3;
        bits10no10000zyo=bits10no10000zyo<<1;
        my_add(temp);
}
 
void seedPermSet(int seedSet[3][3],int r,int c){
       if(r==3){
               int d=1,hashD=0,hashU=0;
               for(int i=0;i<3;i++){
                       for(int j=0;j<3;j++){
                               hashD+=d*(seedSet[i][j]==0);
                               d*=2;
                       }
               }
               d=1;
               for(int i=0;i<3;i++){
                       for(int j=0;j<3;j++){
                               hashU+=d*(seedSet[i][j]==2);
                               d*=2;
                       }
               }
               perm[0][hashD][hashU]++;
       }else if(c==3){
               c=0;
               r++;
               seedPermSet(seedSet,r,c);
       }else{
               int c1 = c+1;
               if(seedSet[r][c]>0){
                       seedPermSet(seedSet,r,c1);
               }else{
                       if(c<2&&seedSet[r][c+1]==0){
                               seedSet[r][c  ]=1;
                               seedSet[r][c+1]=1;
                               seedPermSet(seedSet,r,c1);
                               seedSet[r][c  ]=0;
                               seedSet[r][c+1]=0;
                       }
                       if(r<2&&seedSet[r+1][c]==0){
                               seedSet[r  ][c]=1;
                               seedSet[r+1][c]=1;
                               seedPermSet(seedSet,r,c1);
                               seedSet[r  ][c]=0;
                               seedSet[r+1][c]=0;
                       }
                       seedSet[r][c]=2;
                       seedPermSet(seedSet,r,c1);
                       seedSet[r][c]=0;
                       seedPermSet(seedSet,r,c1);
               }
               
       }
}
 
 
int main(){
       memset(memo,0,sizeof(memo));
       bits10no10000zyo=10;
       for(int i=1;i<10000;i++){
               mult10();
       }
       
       int seedSet[3][3];
       memset(perm,0,sizeof(perm));
       memset(ans,0,sizeof(ans));
       
       
       
       memset(seedSet,0,sizeof(seedSet));
       seedPermSet(seedSet,0,0);
       bool isFirst=true;
       for(int i=0;i<33300;i++){
               printf("２＾%d段を計算",i);
               int next=(i+1)%2;
               int now=i%2;
               if(bits10no10000zyo.test(i)){
                       printf("ok");
                       memset(nextAns,0,sizeof(nextAns));
                       if(isFirst==true){
                               isFirst=false;
                               for(int j=0;j<512;j++){
                                       nextAns[j]=perm[now][0][j];
                               }
                       }else{
                               for(int j=0;j<512;j++){
                                       for(int k=0;k<512;k++){
                                               nextAns[k]+=(ans[j]*perm[now][j][k])%mod;
                                               nextAns[k]%=mod;
                                       }
                               }
                       }
                       memcpy(ans,nextAns,sizeof(nextAns));
               }
               
               memset(perm[next],0,sizeof(perm[next]));
               for(int j=0;j<512;j++){
                       for(int k=0;k<512;k++){
                               if(perm[now][j][k]==0)continue;
                               for(int l=0;l<512;l++){
                                       perm[next][j][l]+=(perm[now][j][k]*perm[now][k][l])%mod;
                                       perm[next][j][l]%=mod;
                               }
                       }
                }
 
                printf("\n");
        }
        std::cout<<"10^10000段の答えは"<<ans[0];
}