This one is wrong even with the triangular storage (including default, i.e. without storage specialization). In this case the shape certainly has to be changed for such a map (or give an error).

Certainly, I agree that producing non-skewsymmetric matrices with skewsymmetric shape and rectangular storage is a bug. It is interesting that diagonal elements are not changed,

map(1,Matrix(2,shape=skewsymmetric,storage=rectangular));

                               Matrix([[0,1],[1,0]])

With other elements, Maple seems to be not able to choose which storage to use, upper triangular, or lower triangular, so it applies map to both of them :) Without any doubt, most of the help pages could be written much better.

With the triangular storage, map seems to work OK with skewsymmetric matrices,

map(1,Matrix(2,shape=skewsymmetric,storage=triangular[upper]));

                    Matrix([[0,1],[-1,0]])

map(1,Matrix(2,shape=skewsymmetric,storage=triangular[lower]));

                    Matrix([[0,-1],[1,0]])

It's not a bug - it is a feature. The elements of matrices with shape=identity and shape=zero can not be changed, i.e. one can't map on them, assign them etc. The same as, say, with shape=triangular[upper] - the elements under the diagonal can not be assigned.

A partial solution of the second problem is given by the following improved code,

Define:=()->
define(subsindets([args],specfunc(anything,args[1]),
y->subsindets(y,indexed,x->convert(op(0,x),function,[op(x)])))[]):

Subs:=proc(a,b) local T,c;
Define(T,T(subs(`&*`=`*`,`&+`=`+`,lhs(a)))=rhs(a));
c:=subsindets(b,subs(map(y->y=op(2,y),
indets(subsindets(lhs(a),indexed,
x->convert(op(0,x),function,[op(x)])),symbol::type)),lhs(a)),T);
eval(c,T=(x->x)) end:

`&/`:=proc(a,b) local S;
S:=selectremove(x->sprintf("%s",x)[-1]="_",
indets(subsindets(lhs(b),indexed,
x->convert(op(0,x),function,[op(x)])),name));
Subs(subs(map(x->x=x::anything,S[1]),
map(x->x=cat(x,`_`)::identical(x),S[2]),lhs(b))=
subs(map(x->nprintf("%s",sprintf("%s",x)[1..-2])=x,S[1]),
map(x->x=cat(x,`_`),S[2]),rhs(b)),a) end:

For example,

(K[a,c]*K[b,c]) &/ (K[i_,k_]&*K[j_,k_] = s(i,j));

                               s(a, b)

That doesn't work with another example though,

(K[a][c]*K[b][c]) &/ (K[i_][k_]&*K[j_][k_] = s(i,j));

                           K[a][c] K[b][c]

The first problem has roots in the way subs works,

subs(a+b=x,a+b);

                                  x

subs(a+b=x,a+b+c);

                              a + b + c

If one needs to do such a replacement, applyrule can be used instead of subs,

applyrule(a+b=x,a+b+c);
                                x + c

Similarly here, I can write ApplyRule procedure working analogously, and incorporate it in &/ instead of Subs. I need some time for that though. Type "sequence" also can be added manually, I think, without big problems. I'll look into that - then double and triple underscores could be added to &/.

Works OK in Maple 10.

With sequences the situation is slightly different, because, as you noticed, there is no such type as an expression sequence. Nevertheless, Subs can be used. In your 1st examples,

Subs(a::specfunc(anything,f)=p('op(a),op(a),op(a)'), f(a,b,c));

                     p(a, b, c, a, b, c, a, b, c)

Subs(a::specfunc(anything,f)=g('0,op(a),0,op(a),0'), f(x,y,z)+f());

              g(0, x, y, z, 0, x, y, z, 0) + g(0, 0, 0)

For 2nd example, one has to use the modified version of Subs (which works OK on the previous examples, too),

Subs:=proc(a,b) local T,c;
define(T,T(subs(`&*`=`*`,`&+`=`+`,lhs(a)))=rhs(a));
c:=subsindets(b,subsindets(lhs(a),symbol::type,curry(op,2)),T);
eval(c,T=(x->x))
end:

Subs(eps(i::anything,j::anything,k::anything) &* 
eps(i::anything,j::anything,l::anything)=delta(k,l), 
eps(a,b,c)*eps(a,b,d));

                             delta(c, d)

Actually, for cases with one underscore (i.e. not for sequences), notation can be made very close to Mathematica notation,

`&/`:=proc(a,b) local S;
S:=selectremove(x->sprintf("%s",x)[-1]="_",indets(lhs(b),name));
Subs(subs(map(x->x=x::anything,S[1]),
map(x->x=cat(x,`_`)::identical(x),S[2]),lhs(b))=
subs(map(x->nprintf("%s",sprintf("%s",x)[1..-2])=x,S[1]),
map(x->x=cat(x,`_`),S[2]),rhs(b)),a) end:

(f([a,b])+f(c)) &/ (f([x_,y_])=p(x+y));

                           p(a + b) + f(c)

[1,x,x^2,x^3] &/ (x^n_=r(n));

                          [1, x, r(2), r(3)]

(eps(a,b,c)*eps(a,b,d)) &/ (eps(i_,j_,k_)&*eps(i_,j_,l_)=delta(k,l));

                             delta(c, d)

For the operator precedence reasons, one has to add parentheses in both sides of &/.

The following procedure does that,

Subs:=proc(a,b) local T;
define(T,T(lhs(a))=rhs(a));
subsindets(b,subsindets(lhs(a),symbol::type,curry(op,2)),T)
end:

For example,

Subs(f([x::anything,y::anything])=p(x+y), f([a,b])+f(c));

                           p(a + b) + f(c)

Subs(a::identical(x)^n::anything=r(n), [1,x,x^2,x^3]);

                          [1, x, r(2), r(3)]

As Joe Riel said, RealDomain can be used,

eqn:=-p*(x^2-a^2)^(3/2)+q*a^3:
RealDomain:-solve(eqn,a):
simplify([%],symbolic);

                     1/3                     1/3
                    p    x                  p    x
            [--------------------, - --------------------]
               (2/3)    (2/3) 1/2      (2/3)    (2/3) 1/2
             (q      + p     )       (q      + p     )

As Joe Riel said, RealDomain can be used,

eqn:=-p*(x^2-a^2)^(3/2)+q*a^3:
RealDomain:-solve(eqn,a):
simplify([%],symbolic);

                     1/3                     1/3
                    p    x                  p    x
            [--------------------, - --------------------]
               (2/3)    (2/3) 1/2      (2/3)    (2/3) 1/2
             (q      + p     )       (q      + p     )

By the way, Sarah and I graduated the same year from UPenn. Another connection with Maple is that Appalachian State Math Dept Chair is Bill Bauldry, a well known Maple expert (they are searching for the new chair starting next year though).

Also, &* can be used. For example,

expand((a+2*b)&*(2*a+3*b));

          2 (a &* a) + 3 (a &* b) + 4 (b &* a) + 6 (b &* b)

`&^`:=(a,n)->expand(`&*`(a$n)):
           
(a+b)&^2;

              (a &* a) + (a &* b) + (b &* a) + (b &* b)

(a+b)&^3;

  &*(a, a, a) + &*(a, a, b) + &*(a, b, a) + &*(a, b, b) + &*(b, a, a)

         + &*(b, a, b) + &*(b, b, a) + &*(b, b, b)

(2*a+3*b)&^2;

          4 (a &* a) + 6 (a &* b) + 6 (b &* a) + 9 (b &* b)

From the plot, one can see that one root is 0 and another one is located between -6 and -5. It can be found using fsolve,

fsolve(x^4 + 5*x^3 + 4*x=sin(x),x=-6..-5);
                             -5.156989371

If I didn't make arithmetical mistakes, it is

2*int(int(int(int(1,z=-(x-w)^2/4/y..1),y=-1..-(x-w)^2/4),w=-x..1),x=-1..1);

                                 23/9

If I didn't make arithmetical mistakes, it is

2*int(int(int(int(1,z=-(x-w)^2/4/y..1),y=-1..-(x-w)^2/4),w=-x..1),x=-1..1);

                                 23/9