完全数を求める

完全数

完全数というのは、自分自身以外の約数の総和が自分自身に等しくなる自然数の事です。例えば6の約数は1, 2, 3でその和は6となりますので完全数です。まずは、この定義通りにプログラムを作り完全数を調べだしてみます。

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

int main( int ac, char *av[] )
{
    int n, min_n, max_n, sum, i;

    /*  コマンドラインから完全数探索範囲を決定する    */
    if( ac < 3 ) {
        printf( "usage : perfect min_n max_n\n" );
        return( 1 );
    }
    min_n = strtol( av[1], NULL, 10 );
    max_n = strtol( av[2], NULL, 10 );

    /*  探索範囲の数を調べる    */
    for( n = min_n; n <= max_n; n++ ) {
        /*  約数の総和を求める  */
        sum = 0;
        for( i = 1; i < n; i++ ) {
            if( n % i == 0 )
                sum += i;
        }

        /*  約数の総和が自身に等しければ完全数である    */
        if( sum == n )
            printf( "%d\n", n );
    }

    return( 0 );
}

まず10000以下の範囲で調べてみた所、見つかった完全数は6、28、496、8128でした。このプログラムは素数の場合と同様、広い範囲の数を探そうとすると膨大な時間が掛かってしまうという欠点があります。

(2012/2/12追記)

Athlon 64 X2で実行したところ、10万以下の値について探索するのに78秒掛かりました。

(2020/4/25追記)

ソースコードを少し整理しました。

約数和配列

完全数を単純な方法で探し出そうとすると、とてつもなく計算コストが掛かります。そこで、効率の良い方法として、エラトステネスのふるいをヒントに次のようなものを考えました。

まず配列を用意します。次に2以上の整数の倍数を求め、その値の配列に元の整数を加えます。整数を配列要素数まで繰り返せば、各配列には添え字の約数の総和が入っていることになります。添え字と約数の和が等しければ、その値は完全数となります。

素数を求める時と同様、この配列を繰り返し使うことによって、大きな範囲を探せるようにしました。

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

#define ARRAY_SIZE_MAX   (100000000)

mpz_t array[ARRAY_SIZE_MAX];

int main( int ac, char *av[] )
{
    int arraysize, i, j;
    mpz_t base, min_n, max_n, to_n, tmp, mpz_i, mpz_j, to_i, to_j;

    /*  コマンドラインから完全数探索範囲を決定する    */
    if( ac < 3 ) {
        printf( "usage : gmpperfectarray min_n max_n\n" );
        return( 1 );
    }
    mpz_init_set_str( min_n, av[1], 10 );
    mpz_init_set_str( max_n, av[2], 10 );

    mpz_init( tmp );
    mpz_init( mpz_i );
    mpz_init( mpz_j );
    mpz_init( to_i );
    mpz_init( to_j );
    mpz_init( to_n );
    for( i = 0; i < ARRAY_SIZE_MAX; i++ ) {
        mpz_init( array[i] );
    }

    /*  ARRAY_SIZE_MAX分ごとの整数区間を調べる   */
    for( mpz_init_set( base, min_n ); mpz_cmp( base, max_n ) < 0; mpz_add_ui( base, base, ARRAY_SIZE_MAX ) ) {
        /*  整数区間配列の大きさを決める    */
        mpz_sub( tmp, max_n, base );
        mpz_add_ui( tmp, tmp, 1 );
        if( mpz_cmp_ui( tmp, ARRAY_SIZE_MAX ) > 0 ) {
            arraysize = ARRAY_SIZE_MAX;
        }
        else {
            arraysize = mpz_get_ui( tmp );
        }
        /*  配列を初期化する    */
        for( i = 0; i < arraysize; i++ ) {
            mpz_set_ui( array[i], 0 );
        }

        /*  約数の和を求める    */
        mpz_add_ui( to_i, base, arraysize );
        mpz_tdiv_q_ui( to_i, to_i, 2 );
        for( mpz_set_ui( mpz_i, 1 ); mpz_cmp( mpz_i, to_i ) < 0; mpz_add_ui( mpz_i, mpz_i, 1 ) ) {
            if( mpz_cmp( base, mpz_i ) <= 0 ) {
                mpz_mul_ui( mpz_j, mpz_i, 2 );
            }
            else if( mpz_mod( tmp, base, mpz_i ), mpz_cmp_ui( tmp, 0 ) == 0 ) {
                mpz_set( mpz_j, base );
            }
            else {
                mpz_tdiv_q( mpz_j, base, mpz_i );
                mpz_add_ui( mpz_j, mpz_j, 1 );
                mpz_mul( mpz_j, mpz_j, mpz_i );
            }
            mpz_add_ui( to_j, base, arraysize );
            for( ; mpz_cmp( mpz_j, to_j ) < 0; mpz_add( mpz_j, mpz_j, mpz_i ) ) {
                mpz_sub( tmp, mpz_j, base );
                j = mpz_get_ui( tmp );
                mpz_add( array[j], array[j], mpz_i );
            }
        }

        /*  約数の和とその数自身が等しければ完全数である    */
        for( i = 0; i < arraysize; i++ ) {
            mpz_add_ui( tmp, base, i );
            if( mpz_cmp( array[i], tmp ) == 0 )
                gmp_printf( "%Zd\n", tmp );
        }
    }

    return( 0 );
}

このプログラムを使用して10億以下の値について探索した結果、更に33550336が完全数であることがさらに判明しました。ここまでの計算時間は約1日ですから、単純な方法に比べて随分と計算効率が良いことが分かります。

(2012/2/12追記)

Athlon 64 X2 5600+(2.9GHz)で実行したところ、10万以下の値について探索するのに0.01秒、10億以下の探索は1時間半ほどで終了しました。やっぱり今時のCPUは速いですね。

(2012/2/16追記)

メモリが潤沢に使用できるようになったので配列のサイズを1億に増やしました。これにより10億以下の値についての探索は約15分で終了するようになりました。

(2020/4/9追記)

GMPを使用するように書き換えてより広い範囲を探索できるようにした版を用いて、100億以下の範囲で探索したところ、8589869056が見つかりました。

(2020/4/25追記)

ソースコードを少し整理しました。

完全数の一覧

これまでにこのページで見つけた完全数の一覧を載せておきます。

世の中では現時点で51個の完全数が見つかっているそうですが、このページでは自力で求めたものだけを掲載することにしていますので、今のところ下記のみとなります。