Implement Math.log1p using port from fdlibm.

Implement Math.log1p using port from fdlibm.

Port contributed by Raymond Toy <rtoy@google.com>;.

R=rtoy@chromium.org
LOG=N
BUG=v8:3481
Committed: https://code.google.com/p/v8/source/detail?r=23082

[Raymond Toy, 2014-08-08 17:32:18]

https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p.js
File test/mjsunit/es6/math-log1p.js (right):

https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p....
test/mjsunit/es6/math-log1p.js:20: 
Although not part of the CL, I noticed this. I think this test should avoid
testing down to 0.1, where log(x+1) and log1p(x) are supposed to be different.
Based on the comments from fdlibm log1p(x) = log(1+x) when x > 2^53.  Then your
delta should be zero, or at least closer to 1.1e-16 (float epsilon).

https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p....
test/mjsunit/es6/math-log1p.js:38: assertEqualsDelta(expected, Math.log1p(x),
expected * 1E-14);
Shouldn't the threshold be closer to 1.1e-16 instead of 1e-14?

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.h
File third_party/fdlibm/fdlibm.h (right):

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.h#...
third_party/fdlibm/fdlibm.h:25: struct TrigonometricConstants {
Should this be renamed? This no longer holds just trig constants.  Maybe
FdlibmConstants instead?

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js
File third_party/fdlibm/fdlibm.js (right):

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:17: // sin, cos, and tan, by Raymond Toy
(rtoy@google.com).
Update to say log1p is also done? Or maybe just remove the comment about sin,
cos, tan?

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:427: const TWO_THIRD = kMath[37];
I'm curious on why you added these constants here for logp1. You didn't do that
for the trig functions. Is there some advantage here (or disadvantage in the
trig functions)?

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:437: var f = 0;
Perhaps initialize f to x here and remove the assignment in line 464.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:440: var u = 0;
Maybe intialize u to x here and remove the assignment in line 480.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:444: if (ax >= 0x3ff00000) {
I would add a comment that |x| >= 1. (Really hard to know that's what ax >=
0x3ff00000 means.)

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:453: if ((TWO54 + x > 0) && (ax < 0x3c900000)) {
The original fdlibm code said two54+x > 0 was meant to raise the inexact flag. I
don't think JS exposes this, so perhaps this test can be removed.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:454: // |x| < 2^-54, so just return x
Probably should update the comment to say

// |x| < 2^-29. For |x| < 2^-54, we can return x. Otherwise use simple 2-term
Taylor series.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:465: hu = 1;
Line 438 initializes hu to 1, so this can probably be removed.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:469: if (hx >= 0x7ff00000) return x + x;
Add comment to say this handles infinities and NaN.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:473: // x < 9.007199254740992e15
Maybe note that 9.007...e15 is actually exactly 2^53.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:483: c = 0;
Since c is initialized to 0 in line 439, perhaps we can remove this.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:491: u = %_ConstructDouble(hu | 0x3fe00000,
%_DoubleLo(u));
fdlibm said this is normalizing u/2, so probably should add that comment.

[Raymond Toy, 2014-08-08 17:40:30]

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js
File third_party/fdlibm/fdlibm.js (right):

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:427: const TWO_THIRD = kMath[37];
On 2014/08/08 17:32:17, Raymond Toy wrote:
> I'm curious on why you added these constants here for logp1. You didn't do
that
> for the trig functions. Is there some advantage here (or disadvantage in the
> trig functions)?

Never mind. You actually did do this for the trig functions.

[Yang, 2014-08-11 07:51:07]

https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p.js
File test/mjsunit/es6/math-log1p.js (right):

https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p....
test/mjsunit/es6/math-log1p.js:20: 
On 2014/08/08 17:32:16, Raymond Toy wrote:
> Although not part of the CL, I noticed this. I think this test should avoid
> testing down to 0.1, where log(x+1) and log1p(x) are supposed to be different.
> Based on the comments from fdlibm log1p(x) = log(1+x) when x > 2^53.  Then
your
> delta should be zero, or at least closer to 1.1e-16 (float epsilon).

Done.

https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p....
test/mjsunit/es6/math-log1p.js:38: assertEqualsDelta(expected, Math.log1p(x),
expected * 1E-14);
On 2014/08/08 17:32:16, Raymond Toy wrote:
> Shouldn't the threshold be closer to 1.1e-16 instead of 1e-14?

Done.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.h
File third_party/fdlibm/fdlibm.h (right):

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.h#...
third_party/fdlibm/fdlibm.h:25: struct TrigonometricConstants {
On 2014/08/08 17:32:16, Raymond Toy wrote:
> Should this be renamed? This no longer holds just trig constants.  Maybe
> FdlibmConstants instead?

Renamed to MathConstants. Fdlibm is implied by namespace.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js
File third_party/fdlibm/fdlibm.js (right):

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:17: // sin, cos, and tan, by Raymond Toy
(rtoy@google.com).
On 2014/08/08 17:32:17, Raymond Toy wrote:
> Update to say log1p is also done? Or maybe just remove the comment about sin,
> cos, tan?

Done.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:437: var f = 0;
On 2014/08/08 17:32:16, Raymond Toy wrote:
> Perhaps initialize f to x here and remove the assignment in line 464.

Done.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:440: var u = 0;
On 2014/08/08 17:32:17, Raymond Toy wrote:
> Maybe intialize u to x here and remove the assignment in line 480.

Done.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:444: if (ax >= 0x3ff00000) {
On 2014/08/08 17:32:17, Raymond Toy wrote:
> I would add a comment that |x| >= 1. (Really hard to know that's what ax >=
> 0x3ff00000 means.)

Done.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:453: if ((TWO54 + x > 0) && (ax < 0x3c900000)) {
On 2014/08/08 17:32:17, Raymond Toy wrote:
> The original fdlibm code said two54+x > 0 was meant to raise the inexact flag.
I
> don't think JS exposes this, so perhaps this test can be removed.

Done.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:454: // |x| < 2^-54, so just return x
On 2014/08/08 17:32:17, Raymond Toy wrote:
> Probably should update the comment to say
> 
> // |x| < 2^-29. For |x| < 2^-54, we can return x. Otherwise use simple 2-term
> Taylor series.

Done.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:465: hu = 1;
On 2014/08/08 17:32:17, Raymond Toy wrote:
> Line 438 initializes hu to 1, so this can probably be removed.

Done.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:469: if (hx >= 0x7ff00000) return x + x;
On 2014/08/08 17:32:17, Raymond Toy wrote:
> Add comment to say this handles infinities and NaN.

-Infinity is actually caught further up. This only handles Infinity and NaN.
Added comment.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:473: // x < 9.007199254740992e15
On 2014/08/08 17:32:18, Raymond Toy wrote:
> Maybe note that 9.007...e15 is actually exactly 2^53.

Done.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:483: c = 0;
On 2014/08/08 17:32:17, Raymond Toy wrote:
> Since c is initialized to 0 in line 439, perhaps we can remove this.

Done.

https://codereview.chromium.org/457643002/diff/1/third_party/fdlibm/fdlibm.js...
third_party/fdlibm/fdlibm.js:491: u = %_ConstructDouble(hu | 0x3fe00000,
%_DoubleLo(u));
On 2014/08/08 17:32:17, Raymond Toy wrote:
> fdlibm said this is normalizing u/2, so probably should add that comment.

Done.

[Raymond Toy, 2014-08-11 16:40:03]

Did you forget to upload a new patch set? I can't see any changes.

[Yang, 2014-08-11 19:03:42]

On 2014/08/11 16:40:03, Raymond Toy wrote:
> Did you forget to upload a new patch set? I can't see any changes.

You are right. Thanks for reminding. I uploaded the new patch set.

[Raymond Toy, 2014-08-11 20:04:35]

Looks good, but I have one last comment.

https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p.js
File test/mjsunit/es6/math-log1p.js (right):

https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p....
test/mjsunit/es6/math-log1p.js:20: 
On 2014/08/11 07:51:06, Yang wrote:
> On 2014/08/08 17:32:16, Raymond Toy wrote:
> > Although not part of the CL, I noticed this. I think this test should avoid
> > testing down to 0.1, where log(x+1) and log1p(x) are supposed to be
different.
> > Based on the comments from fdlibm log1p(x) = log(1+x) when x > 2^53.  Then
> your
> > delta should be zero, or at least closer to 1.1e-16 (float epsilon).
> 
> Done.

Why 1e2 and not 2^53?. It's obviously true for x > 2^53 that 1+x = x so log1p(x)
== log(1+x). In fact, testing a few random values quickly finds, for example,

x = 203.02512450909163d0
log(1+x) = 5.318243145619017d0
log1p(x) = 5.318243145619018d0

(Both log and log1p used fdlibm algorithm)

So, what is the intent of this test? That log1p is close to log(1+x)? If so, you
should say so in a comment, and describe why 1e-14 is a good threshold for this.

[Raymond Toy, 2014-08-12 05:57:28]

lgtm

I don't want to block this CL because of the comment on the test of log1p(x) vs
log(1+x). I think the other parts of the code are correct.

[Yang, 2014-08-12 08:59:52]

On 2014/08/11 20:04:35, Raymond Toy wrote:
> Looks good, but I have one last comment.
> 
>
https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p.js
> File test/mjsunit/es6/math-log1p.js (right):
> 
>
https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p....
> test/mjsunit/es6/math-log1p.js:20: 
> On 2014/08/11 07:51:06, Yang wrote:
> > On 2014/08/08 17:32:16, Raymond Toy wrote:
> > > Although not part of the CL, I noticed this. I think this test should
avoid
> > > testing down to 0.1, where log(x+1) and log1p(x) are supposed to be
> different.
> > > Based on the comments from fdlibm log1p(x) = log(1+x) when x > 2^53.  Then
> > your
> > > delta should be zero, or at least closer to 1.1e-16 (float epsilon).
> > 
> > Done.
> 
> Why 1e2 and not 2^53?. It's obviously true for x > 2^53 that 1+x = x so
log1p(x)
> == log(1+x). In fact, testing a few random values quickly finds, for example,
> 
> x = 203.02512450909163d0
> log(1+x) = 5.318243145619017d0
> log1p(x) = 5.318243145619018d0
> 
> (Both log and log1p used fdlibm algorithm)
> 
> So, what is the intent of this test? That log1p is close to log(1+x)? If so,
you
> should say so in a comment, and describe why 1e-14 is a good threshold for
this.

I see. The original implementation did not have high accuracy as goal. The test
is more of a sanity check: for reasonably large numbers, log1p(x) should stay
reasonably close to log(1+x), with "reasonably" being 1E-14. Obviously, now that
we have a more accurate implementation, that test can be made more strict. I
changed it to x > 1E16.

[Yang, 2014-08-12 13:36:45]

Committed patchset #3 manually as 23082.

[Raymond Toy, 2014-08-12 17:28:01]

On 2014/08/12 08:59:52, Yang wrote:
> On 2014/08/11 20:04:35, Raymond Toy wrote:
> > Looks good, but I have one last comment.
> > 
> >
>
https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p.js
> > File test/mjsunit/es6/math-log1p.js (right):
> > 
> >
>
https://codereview.chromium.org/457643002/diff/1/test/mjsunit/es6/math-log1p....
> > test/mjsunit/es6/math-log1p.js:20: 
> > On 2014/08/11 07:51:06, Yang wrote:
> > > On 2014/08/08 17:32:16, Raymond Toy wrote:
> > > > Although not part of the CL, I noticed this. I think this test should
> avoid
> > > > testing down to 0.1, where log(x+1) and log1p(x) are supposed to be
> > different.
> > > > Based on the comments from fdlibm log1p(x) = log(1+x) when x > 2^53. 
Then
> > > your
> > > > delta should be zero, or at least closer to 1.1e-16 (float epsilon).
> > > 
> > > Done.
> > 
> > Why 1e2 and not 2^53?. It's obviously true for x > 2^53 that 1+x = x so
> log1p(x)
> > == log(1+x). In fact, testing a few random values quickly finds, for
example,
> > 
> > x = 203.02512450909163d0
> > log(1+x) = 5.318243145619017d0
> > log1p(x) = 5.318243145619018d0
> > 
> > (Both log and log1p used fdlibm algorithm)
> > 
> > So, what is the intent of this test? That log1p is close to log(1+x)? If so,
> you
> > should say so in a comment, and describe why 1e-14 is a good threshold for
> this.
> 
> I see. The original implementation did not have high accuracy as goal. The
test
> is more of a sanity check: for reasonably large numbers, log1p(x) should stay
> reasonably close to log(1+x), with "reasonably" being 1E-14. Obviously, now
that
> we have a more accurate implementation, that test can be made more strict. I
> changed it to x > 1E16.

I think you should have added a comment for the test that log1p(x) = log(1+x)
for x > 2^53 (or 1e16). 1e16 now looks like some random number without
explanation of why. 

I'm also curious if log1p(x) really is equal to log(1+x) so that the threshold
of 1e-16 could actually be 0. But that's for another day.