素数生成式を検証する

必ず結果が素数となる公式を見つけることは多くの数学者の夢だそうですが、そのような公式は未だに見つかっていないそうです。とはいえ、これまでいくつかかなりの確率で素数を生成する公式が見つかっていますので、どれくらいの精度で見つけられるのか、検証します。

検証してみるのは次の四つの式です。

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

int isprime( mpz_t n );

int main( int ac, char *av[] )
{
    int   a, b, c;
    mpz_t n, n2, max_n, tmp;

    /*  コマンドラインから係数と探索範囲を決定する    */
    if( ac < 5 ) {
        printf( "usage : primegen a b c max_n\n" );
        return( 1 );
    }
    a = strtol( av[1], NULL, 10 );
    b = strtol( av[2], NULL, 10 );
    c = strtol( av[3], NULL, 10 );
    mpz_init_set_str( max_n, av[4], 10 );

    /*  探索範囲の数を調べる    */
    mpz_init_set_ui( n, 0 );
    mpz_init( n2 );
    mpz_init( tmp );
    while( mpz_cmp( n, max_n ) <= 0 ) {
        /*  a * n^2 + b * n + c を求める   */
        mpz_mul( n2, n, n );
        mpz_mul_ui( n2, n2, a );
        mpz_mul_ui( tmp, n, b );
        mpz_add( n2, n2, tmp );
        mpz_add_ui( n2, n2, c );

        /*  素数かどうか調べる  */
        if( isprime( n2 ) )
            gmp_printf( "%Zd %Zd OK\n", n, n2 );
        else
            gmp_printf( "%Zd %Zd NG\n", n, n2 );

        mpz_add_ui( n, n, 1 );
    }

    mpz_clear( n );
    mpz_clear( n2 );
    mpz_clear( max_n );
    mpz_clear( tmp );

    return( 0 );
}

int isprime( mpz_t n )
{
    mpz_t i, n2;
    int d, ret;

    /*  1は素数ではない */
    if( mpz_cmp_ui( n, 1 ) == 0 )
        return( 0 );

    /*  2,3は素数 */
    if( mpz_cmp_ui( n, 2 ) == 0 || mpz_cmp_ui( n, 3 ) == 0 )
        return( 1 );

    /*  2,3で割り切れたら合成数    */
    if( mpz_even_p( n ) || mpz_divisible_ui_p( n, 3 ) )
        return( 0 );

    /*  sqrt(n)を求める */
    mpz_init( n2 );
    mpz_sqrt( n2, n );

    /*  n2以下の2,3の倍数以外での剰余が0かどうか調べる   */
    d = 2;
    mpz_init_set_ui( i, 5 );
    ret = 1;
    while( mpz_cmp( i, n2 ) <= 0 ) {
        if( mpz_divisible_p( n, i ) ) {
            ret = 0;
            break;
        }
        mpz_add_ui( i, i, d );
        d = ( d == 2 ? 4 : 2 );
    }

    mpz_clear( i );
    mpz_clear( n2 );

    return( ret );
}

(2012/2/15追記)

大きな値にも対応可能なようにGMPで書き直しました。

(2012/2/19追記)

ほとんど同じプログラムを4個書くのは無駄であるということに15年目にしてようやく気付きましたので、係数をコマンドラインから与えられるように修正しました。合わせて素数と同様に、3の倍数による除算を除外するようにして素数判定の処理を高速化しました。

(2020/4/25追記)

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

結果検証

実際に100万以下のnに対して調べてみた所、以下の結果になりました。

nに対する素数の存在
2n+1 n2+n+41 4n2+170n+1847 4n2+4n+59
10以下 7個 (63.6%) 11個 (100.0%) 9個 (81.8%) 11個 (100.0%)
100以下 45個 (44.5%) 87個 (86.1%) 71個 (70.2%) 72個 (71.2%)
1,000以下 302個 (30.1%) 582個 (58.1%) 509個 (50.8%) 460個 (45.9%)
1万以下 2,261個 (22.6%) 4,149個 (41.4%) 3,822個 (38.2%) 3,408個 (34.0%)
10万以下 17,983個 (17.9%) 31,985個 (31.9%) 29,906個 (29.9%) 26,656個 (26.6%)
100万以下 148,932個 (14.8%) 261,081個 (26.1%) 247,425個 (24.7%) 220,611個 (22.0%)

 

式の値に対する素数の存在
2n+1 n2+n+41 4n2+170n+1847 4n2+4n+59
10以下 3個 (60.0%) - - -
100以下 24個 (48.0%) 8個 (100.0%) - 3個 (100.0%)
1,000以下 167個 (33.4%) 31個 (100.0%) - 14個 (93.3%)
1万以下 1,228個 (24.5%) 86個 (86.0%) 23個 (79.3%) 41個 (82.0%)
10万以下 9,591個 (19.1%) 221個 (69.9%) 90個 (65.6%) 102個 (64.5%)
100万以下 78,497個 (15.6%) 581個 (58.1%) 270個 (56.3%) 259個 (51.8%)
1000万以下 - 1,503個 (47.5%) 727個 (46.6%) 691個 (43.7%)
1億以下 - 4,149個 (41.4%) 2,064個 (41.4%) 1,843個 (36.8%)
10億以下 - 11,355個 (35.9%) 5,706個 (36.1%) 5,073個 (32.0%)
100億以下 - 31,985個 (31.9%) 15,971個 (31.9%) 14,288個 (28.5%)
1000億以下 - 90,940個 (28.7%) 45,418個 (28.7%) 40,478個 (25.6%)
1兆以下 - 261,081個 (26.1%) 130,423個 (26.0%) 116,545個 (23.3%)

まず奇数の場合の結果ですが、これはほぼ素数の存在確率そのものを表していることになります。nの値の結果から見ると好成績のようですが、これは式の値自体が他の式に比べると小さいですから直接比較する事は出来ません。素数は小さな数の範囲では存在確率が大きくなることが知られていますから。

他の三つの式はなかなか高い確率で素数を生成していることが分かります。奇数での素数の存在確率が15%程度になってしまう100万以下の範囲でも50%以上の結果で素数を生成しています。また、式の値が1兆以下でも25%程度の結果が素数となっていますから、これは驚異的と言ってもいいでしょう。

ちなみにそれぞれの式に当てはめるnの内、最初から連続してどこまで素数を生成し続けたかを比べてみます。まず2n+1の場合は式の値が0でいきなり素数ではありません。しかしこれは例外として除くと、n=3の場合の7までとなります。次にオイラーのn2+n+41では、n=39の場合の1601までが素数となります。4n2+170n+1847では、n=1でいきなり合成数となってしまうため、n=0の場合の1847までとなります。4n2+4n+59では、n=13の場合の787までが連続して素数となります。

こうしてみると、以上の4式の中ではオイラーの見つけたn2+n+41が僅差で最も優秀な素数生成式と言えるのではないでしょうか。