35 views

Skip to first unread message

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

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

May 18, 2013, 9:36:55 PM5/18/13

to sage-...@googlegroups.com

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
> 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.

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

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
> 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).

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

}}}

instead of

{{{

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

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 !

{{{

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 !

May 19, 2013, 5:52:24 AM5/19/13

to sage-...@googlegroups.com

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.

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).

May 19, 2013, 6:34:11 AM5/19/13

to sage-...@googlegroups.com

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.

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
> 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.

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.

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

May 19, 2013, 11:41:47 AM5/19/13

to sage-...@googlegroups.com

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?
> 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.

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

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

-leif

May 19, 2013, 11:50:46 AM5/19/13

to sage-...@googlegroups.com

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.

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
> 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.

*always* have a good answer. Sometimes the answer is nothing useful

(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.

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

like .n(adaptive=False))

- 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

like .n(adaptive=True))

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).

May 20, 2013, 10:02:18 AM5/20/13

to sage-...@googlegroups.com

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

*always* have a good answer.

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.?

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu