       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:=
>> 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=
>> DirectedInfinity[Pi - (4*I)*Log[-(-1)^(3/4)]]
>>
>> Out=
>> (-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=
>> ((-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
>> is
>> incorrect because it is equal to ComplexInfinity and the limit of 0/x
>> is
>> 0; Out is also equal to ComplexInfinity (the denominator is zero)
>> and
>> therefore wrong as well; finally, Out 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:=
> a = Pi/4; b = I*Log[-(-1)^(3/4)];
>
> In:=
> Simplify[a - b]
>
> Out=
> (1/4)*(Pi - 4*I*Log[-(-1)^(3/4)])
>
> In:=
> FullSimplify[a - b]
>
> Out=
> 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:=
>> Sign[Im[Sqrt[-1 - 10^-18*I]]]
>>
>> Out=
>> 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:=
>> ((1 + I*(-167594143/78256779 + Pi))^4)^(1/2)
>>
>> Out=
>> (1 + I*(-167594143/78256779 + Pi))^2
>>
>> In fact In is equal to -1 and In 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:=
> N[%, 20]
>
> Out=
> 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
is also equal to ComplexInfinity (the denominator is zero) and therefore
wrong as well", and you're talking about In; I hope you can see the
difference between the two. The denominator of Out 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:=
Reduce[Pi*x >
754334734322669483655537561140633/240112203426728897943499903682238*x]

Reduce::ncomp: Unable to decide whether
HoldForm[754334734322669483655537561140633/240112203426728897943499903682238]
and HoldForm[Pi] are equal. Assuming they are.

Out=
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:=
Sign[N[Im[Sqrt[-1-10^-310*I]]]]

Out=
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