Find (prime) factors.

choco1.pl

#!/usr/bin/perl

use strict;
use warnings;
use utf8;

use Math::BigInt lib=>'GMP';
use Math::Primality qw(is_prime);

my $n = Math::BigInt->new('280671392065546467397265294532969672241810318954163887187279320454220348884327');
my $original_n = $n->copy();
#my $n = Math::BigInt->new('481362815814826159');
#my $n = Math::BigInt->new('225');

# is_prime で返ってくる値から文字列に

my @factor_kind = ( 'composite', '(possibly) prime', 'prime' );

# 疑似乱数列を作成する関数
# f(x) = x^2+1 mod n
sub f {
    my $x = shift;

    return $x->copy()->bpow( 2 )->binc()->bmod( $n );
}

# 先に n が素数かどうかを確認する
if ( my $is_prime = is_prime( $n ) ) {
    print "$n is $factor_kind[$is_prime].\n";
    exit 0;
}

# 得られた因数のリスト

my @factors = ();

# Pollard's rho algorithm
while (1) {
    print "n = $n\n";

    my $x = Math::BigInt->new('2');
    my $y = Math::BigInt->new('2');
    my $d = Math::BigInt->new('1');
    while ( $d == 1 ) {
        $x = f($x);#print "$x\n";
        $y = f(f($y));#print "$y\n";
        $d = Math::BigInt::bgcd( (Math::BigInt->babs( $x - $y ), $n) );
        if ( $d == $n ) {
            last;
        }
    }
    if ( $d == $n ) {
        # 疑似乱数列が循環したら、残りが素数かどうかを確認しつつ終了する

        my $is_prime = is_prime( $n );
        print "rest ($factor_kind[$is_prime]) $n\n";
        push @factors, $n;
        last;
    }

    # 因数が見つかった

    push @factors, $d;
    my $is_prime = is_prime( $d );
    print "factor ($factor_kind[$is_prime]) $d\n";

    $n->bdiv( $d );

    # 因数で割った残りが素数(または疑似素数)であれば、終了

    if ( $is_prime = is_prime( $n ) ) {
        print "factor ($factor_kind[$is_prime]) $n\n";
        push @factors, $n;
        last;
    }
}

# 結果を表示
# 結果は、表面積が最小のパターンを探すのに利用する

print "$original_n = \n";
print join ( 'x', @factors ), "\n";

exit 0;

1;