While reviewing code the other day, I came across the following snippet (here converted to a procedure).

Ds := proc(V::set, n::posint, t)
local i,v;
    {seq(seq((D@@i)(v)(t), i=1..n), v in V)};
end proc:

The purpose of this is to generate a set of derivatives at a point of a set of unassigned names. For example

 Ds({x,y},2,0);
                                     (2)           (2)
                {D(x)(0), D(y)(0), (D   )(x)(0), (D   )(y)(0)}

Despite its simplicity, the particular implementation is not the most efficient.  To see that, do

printlevel := 10: 
Ds({x},2,0);
{--> enter Ds, args = {x}, 2, 0
{--> enter @@, args = D, 1
<-- exit @@ (now in Ds) = D}
{--> enter D, args = x
                                   Opr := D

                                     D(x)

<-- exit D (now in Ds) = D(x)}
{--> enter evalapply, args = D(x), [0]
<-- exit evalapply (now in Ds) = D(x)(0)}
{--> enter @@, args = D, 2
                                      (2)
                                     D

<-- exit @@ (now in Ds) = `@@`(D,2)}
{--> enter evalapply, args = `@@`(D,2), [x]
<-- exit evalapply (now in Ds) = `@@`(D,2)(x)}
{--> enter evalapply, args = `@@`(D,2)(x), [0]
<-- exit evalapply (now in Ds) = `@@`(D,2)(x)(0)}
                                        (2)
                            {D(x)(0), (D   )(x)(0)}

<-- exit Ds (now at top level) = {D(x)(0), `@@`(D,2)(x)(0)}}
                                        (2)
                            {D(x)(0), (D   )(x)(0)}
printlevel := 1:

Both D and `@@` are library procedures. Furthermore, applying a Maple function, D(x), to an argument, 0, to get D(x)(0) invokes the library procedure ?evalapply. All of these cost time and memory. Consider

CodeTools:-Usage(Ds({cat(x,1..100)},2,0)):
memory used=2.94MiB, alloc change=2.62MiB, cpu time=0.08s, real time=0.08s

This might not seem like much, and in the context of the actual code it was negligible, however, there are places where such things matter.

So here's the modest challenge.  Rewrite Ds to be faster and use less memory. I'll post my solution in a day or so.

Followup

It's been a couple days.  Because--as with most toy problems---the process of solving is more useful than then the answer, I'll explain how I arrived at one efficient alternative.

A reasonable first step is to figure out what, precisely, is required. We want a set of symbolic derivatives at a point.  In the original procedure, the derivatives are generated by calling Maple library procedures.  For example, the expression (D@@2)(x)(0) calls the Maple library procedures D, `@@`, and evalapply, returning an expression that represents the derivative.  We want to create that same expression, but without invoking the Maple libraries.  Use ?lprint to determine the final expression:

lprint((D@@2)(x)(0));
                       `@@`(D,2)(x)(0)

Generating that particular expression without calling Maple libraries is trivial, just surround it with forward-quotes.  To check that that works, assign a small procedure, P, that evaluates its argument (declared with the ?uneval modifier) while temporarily assigning  ?printlevel to a suitable value.

P := proc(expr :: uneval)
local val;
    try
        printlevel := 20;
        val := eval(expr);
    finally
        printlevel := 1;
    end try;
    print(val);
    return val;
end proc:

Here we use P to evaluate the expression with and without forward-quotes

d1 := P((D@@2)(x)(0)):
{--> enter @@, args = D, 2
<-- exit @@ (now in P) = D@@2}
{--> enter evalapply, args = D@@2, [x]
<-- exit evalapply (now in P) = (D@@2)(x)}
{--> enter evalapply, args = (D@@2)(x), [0]
<-- exit evalapply (now in P) = ((D@@2)(x))(0)}
                         @@(D, 2)(x)(0)
d2 := P('(D@@2)(x)(0)'):
                         @@(D, 2)(x)(0)

While they appear to be the same, how do we know that the values assigned to d1 and d2 are identical? We could use ?evalb, however, a direct usage is not definitive in that the expressions are fully evaluated, which we can see by using P:

P(evalb(d1=d2)):
{--> enter @@, args = D, 2
<-- exit @@ (now in P) = D@@2}
{--> enter evalapply, args = D@@2, [x]
<-- exit evalapply (now in P) = (D@@2)(x)}
{--> enter evalapply, args = (D@@2)(x), [0]
<-- exit evalapply (now in P) = ((D@@2)(x))(0)}
{--> enter @@, args = D, 2
<-- exit @@ (now in P) = D@@2}
{--> enter evalapply, args = D@@2, [x]
<-- exit evalapply (now in P) = (D@@2)(x)}
{--> enter evalapply, args = (D@@2)(x), [0]
<-- exit evalapply (now in P) = ((D@@2)(x))(0)}
                              true

How can we compare the two without fully evaluating them? One way is to wrap the expression in a call to ?eval with the second argument 1, that limits the evaluation of the expression to one level:

P(evalb(eval(d1,1)=eval(d2,1))):
true

Another method, more hackish but definitive, is to use ?disassemble on 'd1' and 'd2' and note that the first address (22279184 below) in their expansions, a pointer to their value, are identical:

disassemble(addressof('d1'));
                  8, 22279184, 18104160, 12644
disassemble(addressof('d2'));
                  8, 22279184, 18104160, 12900

We now know how to generate a given instance of a symbolic derivative in an efficient manner.   How do we extend this to a sequence of derivatives of symbols.  The ?seq function seems a natural choice.   Here is a first attempt

seq('`@@`(D,i)(x)(0)', i=1..2);
                       D(x)(0), @@(D, 2)(x)(0)

That appeared to work, but let's evaluate it in P to see what really happens

P(seq('`@@`(D,i)(x)(0)', i=1..2)):
{--> enter @@, args = D, 1 ...

Clearly that is no good.  Using seq causes the expression to be evaluated, which then calls the libraries. We need to do better.

At this point if you have no idea how to proceed you might look at how Maple's library procedures handle this.  Clearly they have a way of generating the expression without calling themselves, otherwise there would be an infinite recursion.  There is, in fact, some recursion involved, but it is not infinite. Tracing through the operation of the Maple evalapply procedure is instructive.  You can inspect the procedure with

showstat(evalapply);

evalapply := proc(f, t)
   [...]
       ('f')(op(t))
   [...]
end proc

In the above I left only the statement of interest. It applies a function (f) to arguments (op(t)), but does so without calling the function. Applying forward-quotes around f does the trick.  Here we use that technique to partially accomplish our goal

InertApply := f -> 'f'(_rest):
P(seq(InertApply(('D')(v),0), v in {x,y,z})):
{--> enter InertApply, args = D(x), 0
<-- exit InertApply (now in P) = (D(x))(0)}
{--> enter InertApply, args = D(y), 0
<-- exit InertApply (now in P) = (D(y))(0)}
{--> enter InertApply, args = D(z), 0
<-- exit InertApply (now in P) = (D(z))(0)}
                   D(x)(0), D(y)(0), D(z)(0)

Notice that there are no library calls to D or evalapply.  Here is how this technique can be used to rewrite Ds

Ds1 := proc(V::set,n::posint,t)
local i,v;
    {seq(InertApply('D'(v),t), v in V)
     , seq(seq(InertApply(InertApply('`@@`'('D',i),v),t), v in V), i=2..n)
    }:
end proc:
CodeTools:-Usage(Ds1({cat(x,1..100)},2,0)):
memory used=48.57KiB, alloc change=0 bytes, cpu time=0ns, real time=1000.00us

That is a substantial improvement from the original, which was
CodeTools:-Usage(Ds({cat(x,1..100)},2,0)):
memory used=2.94MiB, alloc change=2.62MiB, cpu time=0.08s, real time=0.08s

Can we do better?  One place that can be improved is to eliminate nested call to InertApply in the code that generates derivatives greater than one. Here we use some specialized procedures:

Dxt  := proc(x,t) 'D(x)(t)' end proc:
Dixt := proc(i,x,t) '`@@`(D,i)(x)(t)' end proc:

Ds7 := proc(V::set,n::posint,t)
local i,v;
    { seq(Dxt(v,t), v in V)
      , seq(seq(Dixt(i,v,t), v in V), i=2..n)
    };
end proc:

That gives a small improvement in speed; you'll need to increase the number of symbols in V to see the difference.

There are other ways to approach this.  One is to use the ?subs procedure to substitute into an expression. Because subs does not evaluate the resulting expression, we can avoid the library calls.  Here is one such approach

Ds4 := proc(V::set,n::posint,t)
local i,v,f,ft;
    ft := f(t);
    {seq(subs(f='D'(v), ft), v in V)
     , seq(subs(f='`@@`'('D',i), [seq(subsop(0=f(v), ft), v in V)])[], i=2..n)
    }:
end proc:

It has about the same performance as Ds7.

Postscript

You might wonder why InertApply (see above) works but

proc(x) local y;  y := 3;  'x(y)' end proc(a);
                                a(y)

That is, why doesn't this return 'a(3)'.  How come the value of the parameter x is substituted under the forward-quotes, but the value of the local variable y is not?  The simple answer, sans explanation, is that parameters behave differently than variables. 


Please Wait...