Re: Zero testing

*To*: mathgroup at smc.vnet.net*Subject*: [mg53152] Re: Zero testing*From*: Maxim <ab_def at prontomail.com>*Date*: Sat, 25 Dec 2004 04:01:01 -0500 (EST)*References*: <200412231259.HAA21206@smc.vnet.net> <cqgu81$632$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

On Fri, 24 Dec 2004 11:23:45 +0000 (UTC), Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote: > > 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/ > I suggest that you re-read my post more carefully. I didn't claim that Integrate[1/((x - Pi/4)^2), x] was ComplexInfinity: my post says "Out[3] is also equal to ComplexInfinity (the denominator is zero) and therefore wrong as well", and you're talking about In[3]; I hope you can see the difference between the two. The denominator of Out[3] does equal zero, therefore the whole expression is identical to ComplexInfinity, which is why it is an incorrect result for Integrate[1/((x - Pi/4)^2), x]. On the whole, you pretty much miss the point. Certainly I do not suggest that FullSimplify be applied at every step of the evaluation or anything of the sort. In any case, FullSimplify won't always be able to verify a == b (take a = Derivative[0, 0, 1, 0][Hypergeometric2F1][1, 1, 1, 1 - E] and b = 1/E). When I'm saying that necessary checks aren't performed, in most cases that would mean checks based on significance arithmetic. Like I said, many discontinuous functions such as Floor do that (read the documentation on Floor and think about what should happen when we evaluate, say, Floor[Pi/Rationalize[Pi, 10^-100]]). In particular, Reduce does just that, generating Reduce::ncomp message when needed: In[1]:= Reduce[Pi*x > 754334734322669483655537561140633/240112203426728897943499903682238*x] Reduce::ncomp: Unable to decide whether HoldForm[754334734322669483655537561140633/240112203426728897943499903682238] and HoldForm[Pi] are equal. Assuming they are. Out[1]= False Even though in this particular example the result is wrong, this is a reasonable approach, because if two different but close numbers are incorrectly assumed equal, it should be possible to correct this by increasing $MaxExtraPrecision, while if a equals b but they're incorrectly assumed unequal, there is no way to fix that. That's what happens in the Limit and Integrate examples: Reduce performs necessary checks (numerical checks!), many other functions don't, that's all. I hope you realize that this isn't related to Simplify/FullSimplify in any way. Your suggestion to use Sign[N[Im[Sqrt[-1-10^-18*I]]]] is perhaps the worst possible method of dealing with this; since N doesn't guarantee to give any correct digits, nothing can be concluded from its results: In[2]:= Sign[N[Im[Sqrt[-1-10^-310*I]]]] Out[2]= 1. You will notice that this "easily and reliably" gives incorrect results for small values of the imaginary part (starting roughly on the scale of 10^-Accuracy[0.]). Not that there's any problem with that, it just means that machine-precision N is inapplicable here. Maxim Rytin m.r at inbox.ru

**Follow-Ups**:**Re: Re: Zero testing***From:*Andrzej Kozlowski <akoz@mimuw.edu.pl>

**References**:**Zero testing***From:*Maxim <ab_def@prontomail.com>