New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Accuracy of sinh and complex cos/sin #247

Merged
merged 3 commits into from Dec 14, 2014

Conversation

Projects
None yet
4 participants
@pavpanchekha
Contributor

pavpanchekha commented Dec 12, 2014

The sin and cos function for complex arguments, and the sinh function for real arguments, are inaccurate when the inputs are very small. This is because Math.exp(x) - Math.exp(-x) returns zero for small x, instead of the more accurate 2x.

This patch replaces sinh by a Taylor expansion when the input is small, which increases accuracy.

return x + (x * x * x) / 6 + (x * x * x * x * x) / 120;
} else {
return (Math.exp(x) - Math.exp(-x)) / 2;
}

This comment has been minimized.

@BigFav

BigFav Dec 13, 2014

Collaborator

Style note: You have a return statement in the if, so you don't need the else.

if (Math.abs(x) < 1) {
  return x + (x * x * x) / 6 + (x * x * x * x * x) / 120;
}
return (Math.exp(x) - Math.exp(-x)) / 2;

This comment has been minimized.

@josdejong

josdejong Dec 13, 2014

Owner

I think both versions are ok. @BigFav version is more compact, but on the other hand, an explicit else shows very clearly that the second return is only occurring in an else case (there are rules of thumb that say that every if should have an else prevents cases similar to these tricky fall-though cases of switches). Both are fine with me.

This comment has been minimized.

@pavpanchekha

pavpanchekha Dec 14, 2014

Contributor

I'll leave it as is for now, but would be happy with whatever you decide.

@@ -30,6 +30,7 @@ describe('cos', function() {
approx.equal(cos(pi*7/4), 0.707106781186548);
approx.equal(cos(pi*8/4), 1);
approx.equal(cos(pi/4), math.sqrt(2)/2);
assert.deepEqual(cos(complex('1e-50+1e-50i')), complex(1, -1e-100));

This comment has been minimized.

@BigFav

BigFav Dec 13, 2014

Collaborator

You fail this test case.

1) cos should return the cosine of a number:

AssertionError: {"re":1,"im":-1.0000000000000001e-100} deepEqual {"re":1,"im":-1e-100}

I think you should assert to the value you get (which is 1 - 1.0000000000000001e-100i), and comment that the true value is transcendental.

This comment has been minimized.

@josdejong

josdejong Dec 13, 2014

Owner

Yes, there is a round-off error here. Maybe replace the added assert.deepEqual(...) with

approx.deepEqual(cos(complex('1e-50+1e-50i')), complex(1, -1e-100));

This comment has been minimized.

@BigFav

BigFav Dec 13, 2014

Collaborator

The issue is that approx will pass this test case even if the resulting imaginary value is zero, tricky stuff lol.

This comment has been minimized.

@josdejong

josdejong Dec 13, 2014

Owner

Ah good point. Well then some more precise comparison is needed, or another test value instead of 1e-50+1e-50i, which does not result in a round-off error.

This comment has been minimized.

@pavpanchekha

pavpanchekha Dec 14, 2014

Contributor

This is a tricky situation indeed; on my machine the current test passes. I think there is a hardware issue here: this will pass, or not, depending on whether your copy of Node.js was compiled to use x87 or SSE floating point instructions, and which the JIT decided to use. I'll look for an adequate solution, but I'm a bit stymied without a buggy software / hardware combo. What do you use?

This comment has been minimized.

@pavpanchekha

pavpanchekha Dec 14, 2014

Contributor

I just updated my PR with a simpler solution: test that the imaginary part is nonzero, not its exact value.

@@ -24,6 +24,12 @@ describe('sinh', function() {
approx.equal(sinh(1), 1.1752011936438014);
});
it('should return the sinh of very small numbers', function() {
// If sinh returns 0, that is bad, so we are using assert.equal, not approx.equal

This comment has been minimized.

@BigFav

BigFav Dec 13, 2014

Collaborator

I think you should merge this comment into the test description, or at least be more explicit in the description. Notice how in line 33 it explicitly states (downgrades to number), I mean something along those lines; maybe add a simple (avoid returning zeros).

This comment has been minimized.

@pavpanchekha

pavpanchekha Dec 14, 2014

Contributor

Done.

@BigFav

This comment has been minimized.

Collaborator

BigFav commented Dec 13, 2014

Good job catching this case. Most of my comments are style related (with the only exception of the failed test case), so hopefully this will be merged pretty smoothly.

@josdejong

This comment has been minimized.

Owner

josdejong commented Dec 13, 2014

Thanks @pavpanchekha, nice improvement! Can you please fix this round-off error making the unit tests fail? After that I will merge your PR.

Fix the test case so it passes on all hardware.
The fix is to test that the imaginary part is nonzero,
not that it is exactly 1e-100.
In some software / hardware combos,
it seems to return 1.0000000000000001e-100.
@josdejong

This comment has been minimized.

Owner

josdejong commented Dec 14, 2014

Thanks a lot @pavpanchekha

josdejong added a commit that referenced this pull request Dec 14, 2014

@josdejong josdejong merged commit 77e32bb into josdejong:develop Dec 14, 2014

1 check passed

continuous-integration/travis-ci The Travis CI build passed
Details
@jrus

This comment has been minimized.

jrus commented Oct 19, 2015

This 5th degree taylor expansion of sinh is quite far off when you get close to -1 or 1. For example, for x = 0.99, it’s off by about 1.87e-4. I assume that’s dramatically more error than you want in a math library.

For a graph of the error, here’s Wolfram Alpha: http://www.wolframalpha.com/input/?i=sinh%28x%29+-+%28x+%2B+%281%2F6%29+*+x%5E3+%2B+%281%2F120%29+*+x+%5E+5%29+from+-1+to+1

The GNU Scientific Library uses the degree 17 Taylor expansion for computing sinh in the range [-1, 1]:
https://github.com/ampl/gsl/blob/1a819c2daab24f5f/specfunc/trig.c#L39-L52

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment