# RealField rounding error with negative input

35 views

### vdelecroix

May 18, 2013, 3:09:32 PM5/18/13
to sage-devel
Hi,

In an attempt to improve #14567, I found the following thing (which
for me is a bug). Consider the different real fields with rounding
RNDD or RNDU. Then -pi is systematically wrongly approximated.
The given input
{{{
a = -pi
for rnd in 'RNDD','RNDU','RNDZ':
R = RealField(15,rnd=rnd)
s,m,e = R(a).sign_mantissa_exponent()
print rnd, m
}}}
gives me
{{{
RNDD 25735
RNDU 25736
RNDZ 25735
}}}
Where I do expect
{{{
RNDD 25736
RNDU 25735
RNDZ 25735
}}}

Note that the conversion works very well for rational numbers but also
fails on number field elements (by replacing a above by
NumberField(x^3 - 2, 'a', embedding=1.25).gen()).

Best,
Vincent

### leif

May 18, 2013, 9:36:55 PM5/18/13
vdelecroix wrote:
> In an attempt to improve #14567, I found the following thing (which
> for me is a bug). Consider the different real fields with rounding
> RNDD or RNDU. Then -pi is systematically wrongly approximated.

AFAIK that's a (probably unexpected) feature rather than a bug, because

R(-pi) = R(-1)*R(pi)

(since pi is a symbolic constant, hence -pi a symbolic expression).

-leif

> The given input
> {{{
> a = -pi
> for rnd in 'RNDD','RNDU','RNDZ':
> R = RealField(15,rnd=rnd)
> s,m,e = R(a).sign_mantissa_exponent()
> print rnd, m
> }}}
> gives me
> {{{
> RNDD 25735
> RNDU 25736
> RNDZ 25735
> }}}
> Where I do expect
> {{{
> RNDD 25736
> RNDU 25735
> RNDZ 25735
> }}}
>
> Note that the conversion works very well for rational numbers but also
> fails on number field elements (by replacing a above by
> NumberField(x^3 - 2, 'a', embedding=1.25).gen()).

--
() The ASCII Ribbon Campaign
/\ Help Cure HTML E-Mail

### vdelecroix

May 19, 2013, 5:36:18 AM5/19/13
to sage-devel
On 19 mai, 03:36, leif <not.rea...@online.de> wrote:
> vdelecroix wrote:
> > In an attempt to improve #14567, I found the following thing (which
> > for me is a bug). Consider the different real fields with rounding
> > RNDD or RNDU. Then -pi is systematically wrongly approximated.
>
> AFAIK that's a (probably unexpected) feature rather than a bug, because
>
> R(-pi) = R(-1)*R(pi)
>
> (since pi is a symbolic constant, hence -pi a symbolic expression).

We should not have R(-pi) = R(-1) * R(pi) with rounding mode RNDD or
RNDU because it is how it is documented and how it works for
rationals. I agree that -pi is a symbolic expression but it is not a
reason for having a wrong answer. If we want to round down the
expression "- expr" then we should round up "expr" and then take the
opposite. Note that there are more dramatic problems with the symbolic
ring due to cancellation
{{{
sage: sage: a.n(digits=8)
293.65079
sage: sage: a.n(digits=10)
292.6277805
sage: sage: a.n(digits=12)
292.634618133
sage: sage: a.n(digits=15)
292.634591053247
}}}
What is shown above is that the rounding of ``a`` is *very far* from
being the nearest element of RR...

In my first post I also mentionned that it does not work either for
number field elements
{{{
a = NumberField(x^3-2, 'a', embedding=1.25).gen()
for rnd in 'RNDD','RNDU':
R = RealField(15,rnd=rnd)
s,m,e = R(-a).sign_mantissa_exponent()
}}}
gives
{{{
RNDD 20642
RNDU 20643
}}}
{{{
RNDD 20643
RNDU 20642
}}}
Moreover, the approximation does not work as expected either (the last
digits are always False)
{{{
sage: sage: b = (76504*a - 96389) / (90325 - 71691*a)
sage: b.n(digits=5)
0.096557
sage: b.n(digits=6)
0.0965556
sage: b.n(digits=7)
0.09655549
}}}

Vincent

### vdelecroix

May 19, 2013, 5:49:29 AM5/19/13
to sage-devel
Sorry, the first block of examples should be replaced with
{{{
sage: a = (106*pi - 333) / (355 - 113*pi)
sage: a.n(digits=5)
Traceback (most recent call last)
...
ValueError: power::eval(): division by zero
sage: a.n(digits=6)
289.000
sage: a.n(digits=7)
289.0000
sage: a.n(digits=8)
293.65079
sage: a.n(digits=9)
292.772502
sage: a.n(digits=10)
292.6277805
}}}
which is even worse !

### Volker Braun

May 19, 2013, 5:52:24 AM5/19/13
On Sunday, May 19, 2013 10:36:18 AM UTC+1, vdelecroix wrote:
If we want to round down the
expression "- expr" then we should round up "expr" and then take the
opposite.

If you care about errors then you should be evaluating the symbolic expressions with interval arithmetic. Your proposal is just to do interval arithmetic and then take the upper/lower bound as desired. That would be correct, but there should also be a way to evaluate symbolic expressions *quickly* without doing interval arithmetic.

### vdelecroix

May 19, 2013, 6:17:27 AM5/19/13
to sage-devel
- the specification and the behavior of numerical_approx
- how should be implemented conversion from an exact real constant to
a RealField or a RealIntervalField

In numerical approximation, it is written "Return a numerical
approximation with at least the requested number of bits or digits of
precision". It is False and even critical with the example a = (106*pi
- 333) / (355 - 113*pi) as even the third digit is False with the
option digits=8. On the other hand, I definitely think that a
conversion from an exact real number to a RealField or a
RealIntervalField must be the perfect answer (both are
unsatisfactory) !

@Volker: Do you have examples of a conversion of a symbolic expression
to RR with critical speed ? I found that interval arithmetic prevent
from getting False answer while it does not affect too much the speed
(from a user point of vue).

### Volker Braun

May 19, 2013, 6:34:11 AM5/19/13
For starters, most of the plotting stuff will rely on quickly computing numerical approximations to symbolic expressions. There you certainly don't want to start refining the RIF precision iteratively if you hit a numerically unstable point.

The only bug that I see is that the documentation of n() could be improved to explain that input precision =! output precision for symbolic expressions. But then thats pretty obvious.

### vdelecroix

May 19, 2013, 9:28:29 AM5/19/13
to sage-devel
On 19 mai, 12:34, Volker Braun <vbraun.n...@gmail.com> wrote:
> For starters, most of the plotting stuff will rely on quickly computing
> numerical approximations to symbolic expressions. There you certainly don't
> want to start refining the RIF precision iteratively if you hit a
> numerically unstable point.

I got it but it would be still prefarable to use RIF in doing plots
rather than plotting something False.

More seriously, I do make a difference between making an evaluation of
a function (ie computation of cos(x) for some given approximate real
number x) and a conversion of a constant (ie the conversion of (106*pi
- 333) / (355 - 113*pi) to an approximate field). In particular, when
I write ((106*pi - 333) / (355 - 113*pi)).n(digits=4) I do expect the
four first digits of my real number and not a "ValueError:
power::eval(): division by zero". Even if I am a starter.

> The only bug that I see is that the documentation of n() could be improved
> to explain that input precision =! output precision for symbolic
> expressions. But then thats pretty obvious.

It is not obvious as you do not get in the answer how much digits are
correct (claiming that they are different is not interesting nor
meaningful). Moreover, it is not only about symbolic expression as it
fails for number field elements.

Vincent

### leif

May 19, 2013, 11:41:47 AM5/19/13
vdelecroix wrote:
> In particular, when
> I write ((106*pi - 333) / (355 - 113*pi)).n(digits=4) I do expect the
> four first digits of my real number and not a "ValueError:
> power::eval(): division by zero". Even if I am a starter.

Isn't that pretty educational?

From the message you can immediately guess what you made wrong... ;-)

This would indeed be a good example to add to the docstring.

-leif

### Volker Braun

May 19, 2013, 11:50:46 AM5/19/13
Floating-point and numerically unstable formula just don't mix. Sure, working with RIF makes it much harder to get something that is plainly wrong, but you can still construct quickly-varying functions where it will give you the wrong result. There is no magic bullet, you have to learn about numerics if you want to use floating point numbers effectively.

It would still be nice to have an option .n(adaptive=True) or so to adapt the input precision until a desired output precision in reached. You also need a fuzzy zero in case the actual value is zero.

But you shouldn't write library code with the expectation that numerical_approx will always give you the correct value to the desired precision in face of numerically unstable equations. You should rewrite your algorithm in a numerically stable way and estimate the required input precision that you need to get your answer.

### vdelecroix

May 20, 2013, 9:40:57 AM5/20/13
to sage-devel
On May 19, 4:50 pm, Volker Braun <vbraun.n...@gmail.com> wrote:
> Floating-point and numerically unstable formula just don't mix. Sure,
> working with RIF makes it much harder to get something that is plainly
> wrong, but you can still construct quickly-varying functions where it will
> give you the wrong result. There is no magic bullet, you have to learn
> about numerics if you want to use floating point numbers effectively.

mpfi (on which RIF is based) is implemented in such way that you
(like the interval [-infinity...+infinity]) but it is always such that
the computer answer of f(I) contains the mathematical f(I). I am
curious to see your example and will add it to the documentation if it
exists !

> It would still be nice to have an option .n(adaptive=True) or so to adapt
> the input precision until a desired output precision in reached. You also
> need a fuzzy zero in case the actual value is zero.

Thanks for the suggestion. I opened a ticket in that direction #14620.

I would also like to ask a more general design question. There is a
conversion from SR to RR and RIF implemented respectively through the
method _mpfr_ and _real_mpfi_ (why the second one is explicitly
implemented in the __init__ of RealIntervalFieldElement and not as a
conversion method as in the case of RR ?). Both of them uses the
method _eval_self of a symbolic expression which does what we expect
and does it quick! Now the thing is that by doing
{{{
sage: expr = a = (106*pi - 333) / (355 - 113*pi)
sage: RR(expr)
}}}
I would like to get the floating point that is the nearest to my
expression (perhaps this is not what everybody wants ?). As you
mentioned there is a problem when the result is actually zero and expr
does not know itself it is zero. This might be a problem of SR and I
would be happy if there was a way to work with the transcendental
extension QQ[pi] where I can assert that an expression is zero.

More concretely, once the adaptive algorithm is implemented, it will
be possible to have two conversions from SR to RR:
- one which just works through _eval_self (and behaves
- one which try to check equality to zero, then make some adaptive
evaluation using intervals and if the precision increase too much
raise an error arguing that the answer is too small (behave
My question is what should be the default ?

The same question should be asked for number fields (and the answer
need not be the same).

### Volker Braun

May 20, 2013, 10:02:18 AM5/20/13
On Monday, May 20, 2013 2:40:57 PM UTC+1, vdelecroix wrote:
mpfi (on which RIF is based) is implemented in such way that you

A good answer, but not the mathematically correct answer for the image of the error interval. For example (I'll leave it as an exercise to write a similar unstable function with sin and 1/x:

def unstable(x):
if x == 0.01:
return 1
else:
return x

sage: unstable( RIF(0.0, 0.0) )    # correct
0
sage: unstable(0.01)
1
sage: unstable( RIF(-0.1, 0.1) )    # interval contains 0.01, so value attains 1 on interval
0.?