123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899 |
- #!/usr/bin/perl
- # a(n) is the least integer such that sigma(sigma(k)) >= n*k where sigma is A000203, the sum of divisors.
- # https://oeis.org/A327630
- # New terms:
- # a(23) = 109549440
- # a(24) = 438197760
- # Upper-bounds for other terms:
- # a(25) <= 766846080
- # a(26) <= 3834230400
- # a(27) <= 9081072000
- # a(28) <= 32974381440
- # a(29) <= 147516969600
- # a(30) <= 880887047040
- # a(31) <= 2802822422400
- # a(32) <= 14814918518400
- # a(33) <= 64464915715200
- # a(34) <= 709114072867200
- # a(35) <= 9881550651772800
- # a(36) <= 76648784785372800
- # a(37) <= 2376112328346556800
- # See also:
- # https://oeis.org/A098221
- 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) < 10;
- }
- if ($p == 3) {
- return valuation($n, $p) < 4;
- }
- if ($p == 5) {
- return valuation($n, $p) < 4;
- }
- if ($p == 7) {
- return valuation($n, $p) < 3;
- }
- if ($p == 11) {
- return valuation($n, $p) < 3;
- }
- if ($p == 13) {
- return valuation($n, $p) < 2;
- }
- ($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;
- }
- my $h = smooth_numbers(10**17, primes(57));
- say "\nFound: ", scalar(@$h), " terms";
- @$h = sort { $a <=> $b } @$h;
- my $n = 20;
- foreach my $k (@$h) {
- next if ($k < 360360);
- while (divisor_sum(divisor_sum($k, 1)) >= $k * $n) {
- say "a($n) <= $k";
- ++$n;
- }
- }
|