The correct syntax for the select command is

A:= select(g-> PermOrder(g)=4, H);

Think of it as having this English translation: "Select g such that PermOrder(g)=4 from H."

It's never allowed for an indexed name such as H[i] to be a left-side operand of ->, i.e., a parameter. That was the proximate cause of your error message.

I posted code for extracting the extrema from plots, so you're not "on your own" as Mac Dude says. The code can be found here:

http://www.mapleprimes.com/questions/206660-Maximum-Point-Of-Plot#answer221288

The commands that are (guaranteed to be) thread safe are listed in the help page ?index,threadsafe. Unfortunately, it's a very small list.

Detecting whether a large number is prime is much, much easier computationally than factoring a product of two medium-sized primes. This fact is the basis of most modern cryptography systems. (There is no known theoretical reason for this fact---it may just be the current state of the known algorithms.)

restart:
p:= nextprime(2^99):
q:= nextprime(p):
n:= q*p:
CodeTools:-Usage(ifactor(n));

memory used=313.25MiB, alloc change=164.03MiB, cpu time=9.28s, real time=9.36s
(633825300114114700748351602943) (633825300114114700748351603131)

p:= nextprime(2^999):
%-2^999;

1239

restart:
p:= 2^999+1239:
CodeTools:-Usage(isprime(p));

memory used=3.33MiB, alloc change=8.00MiB, cpu time=32.00ms, real time=39.00ms
true

Note that the second number has five times as many digits as the first!

So, if n is a product of several small, easy-to-find factors and one large prime, then ifactor knows to stop after finding the small factors because it can quickly detect that the remaining factor is prime.

(A side comment on your timing technique: It's more accurate to include the timing commands and the code being timed in a single execution group.)

Regarding Christopher's Answer: There's no need for you to be forced to accept seq's default limitation to integer values of m. Any number of a values of m in the range 1..2 (or in any range) can be accomodated like this:

linspace:= (ab::range, n::posint)-> [seq(op(1,ab)-`-`(op(ab))/(n-1)*_k, _k= 0..n-1)]:plots:-display(seq(plot(sin(m*x), x= -Pi..Pi), m= linspace(1..2, 9)));

The above plot uses 9 values of m.

Regarding Kitonum's Answer: I wouldn't be so quick to reject the idea of a 3-d plot. If contructed properly and viewed from a certain angle, it'll look like the 2-d plot that you want. I think that the continuous color gradation used in 3-d plots improves this particular plot:

plot3d(
     sin(m*x), m= 1..2, x= -Pi..Pi,
     style= wireframe, thickness= 2,
     orientation= [0,90], axes= normal, labels= [``, x, y]
);

The aboves uses plot3d's default of 49 values of m. That can be controlled with the grid option.

Here are six more easy (one-line) methods to plot the ellipse, all starting with the defining expression

f:= (x,y)-> 3*x^2-3*x*y+6*y^2-6*x+7*y-9;

All of these methods produce fine plots without the need for increasing the default number of points plotted.

1. The command algcurves[plot_real_curve]  produces very accurate plots of bivariate polynomial equations (of genus 0).

algcurves[plot_real_curve](f(x,y), x, y);

Note that this command does all the work. In particular, it computes the appropriate ranges for x and y.

2. The ellipse is the intersection of the surfaces z = f(x,y) and z = 0, so the command intersectplot can be used.

plots:-intersectplot(f(x,y), 0, x= -2..4, y= -2..2, orientation= [270,0], axes= normal);

This command shows just the xy-plane because of the orientation option.

3. The ellipse is a level curve of f(x,y), so a contourplot showing just one contour can be used.

plots:-contourplot(f(x,y), x= -2..4, y= -2..2, contours= [0]);

4. A slight variation of the above is contourplot3d, which is faster due to (I believe) some externally compiled processing.

plots:-contourplot3d(f(x,y), x= -2..4, y= -2..2, contours= [0], orientation= [270,0], axes= normal);

The next two methods first perform some algebraic manipulation of the expression before plotting. Neither of these require knowledge or specification of ranges for x and/or y.

5. Make a conversion to polar coordinates.

plot(solve(subs([x,y]=~ r*~[cos,sin](t), f(x,y)), r)[1], t= -Pi..Pi, coords= polar);

The [1] is to select one of the two solutions for r. It doesn't matter which is used.

6. Use a parametrization in Cartesian coordinates.

plot([eval(algcurves[parametrization](f(x,y), x, y, t), t= tan(u))[], u= -Pi/2..Pi/2]);

The tangent substitution is needed because parametrization returns a parametrization with domain t= -infinity..infinity.

Here are some slight improvements. Since this is a fundamental algorithm (see the Wikipedia article "Exponentiation by squaring"), small improvements can add up to significant savings. In this version, by using the three-argument form of irem I avoid the recomputation of the quotient n/2, and by using thisproc for recursion I avoid the name evaluation of using Puis.

Puis:= proc(X, n::integer)
option remember;
local q,r;
     if n = 0 then 1
     elif n = 1 then X
     elif n < 0 then thisproc(1/X, -n)
     else
          r:= irem(n,2,q);
          `if`(r=1,X,1)*thisproc(X*X, q);
     end if
end proc;

Here's an efficient and non-recursive way to do it that still uses the definition, in a way:

fib:= (n::nonnegint)-> (<1,1;1,0>^n)[1,2]:
fib(101);

     573147844013817084101

Any constant-coefficient linear recurrence can be computed this way---by taking a matrix power.

FWIW, here's my rational roots procedure:

RationalRoots:= proc(f::And(polynom(integer), satisfies(f-> nops(indets(f,name))=1)))
uses NT= numtheory;
local x:= indets(f,name)[], p, q;
     select(
          r-> eval(f, x= r) = 0,      
          {0, seq(seq([p/q, -p/q][], p= NT:-divisors(tcoeff(f))), q= NT:-divisors(lcoeff(f)))}
     )
end proc:

f:= x-> piecewise(x<0, x^2-2*x+1, x>0, 1-x, undefined);

Why do you say "Piecewise function in Maple are continuous function"?

FWIW, here's my somewhat minimalist Jacobi method procedure. I simply used the Wikipedia article "Jacobi method" and modified the algorithm to use Maple's vector-mode operations. I use two convergence criteria: the magnitudes of both the error and the residual need to be less than epsilon.

Jacobi:= proc(A::Matrix, b::Vector, x0::Vector, epsilon::positive, m::posint)  
uses LA= LinearAlgebra;  
local
     k,
     x:= Vector(x0),
     p:= Vector(op(1,x)),
     D:= Matrix(A, shape= diagonal),
     R:= A-D,
     B:= Vector(b)
;  
     D:= D^(-1);  B:= D.B;  R:= D.R;
     for k to m while LA:-Norm(A.x-b) > epsilon or LA:-Norm(x-p) > epsilon do
          p:= x;
          x:= B-R.p      
     end do;
     if k > m then WARNING("Max iterations reached.") end if;    
     x  
end proc:

Example of use:

n:= 3:
A:= LinearAlgebra:-RandomMatrix(n):
#Force A to be strictly diagonally dominant:
A:= A^+.A:
A:= Matrix(op(1,A), (i,j)-> `if`(i=j, n, 1)*A[i,j]):
b:= LinearAlgebra:-RandomVector(n):
x:= Jacobi(evalf(A), evalf(b), Vector(n, fill= 1.), 1e-6, 99);

A.x-b;

Here's how to create a random real symmetric float matrix with determinant 2:

macro(LA= LinearAlgebra):
n:= 9:
A:= LA:-RandomMatrix(n, datatype= float[8], shape= symmetric):
det:= LA:-Determinant(A):
A:= signum(det)*(2/abs(det))^(1/n)*A;

The major problem with your procedure is that f isn't defined. You've made f local; it should be a parameter.

A second problem is that a parameter can't be indexed. So change p[0] to p0.

A third problem is that a parameter can't (usually) be assigned to. You attempt to assign to p[0]. If you want to update it, then you need to copy it to a local variable.

omega:= (y::function(name))-> y*diff(y, op(y))/op(y):

omega(y(x));

(In the future, please post small worksheets inline as well as attaching them.)

You had many errors, mostly syntax. I don't have time right now to detail all the syntax errors, but I corrected them all. You should carefully compare what I wrote with what you wrote. Many of the changes I made were simply my personal preference syntax rather than correction of errors.

The glaring logic error that you made is

seq(EnergySpacingGUE(1,1), i= 1..10)

The value of n must be smaller than N, because if there are N values, then there are only N-1 spacings. I changed this to

seq(EnergySpacingGUE(i,1), i= 2..11)

but I'm only guessing your intention!

I guess that you realize that a histogram of only ten data points is ridiculous, and that this was only a proof-of-concept case.

restart:
UseHardwareFloats:= true:
GUE:= proc(N::posint)
uses LA= LinearAlgebra, RT= RandomTools;
local
     R:= RT:-Generate(distribution(Normal(0,1)), makeproc),
     G:= 'Matrix((N,N), R, datatype= float[8])' $ 2,
     C:= G[1] + I*~G[2],
     H:= (C + C^*)/~2
;
     convert(LA:-Eigenvalues(H)/~sqrt(4.*N), list)[]
end proc:

EnergySpacingGUE:= proc(N::posint, n::posint)
local
     Eigenvsorted:= sort(Re~([GUE(N)])),
     Spacing:= abs~(Eigenvsorted[..-2] -~ Eigenvsorted[2..])
;
     Spacing[n]*N/`+`(Spacing[])
end proc:
      
data:= [seq(EnergySpacingGUE(i,1), i= 2..10)];

plots:-display([
     plot(32/Pi^2*x*exp(-Pi/4*x^2), x= 0..2),
     Statistics:-Histogram(data, binwidth= 0.05, range= 0..2)
]);