mpfr sin(pi) != 0

[Vincent Lefevre]

On 2012-03-01 18:51:01 +0100, Jacopo Nespolo wrote:
> Hi,
> 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)