タイリング問題を解いてポイントもらおう

http://q.hatena.ne.jp/1228746034
アニーリングの出番ですよ。エネルギーの設定が肝か。

コード書いてみた。現在の記録の19個は10メガ Stepsかけて冷やすと出てきた。a.out 19 10000000 とかすると19個の配列を見付けてプリントします。20個はまだ出てこない。バリチューンしてあるんで可読性０
２種類のペナルティを使い、１つはピースの形が許されていない物の場合、２つめは各ピースの重なり。はじめはペナルティをゆるくしてランダムに動かし、だんだんペナルティを厳しくしていく。ペナルティが０になればそれが解。

以下Cのコード
必殺 Flat Histogram 法
単に a.out 17 とかすると解をみつけます。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
// box size
#define LX 13
#define LY 19

// Max pieces
#define NN 40
// # of pieces
int nn;

// Total penalty = K1 E1 + K2 E2
// E1=disallowed shape penalty
// E2=overlap penalty
//
// K1, K2 must be integer
#define K1 1
#define K2 1

// flat histogram related
#define E1MAX 100
#define E2MAX 20
#define EMAX (E1MAX*K2+E2MAX*K1)

// decrement of weight per hit
#define PIP 0.01
// energy histogram
int hist[EMAX];
// weight function W(E)
//  Instead of the Boltzmann factor exp(-beta*E), we use exp(W(E))
//   so that the histogram becomes flat, that is, W(E)=-Log N(E) where
//   N(E) is number of state which has energy E.
double ww[EMAX];

// current energy
int enow;
int enow1;
int enow2;

// (cx, cy) center pos of pieces (crossing square)
int cx[NN],cy[NN];
// displacement of horiz/vert bars
int dx[NN],dy[NN];
// dx=0    dx=0     dx=1
// dy=-1   dy=-2    dy=-1 (disallowed)
//  #        #        #
//  #        #        #
//  #        #        #
//#####      #       #####
//  #      #####      #

// buffer to calculate overlap
char pbuf[LX][LY];
char pbuf2[LX][LY];
char dbuf[LX][LY];

#define RMAX 0xfffffff

int ene1i(int i)
{
    int a1=abs(dx[i]);
    int a2=abs(dy[i]);
    if(a2<a1) a1=a2;
    return a1;
}

int ene1()
{
    int res=0;
    int i;
    for(i=0;i<nn;i++)
    {
	res+= ene1i(i);
    }
    return res;
}

int puti(char abuf[LX][LY], int i)
{
    int x0=cx[i];
    int y0=cy[i];
    int ene=0;
    int x,y;

    x=x0+dx[i]-2; y=y0;
    char *ptr = &(abuf[x][y]);

    ene += (*ptr); (*ptr)++;ptr+=LY;
    ene += (*ptr); (*ptr)++;ptr+=LY;
    ene += (*ptr); (*ptr)++;ptr+=LY;
    ene += (*ptr); (*ptr)++;ptr+=LY;
    ene += (*ptr); (*ptr)++;

    x=x0; y=y0+dy[i]-2;
    ptr = &(abuf[x][y]);

    ene += (*ptr); (*ptr)++;ptr++;
    ene += (*ptr); (*ptr)++;ptr++;
    ene += (*ptr); (*ptr)++;ptr++;
    ene += (*ptr); (*ptr)++;ptr++;
    ene += (*ptr); (*ptr)++;ptr++;

    x=x0; y=y0;
    ptr = &(abuf[x][y]);
    (*ptr)--; ene -= *ptr;

    //printf("Put %d,%d + %d,%d\n", x0,y0, dx[i],dy[i]);
    return ene;
}

int deli(char abuf[LX][LY], int i)
{
    int x0=cx[i];
    int y0=cy[i];
    int ene=0;
    int x,y;

    x=x0+dx[i]-2; y=y0;
    char *ptr = &(abuf[x][y]);
    (*ptr)--;  ene -= *ptr; ptr+=LY;
    (*ptr)--;  ene -= *ptr; ptr+=LY;
    (*ptr)--;  ene -= *ptr; ptr+=LY;
    (*ptr)--;  ene -= *ptr; ptr+=LY;
    (*ptr)--;  ene -= *ptr; 

    x=x0; y=y0+dy[i]-2;
    ptr = &(abuf[x][y]);
    (*ptr)--;  ene -= *ptr; ptr++;
    (*ptr)--;  ene -= *ptr; ptr++;
    (*ptr)--;  ene -= *ptr; ptr++;
    (*ptr)--;  ene -= *ptr; ptr++;
    (*ptr)--;  ene -= *ptr;

    x=x0; y=y0;
    ptr = &(abuf[x][y]);
    ene += *ptr; (*ptr)++;

    //printf("Put %d,%d + %d,%d\n", x0,y0, dx[i],dy[i]);
    return ene;
}

void clr(char abuf[LX][LY])
{
    int x,y;
    for(x=0;x<LX;x++)
	for(y=0;y<LY;y++)
	    abuf[x][y]=0;
}

void cpb(char abuf[LX][LY], char abuf2[LX][LY])
{
    memcpy(abuf2, abuf, sizeof(char)*LX*LY);
}

int ene2(char abuf[LX][LY])
{
    int x,y;

    //for(i=0;i<NN;i++) puti(abuf,i);

    int res=0;
    for(x=0;x<LX;x++)
    {
	for(y=0;y<LY;y++)
	{
	    int n=abuf[x][y];
	    res += n*(n-1)/2;
	}
    }

    return res;
}

void init_ww()
{
    int i;
    for(i=0;i<EMAX;i++) {
	hist[i]=0;
	ww[i]=0.0;
    }
}

void upd_ww()
{
    int i;
    FILE *fp=fopen("h","w");

    for(i=0;i<EMAX;i++) {
	fprintf(fp,"%d %d %f\n",i, hist[i], ww[i]);
    }
    fclose(fp);
}

void read_ww(char *fn)
{
    int i;
    char buf[300];
    FILE *fp=fopen(fn,"r");
    for(i=0;i<EMAX;i++) {
	fgets(buf, 222,fp);
        int d1,d2;
        sscanf(buf,"%d %d %lf",&d1,&d2, &ww[i]);
    }
    fclose(fp);
}

int reject(int e, int enew)
{

    if(e>=EMAX || enew>=EMAX)
    {
	if(enew<e) return 0;
	return 1;
    }
    double dw = ww[enew]-ww[e];
    if(dw>0.0) return 0;

    int dice=rand()&RMAX;
    int thresh=(int)(RMAX*exp(dw));
    if(dice>thresh) return 1;
    return 0;
}

//move dx
void move1(int i)
{
    int rr=rand()>>10;
    int ndx = dx[i];
    int ndy = dy[i];
    int dx0=ndx;
    int dy0=ndy;
    const int newdx[2][5]={ {-1,0,1,2,1},{-1,-2,-1,0,1}};
    if(rr&1)
    {
	//move dx
	rr>>=1;
	ndx = newdx[rr&1][ndx+2];
	if( cx[i]+ndx-2 <0 || cx[i]+ndx+2 >=LX)return;
    }
    else
    {
	rr>>=1;
	ndy = newdx[rr&1][ndy+2];
	if(cy[i]+ndy-2 <0 || cy[i]+ndy+2>=LY) return;
    }
    // valid
    cpb(pbuf, pbuf2);
    int de2=deli(pbuf2, i);
    //printf("Move %d %d to %d %d\n", dx[i],dy[i],ndx,ndy);
    int e1now = ene1i(i);
    dx[i]=ndx;
    dy[i]=ndy;
    int e1new = ene1i(i);
    de2 += puti(pbuf2,i);

    int de1 = e1new-e1now;
    int ide = de1*K1 + de2*K2;
    int enew = enow+ide;

    if(reject(enow, enew)==1)
    {
	//reject
	dx[i]=dx0;
	dy[i]=dy0;
	return;
    }
    //accept
    cpb(pbuf2, pbuf);
    enow = enew;
    enow1 += de1;
    enow2 += de2;
}

//move dx
void move2(int i)
{

    int rr=rand()>>10;
    int ncx = cx[i];
    int ncy = cy[i];
    int cx0=ncx;
    int cy0=ncy;
    if(rr&1)
    {
	//move cx
	rr>>=1;
	if(rr&1) ncx++;
	else ncx--;
	if((ncx+dx[i]+2>=LX) || (ncx+dx[i]-2<0)) return;
    }
    else
    {
	//move cy
	rr>>=1;
	if(rr&1) ncy++;
	else ncy--;
	if(ncy+dy[i]+2>=LY || ncy+dy[i]-2<0) return;
    }
    // valid
    cpb(pbuf, pbuf2);
    int de2 = deli(pbuf2, i);
    //printf("Move %d %d to %d %d\n", cx[i],cy[i],ncx,ncy);
    cx[i]=ncx;
    cy[i]=ncy;
    int de22 = puti(pbuf2,i);
    de2 += de22;
    de2*=K2;
    int enew = enow + de2;

    if(reject(enow, enew)==1)
    {
	//reject
	cx[i]=cx0;
	cy[i]=cy0;
	return;
    }
    //printf("Accept de=%e  %d %d\n", de, ene2(pbuf), ene2(pbuf2));
    //accept
    cpb(pbuf2, pbuf);
    enow = enew;
    enow2 += de2;
}

void move3(int i)
{
    int rr=rand()>>10;
    int ndx = dx[i];
    int ndy = dy[i];
    int ncx = cx[i];
    int ncy = cy[i];
    int dx0=ndx;
    int dy0=ndy;
    int cx0=ncx;
    int cy0=ncy;
    const int newdx[2][5]={ {1,-1,-1,-1,-1},{1,1,1,1,-1}};
    if(rr&1)
    {
	//move dx
	rr>>=1;
	rr=rr&1;
	int dd = newdx[rr][ndx+2];
	ndx += dd;
	ncx -= dd;
    }
    else
    {
	rr>>=1;
	rr=rr&1;
	int dd = newdx[rr][ndy+2];
	ndy += dd;
	ncy -= dd;
    }
    // valid
    cpb(pbuf, pbuf2);
    int de2=deli(pbuf2, i);
    //printf("Move %d %d to %d %d\n", dx[i],dy[i],ndx,ndy);
    int e1now = ene1i(i);
    dx[i]=ndx;
    dy[i]=ndy;
    cx[i]=ncx;
    cy[i]=ncy;
    int e1new = ene1i(i);
    de2 += puti(pbuf2,i);

    int de1 = e1new-e1now;
    int ide = de1*K1 + de2*K2;
    int enew = enow + ide;

    if(reject(enow, enew)==1)
    {
	//reject
	dx[i]=dx0;
	dy[i]=dy0;
	cx[i]=cx0;
	cy[i]=cy0;
	return;
    }
    //accept
    cpb(pbuf2, pbuf);
    enow = enew;
    enow1 += de1;
    enow2 += de2;
}

void prt()
{
    int x,y;
    for(y=0;y<LY;y++)
    {
	for(x=0;x<LX;x++)
	{
	    putchar('0'+pbuf[x][y]);
	}
	putchar(10);
    }
    putchar(10);
}


void dump()
{
    int x,y,i;
    for(x=0;x<LX;x++)
    for(y=0;y<LY;y++)dbuf[x][y]='.';

    for(i=0;i<nn;i++)
    {
	int j;
	for(j=-2;j<=2;j++)
	{
	    x=cx[i]+dx[i]+j;
	    y=cy[i];
	    dbuf[x][y] = (pbuf[x][y]!=1) ? '*':'A'+i;
	}
	for(j=-2;j<=2;j++)
	{
	    x=cx[i];
	    y=cy[i]+dy[i]+j;
	    dbuf[x][y] = (pbuf[x][y]!=1) ? '*':'A'+i;
	}
    }
    for(y=0;y<LY;y++)
    {
	putchar('#');

	for(x=0;x<LX;x++)
	    putchar(dbuf[x][y]);
	putchar(10);
    }
    int e1=ene1();
    int e2= ene2(pbuf);
    printf("#Penalty %d %d\n\n", e1,e2);
}

int main(int argc, char **argv)
{
    //int cnt=atoi(argv[2]);
    nn =atoi(argv[1]);

    int i,s;

    clr(pbuf);
    for(i=0;i<nn;i++)
    {
	dx[i]=dy[i]=0;
	cx[i]=2+(rand()>>10)%(LX-4);
	cy[i]=2+(rand()>>10)%(LY-4);
	puti(pbuf,i);
    }

    enow1=ene1();
    enow2=ene2(pbuf);
    enow = enow1*K1+enow2*K2;
    init_ww();
    if(argc==3) read_ww(argv[2]);

    //int hmax =UPD_HIST*nn*3;
    int ebest=enow;
    int stotal=0;

    while(1)
    {
	for(s=0;s<100000;s++)
	{
	    stotal++;

	    for(i=0;i<nn;i++)
	    {
		move1(i);
		hist[enow]++;
		ww[enow] -= PIP;
	        if(enow==0)break;
	    }
	    for(i=0;i<nn;i++)
	    {
		move2(i);
		hist[enow]++;
		ww[enow] -= PIP;
	        if(enow==0)break;
	    }
	    for(i=0;i<nn;i++)
	    {
		move3(i);
		hist[enow]++;
		ww[enow] -= PIP;
	        if(enow==0)break;
	    }

	    if(enow==0)break;
	    if(enow<ebest)
	    {
		ebest=enow;
                dump();
	    }
	}
        if(enow==0)break;
	upd_ww();
    }

    dump();
    printf("# Solution found after %d steps\n", stotal);
    return 0;
}

こちらはもっと基本的だけど遅いアニーリング版
a.out 17 100000 とかすると100000ステップ冷やして見付けます。ステップが足りない（速く冷やしすぎる）と解に到達できません。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define LX 13
#define LY 19

// max piece
#define NN 40
// # of piece
int nn;

// energy coefficient
// 1: disallowed piece shape penalty
// 2: overlap penalty
#define K1 1
#define K2 1

// 1/temperature
double beta;

// annealing schedule
double bmin=0.1;
double bmax=6.0;

#define RMAX 0xfffffff
#define DEMAX 20

// Metropolis table
int thtab[DEMAX];
void init_th(double b)
{
    int i;
    thtab[0]=RMAX/2;
    for(i=1;i<DEMAX;i++) thtab[i]=(int)(RMAX*exp(-i*b));
}

// center and shape of pieces
int cx[NN],cy[NN], dx[NN],dy[NN];
char pbuf[LX][LY];
char pbuf2[LX][LY];
char dbuf[LX][LY];

// disallowed shape penalty
int ene1i(int i)
{
    int a1=abs(dx[i]);
    int a2=abs(dy[i]);
    if(a2<a1) a1=a2;
    return a1;
}

int ene1()
{
    int res=0;
    int i;
    for(i=0;i<nn;i++)
    {
	res+= ene1i(i);
    }
    return res;
}

int puti(char abuf[LX][LY], int i)
{
    int x0=cx[i];
    int y0=cy[i];
    int ene=0;
    int x,y;

    x=x0+dx[i]-2; y=y0;
    char *ptr = &(abuf[x][y]);

    ene += (*ptr); (*ptr)++;ptr+=LY;
    ene += (*ptr); (*ptr)++;ptr+=LY;
    ene += (*ptr); (*ptr)++;ptr+=LY;
    ene += (*ptr); (*ptr)++;ptr+=LY;
    ene += (*ptr); (*ptr)++;

    x=x0; y=y0+dy[i]-2;
    ptr = &(abuf[x][y]);

    ene += (*ptr); (*ptr)++;ptr++;
    ene += (*ptr); (*ptr)++;ptr++;
    ene += (*ptr); (*ptr)++;ptr++;
    ene += (*ptr); (*ptr)++;ptr++;
    ene += (*ptr); (*ptr)++;ptr++;

    x=x0; y=y0;
    ptr = &(abuf[x][y]);
    (*ptr)--; ene -= *ptr;

    //printf("Put %d,%d + %d,%d\n", x0,y0, dx[i],dy[i]);
    return ene;
}

int deli(char abuf[LX][LY], int i)
{
    int x0=cx[i];
    int y0=cy[i];
    int ene=0;
    int x,y;

    x=x0+dx[i]-2; y=y0;
    char *ptr = &(abuf[x][y]);
    (*ptr)--;  ene -= *ptr; ptr+=LY;
    (*ptr)--;  ene -= *ptr; ptr+=LY;
    (*ptr)--;  ene -= *ptr; ptr+=LY;
    (*ptr)--;  ene -= *ptr; ptr+=LY;
    (*ptr)--;  ene -= *ptr; 

    x=x0; y=y0+dy[i]-2;
    ptr = &(abuf[x][y]);
    (*ptr)--;  ene -= *ptr; ptr++;
    (*ptr)--;  ene -= *ptr; ptr++;
    (*ptr)--;  ene -= *ptr; ptr++;
    (*ptr)--;  ene -= *ptr; ptr++;
    (*ptr)--;  ene -= *ptr;

    x=x0; y=y0;
    ptr = &(abuf[x][y]);
    ene += *ptr; (*ptr)++;

    //printf("Put %d,%d + %d,%d\n", x0,y0, dx[i],dy[i]);
    return ene;
}

void clr(char abuf[LX][LY])
{
    int x,y;
    for(x=0;x<LX;x++)
	for(y=0;y<LY;y++)
	    abuf[x][y]=0;
}

void cpb(char abuf[LX][LY], char abuf2[LX][LY])
{
    memcpy(abuf2, abuf, sizeof(char)*LX*LY);
}

// overlap penalty
int ene2(char abuf[LX][LY])
{
    int x,y;

    //for(i=0;i<NN;i++) puti(abuf,i);

    int res=0;
    for(x=0;x<LX;x++)
	for(y=0;y<LY;y++)
	    if(abuf[x][y]>=2) res += abuf[x][y]-1;

    return res;
}

//move dx dy
void move1(int i)
{
    int e1now = ene1i(i);

    int rr=rand()>>10;
    int ndx = dx[i];
    int ndy = dy[i];
    int dx0=ndx;
    int dy0=ndy;
    if(rr&1)
    {
	//move dx
	rr>>=1;
	if(rr&1)
	{
	    ndx++;
	    if(ndx==3 || cx[i]+ndx+2 >=LX)return;
	}
	else
	{
	    ndx--;
	    if(ndx==-3 || cx[i]+ndx-2 <0) return;
	}
    }
    else
    {
	rr>>=1;
	if(rr&1) ndy++;
	else ndy--;
	if(ndy==-3 || ndy==3 || cy[i]+ndy+2>=LY || cy[i]+ndy-2 <0) return;
    }
    // valid
    cpb(pbuf, pbuf2);
    int de2=deli(pbuf2, i);
    //printf("Move %d %d to %d %d\n", dx[i],dy[i],ndx,ndy);
    dx[i]=ndx;
    dy[i]=ndy;
    int e1new = ene1i(i);
    de2 += puti(pbuf2,i);

    int ide = (e1new-e1now)*K1 + de2*K2;

    if(ide>=0)
    {
	int dice=rand()&RMAX;
	if(ide>=DEMAX) ide=DEMAX-1;
	int thresh=thtab[ide];
	if(dice>thresh)
	{
	    //reject
	    dx[i]=dx0;
	    dy[i]=dy0;
	    return;
	}
    }
    //accept
    cpb(pbuf2, pbuf);
}

//move cx cy
void move2(int i)
{

    int rr=rand()>>10;
    int ncx = cx[i];
    int ncy = cy[i];
    int cx0=ncx;
    int cy0=ncy;
    if(rr&1)
    {
	//move cx
	rr>>=1;
	if(rr&1) ncx++;
	else ncx--;
	if((ncx+dx[i]+2>=LX) || (ncx+dx[i]-2<0)) return;
    }
    else
    {
	rr>>=1;
	if(rr&1) ncy++;
	else ncy--;
	if(ncy+dy[i]+2>=LY || ncy+dy[i]-2<0) return;
    }
    // valid
    cpb(pbuf, pbuf2);
    int de2 = deli(pbuf2, i);
    //printf("Move %d %d to %d %d\n", cx[i],cy[i],ncx,ncy);
    cx[i]=ncx;
    cy[i]=ncy;
    int de22 = puti(pbuf2,i);
    de2 += de22;
    de2*=K2;

    if(de2>=0)
    {
	if(de2>=DEMAX) de2=DEMAX-1;
	int dice=rand()&RMAX;
	int thresh=thtab[de2];

	if(dice>thresh)
	{
	    //reject
	    cx[i]=cx0;
	    cy[i]=cy0;
	    return;
	}
    }
    //printf("Accept de=%e  %d %d\n", de, ene2(pbuf), ene2(pbuf2));
    //accept
    cpb(pbuf2, pbuf);
}

void prt()
{
    int x,y;
    for(y=0;y<LY;y++)
    {
	for(x=0;x<LX;x++)
	{
	    putchar('0'+pbuf[x][y]);
	}
	putchar(10);
    }
    putchar(10);
}


void dump()
{
	int x,y,i;
	for(x=0;x<LX;x++)
	    for(y=0;y<LY;y++)dbuf[x][y]='.';

	for(i=0;i<nn;i++)
	{
	    int j;
	    for(j=-2;j<=2;j++)
	    {
		dbuf[cx[i]+dx[i]+j][cy[i]] = 'A'+i;
		if(pbuf[cx[i]+dx[i]+j][cy[i]]!=1) dbuf[cx[i]+dx[i]+j][cy[i]] ='*';
	    }
	    for(j=-2;j<=2;j++)
	    {
		dbuf[cx[i]][cy[i]+dy[i]+j] = 'A'+i;
		if(pbuf[cx[i]][cy[i]+dy[i]+j]!=1) dbuf[cx[i]][cy[i]+dy[i]+j] ='*';
	    }
	}
	for(y=0;y<LY;y++)
	{
	    putchar('#');

	    for(x=0;x<LX;x++)
		putchar(dbuf[x][y]);
	    putchar(10);
	}
	int e1=ene1();
	int e2= ene2(pbuf);
	printf("#Ene %d %d\n", e1,e2);
}


int main(int argc, char **argv)
{
    int cnt=atoi(argv[2]);
    nn =atoi(argv[1]);
    int i,s;
    double db=exp( log(bmax/bmin)/cnt);
    int ee;


    clr(pbuf);
    for(i=0;i<nn;i++)
    {
	dx[i]=dy[i]=0;
	cx[i]=2+(rand()>>10)%(LX-4);
	cy[i]=2+(rand()>>10)%(LY-4);
	puti(pbuf,i);
    }


    while(1)
    {
	beta=bmin;
	for(s=0;s<cnt;s++)
	{
	    init_th(beta);

	    for(i=0;i<nn;i++)
	    {
		move1(i);
		move2(i);
	    }
	    for(i=0;i<nn;i++)
	    {
		move1(i);
		move2(i);
	    }

	    beta *= db;

	    int e1=ene1();
	    int e2= ene2(pbuf);
	    ee=e1+e2;
	    if(s%1000==0)
	    {
		printf("%e %d %d %d\n", beta, ee, e1, e2);
	        if(ee==1)dump();
		fflush(stdout);
	    }
	    if(ee==0)break;
	}
	exit(1);

	dump();
	if(ee==0) break;
    }

    return 0;
}