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>
- Re: Re: Zero testing
- References:
- Zero testing
- From: Maxim <ab_def@prontomail.com>
- Zero testing