123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107 |
- #!/usr/bin/perl
- # Daniel "Trizen" Șuteu
- # Date: 25 July 2017
- # https://github.com/trizen
- # An efficient algorithm for computing sigma_k(n!), where k > 0.
- use 5.020;
- use strict;
- use warnings;
- use experimental qw(signatures);
- use ntheory qw(forprimes vecsum todigits :all);
- use Math::AnyNum qw(ipow factorial);
- prime_precalc(3e7);
- my $t = Math::GMPz::Rmpz_init();
- sub sigma_of_factorial ($n) {
- my $sigma = Math::GMPz::Rmpz_init_set_ui(1);
- forprimes {
- my $p = $_;
- my $k = ($n - vecsum(todigits($n, $p))) / ($p - 1);
- if (log($p) * ($k + 1) > log(2) * 63) {
- #$sigma *= ((ipow($p, ($k + 1)) - 1) / ($p - 1));
- Math::GMPz::Rmpz_ui_pow_ui($t, $p, $k + 1);
- Math::GMPz::Rmpz_sub_ui($t, $t, 1);
- Math::GMPz::Rmpz_divexact_ui($t, $t, $p - 1);
- Math::GMPz::Rmpz_mul($sigma, $sigma, $t);
- }
- else {
- #$sigma *= divint(powint($p, $k+1) - 1, $p-1);
- Math::GMPz::Rmpz_mul_ui($sigma, $sigma, divint(powint($p, $k + 1) - 1, $p - 1));
- }
- } $n;
- return $sigma;
- }
- #say sigma_of_factorial(10, 1); # sigma_1(10!) = 15334088
- #say sigma_of_factorial(10, 2); # sigma_2(10!) = 20993420690550
- #say sigma_of_factorial( 8, 3); # sigma_3( 8!) = 78640578066960
- # 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, 230773);
- 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);
- say "Sanity checking...";
- foreach my $n (0 .. $#array) {
- my $k = $array[$n];
- say "Testing: $k";
- my $s = Math::AnyNum->new(sigma_of_factorial($k));
- my $t = $n * factorial($k);
- $s >= $t or die "error for $k";
- next if $k == 1;
- my $s2 = Math::AnyNum->new(sigma_of_factorial($k - 1));
- my $t2 = $n * factorial($k - 1);
- $s2 >= $t2 and die "too large: $k";
- }
- say "Test passed...";
- #say sigma_of_factorial(75037,1);
- # 1, 1, 3, 5, 9, 14, 23, 43, 79, 149, 263, 461, 823, 1451, 2549, 4483, 7879, 13859, 24247, 42683, 75037
- my $n = 0;
- local $| = '1';
- my $k = 1;
- my $factorial = factorial($k);
- $factorial = Math::GMPz->new("$factorial");
- for (; $k <= 1e6 ; ++$k) {
- while (sigma_of_factorial($k) >= $n * $factorial) {
- print($k, ", ");
- ++$n;
- }
- $factorial *= ($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
|