mpfr sin(pi) != 0
vincent at vinc17.net
Thu Mar 1 23:30:09 CET 2012
On 2012-03-01 18:51:01 +0100, Jacopo Nespolo wrote:
> this is slightly, but not so much, off topic, as it is a problem
> related to mpfr rather than gmp strictly speaking.
It should have been posted in the MPFR mailing-list, or better,
a mailing-list related to programming or floating-point, since
this is not specific to MPFR (MPFR just allows you to compute
with more precision and provides "correct rounding", but doesn't
do exact computations, see below).
> I found myself calculating
> sin(pi * n * x)
> with n in N and x in R. I used mpfr_const_pi for pi.
> Sometimes it happens that n*x = 1, and this screwed the whole
> calculation because i would get 10^-(a lot, according to working
> precision) instead of zero.
This is normal. The number pi cannot be represented exactly, thus
its value is rounded, hence the error in the result. Since pi is
not a rational number, sin(expr) in floating point can be 0 only
if expr evaluates to 0.
If you want an exact computation, you need a computer algebra
system. Note that some packages like iRRAM are able to track
the errors, but unless they are particularly smart (like a
computer algebra system), they won't be able to decide the
exact value (0) of sin(pi).
> This was also a particularly unlucky case, because it happened to be
> in the subdiagonal of a sparse matrix, surrounded by zeros, thus the
> diagonalisation algorithms in Eigen considered important to keep it.
> Do you think it could be worth it to file a bug, at least to warn
> about the issue in the documentation of the trigonometric functions,
> or it is just expected behaviour and I got fool because I'm a newbie?
I don't think there is anything to add in the MPFR manual (BTW,
most functions are affected by rounding errors). Users should just
learn how floating-point works.
Actually, there may be something about this added to the manual in
the future, once mpfr_sinpi (and/or mpfr_sinu from LIA-2), etc. are
implemented. Such functions would avoid the problem you are seeing
with trigonometric functions (if the argument is computed with
enough precision so that it rounds to an integer).
Vincent Lefèvre <vincent at vinc17.net> - Web: <http://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)
More information about the gmp-discuss