I have two comments, unrelated to each other:

1. Regarding debug:

1. The debug command cannot help you find syntax errors that prevent the procedure from even "compiling", which is the case here. I believe that the same is true for any language. To find such errors, pay attention to where your cursor is placed immediately after you get the error message. The error is often (but not always) within the same line.
2. While you are typing code in Maple's GUI, and your cursor is immediately after a bracketing character (parenthesis, square bracket, etc.), a faint "shadow cursor" will appear on the matching bracket. Most errors of unbalanced brackets can be found this way.
3. If you enter code in a Code Edit Region (available from the Insert menu), further syntax problems will be (subtly) highlighted as you type them.
    

2.  Regarding the time complexity of your algorithm:

The isolve command can find your solutions much faster (in general, as P or Q -> infinity) than your brute-force double loop. Your procedure could be:

Fract:= proc(P::posint, Q::posint)  
local p,q;
   ((`/`@op)~)~(  #Convert integer pairs to fractions.
      select(  #Filter solutions to ranges 1..P-1 and 1..Q-1.
         type, 
         #Convert solution pairs to pairs of pairs:   
         map2(eval, [[p,q], [P-p, Q-q]], [isolve((P-p)*q-P*(Q-q) = 1)]),
         [[posint$2]$2]  #a pair of pairs of positive integers
      )
   )[]
end proc;  

Use sum to sum an infinite series, a finite indefinite series, or a finite definite series with a large number terms when it's suspected that  some symbolic tricks will work.

Your sum is 

sum(n/2^k, k= 0..infinity) + 1;

The resut is 2*n+1, not the 2*n-1 that you report. This is true for any n; it needn't be a power of 2.

Note that the 1 in your sum is a separate term; it doesn't fit the pattern of the other terms.

Not many infinite series can be handled by rsolve, but this one is especially simple since a closed form for the partial sums can be found. I'd use

RS:= rsolve({S(0)=1, T(1)= n, T(x)= T(x-1)/2, S(x)= S(x-1)+T(x)}, {T(x), S(x)});

where T is the recurrence for the terms and S is the recurrence for the partial sums. To get the sum of the infinite series, take the limit of S(x) as x approaches infinity:

limit(eval(S(x), RS), x= infinity);

Maple only supports 64-bit hardware floats^[*1]. Software floats can be virtually any length, limited only by availlable memory and the exponent being a 1-word integer.

The following procedure converts between hardware floats and 16-character hex strings. If the argument can be represented as a hardware float but hasn't been, then a conversion to hardware float is performed. When a string is returned as a float, it is returned in a one-element Array. This is to prevent unintentional conversion to software float.

`hf<->hex`:= proc(x::{And(string, 16 &under length), Or(hfloat, realcons)}) 
   option hfloat; 
   `if`(x::string, hfarray(sscanf(x, "%y")), sprintf("%Y", evalhf(x)))
end proc:

#Example of use
Hx:= `hf<->hex`(sin(1));
                    Hx := "3FEAED548F090CEE"
Hf:= `hf<->hex`(Hx);
                   Hf := [0.841470984807897]

A useful excercise would be an additional procedure to extract the mantissa and exponent from the hex string (which'll require subdivision of at least one hexit, perhaps more). If you need, I'll do it. For software floats, this extraction is a trivial application of command op.

^[*1]The term double precision is used colloquially throughout the English-speaking technical world for 64-bit hardware floats, and they are formally described by the IEEE-754 standard. It'd probably be best for you to either avoid or update your technical usage of the adjectives single, double, etc. See Wikipedia IEEE 754.

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.