In this previous post, an example is shown that demonstrates the potential problems that can arise following symbolic conversions such as from sqrt(x^2)  to x^(1/2).

Here x is an unknown symbol. The difficulties include the fact that, while `sqrt` can be smart about simplifying numeric values (eg. integers, rationals) the `^` operator has no such opportunity. Once the conversion from `sqrt` or `root` to `^` has been made there is no usual subsequent opportunity for complete simplification of powers, even when evaluation is done at a numeric value.

In short, eval( x^(1/2), x=4 ) will result in 4^(1/2) and not have the opportunity to simplify to the result of 2, which would be more convenient.

In the previous post, an expression 'q' was shown for which evaluation at a numeric point produced a fractional result whose numerator was zero but whose denominator was a "veiled zero". As was shown, Maple's current behaviour is that zero divided by veiled-zero immediately brings about a result of zero, with no opportunity to detect the exception. The result of zero, for that example, is incorrect.

If, on the other hand,  some smarter resolution of certain `^` terms could instead be done, then a veiled-zero in the denominator might be resolved as being a detectable zero. The final result would be the (much more satisfactory) detection of a division by zero -- a numeric exception event.

Below, I show an experiment toward working around that behaviour. It works by examining the input for one-over-integer powers, when evaluation is to be done at a numeric value. Just prior to actual evaluation, such `^` instances are replaced by effective calls to `root` or `sqrt`.

Note the I have to temporarily redefine `eval` in order to do this. That is seriously changing Maple, because `eval` is used in many hidden and important places. Basically, if `eval` is changed badly then key bits of Maple won't work. Of course, my changes only affect the single Maple session, and are not permanent.

> # The original problematic example
> q := (6*((1/3)*A-1/9))/(36*A-116+12*sqrt(12*A^3-3*A^2-54*A+93))^(1/3);
                                   6 (A/3 - 1/9)
            q := --------------------------------------------------
                                       3      2             1/2 1/3
                 (36 A - 116 + 12 (12 A  - 3 A  - 54 A + 93)   )
 
>
   
> eval(q, A = 1/3); # :(
                                       0
 

> `eval/F` := proc()
>   if nargs=2 and type(1/op(2,args[1]),integer) then
>       root(eval(op(1,args[1]),args[2]),1/(op(2,args[1])));
>   else
>     `^`(eval(op(1,args[1]),args[2]),op(2,args[1]));
>   end if;
> end proc:
>
> unprotect(eval): # living on the edge
> origeval := eval(eval):
>
> eval := proc(x)
> local b, dx, nx;
>   if nargs=0 then
>       error "invalid input: eval expects 1 or 2 arguments, but received 0";
>   elif nargs=1 then Z(x);
>   elif nargs=2 and type(args[2],symbol=numeric) then
>       nx := Z(subsindets(numer(x),`^`,t->'F'(op(t))),args[2..-1]);
>       dx := Z(subsindets(denom(x),`^`,t->'F'(op(t))),args[2..-1]);
>       nx/dx;
>   else
>       Z(args[1],args[2]);
>   end if;
> end proc:
>
> eval:=subs(Z=origeval(origeval),origeval(eval)):
> protect(eval);


> eval(q, A = 1/3); # :)
Error, (in eval) numeric exception: division by zero
>
> f:=(a-2)/(a-8^(1/3));
                                       a - 2
                                f := ----------
                                          (1/3)
                                     a - 8
 
> eval(f,a=2);
Error, (in eval) numeric exception: division by zero
>
>
> ee := (34-sqrt(A^2))^(1/3);
                                         2 1/2 (1/3)
                           ee := (34 - (A )   )
 
> eval(ee,A=7); # want result to be 3
                                       3
 
# some sanity checks
> eval(ee);
                                      2 1/2 (1/3)
                              (34 - (A )   )
 
> eval(ee,A=y*p);
                                     2  2 1/2 (1/3)
                             (34 - (y  p )   )
 
> eval('''ee''',1);
                                     ''ee''
 
> eval(ee,[A=7]); # weakness
                                       1/2 (1/3)
                               (34 - 49   )
 
> eval(A^(11/2),A=4); # weakness
                                         1/2
                                   1024 4
                                                                                
>
# Now revert.
> unprotect(eval);
> eval:=origeval(origeval):
> protect(eval);

Improvements to the (initial, simplistic) methodology would be welcome. This is also expected to be measurably slower than usual `eval`.

acer