123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102 |
- #!/usr/bin/perl
- # Positions of records in A306440.
- # https://oeis.org/A307221
- # Where A306440(n) is defined as:
- # Number of different ways of expressing 2*n as (p - 1)*(q - 1), where p < q are primes
- use 5.020;
- use warnings;
- use experimental qw(signatures);
- use Math::GMPz;
- use ntheory qw(:all);
- sub check_valuation ($n, $p) {
- if ($p == 2) {
- return valuation($n, $p) < 5;
- }
- if ($p == 3) {
- return valuation($n, $p) < 3;
- }
- if ($p == 7) {
- return valuation($n, $p) < 3;
- }
- ($n % $p) != 0;
- }
- sub smooth_numbers ($limit, $primes) {
- my @h = (1);
- foreach my $p (@$primes) {
- say "Prime: $p";
- foreach my $n (@h) {
- if ($n * $p <= $limit) { #and check_valuation($n, $p)) {
- push @h, $n * $p;
- }
- }
- }
- return \@h;
- }
- sub decompositions {
- my ($n) = @_;
- #say "N -> $n";
- my $count = 0;
- #my $max = sqrtint($n);
- foreach my $d(divisors($n)) {
- #last if ($d >= $max);
- $d < $n/$d or last;
- if (is_prime($d+1) and is_prime(($n/$d)+1)) {
- #say "$d, ", $n/$d, ' -> ', $d * ($n/$d);
- ++$count;
- }
- }
- $count;
- }
- my $h = smooth_numbers(275675400, primes(43));
- say "\nFound: ", scalar(@$h), " terms";
- my %table;
- foreach my $n (@$h) {
- $n > 151351200 or next;
- my $p = decompositions(2*$n);
- # $p > 32 or next;
- if (not exists $table{$p}) {
- $table{$p} = $n;
- }
- else {
- if ($table{$p} > $n) {
- $table{$p} = $n;
- }
- }
- #if ($p >= 8) {
- # say "a($p) = $n -> ", join(' * ', map { "$_->[0]^$_->[1]" } factor_exp($n));
- # # push @{$table{$p}}, $n;
- # }
- }
- use Data::Dump qw(pp);
- pp \%table;
|