最適ゴロム定規

総当り法で最適ゴロム定規を探す

最適ゴロム定規(OGR:Optimal Golomb Ruler)とは、ゴロム定規(GR:Golomb Ruler)の内最も短いものの事です。ゴロム定規とはn個の数字からなる数列GR-nにおいて、任意の要素a,b,c,d(ただしa>b、c>d、a≠b、c≠d)がa-b≠c-dを満たすものを言います。nの事をマーク数と呼ぶこともあります。

GR-nは全ての要素を自然数倍しても関係が壊れないため、無限に存在します。その内、数列の要素の最大値が最も小さいものをOGR-nと呼びます。GRは反転しても関係が壊れないため、OGR-nは常に偶数個存在することになります。

#include <stdio.h>
#include <stdlib.h>

int check_ruler( int mark_count, int *ruler );
void print_ruler( int mark_count, int *ruler );


int main( int ac, char *av[] )
{
    int mark, mark_count, *ruler, length, found, nodes;

    /* コマンドラインからマーク数を決定する */
    if( ac < 2 )
        return( 1 );
    mark_count = strtol( av[1], NULL, 10 );
    if( mark_count < 2 )
        mark_count = 2;

    /* ruler用配列を確保する */
    ruler = calloc( mark_count, sizeof(int) );
    if( ruler == NULL )
        return( 1 );

    /* ループ変数を初期化する */
    length = ( mark_count - 1 ) * mark_count / 2;
    found = 0;
    nodes = 0;

    /* OGR探索ループ */
    while( found == 0 ) {
        /* rulerを0,1,2...,lengthで初期化する */
        for( mark = 1; mark < mark_count - 1; mark++ )
            ruler[mark] = mark;
        ruler[mark_count - 1] = length;

        /* 長さがlengthのGRを探す */
        while( 1 ) {
            /* 計算回数をカウント */
            nodes++;

            /* golomb rulerであるか */
            if( check_ruler( mark_count, ruler ) == 1 ) {
                print_ruler( mark_count, ruler );
                found = 1;
            }

            /* markを一つ先に進める */
            for( mark = mark_count - 2; mark > 0; mark-- ) {
                if( ruler[mark] < length - ( mark_count - mark - 1 ) ) {
                    ruler[mark]++;
                    for( mark = mark + 1; mark < mark_count - 1; mark++ )
                        ruler[mark] = ruler[mark - 1] + 1;
                    break;
                }
            }
            /* 進めるマークが無ければ終了 */
            if( mark == 0 )
                break;
        }
        /* lengthを増やす */
        length++;
    }

    /* 計算回数を表示 */
    printf( "%d nodes\n", nodes );

    return( 0 );
}


/* golomb ruler判定 */
int check_ruler( int mark_count, int *ruler )
{
    int n1, n2, n3, n4;
    int d1, d2;

    for( n1 = 0; n1 < mark_count - 1; n1++ ) {
        for( n2 = n1 + 1; n2 < mark_count; n2++ ) {
            d1 = ruler[n2] - ruler[n1];
            for( n3 = n1 + 1; n3 < mark_count - 1; n3++ ) {
                for( n4 = n3 + 1; n4 < mark_count; n4++ ) {
                    d2 = ruler[n4] - ruler[n3];
                    /* 差の等しい要素の組があれば偽 */
                    if( d1 == d2 )
                        return( 0 );
                }
            }
        }
    }

    return( 1 );
}


/* ruler表示 */
void print_ruler( int mark_count, int *ruler )
{
    int mark;

    for( mark = 0; mark < mark_count; mark++ ) {
        printf( "%d", ruler[mark] );
        if( mark != mark_count - 1 )
            printf( "-" );
    }
    printf( "\n" );
}

定義に従い定規を短いものから順に長くしていきながら全ての数列をゴロム定規となるか試しています。当然ながらマーク数が多くなるほど計算に時間がかかり、このプログラムではOGR-11までしか求められていません。その結果は以下のとおりです。

OGR-2 0-1
OGR-3 0-1-3
0-2-3
OGR-4 0-1-4-6
0-2-5-6
OGR-5 0-1-4-9-11
0-2-7-8-11
0-2-7-10-11
0-3-4-9-11
OGR-6 0-1-4-10-12-17
0-1-4-10-15-17
0-1-8-11-13-17
0-1-8-12-14-17
0-2-7-13-16-17
0-3-5-9-16-17
0-4-6-9-16-17
0-5-7-13-16-17
OGR-7 0-1-4-10-18-23-25
0-1-7-11-20-23-25
0-1-11-16-19-23-25
0-2-3-10-16-21-25
0-2-5-14-18-24-25
0-2-6-9-14-24-25
0-2-7-13-21-22-25
0-2-7-15-21-24-25
0-3-4-12-18-23-25
0-4-9-15-22-23-25
OGR-8 0-1-4-9-15-22-32-34
0-2-12-19-25-30-33-34
OGR-9 0-1-5-12-25-27-35-41-44
0-3-9-17-19-32-39-43-44
OGR-10 0-1-6-10-23-26-34-41-53-55
0-2-14-21-29-32-45-49-54-55
OGR-11 0-1-4-13-28-33-47-54-64-70-72
0-1-9-19-24-31-52-56-58-69-72
0-2-8-18-25-39-44-59-68-71-72
0-3-14-16-20-41-48-53-63-71-72

世界では、OGR-23まで既に求められているそうです。また、OGRでは無いかもしれないが、ほぼOGRと思われる近似解はGR-150まで求められているそうです。その一覧はhttp://www.research.ibm.com/people/s/shearer/grtab.htmlにあります。

(2012/2/12追記)

distributed.netによりOGR-26まで求まったそうです。現在、OGR-27の探索が行われています。

(2012/2/21追記)

check_ruler()で余分な比較を行っていたのを修正しました。Athlon 64 X2 5600+(2.9GHz)でのOGR-10の探索が約30分から約25分に短縮されました。

(2012/2/22追記)

これまではOGR-nを求めるのに、まず長さnから徐々に定規を伸ばしていく方法をとっていました。この方法自体に問題はないのですが、スタートの長さがnが短すぎるため、余分な計算をしていたと思われます。なぜなら、OGR-nはOGR-(n-1)よりも必ず長くなるからです。

なぜOGR-nがOGR-(n-1)より長くなるかというと、OGR-nがOGR-(n-1)より短いと仮定します。ゴロム定規は末尾の数字を取去ってもゴロム定規ですのでOGR-nから末尾の数字を取去ってもゴロム定規です。OGR-nはOGR-(n-1)より短いと仮定していますから、OGR-nから末尾の数字を取去ったゴロム定規もOGR-(n-1)よりも短いです。そしてこれもゴロム定規ですから、OGR-(n-1)より短いゴロム定規が存在することになりOGR-(n-1)はOGRではないということになります。これは矛盾しますから、仮定が間違っていたということで、OGR-nは常にOGR-(n-1)より長いということが保証されます。

と思うんですが、ほんとにこれで証明があっているかどうかと問われると自信はないです。数学素人なもんで。おそらく正しいでしょうけれど、上記の証明では不十分かもしれません。ともあれ、OGR-nはOGR-(n-1)より常に長いとするならば、OGR-nを求めるときのスタートの定規の長さはOGR-(n-1)を用いることで計算をだいぶ省力化できるものと思われます。

ということで、長さの初期値をコマンドラインから与えられるように修正してみたのですが、効果は意外とたいしたことはなくて、OGR-10の探索で約1分短縮されただけでした。

(2012/3/9追記)

ゴロム定規にはなり得ない長さから始めて、徐々に長さを長くしていくという方法をとっていますが、この場合始める長さが出来るだけ長い方が計算量が少なくなります。これまではOGR-nを求めるのに長さの初期値をn(0-1-2-3-4...nという定規が初期値)としていましたが、これはあまり効率がよくありません。もう少し効率が良い方法として、上述したようにOGR-(n-1)を用いるという方法もありますが、これはマジックナンバーを要求するということになり、仕様上あまり美しくありません。美しくなくても問題が速く解けるならそれは別に問題はないのですが。他に何かいい手はないかと考えたところ、ゴロム定規の定義上、数列の差が1,2,3,4,..よりも短いゴロム定規は存在しないことが思いつきました。厳密にはこれにも証明が必要なのですが、そこはいつもの言い訳(数学素人ですから)を使わせてもらって、これを用いることにしました。効果のほどはOGR-(n-1)を用いるのと大して変りはないのですが、マジックナンバーを用いなくてよいという点で優れていると思います。

(2012/3/10追記)

check_ruler()を手作業で最適化したところ、数パーセントですが高速化されたので反映しました。コンパイラによる最適化に期待して特に手を入れていなかったのですが、意外と効果があるようです。もちろん、使用しているコンパイラやターゲットCPUによって大きく事情が変わると思います。今回、たまたまFreeBSD-9.0 RELEASE(i386) + gcc4.2.1 + Athlon 64 X2 5600+(2.9GHz)において有意な効果があったということです。ちなみにこの最適化でOGR-10が約22分で完了するようになりました。

(2015/12/17追記)

distributed.netにより2014年2月19日にOGR-27が求まったことを書き忘れていました。現在、OGR-28の探索が行われています。

バックトラック法でOGRを探す

(2012/2/24追記)

総当り法では効率が良くないので、バックトラック法というアルゴリズムを用いてプログラムを書き直してみました。

#include <stdio.h>
#include <stdlib.h>

int set_ruler( int pos, int length, int mark_count, int *ruler );
int check_ruler( int mark_count, int *ruler );
void print_ruler( int mark_count, int *ruler );

int nodes;
int ogrlength[] = { 0, 1, 3, 6, 11, 17, 25, 34, 44, 55, 72, 85, 106, 127, 151, 177, 199, 216, 246, 283, 333, 356, 372, 425, 480, 492 };

int main( int ac, char *av[] )
{
    int mark_count, *ruler, length;
    int i;

    /* コマンドラインからマーク数を決定する */
    if( ac < 2 )
        return( 1 );
    mark_count = strtol( av[1], NULL, 10 );
    if( mark_count < 2 )
        mark_count = 2;

    /* ruler用配列を確保する */
    ruler = calloc( mark_count, sizeof(int) );
    if( ruler == NULL )
        return( 1 );

    /* 変数を初期化する */
    length = ogrlength[mark_count - 1];
    nodes = 0;

    setbuf( stdout, NULL );
    while( 1 ) {
        /* ruler初期化 */
        for( i = 0; i < mark_count; i++ )
            ruler[i] = 0;

        /* ruler生成 */
        if( set_ruler( 1, length, mark_count, ruler ) )
            break;

        /* 長さを1伸ばす */
        length++;
    }

    /* 計算回数を表示 */
    printf( "%d nodes\n", nodes );

    return( 0 );
}


/* ruler生成 */
int set_ruler( int pos, int length, int mark_count, int *ruler )
{
    int i, length2, ret;

    if( pos <= 2 || check_ruler( pos, ruler ) ) {
        /* 指定markのgolomb rulerならOGRである */
        if( pos == mark_count ) {
            print_ruler( mark_count, ruler );
            return( 1 );
        }
    }
    else {
        return( 0 );
    }

    /* 要素を一つ増やして再帰呼び出し */
    ret = 0;
    length2 = length - ogrlength[mark_count - pos - 1];
    for( i = ruler[pos - 1] + 1; i <= length2; i++ ) {
        ruler[pos] = i;
        if( set_ruler( pos + 1, length, mark_count, ruler ) )
            ret = 1;
    }
    return( ret );
}


/* golomb ruler判定 */
int check_ruler( int mark_count, int *ruler )
{
    int n1, n2, n3, n4;
    int d1, d2;

    nodes++;

    for( n1 = 0; n1 < mark_count - 1; n1++ ) {
        for( n2 = n1 + 1; n2 < mark_count; n2++ ) {
            d1 = ruler[n2] - ruler[n1];
            for( n3 = n1 + 1; n3 < mark_count - 1; n3++ ) {
                for( n4 = n3 + 1; n4 < mark_count; n4++ ) {
                    d2 = ruler[n4] - ruler[n3];
                    /* 差の等しい要素の組があれば偽 */
                    if( d1 == d2 )
                        return( 0 );
                }
            }
        }
    }

    return( 1 );
}


/* ruler表示 */
void print_ruler( int mark_count, int *ruler )
{
    int mark;

    for( mark = 0; mark < mark_count; mark++ ) {
        printf( "%d", ruler[mark] );
        if( mark != mark_count - 1 )
            printf( "-" );
    }
    printf( "\n" );
}

このプログラムは大変効率がよく、Athlon 64 X2 5600+(2.9GHz)でOGR-10を求めるのに90秒しか掛かりませんでした。check_ruler()を行った回数も、総当り法では1,499,512,759回だったのがバックトラック法では130,019,067回と大幅に減少しています。同じくAthlon 64 X2 5600+(2.9GHz)でOGR-11を求めたところ、約52分で終了しました。

(2012/2/26追記)

プログラムが高速化されたのでOGR-12を求めてみたところ、Athlon 64 X2 5600+(2.9GHz)で約10時間で求まりました。

OGR-12 0-2-6-24-29-40-43-55-68-75-76-85
0-9-10-17-30-42-45-56-61-79-83-85

(2012/3/9追記)

総当り法と同様に、lengthの初期値を(マーク数-1)*(マーク数)/2に設定しました。また、バックトラックでの後戻りの条件判定を少し変更しました。この変更によりAthlon 64 X2 5600+(2.9GHz)でOGR-10を求めるのに24秒と大幅に短縮されました。check_ruler()を行った回数も34,543,691回と更に減少しました。OGR-11を求めるのには約15分、OGR-12では135分となりました。

(2012/3/10追記)

総当り法と同様に、check_ruler()を修正しました。ただ、こちらではcheck_ruler()を行う回数が少ないため、全体としての効果は総当り法ほどではありませんでした。

自明なcheck_ruler()を行わないように修正しました。具体的にはマーク数2の場合はどんな数列でもゴロム定規になりますのでcheck_ruler()を行う必要はありません。この修正により、ほんの少しだけ高速化しました。

(2012/3/11追記)

上のほうで「マジックナンバーは使いたくない」と書いておきながらなんですが、やっぱり処理速度を上げるには使った方が有利なので、既知のOGRの長さをプログラム上に固定で持つように修正しました。既知の長さを全て持っていますが、OGR-nを求めるときにはOGR-(n--1)の長さまでしか使用しませんので、短いものから順番に求めていっていれば、問題はないと思います。この修正の効果は顕著で、OGR-10がわずか5秒で求まります。check_ruler()の回数も6,355,224回と更に激減しました。OGR-11は2分半、OGR-12は18分で求まりました。

(2012/3/12追記)

3度目のチャレンジでようやくOGR-13が求まりました。計算時間は約10時間で、OGR-12から大幅に増加しました。OGR-14も求めたいと思っていましたが、この様子では困難かもしれません。

OGR-13 0-2-5-25-37-43-59-70-85-89-98-99-106
0-7-8-17-21-36-47-63-69-81-101-104-106

(2012/3/18追記)

これまで動作環境としてFreeBSD(i386)を使用していたのですが、よく考えると64bit CPUなのですからFreeBSD(amd64)が使用できます。ということで環境を変えてみたところ、OGR-13が約5時間と倍近い速度になりました。これなら可能かもとOGR-14も計算してみたところ、約92時間で求まりました。

OGR-14 0-4-6-20-35-52-59-77-78-86-89-99-122-127
0-5-28-38-41-49-50-68-75-92-107-121-123-127

あおやぎのさいと2.0初歩の整数論プログラミング