倍積完全数

倍積完全数

約数和配列を利用して、約数の和が元の数のn倍になるような数を探してみます。約数の和がその数自身より小さいものを不足数、大きいものを過剰数と呼ぶそうです。という事は全ての自然数は不足数、完全数、過剰数のどれかに分類されることになります。倍積完全数は過剰数の特別なものになります。


#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, q, r, 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 );
    mpz_init( q );
    mpz_init( r );
    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 );
            }
        }

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

    return( 0 );
}

1億以下の数について10倍完全数までを探索したところ、3倍完全数として120、672、523776が、4倍完全数として30240、32760、2178540、23569920、45532800が見つかりました。5倍以上の完全数は見つかっていません。直感的にも倍数の大きな完全数は存在しにくいだろうと思われますので、この結果にも納得がいきます。

(2012/2/16追記)

以前ここでは「n倍完全数」と書いていましたが「倍積完全数」と呼ぶのが正しいようです。また、n倍完全数とはn自身も加えて計算するそうなので、一般の完全数は2倍完全数となり、これまでこのページで呼んでいたn倍完全数はn+1倍完全数と呼ぶのが正しいようです。これらの点を修正しました。

また、メモリが潤沢に使用できるようになったので配列のサイズを1億に増やしました。探索範囲を10億以下に広げたところ、3倍完全数として459818240が、4倍完全数として142990848が見つかりました。

(2020/4/9追記)

GMPを用いてより広い範囲で探索できる版を作成し、100億以下を探索したところ3倍完全数として1476304896が、4倍完全数として1379454720が見つかりました。

(2020/4/25追記)

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