123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150 |
- #!/usr/bin/perl
- # Daniel "Trizen" Șuteu
- # Date: 27 June 2019
- # https://github.com/trizen
- # Find highly composite numbers H(n) such that H(n)-1 and H(n)+1 are both primes.
- # In other words, find highly composite numbers sandwiched between twin primes.
- # The sequence of indices of highly composite numbers A002182 which are between a twin prime pair, begins as:
- # 3, 4, 5, 9, 11, 12, 20, 28, 30, 84, 108, 118, 143, 149, 208, 330, 362, 1002, 2395, 3160, 10535
- # Example:
- # H(10535) = A108951(52900585920)
- # = 14742217487368791965347653720647452690286549052234444179664342042930370966727413549068727214664401976854238590421417268673037399536054005777393104248210539172848500736334237168727231561710827753972114334247396552090671649834020135652920430241738510495400044737265204738821393451152066370913670083496651044937158497896720493198891148968218874744806522767468280764179516341996273430700779982929787918221844760577694188288275419541410142336911631623319041967633591283303769044016192030492715535641753600000
- # where H(10535)-1 and H(10535)+1 are both prime numbers.
- # See also:
- # https://oeis.org/A321995 -- Indices of highly composite numbers A002182 which are between a twin prime pair.
- # https://oeis.org/A108951 -- Completely multiplicative with a(p) = p# for prime p, where x# is the primorial A034386(x).
- # https://oeis.org/A002182 -- Highly composite numbers, definition (1): where d(n), the number of divisors of n (A000005), increases to a record.
- use 5.020;
- use strict;
- use warnings;
- use Math::GMPz;
- use POSIX qw(ULONG_MAX);
- use ntheory qw(:all);
- use IO::Uncompress::Bunzip2;
- use experimental qw(signatures);
- local $| = 1;
- prime_precalc(1e7);
- sub primality_pretest ($n) {
- # Must be positive
- (Math::GMPz::Rmpz_sgn($n) > 0) || return;
- # Check for divisibilty by 2
- if (Math::GMPz::Rmpz_even_p($n)) {
- return (Math::GMPz::Rmpz_cmp_ui($n, 2) == 0);
- }
- # Return early if n is too small
- Math::GMPz::Rmpz_cmp_ui($n, 101) > 0 or return 1;
- # Check for very small factors
- if (ULONG_MAX >= 18446744073709551615) {
- Math::GMPz::Rmpz_gcd_ui($Math::GMPz::NULL, $n, 16294579238595022365) == 1 or return 0;
- Math::GMPz::Rmpz_gcd_ui($Math::GMPz::NULL, $n, 7145393598349078859) == 1 or return 0;
- }
- else {
- Math::GMPz::Rmpz_gcd_ui($Math::GMPz::NULL, $n, 3234846615) == 1 or return 0;
- }
- # Size of n in base-2
- my $size = Math::GMPz::Rmpz_sizeinbase($n, 2);
- # When n is large enough, try to find a small factor (up to 10^8)
- if ($size > 10_000) {
- state %cache;
- state $g = Math::GMPz::Rmpz_init_nobless();
- my @checks = (1e4);
- push(@checks, 1e6) if ($size > 15_000);
- push(@checks, 1e7) if ($size > 20_000);
- push(@checks, 1e8) if ($size > 30_000);
- my $prev;
- foreach my $k (@checks) {
- my $primorial = (
- $cache{$k} //= do {
- my $z = Math::GMPz::Rmpz_init_nobless();
- Math::GMPz::Rmpz_primorial_ui($z, $k);
- Math::GMPz::Rmpz_divexact($z, $z, $prev) if defined($prev);
- $z;
- }
- );
- Math::GMPz::Rmpz_gcd($g, $primorial, $n);
- if (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0) {
- return 0;
- }
- $prev = $primorial;
- }
- }
- return 1;
- }
- # "HCN.bz2" was generated by Achim Flammenkamp, and is available at:
- # http://wwwhomes.uni-bielefeld.de/achim/HCN.bz2
- my $z = IO::Uncompress::Bunzip2->new("HCN.bz2");
- my $tmp = Math::GMPz->new(1);
- while (defined(my $line = $z->getline())) {
- my @fields = split(' ', $line);
- my $len = shift(@fields);
- if ($len == 0) {
- next;
- }
- my @primes = @{primes(nth_prime($len))};
- my $prod = Math::GMPz->new(1);
- while (@primes) {
- my $k = shift(@fields) // die "error: $line";
- my $e = 1;
- if ($k =~ /^(\d+)\^(\d+)\z/) {
- $k = $1;
- $e = $2;
- }
- for (1 .. $e) {
- my $p = shift(@primes);
- if ($k == 1) {
- Math::GMPz::Rmpz_mul_ui($prod, $prod, $p);
- }
- elsif ($p**$k < ULONG_MAX) {
- Math::GMPz::Rmpz_mul_ui($prod, $prod, powint($p, $k));
- }
- else {
- Math::GMPz::Rmpz_ui_pow_ui($tmp, $p, $k);
- Math::GMPz::Rmpz_mul($prod, $prod, $tmp);
- }
- }
- }
- if ( primality_pretest($prod - 1)
- and primality_pretest($prod + 1)
- and is_prob_prime($prod + 1)
- and is_prob_prime($prod - 1)) {
- print $., ", ";
- }
- }
|