I think that you're unduly concerned due to the "memory used" reported by CodeTools:-Usage. This amount is not cumulative; it includes memory that has been used temporarily and returned to the system to be garbage collected.

Using envelope is much slower than not using it (about 1000 times slower in this case). But you say that you must use it. You can vary the range based on the standard deviation. In this worksheet, I do that. The plots show that

• The growth of CPU time is linear wrt the number of samples.
• The samples' standard deviations have reasonable correspondence with their respective disistributions' standard deviations.

Download envelopeSample.mw

@FrancescoLR By "multiple disk drives", I mean more than one disk drive. A CPU can output data much faster than a disk drive can process it. So, it seems unlikely to me that one could derive any benefit from multi-processing for a job that involves a lot of disk output unless there's at least one disk drive per CPU. I'm not saying that you can't do it (you can), nor am I saying that this has anything to do with your errors (it doesn't); I'm just saying that you won't get any benefit from doing it. The only benefit that matters that I can imagine is a reduction in the overall time for the task. That is the reason that you want to do this, right?

Maple has four packages (that I know of) for multi-processing: Threads, Grid, Grid (Computing Toolbox), and process (deprecated).

The most widely known is I think Threads, which you've been trying to use. In this, the simultaneousy executing tasks share all global variables^[*1], which is a severe restriction in Maple because older Maple library code---still very much in use in the background---uses global variables extensively. Processes that inadvertently overwrite each others' global variables will not work: Sometimes they'll return difficult-to-understand errors, and sometimes they'll just return nonsense results; it's unpredictable. Consequenty, this package can only be used naively for fairly simple jobs; by using sophisticated access control to the global (and  especially lexical^[*1]) variables, it can also be used for a few more jobs of middling complexity. Most Maple commands can't be used with Threads. The tiny fraction that can be are listed at ?index,threadsafe. So, you're never going to be able to do code-generation tasks with Threads.

The lesser-known (curiously and unfortunately) package Grid is much more likely to be useful for most Maple multi-processing applications. Using this package, tasks execute as separate processes, each in its own "kernel" or memory area, so they can't accidentally overwrite each others' globals. There is significant overhead to setup these kernels, so it's often not worth it unless each task takes about 2 seconds or more. And, obviously, using Grid requires much more memory than using Threads.

The package Grid Computing Toolbox is for distributed processing. It is very similar to Grid, but tasks can also be executed on remote computers (that have their own copy of Maple). It needs to be purchased separately as a Maple add-on.

The package process is old and deprecated, meaning that it has been superceded by the others and shouldn't be used anymore.

^[*1]For brevity in this article, I am including lexical variables with the globals. Lexical variables are those that are local to a parent procedure and used in a child procedure. Strictly speaking, the term global only applies to variables that are visible to all procedures.

For your first computation, Maple does give me 1, as shown by the following code:

restart:
p:= 3:
firred:= T^2+2*T+2:
d:= degree(firred):
q:= p^d:
alias(tau = RootOf(firred)):
Expand((tau+1)*(2*tau+2)) mod 3;

The Expand could also be Normal in this example. 

Your error is caused by using 2D Input (not your fault). In 2D Input, (tau+1)(2tau+2) is not interpretted as a product of two binomials. Unfortunately, it's interpretted as a "function" (tau+1) being applied to an argument (2tau+2). The degenerate "function" is a constant function, so the result will be tau+1 regardless of what's in the second pair of parentheses. There are two ways to correct this:

1. Using 2D Input, put a space between the two pairs of parentheses: (tau+1) (2tau+2). In all cases where there's an ambiguity between implied multiplication and function application, a space is needed to force it to be implied multiplication (no matter how stupid or degenerate the "function" may be). For example 2(x+3) needs to be 2 (x+3).
2. In any form of input, you can use * for multiplication. You just need to get out of the habit of using juxtaposition as implied multiplication.

Regarding your second computation, I have done it meticulously with pen and paper, and I get what Maple gets: Rem(x^q, f, x) mod 3 is x, not (tau+1)*x. You'll need to show me the computation that leads you to believe your answer. Indeed, my pen-and-paper computation confirms all 9 entries in this table:

f:= x^2 + tau*x:
interface(rtablesize= 2+q):
<
   <` k ` | x^` k`>,
   <`===` | `=====`>,
   < 
      <seq(1..q)> | 
      Vector(q, k-> collect(Rem(x^k, f, x) mod p, x))
   >
>;

Using anything may "work"---if by "work" you simply mean "returns true"---but that's not achieving the effect that you want, so I wouldn't say that it "works correctly". The coefficients of g aren't simply tau; they are themselves polynomials in tau with integer coefficients. So, what you want is

type(g,  polynom(polynom(integer, tau), x));

I'm not sure exactly what you mean, but I think that what you want is given by the command simplex[ratio]. Certainly, this command returns a list of multipliers. See ?simplex and ?simplex,ratio.

The command select can be used to link your two pieces of code. The general form is select(P, L), where P is a boolean procedure (i.e., one that returns true or false) acting on the elements of a container L (such as a list or set). 

I wouldn't use your code above, because both pieces are extremely inefficient (although they do perform the computation correctly, which is the most important thing, and they're both fairly clean stylistically), but if I was going to use it with select, I'd do

select(p-> findOrderOf(2,p) = p-1 and findOrderOf(3,p) = p-1, primes),

which can be improved with andmap to 

select(p-> andmap(g-> findOrderOf(g,p) = p-1, [2,3]), primes).

Your order-finding procedure is inefficient for a mathematical reason: It doesn't make use of the fact that the order must be a divisor of the totient^[*1] of the modulus, which means that there are only a few exponents that actually need to be checked.

Your prime-listing code is inefficient due to an idiosyncracy of Maple rather than a mathematical reason: Lists (and sequences and sets) are effectively read-only; they cannot be added to (or even changed in any way). Any attempt to do it actually creates a fresh list, copying the parts of the original that are needed, no matter how many elements need to be copied. That's done in the background, transparently to the user, who often doesn't realize that something rather inefficient is occuring. Some other Maple containers such as Vectors and tables do not have this read-only behavior.^[*2]

So here's some efficient code: If you're simply interested in answering the given homework problem, you could do

select(
   p-> isprime(p) and NumberTheory:-MultiplicativeOrder~({2,3}, p) = {p-1}, 
   [seq(2001..2999, 2)] #only need to check odds
);  

If you'd rather not rely on Maple's prewritten MultiplicativeOrder, you could do

IsGenerator:= (x::integer, p::prime)->
   not ormap((d,g)-> g &^ d mod p in {0,1}, (p-1) /~ op~(1, ifactors(p-1)[2]), irem(x,p))
:
select(p-> isprime(p) and andmap(g-> IsGenerator(g,p), [2,3]), [seq(2001..2999, 2)]);

My procedure IsGenerator makes use of a special efficiency that applies only when the modulus p is prime: that the totient in that case is simply p-1. There's also an efficiency (regardless of whether the modulus is prime) that results from if the order is not the totient, then we know that x is not a generator without needing to actually compute its order.

The special operator &^ (as in g &^ d mod p) computes g^d mod p without first computing g^d. Thus, it can be used efficiently even when the exponent is extremely large. (The operator itself is inert: g &^ d by itself, without mod, does nothing.)

^[*1]As you probably know, the totient of a positive integer n is the count of positive integers x <= n such that gcd(x,n) = 1, i.e.,  x and n are relatively prime. Thus, the totient of n is the order of the multiplicative group mod n. The totient is also called "Euler's phi function".

^[*2]In Maple documentation, structures that are not read-only (in the way that I described above) are called mutable, and the others are called immutable (or perhaps non-mutable). I think that the term read-only is more widely known and the concept more widely understood.

There are two problems with the worksheet that you posted, one typographical and one logical. For the typographical problem, look closely at the line that begins consd1:=. Notice how that line is printed in a different way than the line immediately above it? I think that that means that that line is considered textual commentary, and so it isn't being executed. I suggest that you simply retype that line. You can ensure that it gets placed in an execution group by using Ctrl-J (while your cursor is on that line) and retyping it at the prompt that appears.

The logical problem is that the variables x[2,1], x[3,1], x[4,1], and x[5,1] do not have any constraints. Hence they can take any negative values, making the problem unbounded towards negative infinity.

Your line

cons[k] := add(x[i, k], i = N)+add(x[k, j], N) = 2

should be

cons[k] := add(x[i, k], i = N)+add(x[k, j], i= N) = 2​​​​​​​

As a workaround, the original 33x32 system can be easily and quickly solved by the direct use of the LinearAlgebra package:

macro(LA= LinearAlgebra):
(A,B):= LA:-GenerateMatrix(liseq1, lisvar1):
(P,L,U):= LA:-LUDecomposition(A):
Sol:= lisvar1 =~ seq(LA:-BackwardSubstitute(U, LA:-ForwardSubstitute(L, P^+.B))):
andmap(evalb, eval(liseq1, Sol)); #verify correctness

A small adjustment to Kitonum's second method, the one that uses subsindets, will make it safer (more robust) in that the conversion of  `*` to `.` will only occur when intended:

subsindets['flat'](A.B, 'specop'(rtable, `*`), `.`@op);

If your Vectors have purely numeric entries, this won't make any difference; but if they themsleves contain products, it might.

Maple's assignment operator is two ordinary ASCII characters: a colon (:) followed by an equals (=). The character that you have after X is some exotic non-ASCII character that looks like :=, but isn't.

If you make this small change, then your example will work as expected.

Subscripting is not quite the same thing as indexing, and the underscore _ is not an indexing operator in Maple (nor is it a subscripting operator, though double underscore __ is). 

There are many ways to do what you want:

Most-straightforward way:

seq([X[k], Y[k]], k= 1..nops(Y));

Clever, minimal-code, way:

`[]`~(X,Y)[]:

Perhaps the most efficient (although I haven't checked):

convert(<<X> | <Y>>, list, 'nested')[];

Most versatile: Just make it into a two-column Matrix:

XY:= <<X> | <Y>>:

This last form is accepted by most Statistics commands and plot; indeed, they are more efficient when the input data is presented thus.

As you can see, of the four methods, only one uses an index. I tend to prefer the clean syntax of these index-free methods over those that require introducing a bound variable for the index; however, these latter methods are sometimes a little more efficient.

One way is to use Rem with mod. Here's an example over GF(3,2):

restart:
alias(alpha= RootOf(Randprime(2,z) mod 3)):
f:= x*(x-1):
Rem(alpha^5, f, x) mod 3;

                          2*alpha

In addition, you could redefine your ring operations `+`, `*`, `-` to perform this Rem operation; or you could give them other names, which would still be simpler than frequently typing Rem(...., f, x) mod 3. If the ring has only a few hundred elements or fewer, remember tables for the operations are very effective for high-speed computations. 

My preference is 

n:= 10:
cat(seq("TH"[x], x= rtable(1..n, random(1..2))));

I realize that this probably looks weird to most readers. I like it because I have a suspicion (not confirmed yet) that rtable with a random or frandom initializer is very efficient, more efficient than a sequence of calls to a procedure produced by rand.