Use

simplify(evalc([(()->"Re")=Re, (()->"Im")=Im](Ec))) assuming positive;

Hint: Substitute epsilon= 1/upsilon, then simplify, then consider "by hand" the limit as upsilon approaches +infinity. I think that this is easier to see than the limit as epsilon approaches 0 from the right. Likewise, limit(..., epsilon= 0, left) is equivalent to limit(..., upsilon= -infinity).

And, of course, to make your problem more understandable to Maple, you must replace [ ] with ( ) and 0.5 with 1/2 (as VV said).

Here is an extendable template that should let you answer (in a natural way if you know statistical language) any such question for any discrete distribution with integer support. All results come directly from applying simplify (with assuming) and sum to the pmfs of the distribution. Each distribution is represented by a Maple container structure called a Record, which contains fields for the pmf, support, and assumptions on the parameters.

All parameters can be specified as any algebraic expressions (including constants); they need not be simple names. The `if` expressions in the Records adjust the assumptions appropriately.

The table (actually a DataFrame) at the end actually looks like a table when viewed in a worksheet. The MaplePrimes renderer squashes it here.
 

Download DiscreteDist.mw

Joe has pointed out, essentially, that's it's often erroneous and usually suspicious to refer to a for loop's index variable outside the loop. The practice does have a few valid uses, such as determining whether the loop ended due to exceeding the maximum index or due to some other reason.

You have another major error: "Product of the eigenvalues" does not mean the dot product of the vector of eigenvalues with itself; rather, it means the (scalar) product of the three eigenvalues. You can do this with command mul, like this:

mul(N[i])

This is something like what your professor has in mind:

with(LinearAlgebra):
N:= 8:  n:= 3:  r:= true:
to N do
   M:= RandomMatrix(n,n);
   if Determinant(M) <> simplify(expand(mul(Eigenvalues(M)))) then 
      r:= false;
      break
   end if
end do: 
r;

I prefer this:

macro(LA=LinearAlgebra):
N:= 8:  n:= 3:  r:= true:
to N do
   if (radnormal@(LA:-Determinant - mul@LA:-Eigenvalues))(LA:-RandomMatrix((n,n))) <> 0 then
      r:= false;
      break
   fi
od:
r; 

I prefer it primarily because I trust radnormal when applied to a zero-equivalent radical expression more than I trust simplify@expand when applied to an integer-equivalent radical expression.

If F is list (or set or sequence) of functions (or procedures) and S is a sequence of arguments that would work for any of them, then simply F(S) is what you want. Note that this uses neither map nor an elementwise operator with ~. Example:

[iquo, irem](7,5)

or

(iquo, irem)(7,5)

I often use this as (min, max)(S). So, the number of arguments doesn't matter.  In particular, it could be one argument.

If X and Y are lists of the same length, where X are arguments for the first position and Y for the second, then do

F~(X, Y)

Unlike zip, this idea can be extended to any number of arguments. Any of the argument lists can be replaced by a single element, as long as that element is not itself a list, set, table, or rtable. The take-away lesson from all this is that there's no syntactic difference between a single function (or procedure) and a list, set, or sequence of them (an explicit sequence needs to be in parentheses).

Your code

map(eval~, [f(x), g(x)], x=~ [p,q,t])

can be (and should be) replaced by 

map(map, [f, g], [p, q, t])

because it's not necessarily possible to evaluate every function at symbolic x, and even if it's possible, it's not necessarily efficient.

The following procedure detects and reports cycles of exact algebraic numbers:

CycleDetectAlgnum:= proc(f, x0::radalgnum, Max::posint:= 999)
local x:= evala(Normal(x0)), R_inds:= table([x=0]), R_vals:= table([0=x]), k, j; 
   for k to Max do
      x:= evala(Normal(f(x)));
      if assigned(R_inds[x]) then return seq(R_vals[j], j= R_inds[x]..k-1) fi; 
      R_inds[x]:= k;  R_vals[k]:= x
   od;
   return
end proc
:   

It's usage is exactly like the other's:

CycleDetectAlgnum(y-> 4*y*(1-y), (5-sqrt(5))/8);

If you're only interested in detecting and reporting cycles rather than analyzing algebraically why they happen, then the following code---which does it all in hardware complex floats---will be much faster than either iterating the map or using a closed form from rsolve.

CycleDetectInner:= proc(
   C::Array(datatype= complex[8], order= C_order),
   k::integer[4],
   x::complex[8],
   eps::float[8]
)::integer[4];
option autocompile;
local j::integer[4];
   for j to k do
      if abs((C[j]-x)/`if`(abs(x) <= eps, 1, x)) <= eps then return j fi
   od;
   0
end proc
:
CycleDetect:= proc(
   f::procedure, 
   x0::complexcons, 
   Max::posint:= 999, 
   eps::{float,positive}:= evalhf(DBL_EPSILON)
)
option hfloat;
local x:= evalhf(x0), C:= Array(1..1, [x], datatype= complex[8], order= C_order), k, p; 
   for k to Max do
      x:= f(x);
      p:= CycleDetectInner(C, k, x, eps);
      if p > 0 then return C[p..] fi;
      C(k+1):= x
   od;
   return
end proc
:   

Example usage:

CycleDetect(y-> 4*y*(1-y), (5-sqrt(5))/8);

You may want to apply Re~ to the results to remove the zero imaginary parts. Any returned result is a 1-D Array containing a detected cycle. A return of NULL means that Max iterations were performed without a repeated value. 

The reason that add doesn't work as expected is that Tolerances doesn't overload it, whereas it does overload `+`. So, if you replace add with (`+`@op), you'll get your expected results.

Another option is post-processing of the partially evaluated results returned by add:

subsindets(Z, :-`+`, `+`);

The second `+` in that command is properly Tolerances:-`+`, but under the auspices of with(Tolerances), simply `+`  is sufficient.

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))
   >
>;