1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192 |
- #!/usr/bin/perl
- # The number of terms of A354558 that are <= 10^n.
- # https://oeis.org/A354559
- # Known terms:
- # 1, 2, 5, 13, 28, 79, 204, 549, 1509, 4231, 12072, 36426, 112589
- use 5.020;
- use strict;
- use warnings;
- use ntheory qw(:all);
- use experimental qw(signatures);
- sub smooth_numbers ($limit, $primes) {
- if ($limit <= $primes->[-1]) {
- return [1 .. $limit];
- }
- if ($limit <= 5e4) {
- my @list;
- my $B = $primes->[-1];
- foreach my $k (1 .. $limit) {
- if (is_smooth($k, $B)) {
- push @list, $k;
- }
- }
- return \@list;
- }
- my @h = (1);
- foreach my $p (@$primes) {
- foreach my $n (@h) {
- if ($n * $p <= $limit) {
- push @h, $n * $p;
- }
- }
- }
- return \@h;
- }
- sub upto ($k, $j) {
- my $limit = rootint($k, $j);
- my @smooth;
- my @primes;
- my $i = 0;
- my $pi = prime_count($limit);
- foreach my $p (@{primes($limit)}) {
- ++$i;
- say "[$i / $pi] Processing prime $p";
- my $pj = powint($p, $j);
- push @primes, $p;
- push @smooth, grep {
- my $m = addint($_, 1);
- valuation($m, (factor($m))[-1]) >= $j;
- } map { mulint($_, $pj) } @{smooth_numbers(divint($k, $pj), \@primes)};
- }
- return sort { $a <=> $b } @smooth;
- }
- #~ my $n = 7;
- #~ say join(', ', upto($n, 2));
- #~ __END__
- my $n = 13;
- my $i = 0;
- open my $fh, '>', 'bfile.txt';
- foreach my $k (upto(powint(10, $n), 2)) {
- my $row = sprintf("%s %s\n", ++$i, $k);
- print $row;
- print $fh $row;
- }
- close $fh;
|