The function that determines the sign of a^b is full of bugs; here is a quick example:
(%i1) assume(x<0);
(%o1) [x<0]
(%i2) csign(x^x);
(%o2) neg
But, of course, csign((-2)^(-2)) = csign(1/4) = pos.
The function mexpt-sign has to handle many cases--since csign returns, I think, 8 possible values (neg, nz, zero, pz, pos, pn,
imaginary, complex), if csign(a^b) only examines the signs of a and b, it has to handle at least 64 cases. Additionally there are
special cases of even and odd powers and roots, so mexpt-sign has to handle 100s of cases. Maybe a (gaint) nested conditional
isn't the best approach?
Here is a bit of code that tests csign(a^b)--it shows many related bugs:
(%i1) assume(n < 0, nz <=0, pz >=0, p >0, equal(z,0), notequal(pn,0))$
(%i2) l : [n,nz,pz,p,z,pn]$
(%i3) outermap(lambda([a,b], 'csign(a^b) = block([w : errcatch(csign(a^b))], if w = [] then false else first(w))),l,l)$
sign: division by zero in z^n
sign: division by zero in z^nz
sign: division by zero in z^z
(%i4) funmake('matrix,%);
(%o4) matrix([csign(n^n)=neg,csign(n^nz)=pn,csign(n^pz)=pn,csign(n^p)=neg,csign(n^z)=pos,csign(n^pn)=neg],[csign(nz^n)=neg,csign(nz^nz)=neg,csign(nz^pz)=pn,csign(nz^p)=neg,csign(nz^z)=pos,csign(nz^pn)=neg],[csign(pz^n)=pos,csign(pz^nz)=pos,csign(pz^pz)=pos,csign(pz^p)=pos,csign(pz^z)=pos,csign(pz^pn)=pos],[csign(p^n)=pos,csign(p^nz)=pos,csign(p^pz)=pos,csign(p^p)=pos,csign(p^z)=pos,csign(p^pn)=pos],[csign(z^n)=false,csign(z^nz)=false,csign(z^pz)=zero,csign(z^p)=zero,csign(z^z)=false,csign(z^pn)=zero],[csign(pn^n)=pos,csign(pn^nz)=pos,csign(pn^pz)=pos,csign(pn^p)=pos,csign(pn^z)=pos,csign(pn^pn)=pos])