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

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

No, it looks like the diagonal entries cannot be replaced by anything but zero while still maintaining the option shape=antisymmetric. > generators[1,1]:=Matrix(4,4); Error, attempt to assign non-zero to antisymmetric diagonal This seems wrong to me. Why couldn't the restriction on the diagonal elements instead be that a[i,j]=-a[i,j] ? That way, having rtables on the diagonal might be possible. Well probably still not, since evalb(a=-a) doesn't return true for zero-matrix `a`. But it would have been nice if you could have done, directly, 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=antisymmetric); I don't see a `index/antisymmetric` or `index/skewsymmetric` routine, so the antisymmetric indexing function looks built-in. You could try creating your own Antisymmetric indexing function with something more sophisticated for handling the diagonal's qualities. acer

I second the idea of a longer (or expandable) list of "Recent comments". An exapandable "Active forum topics" would also be nice. And I don't mean like its current "more", which immediately leads to too many choices to scan. I mean just like the sidebar, but longer. A problem with the MaplePrimes Forums's pages is that it's only possible to see the whole thread, and it's not possible to get to just the last entry posted in each. So once an entry gets pushed off the main page and sidebars, it's difficult to access easily. Also, could someone please improve the tracking pages of each individual. Instead of just showing a single entry to each forum topic, perhaps it could list each entry separately. That way, one could get to the users posts directly. acer

Trying to be tricky. unprotect(LinearAlgebra:-Multiply); Z:=copy(LinearAlgebra:-Multiply): LinearAlgebra:-Multiply:=proc() if args[1]=`#mo("∇")` or args[1]=`#mo("∇")` then args[1](args[2]); else Z(args); end if; end proc: `#mo("∇")`:=proc(V::Vector) map(sin,V); end proc: `#mo("∇")`:=`#mo("∇")`: Now insert the Del symbol from the palette, and toggle it as Atomic Identifier, and it can be used without the brackets. That is to say, I was able to do 2D things like, V V V I didn't check to see whether I'd broken LinearAlgebra:-Multiply. I could still do, V V There may be some simple way to customize it, which someone else can suggest. It might also be possible to use instead a module for the redefined multiplication, but I didn't try. acer