MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • References:
  • Prev by Date: Help with a summation
  • Next by Date: Re: Re: Re: Mathematica language issues
  • Previous by thread: Re: Zero testing
  • Next by thread: Re: Re: Zero testing