series(q,a=1/3,2);

                       1     (2/3)   (1/3)    /    1\
                       -- 108      13      + O|a - -|
                       54                     \    3/
map(evalf,%);
                                           /    1\
                           0.9874986897 + O|a - -|
                                           \    3/

limit(q,a=1/3);
                                  undefined
MultiSeries:-limit(q,a=1/3);
                                  undefined

MultiSeries:-series(q,a=1/3,2);
map(radnormal,%);
map(evalf,%);


                                 2                   /    1\
                 - ------------------------------ + O|a - -|
                                            (1/3)    \    3/
                     /   2     (1/2)  (1/2)\                
                   3 |- --- 676      9     |                
                     \  507                /                
                       1     (2/3)   (1/3)    /    1\
                       - (-4)      13      + O|a - -|
                       6                      \    3/
                                                    /    1\
                  -0.4937493449 + 0.8551989515 I + O|a - -|
                                                    \    3/

This MultiSeries is the latest algolib version on Maple 12.

PS Interesting: This plot 

dq:=denom(q);
plot([Re(dq),Im(dq)],a=0..1/2);

may explain why different series expansions produce different results, and the limits give "undefined" if they use a L'Hospital rule or something equivalent, as the derivative of the denominator change at a=1/3.

PS2 These other plots seem to show that the real and imaginary parts of the denominator on the complex plane are piecewise linear, with the planes meeting in a kind of "origami" pattern at a=1/3:

plots:-complexplot3d([Re(dq),.4],a=0-0.2*I..0.5+0.2*I,axes=boxed); plots:-complexplot3d([Im(dq),.2],a=0-0.2*I..0.5+0.2*I,axes=boxed);

Not an error due to documentation, which says:

  "Since eval does pointwise evaluation, eval cannot be used to
   evaluate an expression at a singularity. Use limit instead."

That's an ugly trap. In which sense ever one wants to consider q
(as expression to comute a number or in some algebraic structure).

So it meets the specification - it does not meet user's needs.


To have it simplier I used eval(q,a = a/3); eval(%,a=a+1)/2*3;
to have a=0 as the problematic point:

  Q := a/(12*a-104+4*(4*a^3+9*a^2-156*a+676)^(1/2))^(1/3);

The singularity is easy to find using singular: it is just a=0.

And the limit is seen not to exist (@Jakubi: take a wider range
-16 ... 10 to see, that Re or Im(denom) is not piecewise linear).


I am a fan of save results (i.e. Maple should check, whether eval
can be used) - but they may last long, so some 'lame and fast' eval
then would be helpfull.

NB: for the above subs(a=0,Q) or subs(a=0.0,Q) behaves similar.

Funny, that rationalize(Q) gives it as algebraic/a^2 and thus it
gives no problems.


Besides Alec's nice example I have some, where eval is not allowed
due to the documentation, but I see no easy way to find it out:

  J:= x -> Int( ln( (x-1)^2 + 4*x*sin(theta)^2 ), theta=0..Pi/2 );

  theQ:= x/J(x);

Now I want to check singularities by singular(theQ,x), which gives
{x = infinity}, {x = -infinity}. So the following should be correct
eval(theQ,x=0), eval(theQ,x=0.0); subs(x=0,theQ), subs(x=0.0,theQ);
They all return zeros.

However the following shows J to be identical 0 for -1 <= x <= 1:

  op(J(x)), method = _d01ajc: Int(%);
  plot([Re(J(x)), Im(J(x))], x= - 2..2, numpoints=5, thickness=2,
    color=[red,blue]);

I filed that function once and IIRC it was given by Robert Israel.

It also shows how dangerous it may be to use functions from some
expressions without thinking about the domain of definitions: for
the q posted here one should have in mind, that 3rd root is 'only'
a function on positive reals (which Alex is aware of, for sure;
Edited: however what is if I do not want to see that as function,
but as rational polynomial, say for x/x ?).

Just my 2 cents ... and thx - it is interesting!


BTW: how all that is done in other CAS (if comparable)?

SAGE gives an error returned from Maxima,

sage: a=var('a')
sage: q = (6*((1/3)*a-1/9))/(36*a-116+12*sqrt(12*a^3-3*a^2-54*a+93))^(1/3)
sage: q(1/3)
---------------------------------------------------------------------------
             Traceback (most recent call last)
... [long output skipped]
Error executing code in Maxima
CODE:
	(0/1) / (((-104/1) + ((12) * (sqrt(676/9)))) ^ (1/3));
Maxima ERROR:
	
Division by 0

Alec

If one uses the root command instead of 1/3 as the exponent to guarantee that you get the principle root, you get an immediate simplification for

F:= (a-2)/(a-root(8,3));

to 1.  The function form

f2 := a -> (a-2)/(a-root(8,3));

works perfectly at 2, but gives 1 for any variable input.

With the first example, if one uses simplify on the denominator for a=1/3 one gets 0.  Shouldn't a check for a zero denominator always be done?

q := (6*((1/3)*a-1/9))/(36*a-116+12*sqrt(12*a^3-3*a^2-54*a+93))^(1/3);
simplify(eval(denom(q),a=1/3));

Finally, a question:  Is the basic problem in Maple's evaluation of this expressions that it does not check for a 0 denominator?
 

I said that I would give my opinion on this. So, what follows is that: just my own personal opinion.

Let's call the functionality evalat, just to distinguish it from 1-argument eval.

I think that the behaviour in the posted example shows something that could (and should) be improved.

One scenario that concerns me here is that in which the calling of evalat is made deep inside some routine in the shipped Maple Library. That's quite different from cases where the user calls it interactively (directly) at the top-level. I'm thinking of cases where there is no chance for the user to intercede and choose to do some different analysis of the expression at or near the point.

And, as I only found out only a few days ago, this very example was such a case. The user received an answer for a complicated computation that turned out to be wrong. This (longstanding, expert) user was fortunately able to recognize that the result was wrong, and was able to trace the problem back to this particular evaluation at a point, which was being done within some Library routine. So there was no clean and simple workaround -- no immediate way to make the Library routine take some other approach.

As a Maple user, I would be prepared to take a small performance penalty if evalat could be improved here. If one wanted to use a super fast but mathematically unsafe evaluator then one could make do with subs, which is more intended for substitution into a data structure. But eval is supposed to be the smarter, mathematical safer evaluator.

And the improvement should be centralized. I don't buy a line of argument which says that individual Library routines all ought to be more careful in their use of evalat. The evalat functionality is used very widely within the Maple Library. If Maple's going to be improved to be stronger in its zero detection when evaluating expressions at points, for mathematical computation, then that should be done centrally.

Maple is not finished. No computer algebra system will ever be perfect. But the competition will improve, and thus so must Maple. I think that we ought to be able to seriously revisit this kind of functionality at least once every few years. It's in the class of issues for which -- while perfection is unattainable -- incremental, practical improvements are good.

This problem should not be impossible to cover adequately. Perhaps radicals could be dealt with at object simplification/uniquification time. Eg. 27^(1/2) -> 3*3^(1/2) automatically. Having evalat revert to substituting `^` by root and sqrt (as shown here) would be prohibitively expensive.

An extensible, user-customizable solution might be cheapest, if the default behaviour were not changed. That's because the system would only have to check something simple like a flag. And the user could choose how hard the system ought to work at zero testing. But then, how helpful would it be, if the default behaviour were left as it is?

Dave Linder
Mathematical Software, Maplesoft