@vv TIL (an Internet-forum abbreviation for Today I Learned) that numboccur is significantly faster than add for counting 1 bits (a.k.a. the Hamming weight). So, you might as well go back to numboccur.

@2cUniverse If you're using Maple 2021 or later and 1D input, and you change $b+1 to $1..b+1, then I think that it'll work for you. You must use 1D input. There's no way that I'm going to modify that code to work in 2D Input.

To get a gray code box in MaplePrimes, click the symbol < > on the toolbar; it's the 4th icon from the right end of the 2nd toolbar. Then copy and paste your code into the box (which is initially white), and click OK. Ignore any awkward spacing of your code in the white box; it'll correct when you click OK and it turns into the grey box. Although you can edit in the white box, I recommend that you don't; it's much easier to edit in the grey box. You need to put something in the white though or it'll disappear when you click OK. Note that the colors and character styles that you see in my gray boxes are the result of my hand editing within the gray boxes. 

You mentioned the efficiency of log2 and Bits:-Split. Note that the efficient command is ilog2, not log2, which is essentially the same as ln(...)/ln(2) and is thus not efficient.

There's no good reason to ever use evalm. It's a very old command that was replaced over 20 years ago. If you want to multiply Matrices A, B, and C (for example), all that you need is

MM:= A.B.C;

I don't know if this will completely solve your problem, but, regardless, you should still eliminate all uses of evalm, because they're certainly not doing anything useful.

@sursumCorda Using the Maple code given for A160414  is essentially circular because that code essentially computes the sequence from this Question, multiplies it by 4, and adds 1. This kind of coding leads, over the years, to inefficient libraries full of crufty redundant code (which nonetheless may produce accurate results).

@vv Okay, you're right that it's better here to avoid recursion and a remember table and instead use seq[scan= `+`] (which is an excellent recent innovation from MapleSoft). So I incorporated that with the Bits:-Split, which makes a major efficiency improvement because it's externally compiled. I also removed some redundant counting and powering of 2s and 3s. The overall efficiency improvement is a factor of 30 to get 2^16 terms.

#VV's procedure:
restart:
A219954list := nmax -> [ 
seq[scan=`+`](
    ifelse(n=1, 0, 3^numboccur(convert(n-1, 'base', 2), 1) - ifelse(2^ilog2(n)=n, n/2, 0)),
    n = 1..nmax
)]:
A219954list(33);
 [0, 2, 5, 12, 15, 24, 33, 56, 59, 68, 77, 104, 113, 140, 167, 240, 243, 252, 
   261, 288, 297, 324, 351, 432, 441, 468, 495, 576, 603, 684, 765, 992, 995]

CodeTools:-Usage(A219954list(2^16)): %[-1];
memory used=0.64GiB, alloc change=-4.00MiB, 
cpu time=4.61s, real time=4.29s, gc time=718.75ms
                           4294901760
#Carl's procedure:
restart:
OEIS_A219954_B:= proc(n::posint)
uses S= Bits:-Split, C= numboccur;
local b:= ilog2(n), k, P:= 2^~[$0..b], Q:= 3^~[$b+1], L:= 1, c;
    [0, seq['scan'= `+`](Q[(c:= C(S(k-1), 1))] - `if`(c=L, P[L++], 0), k= 2..n)]
end proc
:
OEIS_A219954_B(33);
 [0, 2, 5, 12, 15, 24, 33, 56, 59, 68, 77, 104, 113, 140, 167, 240, 243, 252,
   261, 288, 297, 324, 351, 432, 441, 468, 495, 576, 603, 684, 765, 992, 995]

CodeTools:-Usage(OEIS_A219954_B(2^16)): %[-1];
memory used=17.07MiB, alloc change=2.00MiB, 
cpu time=125.00ms, real time=140.00ms, gc time=0ns
                           4294901760

I give you a Vote Up for using seq[scan= `+`], which I think is the best possible thing for this sequence, or indeed any sequence where each term is computed by adding something to the previous term. I learned that today, while coding for this Question; so, thank you @vv.

@2cUniverse The error is due to either using 2D Input, using an older version of Maple, or both. Please use the edited code in my Answer, which'll work in any input mode and in any version from the last 10 years or so.

My originally posted procedure used two relatively new syntax enhancements:

1. Arrow procedures can now have local declarations.
2. Embedded assignments: Assignments with := (and other infix assignment operators) can now appear within any expression in any position where their output alone (i.e., the right-side operand of the assignment) would be syntactically allowed, provided that the entire assignment is enclosed in parentheses (i.e., the operand and both operands within the same pair of parentheses). Previously, those assignments were required to be stand-alone statements.

Almost all syntax enhancements (other than changes in command-procedure arguments) take many years to be incorporated into 2D Input.

@dharr If one wants a consistent aspect ratio other than the one given by scaling= constrained, it might be possible with the size option. For example, plot(..., size= [800, 500]) might preserve an aspect ratio of 800/500 (horizontal/vertical). I don't know how to test this, but presumably you do.

Even a Google search for that package turns up nothing.

@C_R That's a good point that the number of gridlines displayed (the radial spokes in this case) may be less than the number used for computation. I don't know if the reduction factor is documented, but requesting 49 will display (49-1)/2 = 24. ((I can see that you know this; I'm pointing it out for other readers): We use an odd number (49, 25, etc.) in the grid specification because of "fenceposting": That is, the number of real intervals is 1 less than the number of evaluation points, which are the endpoints of the intervals.) The extra evaluation points round out that spider-web effect. 

I've also given the featured surface a bit of transparency with transparency= 0.15 to allow the underlying polar grid to show through. 

R:= 2:  #max r of featured surface 
spokes:= 24:  #number of radial gridlines in base

plots:-display(
    plot3d(  #the featured surface
        [r, theta, (csc(theta)/r)^2], r= 0..R, theta= -Pi..Pi,
        coords= cylindrical, style= surface, shading= zhue, transparency= 0.15
    ),
    plot3d(  #the base
        [r, theta, 0], r= 0..3/2*R, theta= -Pi..Pi, coords= cylindrical,
        grid= [7, 2*spokes+1], color= COLOR(RGB, 0.97$3), thickness= 2.5, glossiness= 0
    ),
    plots:-textplot3d(  #the spoke labels
        [seq](
            [7/4*R, i, 0, i],
            i= (2*Pi/spokes)*~([$0..spokes-1] -~ iquo(spokes,2))
        ), 
        coords= cylindrical, font= ["times", 10]
    ),
    view= [map(`*`, -1..1, 2*R)$2, 0..20], axes= framed, 
    orientation= [40,40,0], labels= [r$2, z], labelfont= ["times", 16], 
    axis[1,2]= [tickmarks= [-R, R]]
);

@mmcdara Vote up. Although I think that the OP still needs to provide more detail, you've managed to make a remarkable elaboration upon a minimally informative Question.

Obviously, you need to provide more context. Do you have any Maple code, or perhaps even a journal article, to accompany this request?

What happens if you use Array instrad of array?

@Art Kalb What kind of optimization do you expect for a procedure that only does abs of a complex argument? I don't think any optimization is possible. 

@C_R You can control the polar "base" grid size like this:

plots:-display(
    plot3d(
        [r, theta, (csc(theta)/r)^2], r= 0..2, theta= -Pi..Pi,
        coords= cylindrical, style= surface, shading= zhue
    ),
    plot3d(
        [r, theta, 0], r= 0..3, theta= -Pi..Pi, coords= cylindrical,
        grid= [6, 33], color= white, thickness= 2, glossiness= 0,
        shading= none
    ),
    view= [-3..3, -3..3, 0..99], axes= none
);

@C_R What you're calling the "z" axis is showing what you call "theta", and what you're calling "r" is the angular coordinate. This is more clear if you look at the axes with numeric tickmarks. So, the "r" and "theta" in this plot are not the radial and angular coordinates, and thus they can't correspond with the r and theta as used in the Question (despite the OP's apparent approval of this plot). In other words, taking the variables to mean what they usually mean in cylindrical coordinates, you've plotted
r = 1/theta^2/sin(z)^2.

If you plot in cylindrical coordinates without giving the three coordinates parametrically, then the 1st argument to plot3d is r, the 2nd is the range of theta, and the 3rd is the range of z. This is determined by the argument order; the actual variable names that you use are irrelevant. This seems weird to me, but that's the way that it's been for as long as I can remember. Like you, I find it more natural to think of z as function of r and theta.