From bf27525b9e6f6e315b1d5d7bdefdb406b6e7ad69 Mon Sep 17 00:00:00 2001
From: Abigail <abigail@abigail.be>
Date: Wed, 14 Oct 2020 13:55:17 +0200
Subject: Improve the performance of week 082/part 1.

Improved the running time of the algorithm from O (sqrt (gcd ($M, $N)))
to O (sqrt (p) + |f| log |f|) where p is the largest prime factor
of gcd ($M, $N), and |f| is total number of factors.
This is not an improvement if the gcd is a prime, but a dramatic
improvement if the gcd contains lots of prime factors.
---
 challenge-082/abigail/perl/ch-1.pl | 77 ++++++++++++++++++++++++++++++++------
 1 file changed, 65 insertions(+), 12 deletions(-)

diff --git a/challenge-082/abigail/perl/ch-1.pl b/challenge-082/abigail/perl/ch-1.pl
index 34db614e9b..76ec111012 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,72 @@ 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, we try *all* numbers starting from 2,
+# not just primes. But any composite number will have a prime factor
+# smaller than itself, and this will already be divided out of $gcd,
+# so for any composite number $n, $n will never divide $gcd evenly.
 #
-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 ++;
 }
-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__
-- 
cgit