プロジェクトオイラー問520
ｎ桁の数字の中の
偶数個の偶数の数、奇数個の偶数の数、奇数個の奇数の数、偶数個の奇数の数
を主キーに漸化式を立て。
後は一桁と一桁で２桁。
２桁と２桁で４桁
４桁と４桁で８桁。
、、、
と
一桁を作り
1桁と２桁で３桁
３桁と２桁を１桁に足す。
１桁+２桁+３桁と4桁の直積で5,6,7桁
これに４桁を足して4,5、6,7、これに1,2,3桁を足せば７ケタ。
以下同様。
計算時間１０分。
問題の
#include<stdio.h>
#include<iostream>
#include<string.h>
#include<map>
const int MOD=1000000123;
struct E{
int j,k,l,m;
bool operator<(const E& e)const{
if(j!=e.j)return j<e.j;
if(k!=e.k)return k<e.k;
if(l!=e.l)return l<e.l;
return m<e.m;
}
void set(int j1,int k1,int l1,int m1){
j=j1;
k=k1;
l=l1;
m=m1;
}
void print(){
std::cout<<"("<<j<<" "<<k<<" "<<l<<" "<<m<<")";
}
};
int firstSet(E e1){
int j=e1.j;
int k=e1.k;
int l=e1.l;
int m=e1.m;
if((j==0)&&(k==1)&&(l==0)&&(m==0))return 4;
if((j==0)&&(k==0)&&(l==1)&&(m==0))return 5;
return 0;
}
bool secondSet(E e1){
int j=e1.j;
int k=e1.k;
int l=e1.l;
int m=e1.m;
if((k>0)||(m>0))return false;
if((j>0)||(l>0))return true;
return false;
}
int main(){
//一桁の計算量が411*6
//Q(100)
__int64 as[7][7][7][7],as2[7][7][7][7];//[偶数個の偶数][奇数個の偶数の数][奇数個の奇数の数][偶数 個の奇数の数]
__int64 ans=5;
memset(as,0,sizeof(as));
as[0][1][0][0]=4;
as[0][0][1][0]=5;
std::map<E,std::map<E,__int64> > count,sums,nori;
std::map<E,std::map<E,__int64> > allCount;
//まずは１桁と１桁で２桁を作り、その過程で漸化式の構造をcount変数で把握する
//２段と２段で４段、４段と４段で
//以後は１段と２段で３段を作る
//４段と（３，２，１）段で（7,6,5）+4+（３，２，１）段を求め
//この調子で大きな段を計算して行く
E e1;
e1.set(0,1,0,0);
allCount[e1][e1]=1;
e1.set(0,0,1,0);
allCount[e1][e1]=1;
for(int i=2;i<=2;i++){
memset(as2,0,sizeof(as2));
for(int j=0;j<6;j++){
for(int k=0;k+j<6;k++){
for(int l=0;l<6;l++){
for(int m=0;m+l<6;m++){
//偶数個の偶数を一つ増やす
__int64 c=as[j][k][l][m];
E e1,e2;
e1.set(j,k,l,m);
if(k<5){
if(j>0){
as2[j-1][k+1][l][m]=(as2[j-1][k+1][l][ m]+c*j)%MOD;
e2.set(j-1,k+1,l,m);
count[e2][e1]+=j;
nori[e2][e1]+=j;
}
if(j+k<5){
as2[j][k+1][l][m]=(as2[j][k+1][l][m]+c*(5-j-k ))%MOD;
e2.set(j,k+1,l,m);
count[e2][e1]+=(5-j-k);
nori[e2][e1]+=(5-j-k);
}
}
//奇数個の偶数を一つ増やす
if((0<k)&&(j<5)){
as2[j+1][k-1][l][m]=(as2[j+1][k-1][l][m]+c*k)%MOD;
e2.set(j+1,k-1,l,m);
count[e2][e1]+=k;
nori[e2][e1]+=k;
}
//奇数個の奇数を一つ増やす
if((l>0)&&(m<5)){
as2[j][k][l-1][m+1]=(as2[j][k][l-1][m+1]+c*l)%MOD;
e2.set(j,k,l-1,m+1);
count[e2][e1]+=l;
nori[e2][e1]+=l;
}
//偶数個の奇数を一つ増やす。
if(l<5){
if(m>0){
as2[j][k][l+1][m-1]=(as2[j][k][l+1][m-1 ]+c*m)%MOD;
e2.set(j,k,l+1,m-1);
count[e2][e1]+=m;
nori[e2][e1]+=m;
}
if(l+m<5){
as2[j][k][l+1][m]=(as2[j][k][l+1][m]+c*(5-l-m ))%MOD;
e2.set(j,k,l+1,m);
count[e2][e1]+=(5-l-m);
nori[e2][e1]+=(5-l-m);
}
}
}
}
}
}
memcpy(as,as2,sizeof(as2));
__int64 add=0;
for(int i=0;i<6;i++){
for(int j=0;j<6;j++){
add=(add+as[i][0][j][0])%MOD;
}
}
}
ans=29;
for(__int64 i=4;i<=pow(2,39);i*=2){
//まずは２倍
std::map<E,std::map<E,__int64> > count2,allCount2;
std::map<E,std::map<E,__int64> >::iterator itA,itNori;
std::map<E,__int64>::iterator itA2,itNori2,itB;
E e1,e2,e3,e4,e5;
__int64 p1,p2,p3,pAll,cc=0;
//一回目は1段+２段＝３段を計算
for(itA=count.begin();itA!=count.end();itA++){
e1=(*itA).first;
for(itA2=(*itA).second.begin();itA2!=(*itA).second.end();itA2++){
e2=(*itA2).first;
p1=(*itA2).second;
for(itNori2=nori[e2].begin();itNori2!=nori[e2].end();itNori2++){
e3=(*itNori2).first;
p2=(*itNori2).second;
for(itB=allCount[e3].begin();itB!=allCount[e3].end();itB++){
e4=(*itB).first;
p3=(*itB).second;
pAll=(p1*p2)%MOD;
pAll=(pAll*p3)%MOD;
allCount2[e1][e4]+=pAll;
allCount2[e1][e4]%=MOD;
cc++;
}
}
}
}
//最初は３段を追加
for(itA=allCount2.begin();itA!=allCount2.end();itA++){
E e1=(*itA).first;
for(itA2=(*itA).second.begin();itA2!=(*itA).second.end();itA2++){
E e2=(*itA2).first;
allCount[e1][e2]+=(*itA2).second;
allCount[e1][e2]%=MOD;
}
}
//２＾i段を追加
for(itA=count.begin();itA!=count.end();itA++){
e1=(*itA).first;
for(itA2=(*itA).second.begin();itA2!=(*itA).second.end();itA2++){
e2=(*itA2).first;
allCount[e1][e2]+=(*itA2).second;
allCount[e1][e2]%=MOD;
}
}
__int64 add2=0;
for(itA=allCount.begin();itA!=allCount.end();itA++){
E e1=(*itA).first;
for(itA2=(*itA).second.begin();itA2!=(*itA).second.end();itA2++){
E e2=(*itA2).first;
cc++;
if((firstSet(e2)>0)&&(secondSet(e1))){
add2+=firstSet(e2)*allCount[e1][e2];
add2%=MOD;
}
}
}
std::cout<<i<<" "<<add2<<" ";
//2^n段だけを計算
for(itA=count.begin();itA!=count.end();itA++){
e1=(*itA).first;
for(itA2=(*itA).second.begin();itA2!=(*itA).second.end();itA2++){
e2=(*itA2).first;
p1=(*itA2).second;
for(itNori2=nori[e2].begin();itNori2!=nori[e2].end();itNori2++){
e3=(*itNori2).first;
p2=(*itNori2).second;
for(itB=count[e3].begin();itB!=count[e3].end();itB++){
e4=(*itB).first;
p3=(*itB).second;
pAll=(p1*p2)%MOD;
pAll=(pAll*p3)%MOD;
count2[e1][e4]+=pAll;
count2[e1][e4]%=MOD;
cc++;
}
}
}
}
count.clear();
count.insert(count2.begin(),count2.end());
__int64 add=0;
for(itA=count2.begin();itA!=count2.end();itA++){
E e1=(*itA).first;
for(itA2=(*itA).second.begin();itA2!=(*itA).second.end();itA2++){
E e2=(*itA2).first;
if((firstSet(e2)>0)&&(secondSet(e1))){
add+=firstSet(e2)*count[e1][e2];
add%=MOD;
}
}
}
std::cout<<i<<" cc="<<cc<<" add="<<add<<" all="<<add+add2<<"\n";
ans+=add+add2;
ans%=MOD;
}
std::cout<<"Σi=1～39 Q(2^i)="<<ans;
}
最終更新：2015年07月21日 08:21
