@Adam Ledger Never mind cleaning up that code or explaining the algorithm: I got it figured out. It was easier to figure out once I saw that the vast majority of your aliases and packages were unused.

Your procedure says, in English, "k is prime if the number of positive divisors of its positive divisors (counted with multiplicity) is 3." Indeed, those three divisors are [1, 1, k], not the [k, 1, 0] that you claimed.

Written in more-normal Maple (yet maintaining your algorithm) your procedure is equivalent to

P:= proc(N::posint)
local k;
uses numtheory;   #For `divisors`
     {seq(k*`if`(nops((op@[op]@divisors)~([divisors(k)[]])) = 3, 1, 0), k= 1..N)} minus {0}
end proc;

(An astute reader may notice that the [op] is superfluous. I only included it because it corresponds to the outer convert(..., list) in the OP's code.)

Of course that English translation is just a bloated way of saying "k is prime if the number of its positive divisors is 2." A simple Maple translation of that is

P:= N-> select(k-> nops(numtheory:-divisors(k))=2, {$1..N});

As far as I know, Maple's numeric PDE solver doesn't support the use of imaginary coefficients. Often, the pdsolve command itself will complete without error, but any practical attempt to use the solution module will result in error. I think that has happened in this case. Certainly, the direct cause of the problem is the imaginary coefficients.

@Adam Ledger Yes, I'd appreciate it if you made some modifications for readability and re-uploaded it. For my taste, the readability would be greatly improved if there was no use of alias or with whatsoever.

Regarding Statistics:-Count: The number of elements of a set or a list can be found with the top-level commands nops or numelems. These have the benefit of not requiring set-to-list conversion.

Regarding efficiency: You can use CodeTools:-Usage to make some measurements. There is no prepackaged prime-subset command that I'm aware of. But it's trivial to build one from the pre-existing commands. The most trivial I can come up with is

AllPrimes:= (n::posint)-> select(isprime, {$1..n}):

This is also surprisingly efficient. Comparing it to yours:

CodeTools:-Usage(`ℙ`(2^14)):   # `ℙ` is your procedure's name in ASCII.
memory used=0.57GiB, alloc change=0 bytes, cpu time=8.14s, real time=8.18s, gc time=390.62ms

CodeTools:-Usage(AllPrimes(2^14)):
memory used=6.74MiB, alloc change=0 bytes, cpu time=63.00ms, real time=52.00ms, gc time=0ns

(s = seconds, ms = milliseconds)

I have a highly polished Sieve of Eratosthenes procedure that's a little bit faster than AllPrimes. But it's on a different computer, and it'll take me a while to get it.

[A while passes.]

EratosthenesSieve:= proc(N::posint)
# Returns list of primes <= N
   local
      # Only record data for the odds >= 3
      M:= iquo(N-1,2)
     ,prime:= rtable(1..M, fill= true)
     ,p, P, k
   ;
   for p to trunc(evalf(sqrt(N)/2)) do
      if prime[p] then
         P:= 2*p+1;
         for k from (P+1)*p by P to M do  prime[k]:= false  end do
      end if
   end do;
   [`if`(N=1,[][],2), seq(`if`(prime[p], 2*p+1, [][]), p= 1..M)]
end proc:
     
CodeTools:-Usage(EratosthenesSieve(2^14)):
memory used=0.55MiB, alloc change=0 bytes, cpu time=16.00ms, real time=8.00ms, gc time=0ns

The Sieve procedure is written in "pure, native" Maple. I may be able to make it faster with evalhf or compilation.

I combined your three posts into one. 

It's an interesting algorithm that I've never seen before. Some very quick and preliminary tests show that the prime-generation part is correct. The counting part is off by one. 

Your code is exhausting to read due to the overuse of alias (especially) and also overuse of with. I think that there are some packages and aliases that you declare but don't use. Please consider rewriting this with the aliases and package names expanded.

List-to-set conversion of L can be done by {L[]}, and set-to-list conversion of S can be done by [S[]]. I find these easier to read than your single-letter operators. (I'm sure that some others would disagree with that.)

@Kitonum Thanks. While you were voting, I added a significant amount to the Reply, so you'll probably want to reread it.

@Kitonum 

piecewise uses assumptions when evaluating conditions with parameters; `if` doesn't. So, if n isn't manifestly an integer, then the n::odd in the `if` evaluates to false, regardless of assumptions. There is a distinction between types (which are checked with type) and properties (which are checked with is). Many types have the same name as an equivalent property; odd, even, and integer are examples of this. The `::` operator can be used to do a type check or a property check. Which is used depends on the context. piecewise is an is context; `if` is a type context.

The expression piecewise(cond, a, b) is almost equivalent to

(()-> `if`(is(cond)::truefalse, `if`(is(cond), a, b), 'procname'()))()

In other words, is forces a decision to true, false, or FAIL; piecewise just remains unevaluated rather than going to FAIL when the condition can't be decided under the current assumptions.

I say almost equivalent because the call to is from piecewise is moderated through a call to the undocumented procedure PiecewiseTools:-Is, and I don't know the exact details of that. In particular, there are many common circumstances where is(cond) evaluates false, but the piecewise remains unevaluated rather than simplifying to b. The code of PiecewiseTools:-Is is fairly easy to read, easier than the code of is.

@tomleslie 

Yes, I agree that the length check for the purposes of showing the excessive-length message should take into account whether the output is already elided. But that's not the way that it's implemented.

Also, the excessive-length message should include "This is not an error. The result has been computed."

@tomleslie 

Give the command

interface(verboseproc= 3);

Then the procedure definitions returned by dsolve won't be elided. Now compare the output of the following two dsolve commands, and you'll see why the length increases:

sol1:= dsolve(dsys3, numeric, initmesh= 8);
sol2:= dsolve(dsys3, numeric, initmesh= 16);

The OP hasn't provided a reason for such an unusually large value of initmesh, and I doubt that there is any good reason. Inceasing maxmesh to a large value often makes sense. I'd only increase initmesh if the problem seems to have trouble "getting started" (as indicated by an error message), and then I'd only increase it gradually.

@roman_pearce 

This is all-the-more reason why one should be able to set kernelopts(numcpus) before a call to Threads. The current situation---where one can only effectively set it immediately after a restart---is unacceptable.

If I simply change my test list to a million Maple floats rather than 32 million integers, then using any number of threads greater than one leads to a substantial increase in real time! It's about 5 times higher when using two threads and nearly 2 times higher when using eight threads. To run this test, simply change the command the creates the input data to

L:= evalf(RandomTools:-Generate(list(integer, 2^20))):

In another test, I tried to defeat the caching by adding elements in a random order. I changed the input to

L:= evalf(RandomTools:-Generate(list(integer, 2^20))):
J:= RandomTools:-Generate(list(integer(range= 1..2^20), 2^20)):
save L, J, "ThreadsTestData.m":

And I changed the test execution (both places!) to

Threads:-Add(L[k], k= J)

@roman_pearce Thank you for the detailed and very useful information, and thank you for running my test on other CPUs. I think that you've mostly explained the phenomena that I observed. So I should probably ignore CPU times for parallel processes on my CPU. Still, the real-time reduction going from four to eight threads on my CPU is very small.

Something else to consider is the sharing of the caches between threads on the same core.

@roman_pearce 

Intel Core i7-3620QM @ 2.4 GHZ running Windows 8.1 64-bit. Intel's website claims that it has four "cores" and eight "threads".

Please let me know if you can achieve substantially different results on any processor.

@Christopher2222 

Thank you for your comment and your interest in the post.

You wrote:

based on performance gains it's hardly worth devoting more than 4 of 5 cpu's to a task, unless for some reason a deadline looms and calculations are required to be completed before a specific time.

Yes, that was my conclusion exactly. I wonder to what extent this is a hardware issue. Perhaps hyperthreading two threads per core is not effective, which is why my "sweet spot" is at four threads---my number of cores.

I wonder if two dual core processors (two physical processors with two independent cores each) would outperform a quad core processor (single processor with four independent cores).

@sand15 This discussion prompted me to make precise timing measurements using different numbers of cores and analyze the results. It took me several days because there was a big hurdle trying to figure out how to adjust the number of cores when using the Threads package. This was covered in my Question "How to limit the number of processors used by Threads?" The result was that the test needed to be run in a different session for each number of cores, which complicated the programming. But, it's done, and the results are now posted as "Diminishing Returns from Parallel Processing: Is it worth using more than four processors with Threads?" My conclusion was no, but perhaps this is a hardware issue. I need to see results from a machine capable of 16 or more threads.

@Milos Bilik You may be thinking of a continuous Markov chain as being "really hard to compute." This is a discrete Markov chain---I was taught how to do the computations by hand in high school.

I'm not totally convinced that this model is the best for your situation. But, it's trivial to compute, so at the very least it's a good start.

By reading the data with Maple, it would be trivial to compute, empirically, the following conditional probabilities: Let P[i,j] be the probability that j consultations are performed the next business day given that i consultations are performed today. Then P is the transition matrix of a Markov chain.