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>