Re: Zero testing
- To: mathgroup at smc.vnet.net
- Subject: [mg53133] Re: [mg53106] Zero testing
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Fri, 24 Dec 2004 05:59:39 -0500 (EST)
- References: <200412231259.HAA21206@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
On 23 Dec 2004, at 21:59, Maxim wrote:
> By 'zero testing' I mean the problem of deciding whether a given
> numerical
> quantity is exactly equal to zero. Often Mathematica doesn't make
> necessary checks at all:
>
> In[1]:=
> a = Pi/4; b = I*Log[-(-1)^(3/4)];
> Limit[(a - b)/x, x -> 0]
> Integrate[1/((x - a)*(x - b)), x]
> Integrate[1/(x - b), {x, 0, 1}]
>
> Out[2]=
> DirectedInfinity[Pi - (4*I)*Log[-(-1)^(3/4)]]
>
> Out[3]=
> (-4*((-2*I)*ArcTan[Log[-(-1)^(3/4)]/x] - 2*Log[Pi - 4*x] + Log[16*(x^2
> + Log[-(-1)^(3/4)]^2)]))/(2*Pi - (8*I)*Log[-(-1)^(3/4)])
>
> Out[4]=
> ((-I)*Pi - (2*I)*ArcTan[Log[-(-1)^(3/4)]] - Log[Log[-(-1)^(3/4)]^2]
> + Log[1 + Log[-(-1)^(3/4)]^2])/2
>
> The constants a and b are equal to each other. We can see that Out[2]
> is
> incorrect because it is equal to ComplexInfinity and the limit of 0/x
> is
> 0; Out[3] is also equal to ComplexInfinity (the denominator is zero)
> and
> therefore wrong as well; finally, Out[4] is incorrect because the
> integral
> doesn't converge.
This problem amounts to the difference between Simplify and
FullSimplify. FullSimplify is needed to determine that a and b are
equal, just Simplify cannot achieve this:
In[1]:=
a = Pi/4; b = I*Log[-(-1)^(3/4)];
In[2]:=
Simplify[a - b]
Out[2]=
(1/4)*(Pi - 4*I*Log[-(-1)^(3/4)])
In[3]:=
FullSimplify[a - b]
Out[3]=
0
Functions like integrate use Simplify automatically but not
FullSimplify. Why not? The answer you provided yourself (well, almost)
in an earlier posting: FullSimplify generally takes far too long to be
used automatically. Even using Simplify leads to problems, as you have
noted. Are you then arguing that FullSimplify ought to be used
automatically? I know it ould b enice if one were able to eat the cake
and still have it ...
If you apply FullSimplify to the argument in the above examples you
will get the answer you stated as correct. In all cases, that is,
except in one when it is you and not Mathematica that blundered. The
case in question is:
> Integrate[1/((x - a)*(x - b)), x]
Are you seriously claiming that
Integrate[1/((x - Pi/4)^2), x]
is ComplexInfinity? At least Mathematica knows better than that:
Integrate[FullSimplify[1/((x - a)*(x - b))], x]
4/(Pi - 4*x)
>
> A variation of the same problem is when we're numerically evaluating a
> discontinuous function for an argument close to a discontinuity point.
> In
> those cases Mathematica does perform necessary checks but handles more
> complicated cases incorrectly, especially when the branch cuts are
> concerned:
>
> In[5]:=
> Sign[Im[Sqrt[-1 - 10^-18*I]]]
>
> Out[5]=
> 1
The problem here is that Mathematica has great difficulty in
determining Im[Sqrt[-1 - 10^-18*I] is -1. Algebraically. You ned ot use
Sign[FullSimplify[ComplexExpand[Im[Sqrt[-1-10^-18*I]]]]]
-1
This I agree to be a problem that ought to be fixed. It seems to me
that numerical methods should be used since:
Sign[N[Im[Sqrt[-1-10^-18*I]]]]
-1
works easily and reliably. On the other hand it could be argued that
this sort of thing should be approached numerically anyway.
>
> In[6]:=
> ((1 + I*(-167594143/78256779 + Pi))^4)^(1/2)
>
> Out[6]=
> (1 + I*(-167594143/78256779 + Pi))^2
>
> In fact In[5] is equal to -1 and In[6] is equal to -(1
> + I*(-167594143/78256779 + Pi))^2. In the last case we do not even
> apply
> any functions, the error is in the automatic simplifications. Sometimes
> setting a high value for $MinPrecision helps, but here it doesn't have
> any
> effect either.
>
This one is indeed a problem. Even this
FullSimplify[ComplexExpand[
(1 + I*(-167594143/
78256779 + Pi))^4]]^
(1/2)
-(((167594143 + 78256779*I) -
78256779*Pi)^2/
6124123459454841)
gives the wrong answer ( the sign is wrong). To get the right one you
need to use high precision arithmetic and Unevaluated:
N[Unevaluated[
((1 + I*(-167594143/
78256779 + Pi))^4)^
(1/2)], 20]
1.5635873239815087691232`4.043\
607143112669*^-16 -
2.00000000000000015635873239\
81508708003`20.15051499783199*I
The only way that I have found of getting the correct exact answer is:
FullSimplify[ComplexExpand[((1 + I*(-u + Pi))^4)^
(1/2)] /. u -> 167594143/78256779]
((167594143 + 78256779*I) - 78256779*Pi)^2/
6124123459454841
In[54]:=
N[%, 20]
Out[54]=
1.5635873239815087691362`4.043607143112669*^-16 -
2.00000000000000015635873239815087005589`20.15051499783\
1987*I
It is curious that this method works while the one using Unevaluated
does not:
FullSimplify[ComplexExpand[((1 + I*(-u + Pi))^4)^
(1/2)] /. u -> 167594143/78256779]
FullSimplify[ComplexExpand[Unevaluated[
((1 + I*(-167594143/78256779 + Pi))^4)^(1/2)]]]
((167594143 + 78256779*I) - 78256779*Pi)^2/
6124123459454841
-(((167594143 + 78256779*I) - 78256779*Pi)^2/
6124123459454841)
While I think this Mathematica's dealing wiht this situation could be
improved (particularly the last example) I find it difficult to see
how such problems could be avoided in general. The problem seems to be
very subtle. On the one hand, one can't expect the Mathematica
evaluator to automatically apply powerful transformation s such as
ComplexExpand or FullSimplify. On the other hand, certain
transformations have to be automatically performed by the evaluator
otherwise Mathematica would become unbearably tedious to use (one could
leave every expression unevaluated until Simplify or FullSimplify or
possibly a weaker function were applied to it but it would still
require the user to judge when to apply which).
As far as I can tell all this means that the user has to remain
vigilant at all times and always think carefully about the answers
returned by the program.
Andrzej Kozlowski
Chiba, Japan
http://www.akikoz.net/~andrzej/
http://www.mimuw.edu.pl/~akoz/
- References:
- Zero testing
- From: Maxim <ab_def@prontomail.com>
- Zero testing