Merge pull request #2524 from Abigail/abigail/week-082

Improve the performance of week 082/part 1.

1 files changed, 70 insertions, 12 deletions

diff --git a/challenge-082/abigail/perl/ch-1.pl b/challenge-082/abigail/perl/ch-1.pl
index 34db614e9b..6788e1ae45 100644
--- a/challenge-082/abigail/perl/ch-1.pl
+++ b/challenge-082/abigail/perl/ch-1.pl
@@ -37,7 +37,6 @@ use experimental 'lexical_subs';
 chomp (my $M = <>);
 chomp (my $N = <>);
 
-
 #
 # Find the GCD, using Stein's algorithm
 #    (https://en.wikipedia.org/wiki/Binary_GCD_algorithm)
@@ -59,18 +58,77 @@ sub stein ($u, $v) {
 my $gcd = stein $M, $N;
 
 #
-# Brute force finding all factors
+# Now we want to calculate the list of factors of $gcd. 
+#
+# We do this by finding prime factors $p of $gcd. For each such
+# $p, we find how often ($t) it divides $gcd. That is,
+# $p ^ $t divides $gcd, but $p ^ ($t + 1) does not.
+#
+# We also keep a running list of factors (@f), containing all factors
+# of the form q_1^a * .. * q_i^b, with q_1 < ... < q_i < p. This list
+# is initialized to 1.
+#
+# For each new pair ($p, $t) found, and each $f in @f, we add
+# $f * $p ^ $i, 1 <= $i <= $t to @f. We can then also divide
+# $gcd by $p ^ $t.
+#
+# Note that in the loop below, ideally, we want to check only
+# prime numbers (any composite number will have a prime factor
+# smaller than itself). Since, except for 2 and 3, all prime
+# numbers are of the form 6 * k +/- 1, the loop below checks
+# 2, 3, and then numbers of the form 6 * k +/- 1.
 #
-my $max = int sqrt $gcd;
-my @small;   # All factors <= sqrt ($gcd);
-my @large;   # All factors >  sqrt ($gcd);
-for (my $i = 1; $i <= $max; $i ++) {
-    next if $gcd % $i;
-    push @small => $i;
-    # If $gcd is a perfect square, we should not report
-    # its square root twice.
-    push @large => $gcd / $i unless $gcd / $i == $i;
+# Also note that we stop if the current number is larger than the
+# square of $gcd. This means that, at the end of the loop, $gcd may
+# be a prime $r. In that case, for each found factor $f so far, we
+# must add $r * $f to the list of factors.
+# 
+#
+# The run time complexity is O (sqrt (p) + |f| log |f|), where
+# p is largest prime factor, and |f| is the number of factors found.
+#
+
+my @factors = (1);
+my $p = 2;
+while ($p ** 2 <= $gcd) {
+    #
+    # Find out how often $p is a factor of $gcd
+    #
+    my $t = 0;
+    until ($gcd % $p) {
+        #
+        # If we are here, $p cannot be composite.
+        #
+        $t ++;
+        $gcd /= $p;
+    }
+
+    #
+    # Push all new factors. For each 1 <= $i <= $t, and each $f already
+    # in the list of factors, add $f * $p ^ $i to the list of factors.
+    #
+    # The test for $t isn't needed for correctness, just for performance.
+    #
+    if ($t) {
+        push @factors => map {my $f = $_; map {$f * $p ** $_} 1 .. $t} @factors;
+    }
+
+    $p += $p     == 2 ? 1
+        : $p     == 3 ? 2
+        : $p % 6 == 1 ? 4
+        :               2;
+
 }
-say for @small, reverse @large;
+
+#
+# Case where $gcd is left to be a prime.
+#
+push @factors => map {$_ * $gcd} @factors if $gcd > 1;
+
+#
+# Print the list of factors, after sorting.
+#
+say for sort {$a <=> $b} @factors;
+
 
 __END__