Re: Re: Zero testing
- To: mathgroup at smc.vnet.net
- Subject: [mg53155] Re: [mg53152] Re: Zero testing
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Mon, 27 Dec 2004 06:41:24 -0500 (EST)
- References: <200412231259.HAA21206@smc.vnet.net> <cqgu81$632$1@smc.vnet.net> <200412250901.EAA18578@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
On 25 Dec 2004, at 18:01, Maxim wrote:
> 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].
Yes, I agree, I read your post too quickly; I apologize for this. You
did not say that ComplexInfinity was the right answer.
> 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.
There is a combination of different issues involved here.First, Reduce
is a much newer function, if I remember correctly it was completely
re-written in version 5. But also, it is not a good example, because it
is a function that suffers form exactly the same problem as
FullSimplify: it very frequently ends up being aborted since it is
unable to return any answer in reasonable time. This is O.K. because it
is relatively rarely used. Not so Integrate in which case you might
have noticed, if you have been reading this list, complaints about the
loss of speed compared with earlier versions are growing more frequent.
Using significance arithmetic to make verifications in all cases is
not a reasonable choice (as Daniel Lichtblau pointed out recently in
another context). Mathematica has to be able to identify the
situations when a numerical verification is needed and this is of
course the difficulty. Such verifications are not performed when
Mathematica "believes" it has returned the correct answer.
>
> 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:
Well, yes I agree, again I was too quick, I should have used N with a
higher precision setting (though it was quite obvious that this is what
I meant from the rest of my posting-- cf. my reply to your last
example). But the point I was making was actually not that different
form yours, which is that you can't really on being able to determine
such questions symbolically and numerical methods have to be used. The
issue, and the main difference between our views, seems to be whether
such things should be left to the user or performed by Mathematica
itself. For example in the Sign case, in my opinion Mathematica should
return
Sign[Im[Sqrt[-1-10^-18*I]]]
unevaluated, and then it ought to be up to the user to decide whether
to use FullSimplify[ComplexExpand[]] or N[], which may actually be
enough, or N[ ,20] (as I use in the last example) or possibly
sometimes even higher precision. In all cases human judgement will have
to be used for the final decision, which is another reason why I do not
think your demands for automatic checking in all cases are reasonable
(of course I may be misunderstanding you again but to tell the truth i
have never understood what it is that you actually want).
Andrzej Kozlowski
Chiba, Japan
http://www.akikoz.net/~andrzej/
http://www.mimuw.edu.pl/~akoz/
- References:
- Zero testing
- From: Maxim <ab_def@prontomail.com>
- Re: Zero testing
- From: Maxim <ab_def@prontomail.com>
- Zero testing