Merge pull request #5880 from drbaggy/master

First pass at 159

5 files changed, 170 insertions, 86 deletions

diff --git a/challenge-158/james-smith/README.md b/challenge-158/james-smith/README.md
index 342202ccbf..a81bb61a28 100644
--- a/challenge-158/james-smith/README.md
+++ b/challenge-158/james-smith/README.md
@@ -94,7 +94,7 @@ sub additive_primes_split_sum0 {
  * `split`            -  75%
  * `sum0 split`       -  45%
 
-# Challenge 2 - First series buban primes
+# Challenge 2 - First series cuban primes
 
 ***Write a script to compute first series cuban primes <= 1000. (First series cuban primes have the form `((x+1)^3-x^3)/(x+1-x)` = `3x^2+3x+1`)***
 
diff --git a/challenge-159/james-smith/README.md b/challenge-159/james-smith/README.md
index 342202ccbf..31d01d0ee1 100644
--- a/challenge-159/james-smith/README.md
+++ b/challenge-159/james-smith/README.md
@@ -1,6 +1,6 @@
-[< Previous 157](https://github.com/drbaggy/perlweeklychallenge-club/tree/master/challenge-157/james-smith) |
-[Next 159 >](https://github.com/drbaggy/perlweeklychallenge-club/tree/master/challenge-159/james-smith)
-# The Weekly Challenge 158
+[< Previous 156](https://github.com/drbaggy/perlweeklychallenge-club/tree/master/challenge-156/james-smith) |
+[Next 160 >](https://github.com/drbaggy/perlweeklychallenge-club/tree/master/challenge-160/james-smith)
+# The Weekly Challenge 159
 
 You can find more information about this weeks, and previous weeks challenges at:
 
@@ -12,120 +12,104 @@ submit solutions in whichever language you feel comfortable with.
 
 You can find the solutions here on github at:
 
-https://github.com/drbaggy/perlweeklychallenge-club/tree/master/challenge-158/james-smith
+https://github.com/drbaggy/perlweeklychallenge-club/tree/master/challenge-159/james-smith
 
-# Challenge 1 - Additive primes
+# Challenge 1 - Farey Sequence
 
-***Additive primes are prime numbers for which the sum of their decimal digits are also primes.***
+***You are given a positive number, `$n`. Write a script to compute Farey Sequence of the order `$n`. All rational numbers less than with with denominators less than `$n`***
 
 ## The solution
 
-We loop through each prime p, work out the digit sum (by repeated modulo 10/divide by 10) and check that it is prime.
-We craft this as two loops - an outer `for` loop and an inner `do {} while`.
+This is a relatively simple piece of code - we loop through the denominators `1`..`$n` and compute all the fractions greater than 0 and less than or equal to 1.
+We use the value as the key to the hash and "`$_/$d`" as the values.
+We start with the lower denominators so if we find the same fraction twice (e.g. `1/2` & `2/4`) the value we always store has the lowest denominator and so the fraction is in its simplest form. A simple (key based) sort at the end and returning the list of strings of fractions...
 
 ```perl
-use Math::Prime::Util qw(next_prime is_prime);
-
-sub additive_primes {
-  my @res;
-  for(
-    my $p = 2                                             ;
-    my $s = 0,                        ( my $q = $p ) <= $N;
-    (is_prime $s) && (push @res, $p), $p = next_prime $p
-  ) {
-    do { $s += $q % 10 } while $q = int $q / 10;
+sub farley {
+  my($n,%v) = shift;
+  for my $d (1..$n) {
+    $v{$_/$d}||="$_/$d" for 1..$d;
   }
-  @res;
+  map{$v{$_}}sort{$a<=>$b}keys%v;
 }
 ```
-## Output
-```
-2, 3, 5, 7, 11, 23, 29, 41, 43, 47, 61, 67, 83, 89
-```
-## Notes
- * We use a C-style `for` loop, with it's three parts *initialization*, *condition* and *increment*.
- * The *condition* and *increment* statements are complicated, each with two parts separated by `,`s.
- * The *condition* block is called at the start of each loop, and so we use it to initialise the variables as the loop, as well as checking the condition.
- * The *increment* block is called at the end of the loop and so stores the value of `$p` if it is an additive prime, then it increments the loop with the next prime.
- * Rather than doing a split and sum we use repeated dividing and summing, as it is more efficient around 20-30% more efficient. The increased performance is probably due to avoiding the "duality" of perl variables storing numbers as numbers/strings.
 
-## Extra code
+# Challenge 2 - Moebius Number
 
-Rewritten with single line `for` ... (the original version of the code)
+***You are given a positive number `$n`. Write a script to generate the Moebius Number for the given number (see defn below)***
 
-```perl
-sub additive_primes_div {
-  my @res;
-  for(my$p=2; my$s=0,(my $q=$p)<=$N;(is_prime$s)&&(push@res,$p),$p=next_prime$p) {
-    do { $s += $q % 10 } while $q = int $q / 10;
-  }
-  @res;
-}
-```
+## Definition
 
-Alternative form with `for split`:
+The value of the Moebius number is:
 
-```perl
-sub additive_primes_split {
-  my @res;
-  for( my $p = 2; my $s = 0, $p <= $N ; $p = next_prime $p ) {
-    $s+=$_ for split //, $p;
-    push @res, $p if is_prime $s;
-  }
-  @res;
-}
-```
+ * `0` if `$n` has a prime factor
+ * `1` if the number has an even number of prime factors
+ * `-1` otherwise
 
-And an alternate using `sum0`...
+## First pass... *naive* approach
 
 ```perl
-sub additive_primes_split_sum0 {
-  my @res;
-  for( my $p = 2; $p <= $N ; $p = next_prime $p ) {
-    push @res, $p if is_prime sum0 split //, $p;
-  }
-  @res;
+sub moebius {
+  my ($n,$p,$r) = (shift,1,1);
+  return -1 if is_prime($n);
+  $n%($p**2) ? ($n%$p && ($r=-$r)) : return 0 while ($p = next_prime $p) < $n;
+  $r;
 }
 ```
 
-### Relative performance
+Here we just look for prime factors less than `$n` - for large `$n` - the number of primes checked is `$n/log $n`.
 
- * `div`/`mod` method - 100%
- * `split`            -  75%
- * `sum0 split`       -  45%
+## Second pass... *first-pass* optimization
 
-# Challenge 2 - First series buban primes
+We can reduce the number of primes checked by reducing `$n` each time a prime is found, by dividing `$n` by the prime factor.
 
-***Write a script to compute first series cuban primes <= 1000. (First series cuban primes have the form `((x+1)^3-x^3)/(x+1-x)` = `3x^2+3x+1`)***
+```perl
+sub moebius_div {
+  my ($n,$p,$r) = (shift,1,1);
+  return -1 if is_prime $n;
+  $n%($p**2) ? ( $n%$p || ($r=-$r,$n/=$p)) : return 0 while ($p = next_prime $p) && $n>1;
+  $r;
+}
+```
 
-## The solution
+## Third pass... *improved* performance
 
-The solution is rather short. We try each of the value of the sequence `3*$x**2+3*$x+1` in turn up to the value of 1000 (`$N`).
-We output each value which in turn is prime.
+Finally we ideally want to restrict the number of primes checked further. We know that if `$n/$p` is prime we can stop there when looking for factors. It then also reduces the number primes to check to the square root of `$n`, so the maximum number of primes needed to check is `sqrt($n)/log $n/2`.
 
 ```perl
-(is_prime $_) && say while $N >= ($_ = 3*$x*++$x+1);
-```
-## Output
-```
-7, 19, 37, 61, 127, 271, 331, 397, 547, 631, 919
+sub moebius_div_opt {
+  my ($n,$p,$r) = (shift,1,1);
+  return -1 if is_prime $n;
+  $n%$p || ($n/=$p) && $n%$p ? ( is_prime $n ? (return $r) : ($r=-$r) ) : return 0 while ($p = next_prime $p)**2 <= $n;
+  $r;
+}
 ```
-## Notes
- * As we use `$_` as our temporary variable we can use `say` by itself to output it.
- * We use our common trick of `(condition) && (command)` to work as an `command if(condition)` which can be embedded in a postfix loop.
- * There is a "trick" as we increment `$x` half way through the calculation of `$_`.
-   * `$_ = 3*$x**2 + 3*$x + 1` => `$_ = 3 * $x * ($x+1) + 1` => `$_ = 3 * $x * ++$x + 1`
-   * We can replace the `$x+1` by the pre increment operator `++$x` so this becomes `3 * $x * ++$x + 1`.
 
-## Aside
+## Relative performance
 
-Second series cuban primes have the form `((x+2)^3-x^3)/(x+2-x)` = `3x^2+3.2x+4`. We can tweak the code to give:
+On our extended test set the relative performance of the three methods is:
 
-```perl
-(is_prime $_) && say while $N >= ($_ = 3 * $x * (2 + $x++) + 4);
-```
+| Method  | Calls/second | Relative performance |
+| :-----: | -----------: | -------------------: |
+| simple  |     0.0025/s |                  1 x |
+| divide  |        3.0/s |              1,200 x |
+| div opt |     38,000/s |         16,000,000 x |
+
+## Expanded version of optimal solution..
+
+The format of the optimized code is compact - here is an expanded version of the code to show the logic.
 
-which outputs:
 ```
-13, 109, 193, 433, 769
+sub moebius {
+  my ($n,$p,$r) = (shift,1,1);
+  return -1 if is_prime $n;
+  while( ($p = next_prime $p )**2 <= $n ) {
+    next if $n%$p;
+    $n/=$p;
+    return 0 unless $n%$p;
+    return $r if is_prime $n;
+    $r=-$r;
+  }
+  $r;
+}
 ```
diff --git a/challenge-159/james-smith/blog.txt b/challenge-159/james-smith/blog.txt
new file mode 100644
index 0000000000..e631de6584
--- /dev/null
+++ b/challenge-159/james-smith/blog.txt
@@ -0,0 +1 @@
+https://github.com/drbaggy/perlweeklychallenge-club/tree/master/challenge-159/james-smith
diff --git a/challenge-159/james-smith/perl/ch-1.pl b/challenge-159/james-smith/perl/ch-1.pl
new file mode 100644
index 0000000000..5f6402ec41
--- /dev/null
+++ b/challenge-159/james-smith/perl/ch-1.pl
@@ -0,0 +1,20 @@
+#!/usr/local/bin/perl
+
+use strict;
+
+use warnings;
+use feature qw(say);
+use Test::More;
+use Benchmark qw(cmpthese timethis);
+use Data::Dumper qw(Dumper);
+
+say join ', ', farley( $_ ) for 1..10;
+
+sub farley {
+  my($n,%v) = shift;
+  for my $d (1..$n) {
+    $v{$_/$d}||="$_/$d" for 1..$d;
+  }
+  map{$v{$_}}sort{$a<=>$b}keys%v;
+}
+
diff --git a/challenge-159/james-smith/perl/ch-2.pl b/challenge-159/james-smith/perl/ch-2.pl
new file mode 100644
index 0000000000..fb9ee8a893
--- /dev/null
+++ b/challenge-159/james-smith/perl/ch-2.pl
@@ -0,0 +1,79 @@
+#!/usr/local/bin/perl
+
+use strict;
+
+use warnings;
+use feature qw(say);
+use Test::More;
+use Benchmark qw(cmpthese timethis);
+use Data::Dumper qw(Dumper);
+use Math::Prime::Util qw(next_prime is_prime);
+
+my @TESTS = (
+  [  5, -1 ],
+  [ 10,  1 ],
+  [ 20,  0 ],
+## The following are good examples where the second algorithm works well!
+  [ 2,   -1 ],
+  [ 6,   1 ],
+  [ 30,   -1 ],
+  [ 210,   1 ],
+  [ 2310,   -1 ],
+  [ 30030,   1 ],
+  [ 510510,  -1 ],
+  [ 9699690, 1 ],
+  [ 6895531,   -1 ],
+  [ 206865930,  1 ],
+  [ 223092870, -1 ],
+  [ 6469693230, 1 ],
+);
+
+warn "exp";
+is( moebius_exp(     $_->[0] ), $_->[1] ) for @TESTS;
+warn "opt";
+is( moebius_div_opt( $_->[0] ), $_->[1] ) for @TESTS;
+#timethis( 500_000, sub { moebius_div_opt( $_->[0] ) for @TESTS } );
+exit;
+warn "div";
+is( moebius_div(     $_->[0] ), $_->[1] ) for @TESTS;
+timethis(     100, sub { moebius_div(     $_->[0] ) for @TESTS } );
+warn "dumb";
+is( moebius(         $_->[0] ), $_->[1] ) for @TESTS;
+timethis(       5, sub { moebius(         $_->[0] ) for @TESTS } );
+
+done_testing();
+
+sub moebius {
+  my ($n,$p,$r) = (shift,1,1);
+  return -1 if is_prime($n);
+  $n%($p**2) ? ($n%$p || ($r=-$r)) : return 0 while ($p = next_prime $p) <= $n;
+  $r;
+}
+
+sub moebius_div {
+  my ($n,$p,$r) = (shift,1,1);
+  return -1 if is_prime $n;
+  $n%($p**2) ? ( $n%$p || ($r=-$r,$n/=$p)) : return 0 while ($p = next_prime $p) && $n>1;
+  $r;
+}
+
+sub moebius_div_opt {
+  my ($n,$p,$r) = (shift,1,1);
+  return -1 if is_prime $n;
+  $n%$p || ($n/=$p) && $n%$p ? ( is_prime $n ? (return $r) : ($r=-$r) ) : return 0 while ($p = next_prime $p)**2 <= $n;
+  $r;
+}
+
+sub moebius_exp {
+  my ($n,$p,$r) = (shift,1,1);
+  return -1 if is_prime $n;
+  while( ($p = next_prime $p )**2 <= $n ) {
+    next if $n%$p;
+    $r=-$r;
+    $n/=$p;
+    return 0 unless $n%$p;
+    return -$r if is_prime $n;
+  }
+  $r;
+}
+