123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110 |
- #!/usr/bin/perl
- # Least k > 0 such that sigma(k!) >= n*k!, with a(0) = a(1) = 1.
- # https://oeis.org/A061556
- # Known terms:
- # 1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, 1451, 2549, 4483, 7879, 13859, 24247, 42683, 75037
- use 5.020;
- use strict;
- use warnings;
- use experimental qw(signatures);
- use ntheory qw(forprimes vecsum todigits :all);
- use Math::AnyNum qw(ipow factorial lgamma log );
- sub factorial_power ($n, $p) {
- ($n - vecsum(todigits($n, $p))) / ($p - 1);
- }
- my $t = Math::GMPz::Rmpz_init();
- sub sigma_of_factorial ($n) {
- #my $sigma = Math::GMPz::Rmpz_init_set_ui(1);
- my $sigma = 0;
- forprimes {
- my $p = $_;
- my $k = factorial_power($n, $p);
- # (p^(k+1) - 1) / (p-1)
- $sigma += (log($p) * ($k + 1)) - log($p - 1);
- }
- $n;
- return $sigma;
- }
- say sigma_of_factorial(42683);
- say sigma_of_factorial(75037);
- say sigma_of_factorial(126000);
- # Least k > 0 such that sigma(k!) >= n*k!.
- my @array = (1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, 1451, 2549, 4483, 7879, 13859, 24247, 42683, 75037, 131707);
- #my @array = (1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, 1451, 2549, 4483, 7879, 13859, 24247, 42683, 75037, 131707, 230773, 405401, 710569, 1246379, 2185021, 3831913, 6720059, 11781551, 20657677);
- #my @array = (0..12, 1451, 2549, 4483, 7879, 13859, 24247, 42683, 75037, 131707, 230773, 405401, 710569, 1246379, 2185021, 3831913, 6720059, 11781551, 20657677, 36221753, 63503639, 111333529, 195199289, 342243479, 600036989);
- say "Sanity checking...";
- foreach my $n (0 .. $#array) {
- my $k = $array[$n];
- say "[$n] Testing: $k";
- next if $k < 1000;
- next if $k == 1;
- my $s = Math::AnyNum->new(sigma_of_factorial($k));
- my $t = log($n) + lgamma($k + 1);
- $s >= $t or die "error for $k";
- next if $k < 1000;
- my $s2 = sigma_of_factorial($k - 1);
- my $t2 = log($n) + lgamma($k);
- $s2 >= $t2 and die "too large: $k -- ($s2 >= $t2)";
- }
- say "Test passed...";
- my $n = 21;
- local $| = '1';
- my $k = 131707 - 100;
- my $factorial = lgamma($k + 1);
- #$factorial = Math::GMPz->new("$factorial");
- for (; $k <= 1e6 ; ++$k) {
- say "Testing: $k" if ($k % 1000 == 0);
- while (sigma_of_factorial($k) >= log($n) + $factorial) {
- print($k, ", ");
- if ($k > 131707) {
- die "Found: $k";
- }
- ++$n;
- }
- #$factorial *= ($k + 1);
- #Math::GMPz::Rmpz_mul_ui($factorial, $factorial, $k+1);
- $factorial += log($k + 1);
- }
- __END__
- 1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, perl r.pl 2.13s user 0.01s system 99% cpu 2.140 total
- 1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, perl r.pl 0.77s user 0.02s system 99% cpu 0.794 total
- 1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, perl r.pl 0.77s user 0.02s system 99% cpu 0.788 total
- 1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, perl r.pl 0.40s user 0.02s system 99% cpu 0.423 total
- 1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, perl r.pl 0.19s user 0.01s system 99% cpu 0.203 total
- 1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, 1451, 2549, 4483, perl r.pl 2.94s user 0.01s system 99% cpu 2.967 total
- 1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, 1451, 2549, 4483, perl r.pl 2.95s user 0.02s system 99% cpu 2.974 total
- 1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, 1451, 2549, 4483, perl r.pl 2.93s user 0.02s system 98% cpu 3.003 total
|