123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130 |
- #!/usr/bin/perl
- use 5.020;
- use strict;
- use warnings;
- use experimental qw(signatures);
- use Math::GMPz;
- use ntheory qw(is_prime kronecker);
- sub is_fermat_fibonacci_prime($n) {
- if (ref($n) eq 'Math::GMPz') {
- $n = Math::GMPz::Rmpz_init_set($n);
- }
- else {
- $n = Math::GMPz->new("$n");
- }
- if (Math::GMPz::Rmpz_even_p($n)) {
- return ((Math::GMPz::Rmpz_cmp_ui($n, 2) == 0) ? 1 : 0);
- }
- if (Math::GMPz::Rmpz_cmp_ui($n, 5) == 0) {
- return 1;
- }
- my $t = Math::GMPz::Rmpz_init_set_ui(2);
- my $u = Math::GMPz::Rmpz_init();
- Math::GMPz::Rmpz_sub_ui($u, $n, 1);
- Math::GMPz::Rmpz_powm($u, $t, $u, $n);
- Math::GMPz::Rmpz_cmp_ui($u, 1) == 0 or return 0;
- my $v = Math::GMPz::Rmpz_init();
- my $w = Math::GMPz::Rmpz_init();
- Math::GMPz::Rmpz_sub_ui($w, $n, 1);
- my $m = Math::GMPz::Rmpz_init_set($n);
- my $f = Math::GMPz::Rmpz_init_set_ui(0);
- my $g = Math::GMPz::Rmpz_init_set_ui(1);
- my $a = Math::GMPz::Rmpz_init_set_ui(0);
- my $b = Math::GMPz::Rmpz_init_set_ui(1);
- my $c = Math::GMPz::Rmpz_init_set_ui(1);
- my $d = Math::GMPz::Rmpz_init_set_ui(1);
- for (; ;) {
- if (Math::GMPz::Rmpz_odd_p($n)) {
- Math::GMPz::Rmpz_set($t, $f);
- Math::GMPz::Rmpz_mul($f, $f, $a);
- Math::GMPz::Rmpz_addmul($f, $g, $c);
- Math::GMPz::Rmpz_mod($f, $f, $m);
- Math::GMPz::Rmpz_mul($g, $g, $d);
- Math::GMPz::Rmpz_addmul($g, $t, $b);
- Math::GMPz::Rmpz_mod($g, $g, $m);
- }
- Math::GMPz::Rmpz_div_2exp($n, $n, 1);
- Math::GMPz::Rmpz_sgn($n) || last;
- Math::GMPz::Rmpz_set($t, $a);
- Math::GMPz::Rmpz_mul($a, $a, $a);
- Math::GMPz::Rmpz_addmul($a, $b, $c);
- Math::GMPz::Rmpz_mod($a, $a, $m);
- Math::GMPz::Rmpz_set($u, $b);
- Math::GMPz::Rmpz_mul($b, $b, $d);
- Math::GMPz::Rmpz_addmul($b, $u, $t);
- Math::GMPz::Rmpz_mod($b, $b, $m);
- Math::GMPz::Rmpz_set($v, $c);
- Math::GMPz::Rmpz_mul($c, $c, $t);
- Math::GMPz::Rmpz_addmul($c, $v, $d);
- Math::GMPz::Rmpz_mod($c, $c, $m);
- Math::GMPz::Rmpz_mul($d, $d, $d);
- Math::GMPz::Rmpz_addmul($d, $v, $u);
- Math::GMPz::Rmpz_mod($d, $d, $m);
- }
- if ( Math::GMPz::Rmpz_cmp_ui($f, 1) == 0
- or Math::GMPz::Rmpz_cmp($f, $w) == 0) {
- if (not is_prime($w + 1)) {
- say "error for ", $w + 1, " -> with $f";
- }
- return 1;
- }
- return 0;
- }
- my %seen;
- my @nums;
- while (<>) {
- next if /^\h*#/;
- /\S/ or next;
- my $n = (split(' ', $_))[-1];
- $n || next;
- next if $seen{$n}++;
- #say $. if $. % 10000 == 0;
- #if ($n >= (~0 >> 1)) {
- $n = Math::GMPz->new("$n");
- #}
- push @nums, $n;
- #ntheory::is_provable_prime($n) && die "error: $n\n";
- #ntheory::is_prime($n) && die "error: $n\n";
- #ntheory::is_aks_prime($n) && die "error: $n\n";
- #ntheory::miller_rabin_random($n, 7) && die "error: $n\n";
- }
- @nums = sort { $a <=> $b } @nums;
- foreach my $n (@nums) {
- is_fermat_fibonacci_prime($n) && do {
- warn "error: $n\n";
- last;
- };
- }
|