The (maple kernel) built-in routine trunc() may be faster than floor() for you, where you know that all your values will be greater than zero. But the improvement due to using trunc() instead of floor() may be small. There may also be some other, faster, ways to approach the whole problem. But these too might have their own bottlenecks to figure out. Here's an example idea, which I haven't tested. My guess is that done well it might perform faster, even if a first crack at trying it doesn't seem better. Suppose that, right up front and just the once, you created a very long 1-D Array or Vector from the distribution. Say with g * 10 entries. This is by calling Sample just once, like, A := Statistics:-Sample(Statistics:-RandomVariable((':-Poisson')(1)), g*10); Then you could create some new, inner, proc which can walk A. As it walks A, it keeps a counter. It take a value from A, increments the counter by 1. Take the floor of that, and let that number be assigned to C. Then walk and add up the next C entries, and take the floor and assign to y[i] or do what you will with it. Increment the counter by C. And repeat. The point is that, when calling the routine S returned by Sample(), a new Vector for the results is called each time. These Vectors make Maple garbage, and garbage collection is part of what's taking up time. Such time might be avoided by creating just a single Vector and incurring the Vector creation overhead just once. You can choose the length of Array A so that it is very, very likely to contain enough values from the distribution for your full purpose. You might be able to use add() to add up the C entries from current-counter to current-counter+C. Or you might be able to use for-loops and perhaps even use Maple's Compiler on this new inner procedure. acer

For a Matrix whose entries are all explicitly stored, one practical limit is the amount of available memory. If the amount of available memory is the effective limit, then using a hardware datatype (integer[1], float[8], etc) can allow one to create and use larger Matrices and Vectors. If, on the other hand, not all entries must be stored, then a hard limit in each dimension may be something like 2^(kernelopts(wordsize)-1)-1. This can be relevant when the Matrix is set up to have either so-called empty storage with unstored entries specified only by an indexing function, or sparse storage with a only a specific number of entries' storage allocated up front. For example, on 64bit Linux, M:=Matrix(2^63-1,2^63-1,storage=sparse[1]); and on 32bit Windows, M:=Matrix(2^31-1,2^31-1,storage=sparse[1]); One way to conceptually extend even that limitation could be to use multi-dimensional Arrays instead of Matrices or Vectors. An Array can have from 0 to 63 dimensions (as opposed to Matrices which have just 2 and Vectors which have just 1). Another possible extension might be to have Matrices whose entries are themselves Matrices. Eg, M:=Matrix(2^31-1,2^31-1,storage=sparse[1]); M[1,1]:=Matrix(2^31-1,2^31-1,storage=sparse[1]); acer

Suppose that the variable `sol` were assigned something like this, sol:= fsolve({Eq16, Eq15}, {beta = 4 .. 5, N_bar = 0 .. 1000}); Then, 2-argument eval() should do the trick, to get at the separate values. sol := {N_bar = 18.04394737 + 0.*I, beta = 4.474743382}; eval(N_bar,sol); eval(beta,sol); acer

Suppose that the variable `sol` were assigned something like this, sol:= fsolve({Eq16, Eq15}, {beta = 4 .. 5, N_bar = 0 .. 1000}); Then, 2-argument eval() should do the trick, to get at the separate values. sol := {N_bar = 18.04394737 + 0.*I, beta = 4.474743382}; eval(N_bar,sol); eval(beta,sol); acer

Don't call Sample() and RandomVariable() each time through the loop. Also, when forming x[i], since you're apparently only taking a single (1) value from the distribution then you don't need to call the heavyweight convert(..,`+`) on it. FFFF_2 := proc(g) local i, x, y, S, Samp; Samp:=Statistics:-Sample(Statistics:-RandomVariable(':-Poisson'(1))); for i to g do x[i] := floor(Samp(1)[1]); if 0 < x[i] then y[i] := floor(convert(Samp(x[i]),`+`)); else y[i] := 0 end if end do; S := [seq(y[i], i = 1 .. g)]; end proc; Also, why use local x as a table, when a scalar x would suffice (since you don't use the x[i] for anything else? If that's the case, then replace x[i] with simply x, and save some more. Lastly, don't load the Statistics package outside the proc definition. It's not good programming style. Either explicitly write out the commands, in the proc, as I've done or utilize the `use` functionality. acer

Can LPSolve be used for the entire problem? One might look at the problem as one of finding merely a feasible point. A linear programming problem with only equality constraints. Using a constant objective should allow it to return the first feasible point found. Optimization:-LPSolve(Vector(4),[NoUserValue,NoUserValue,Matrix([[1,2,7,5], [1,1,2,1]]),Vector(2,[3,1])],assume=nonnegative); If so, then could this help with some of the potential numerical difficulties? The help-page ?Optimization,LPSolveMatrixForm shows that this routine has extra options such as feasibilitytolerance and integertolerance. acer

Can LPSolve be used for the entire problem? One might look at the problem as one of finding merely a feasible point. A linear programming problem with only equality constraints. Using a constant objective should allow it to return the first feasible point found. Optimization:-LPSolve(Vector(4),[NoUserValue,NoUserValue,Matrix([[1,2,7,5], [1,1,2,1]]),Vector(2,[3,1])],assume=nonnegative); If so, then could this help with some of the potential numerical difficulties? The help-page ?Optimization,LPSolveMatrixForm shows that this routine has extra options such as feasibilitytolerance and integertolerance. acer

I get amusement from taking these sorts of examples as linear programming problems with a constant objective. The goal is then to find a feasible point. eq1:= q + d = 8: eq2:= 0.25*q + 0.10*d = 1.25: Optimization[LPSolve](1,{eq1,eq2},assume=nonnegint); More seriously, I don't think that one should shy away from suggesting use of isolve() for such examples. The fact that only whole coins may be possible is a genuine part of the problem. Consider cases where more than one answer is possible. For them, the result from isolve({eq2}) is better than that from solve({eq2}). acer

There's no object "sqrt(b^2)". Instead, the call sqrt(b^2) passes b^2 to a somewhat smart routine. The output is the object (b^2)^(1/2). So this is the sort of object at hand, this (b^2)^(1/2) thing. It's a structure, or an expression. It's a power thing. It's not a function, or an unevaluated function call. If you evaluate it at a specific b, then doing so won't do anything clever. It won't call sqrt() again. On the other hand, sqrt() is a function. And it is somewhat clever. So sqrt(2^2) produces 2 right away. Notice that (b^2)^(1/2) and sqrt(b^2) are not the same thing. The first is a structure, or expression. The second is a function call, which happens to return the first (under no assumptions on b). The sqrt() function is not the only mechanism that can construct (b^2)^(1/2). It could be constructed directly, typing it in just as it is. (That's a partial rationale for why evaluating the expression at a particular b doesn't call sqrt() again.) There's nothing clever in (b^2)^(1/2), that can take advantage of b>0 say. The cleverness resided in sqrt(), the function. Evaluating (b^2)^(1/2) at b= will not call sqrt(^2) once more. Once sqrt() returns its result, the opportunity for sqrt() to be clever has gone. In fact, you don't need to create the "special procedure " that you mentioned, because sqrt() is such a thing already. The issue that you are seeing is also present with your own procedure, a:=proc(b); return sqrt(b^2); end proc; To see that, just try, a(b); eval(%,b=2); Just like for sqrt(), once your a() gets called it returns an expression. Evaluating the expression doesn't cause anything clever to re-evaluate or simplify it. That's why a suggestion to hit it with a big hammer like simplify() was made. Also, there is no automatic simplification of ^(1/2) for this positive sort of . Depending on whom you ask, you may get different opinions about how much merit there is in that. A slightly smaller hammer is normal(). That is, 4^(1/2); normal(%); acer

There's no object "sqrt(b^2)". Instead, the call sqrt(b^2) passes b^2 to a somewhat smart routine. The output is the object (b^2)^(1/2). So this is the sort of object at hand, this (b^2)^(1/2) thing. It's a structure, or an expression. It's a power thing. It's not a function, or an unevaluated function call. If you evaluate it at a specific b, then doing so won't do anything clever. It won't call sqrt() again. On the other hand, sqrt() is a function. And it is somewhat clever. So sqrt(2^2) produces 2 right away. Notice that (b^2)^(1/2) and sqrt(b^2) are not the same thing. The first is a structure, or expression. The second is a function call, which happens to return the first (under no assumptions on b). The sqrt() function is not the only mechanism that can construct (b^2)^(1/2). It could be constructed directly, typing it in just as it is. (That's a partial rationale for why evaluating the expression at a particular b doesn't call sqrt() again.) There's nothing clever in (b^2)^(1/2), that can take advantage of b>0 say. The cleverness resided in sqrt(), the function. Evaluating (b^2)^(1/2) at b= will not call sqrt(^2) once more. Once sqrt() returns its result, the opportunity for sqrt() to be clever has gone. In fact, you don't need to create the "special procedure " that you mentioned, because sqrt() is such a thing already. The issue that you are seeing is also present with your own procedure, a:=proc(b); return sqrt(b^2); end proc; To see that, just try, a(b); eval(%,b=2); Just like for sqrt(), once your a() gets called it returns an expression. Evaluating the expression doesn't cause anything clever to re-evaluate or simplify it. That's why a suggestion to hit it with a big hammer like simplify() was made. Also, there is no automatic simplification of ^(1/2) for this positive sort of . Depending on whom you ask, you may get different opinions about how much merit there is in that. A slightly smaller hammer is normal(). That is, 4^(1/2); normal(%); acer

It doesn't handle series very well. It misses the variable name(s). acer

To temporarily write displayed output to a file, t := [2,3]; writeto("myfile"); dismantle(t); writeto(terminal); The ToInert() and FromInert() routines allow one to produce a structured breakdown of an expression and subsequently reconstruct it, with programmatic manipulation between the two. The help-page ?dismantle could usefully get a cross-reference to ToInert in its "See Also" section at least. acer

I forgot to mention something that I was thinking about yesterday. The entries on the main diagonal don't have to be either zero or a zero Matrix. They could themselves be antisymmetric Matrices which do not happen to be all zero. Some sort of recursive type-check might work, to establish or allow the main diagonal entries to be antisymmetric Matrices (whose main diagonal entries in turn might be scalar zeros or antisymmetric Matrices... and so on). acer

Thanks for the encouragement, John. How's this?

`index/AntiSym` := proc(inds::[posint, posint], comps::Matrix, entry::list)
local row, col, sign;
    row, col, sign := inds[], 1;
    if nargs = 2 then
        return sign*comps[row, col];
    else
        if row = col then
            if type(procname, 'indexed') and type(op(1, procname), posint) then
               comps[row, col] := Matrix(op(1, procname));
            else
               comps[row, col] := 0;
            end if;
        else
           comps[row, col] := sign*op(entry);
        end if;
    end if;
end proc:

metric := Matrix(4,4,Vector([-1,1,1,1]),shape=diagonal):

generators := Matrix(4,4,
 (a,b) -> `if`(a=b,Matrix(4,4),
               metric . Matrix(4,4,
     (c,d) -> metric[a,c]*metric[b,d] -
              metric[b,c]*metric[a,d])),
 shape=AntiSym[4]);

I don't quite understand why Matrix(generators) doesn't automatically flatten it out. But here it is, flattened and with no indexing function.

Matrix(16,16,[seq([seq(generators[i,j],j=1..4)],i=1..4)]);

acer

I did something similar, wrapping an unapplied `ratio_f0` inside a proc. I then plotted x->evalf[10](ratio_f0(x)) . The purpose of this was to render irrelevant any evalhf use by the plotting routing. The resulting plot, which took a long time to produce, had a jaggedness, more pronounced at smaller f0 value, indicative of the expected roundoff error. acer