123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124 |
- #!/usr/bin/perl
- # Unitary harmonic numbers (A006086) that are not unitary arithmetic numbers (A103826).
- # https://oeis.org/A353038
- # Known terms:
- # 90, 40682250, 81364500, 105773850, 423095400, 1798155450, 14385243600
- # No other terms are known.
- 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) < 7;
- }
- if ($p == 3) {
- return valuation($n, $p) < 7;
- }
- if ($p == 5) {
- return valuation($n, $p) < 6;
- }
- if ($p == 7) {
- return valuation($n, $p) < 5;
- }
- #~ if ($p == 11) {
- #~ return valuation($n, $p) < 2;
- #~ }
- if ($p == 13) {
- return valuation($n, $p) < 2;
- }
- if ($p == 17) {
- return valuation($n, $p) < 2;
- }
- if ($p == 41) {
- return valuation($n, $p) < 2;
- }
- ($n % $p) != 0;
- }
- sub smooth_numbers ($limit, $primes) {
- my @h = (Math::GMPz->new(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;
- }
- #use Math::Sidef qw(usigma usigma0);
- sub isok ($n) {
- #is_power(Math::GMPz->new(divisor_sum($n)) * euler_phi($n));
- # is(n)=my(f=factor(n)); prod(i=1, #f~, f[i, 1]^f[i, 2]+1)%2^#f~==0
- my @f = factor_exp($n);
- my $usigma = vecprod(map { powint($_->[0], $_->[1]) + 1 } @f);
- my $usigma0 = powint(2, scalar(@f));
- modint($usigma, $usigma0) == 0
- and return;
- modint(mulint($n, $usigma0), $usigma) == 0;
- #(usigma0($n)*$n) % usigma($n) == 0;
- }
- my @smooth_primes;
- foreach my $p (@{primes(5000)}) {
- if ($p == 2) {
- push @smooth_primes, $p;
- next;
- }
- my @f1 = factor($p - 1);
- my @f2 = factor($p + 1);
- if ($f1[-1] <= 5 and $f2[-1] <= 7) {
- push @smooth_primes, $p;
- }
- }
- my $h = smooth_numbers(Math::GMPz->new(10)**30, \@smooth_primes);
- say "\nFound: ", scalar(@$h), " terms";
- my %table;
- foreach my $n (@$h) {
- $n > 1e12 or next;
- $n % 2 == 0 or next;
- if (isok($n)) {
- say $n;
- }
- }
|