12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364 |
- #!/usr/bin/perl
- # Daniel "Trizen" Șuteu
- # Date: 27 June 2019
- # https://github.com/trizen
- # Based on the identity:
- # Sum_{d|n} 2^omega(d) = sigma_0(n^2)
- # The algorithm iterates over each number k in 2..N,
- # and computes sigma_0(k^2) = Prod_{p^e | k} (2*e + 1).
- # By keeping track of the partial products, we find sigma_0(k!^2).
- # Runtime: 44.576s
- # https://projecteuler.net/problem=675
- use 5.020;
- use strict;
- use warnings;
- use ntheory qw(:all);
- use experimental qw(signatures declared_refs);
- sub F ($N, $mod) {
- my @table;
- my $total = 0;
- my $partial = 1;
- foreach my $k (2 .. $N) {
- my $old = 1; # old product
- foreach my \@pp (factor_exp($k)) {
- my ($p, $e) = @pp;
- $old = mulmod($old, 2 * $table[$p] + 1, $mod) if defined($table[$p]);
- $table[$p] += $e;
- $partial = mulmod($partial, 2 * $table[$p] + 1, $mod);
- }
- $partial = divmod($partial, $old, $mod); # remove the old product
- $total += $partial;
- $total %= $mod;
- }
- return $total;
- }
- my $MOD = 1000000087;
- foreach my $n (1 .. 7) {
- printf("F(10^%d) = %*s (mod %s)\n", $n, length($MOD) - 1, F(10**$n, $MOD), $MOD);
- }
- __END__
- F(10^1) = 4821 (mod 1000000087)
- F(10^2) = 930751395 (mod 1000000087)
- F(10^3) = 822391759 (mod 1000000087)
- F(10^4) = 979435692 (mod 1000000087)
- F(10^5) = 476618093 (mod 1000000087)
- F(10^6) = 420600586 (mod 1000000087)
|