Could you guide me to the menu on which it's greyed out?

@tomleslie The inert commands Normal and Eval become active when used with mod and a numeric modulus. There are about 45 such commands predefined, and user-defined ones are easy to add. See ?mod for the 45.

In order for the OP's calculation to make any sense, q must be defined as a RootOf an irreducible mod p integer-coefficient polynomial for a prime p that's explicitly given. Otherwise, the whole expression remains inert, as seen in your code above.

As you probably know, and as can be seen above, cancellation of identical factors from a numerator and denominator is an "automatic simplification" that occurs even in inert expressions (unless the arithmetic operators are also made inert, such as %/).

@Christian Wolinski 

You're correct, and the precise rules regarding this sort of thing are on the help page
?operators,precedence. Given any two binary infix operators op1 and op2 (possibly the same) and the expression

a op1 b op2 c

exactly one of the following is always true:

1. the expression evaluates as (a op1 b) op2 c;
2. the expression evaluates as a op1 (b op2 c);
3. it's a syntax error.

Here is EAdd as displayed by showstat(EAdd), except that the lines in red were added by me to show what's actually happening (not what should be happening). In other words, your code acts as if my added lines were already there. 

I used showstat because the indentation style makes it easy to see that there's a NULL logic path. You should adopt a consistent indentation style for your own code writing.

EAdd := proc(f, x, p, P1, P2, Q1, Q2)
local a, mm, R1, R2;
   1   a := coeff(f,x,1);
   2   if P1 = infinity then
   3       return [Q1, Q2]
       elif Q1 = infinity then
   4       return [P1, P2]
       elif Q1 = P1 mod p then
   5       if Q2 = -P2 mod p then
   6           return [infinity, infinity]
           elif Q2 = P2 then
   7           mm := 1/2*(3*P1^2+a)/P2 mod p;
   8           R1 := Normal(mm^2-P1-Q1) mod p;
   9           R2 := Normal(mm*(P1-R1)-P2) mod p;
  10           return [R1, R2] mod p
           else
               return NULL
           end if
       else
  11       mm := (Q2-P2)/(Q1-P1) mod p;
  12       R1 := Normal(mm^2-P1-Q1) mod p;
  13       R2 := Normal(mm*(P1-R1)-P2) mod p;
  14       return [R1, R2] mod p
       end if
end proc

When NULL is returned and you try to index it, you get the error invalid subscript selector. Note that all procedures have a return value, possibly NULL, regardless of whether a return value is explicitly specified.

@mmcdara I've always been fascinated visually by the hurky-jerky way a container half full of water slides along a floor when kicked. A simulation of that would be fascinating also. If the box on Rouben's conveyor belt were half full of water, it'd be even better. 

@epostma Vote up.

Yes, I'd consider two-asterisk-flagged anomalies that disappear upon rerunning with a different randomize seed to have "passed" the test. Indeed, my criteria for the asterisk counts were

*     1/20   <  significance level (alpha) < 1/200  (alpha = 1/20 is widely used historical value)
**    1/200  <                   "         < 1/2000
***   1/2000 <                   "         < 1/2e4

etc. Since my test run posted as an Answer shows 336 p-values, under the null hypothesis of "true randomness" we'd expect about 17 one-asterisk flags and about 2 two-asterisk flags. So, your new sampling code is right on target.

@Kitonum It's a good point in general. But in this case we're dealing with polynomials. 

@666 jvbasha 

eval(diff(f, x$2), p union {x=0});

@Shah Fahd As far as I can tell, there's no way to use ifactor "off-the-shelf" for multiprocessing the factorization of a single large integer. However, I think it'd be possible to pull the Maple code out from the ifactor subprocedures and refactor it for multiprocessing.

I'd be interested in coding Number Field Sieve in Maple if you can provide detailed exposition (with or without pseudocode). From my reading just now, I understand that there are numerous heuristic choices to be made at various steps in the algorithm.

I don't recall ever seeing an error message that gave the exact position and exact library entry where the error occurred. That's very useful. Does anyone know if that information is generally available (perhaps from a debugopts command)?

@Yaongyaong 

Did you not understand the part of my Answer under the plot where I computed the surface area? I'd be happy to explain it in more detail. Do you know cylindrical coordinates: x = r*cos(theta), y = r*sin(theta), z = z?

@Carl Love I got your integral to work in Maple Companion by entering

int ln(x)/x*BesselJ(1,x)dx

where

int is the indefinite integral from the functions keyboard;
ln( is ln from the functions keyboard (the ( comes with it);
x is the variable from the main keyboard (not the alphabetic keyboard);
) is from the main keyboard;
/ is division from the main keyboard;
x ditto;
* is multiplication from the main keyboard;
BesselJ is (carefully) typed from the alphabetic keyboard;
( is from the main keyboard;
1 ditto;
, is from the main keyboard, but you need to hold down the . and swipe up and to the right to get it (the gesture required a little practice to get right, but I've mastered it now);
) is from the main keyboard;
dx comes "for free" when you enter the integral sign, and it's probably best if you ignore it.

The computation is then activated by the blue-highlighted check mark from the main keyboard. The result is
 

Download MC_BesselInt.mw

That result can be substantially simplified, as I showed in previous Replies above.

@Karishma I got BesselJ to work now. I don't know what I did wrong before. Thanks for asking, though.

@Mac Dude This happened for me also. It's frustrating when you've composed an Answer, and, poof, it disappears into the ether. No way to retrieve it.

@HS A warning is not an error, and you say that it "does not seem to work" while I see it working just fine despite the warning.

The reason for the warning is that you have a procedure which is changing the status of an undeclared global name, q. This is usually something worth warning about. A better way to code it is to pass in the name that will be aliased:

restart:
`mod/field`:= proc(k::And(posint, Not(1)), q::name, p::prime)    
local z;
   alias(q= RootOf(Randprime(k,z) mod p, z));
   q
end proc:

field(2, ':-q') mod 13;
field(3, ':-q') mod 13;

The purpose of the ':-q' rather than just q is so that any previous value or alias that q may have will be ignored, and it'll be passed as just an unassigned name.

Now---unlike your original---you can use this procedure to create multiple extensions with different names:

field(2, ':-q1') mod 13;
field(3, ':-q2') mod 13;

You can also return the irreducible polynomial if you insist; however, I find that redundant because that polynomial can always be accessed as op(q), op(q1), etc.