# Numerical approximation error involving cosh and sinh

Hi. I'm new to Sage and to computational mathematics in general. I have been trying to code something, but I am struggling to understand some numerical approximation errors. My problem seems to boil down to the following phenomenon:

sage: x = exp(12*pi*i*i)
sage: x.n()
0.000000000000000


Which is as expected, $x=e^{-12\pi }$ is "very close" to 0. However, I do not understand why Sage gives me

sage: (5+x).n()
4.00000000000000


If I give it more bits of precision, then it works:

sage: (5+x).n(200)
5.0000000000000000424115118301607754401746440550887456940681


But I can't understand why it is making such grotesque error. Sage knows $x$ is very small (it returned 0.000000000000000). So why did it approximate $5+x$ to 4?

What is also weird, is that if I repeat the same computation with 4+x or 6+x instead, then it works just fine:

sage: (6+x).n()
6.00000000000000


My apologies if this is something trivial in computational maths, but can someone clarify what is going on here?

In a situation like this, is there a way for me to know if I have given Sage enough bits of precision to get at least the first few digits correct?

(Remark: I wrote exp(12*pi*i*i) instead of exp(-12*pi) to make Sage use sinh and cosh.)

edit retag close merge delete

( 2021-05-20 17:10:03 +0200 )edit

Congratulations! You found a bug!

( 2021-05-20 17:10:13 +0200 )edit

Thanks for the report. This reminds me of

although it seems quite distinct.

( 2021-05-20 17:13:29 +0200 )edit

So this is a bug? I was trying to do some computations with modular forms, but was getting some weird errors when evaluating them at some specific points of the upper halfplane. After a long time trying to understand what was going on, it boiled down to this issue here. So, should I open a report ticket? Or does this post already counts as one? Many thanks!!

( 2021-05-20 18:58:02 +0200 )edit

If we accept that exp(12*pi*i*i) is transformed into cosh(12*pi) - sinh(12*pi), then the observed behaviour is not a bug.

Indeed, evaluating cosh(12*pi) - sinh(12*pi) is done by evaluating pi to some precision, then computing 12*pi, then cosh(12*pi) and sinh(12*pi), then subtracting, and the loss of precision is such that no information on the value of the expression is left.

However it could be argued that transforming exp(12*pi*i*i) into cosh(12*pi) - sinh(12*pi) is in itself a bug, as it can lead to the loss of information observed in the question.

( 2021-07-05 11:08:22 +0200 )edit

Sort by » oldest newest most voted

A Sage real number $a$ can be written in the form $a=s\cdot m\cdot 2^e$, where $s$, $m$ and $e$ are its sign, mantissa and exponent, given by the method .sign_mantissa_exponent. For example,

sage: a = -12.345
sage: s, m, e = a.sign_mantissa_exponent(); s, m, e
(-1, 6949617174986097, -49)
sage: s*m*2^e.n()
-12.3450000000000


Let us see what is $x=\exp(-12\pi)=\exp(12\pi\cdot i\cdot i)$:

sage: x = exp(12*pi*i*i); x
cosh(12*pi) - sinh(12*pi)


We next compute the mantissa and exponent of $5+\cosh(12\pi)$ and $\sinh(12\pi)$:

sage: (5 + cosh(12*pi)).n().sign_mantissa_exponent()
(1, 5894625992139550, 1)
sage: sinh(12*pi).n().sign_mantissa_exponent()
(1, 5894625992139548, 1)


So SageMath gets $5+x$ as follows:

\begin{aligned}5+x&=5+\cosh(12\pi)-\sinh(12\pi) \\ &\approx 5894625992139550\cdot 2^1 - 5894625992139548\cdot 2^1=2\cdot 2=4.\end{aligned}

This justifies the result you get. You can repeat computations to check what happens with $6+x$:

sage: (6 + cosh(12*pi)).n().sign_mantissa_exponent()
(1, 5894625992139551, 1)


Hence: \begin{aligned}6+x&=6+\cosh(12\pi)-\sinh(12\pi) \\ &\approx 5894625992139551\cdot 2^1 - 5894625992139548\cdot 2^1=3\cdot 2=6.\end{aligned}

The problem here is that you are reaching the limits of the double precision floating point arithmetic. Just use more digits:

sage: x = exp(12*pi*i*i)
sage: R = RealField(100)
sage: R(5 + x)
5.0000000000000000000000000000

more

Thanks you so much!! This answer was incredibly helpful. So, it was not a bug!! Though, is there a way for me to know if I am using enough bits of precision, so that the first few digits of a given calculation are correct? I guess I could be using .sign_mantissa_exponent to check if the exponents in my calculations are small, right?

For example, having exponent 1 means the computer can only see things with precision $2^1$. So 4 and 5 are equal in this case. However, if I give Sage enough bits to make all the exponents in my calculations lower than $-10$, then I would have precision of $2^{-10}$. Is my reasoning correct?

( 2021-05-21 15:48:18 +0200 )edit

Now, with this in mind, I think I could code some kind of precision handling into many parts of my code. But I wonder, doesn't Python/Sage have some sort of command to deal with this kind of issue automatically?

Say, something that I could initiate at the start of my program. And then, if one of my numerical calculations gets exponents too big, it gives me a warning massage. Or something similar, to make sure I have used enough precision. Once again, many thanks!!

( 2021-05-21 15:59:49 +0200 )edit

Have a look at real or complex ball fields, which keep track of approximation error.

sage: x = exp(12*pi*i*i)
sage: RBF(x)
[+/- 67.1]
sage: 5 + RBF(x)  # not accurate enough
[+/- 72.1]
sage: 5 + RealBallField(80)(x)  # more precision
[5.00000 +/- 1.75e-6]
more

That is what I needed. Thanks!

( 2021-05-22 15:06:15 +0200 )edit