123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123 |
- #!/usr/bin/perl
- # Numbers that are both infinitary and noninfinitary harmonic numbers
- # https://oeis.org/A348922
- # Known terms:
- # 45, 60, 54600, 257040, 1801800, 2789640, 4299750, 47297250, 1707259680, 4093362000
- # a(7) > 10^10, if it exists. - Amiram Eldar, Nov 04 2021
- # Equivalently, numbers n such that isigma(n) divides n*isigma_0(n) and nisigma(n) divides n*nisigma_0(n).
- # Non-squarefree numbers n such that A049417(n) divides n*A037445(n) and A348271(n) divides n*A348341(n).
- # See also:
- # https://oeis.org/A247077
- # The sequence also includes:
- # 18779856480
- # 44425017000
- # 13594055202000
- # 27188110404000
- # 299069214444000
- # 6824215711404000
- use 5.020;
- use warnings;
- use experimental qw(signatures);
- use Math::GMPz;
- use ntheory qw(:all);
- use Math::Sidef qw(isigma isigma0 nisigma nisigma0);
- sub check_valuation ($n, $p) {
- if ($p == 2) {
- return valuation($n, $p) < 10;
- }
- if ($p == 3) {
- return valuation($n, $p) < 12;
- }
- if ($p == 5) {
- return valuation($n, $p) < 7;
- }
- if ($p == 7) {
- return valuation($n, $p) < 7;
- }
- if ($p == 11) {
- return valuation($n, $p) < 3;
- }
- ($n % $p) != 0;
- #valuation($n, $p) < 2;
- }
- 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 usigma($n) {
- vecprod(map { addint(powint($_->[0], $_->[1]), 1) } factor_exp($n));
- }
- sub isok ($n) {
- is_square_free($n) && return;
- (($n*isigma0($n)) % isigma($n)) == 0 or return;
- my $t = nisigma($n);
- $t == 0 and return;
- (($n*nisigma0($n)) % $t) == 0 or return;
- return 1;
- }
- my @smooth_primes;
- foreach my $p (@{primes(4801)}) {
- if ($p == 2) {
- push @smooth_primes, $p;
- next;
- }
- if (
- is_smooth($p-1, 5) and
- is_smooth($p+1, 7)
- ) {
- push @smooth_primes, $p;
- }
- }
- my $h = smooth_numbers(~0, \@smooth_primes);
- say "\nGenerated: ", scalar(@$h), " numbers";
- foreach my $n (@$h) {
- if ($n > 1e10 and isok($n)) {
- #if (isok($n)) {
- say $n;
- }
- }
|