969203 - Float32: Take out all Math functions that don't properly commute

Status: RESOLVED FIXED

(Core :: JavaScript Engine: JIT, defect)

Description

[Gary Kwong [:gkw] [:nth10sd] (NOT official MoCo now)]

+++ This bug was initially created as a clone of Bug #958381 +++

f = (function(x) {
    return Math.tan((Math.fround(Math.log(Math.fround(Math.exp(x))))), x)
})
print(f(1))
print(f(1))

$ ./js --ion-parallel-compile=off --ion-eager w98-cj-in.js
1.5574075204780886
1.5574077246549023

$ ./js --ion-parallel-compile=off w98-cj-in.js
1.5574075204780886
1.5574075204780886

Tested on js 64-bit opt threadsafe shell on m-c changeset 1e9f169c9715 on Linux.

My configure flags are:

AR=ar sh ./configure --enable-optimize --disable-debug --enable-profiling --enable-gczeal --enable-debug-symbols --enable-methodjit --enable-type-inference --disable-tests --enable-more-deterministic --disable-exact-rooting --with-ccache --enable-threadsafe <other NSPR options>

Comment 1

[Gary Kwong [:gkw] [:nth10sd] (NOT official MoCo now)]

Attachment: regression window

There's a large regression window due to brokenness, but the following changeset stands out:

changeset:   http://hg.mozilla.org/mozilla-central/rev/e1226725f674
user:        Benjamin Bouvier
date:        Thu Oct 17 08:50:56 2013 +0200
summary:     Bug 918613 - Specialize some Maths function calls for Float32 in Ion. r=sstangl

Benjamin, is bug 918613 a potential regressor?

Comment 2

[Gary Kwong [:gkw] [:nth10sd] (NOT official MoCo now)]

Wrong bug number in changeset message, it should be bug 918163.

Comment 3

[Luke Wagner [:luke]]

For the trig functions, the spec specifically allows for an implementation-defined best effort:
  http://people.mozilla.org/~jorendorff/es6-draft.html#sec-function-properties-of-the-math-object
so this seems fine given that the different here is < 1e-7 (according to bbouvier on irc).  I expect we also get a different result from V8 (since they embed their own impl of various trig functions now).  I'd be interested to hear jorendorff's opinion on this, though, since ISTR he spent some time thinking about trig function precision.  The one unfortunate thing is that the interpreter/baseline will always call the double version; only Ion will call the float version, so the different in execution mode is observable.  (Ion is also observable via Date.now() and recursion limits, so this isn't the first.)

Comment 4

[Benjamin Bouvier [:bbouvier] (inactive)]

So, to be complete, the observed difference is due to the fact that in the interpreter, libc's log is used, and in ion, libc's logf is used. Log doesn't verify the float32 property. The difference between the two log calls is < 1e-7, then the error is propagated with the wrapping Math.tan call and it raises up to 2e-7, which is still < 1e-6.

Comment 5

[Jesse Ruderman]

> (Ion is also observable via Date.now() and recursion limits, so this isn't the first.)

jsfunfuzz already avoids comparing output for scripts that use Date or that hit a recursion limit.

We could do the same for the approximate-float math functions for this issue. That would take care of the hard cases (this bug and bug 948321), but that would make me slightly sad, because we've caught some real bugs involving math functions (bug 941381, bug 942550, bug 946969).

Comment 6

[Jason Orendorff [:jorendorff]]

This doesn't seem to reproduce for me.

Right now I'm thinking this is just a busted optimization. I can only guess at real-world code's requirements for Math functions, but it doesn't seem too much to ask that they return the same number if you call them twice.

The spirit of the Math API is as different from Date.now() as it could possibly be.

Comment 7

[Luke Wagner [:luke]]

(In reply to Jason Orendorff [:jorendorff] from comment #6)
Right, I'm starting to think the same thing.  (To be precise, the only thing that I think that is "busted" is that the observable results differ based on the JIT; I think using float precision for trig ops with float uses/defs would be fine if it were regular across all execution modes.)

Perhaps, then, we should rip out all the float optimizations of math stdlib functions that don't *precisely* commute.  Fortunately, the hottest trig ops seem to be sin/cos and bug 967709 should make these way faster than even now and also regular across all platforms and execution modes.

Comment 8

[Benjamin Bouvier [:bbouvier] (inactive)]

I just made a plain C program that generates random numbers and checks whether a math function commutes (it generates a million random float numbers and directly tests the property upon these inputs -- and it finds counter examples after only a few tries (<10), so this empirical method sounds acceptable).

The functions that can precisely commute are cos, sin, floor, ceil, round (the former is not specialized in Ion at the moment, see also bug 930477). The ones that don't precisely commute and that we do optimize for now are log, tan, atan, asin, acos. pow isn't optimized and doesn't commute precisely neither, but it was planned at some point (see also bug 930477 once again).

Removing these functions would imply that we can't use them in Odin too (otherwise there would be a difference of behaviour between Odin and no-asm modes). This is a shame as this added more expressivity to asm.js, since C code compiled with emscripten using float functions (as logf, for instance), could behave the same in JS and in native code, at the same speed. If we remove optimizations for these functions and still want to have them available in Odin, we'd have to: convert the float input to a double, call the double math function, then convert it back to a float, which is worse in terms of performance. But precision seems more important here, so that's clearly a question of tradeoff. What do you think?

Comment 9

[Luke Wagner [:luke]]

(In reply to Benjamin Bouvier [:bbouvier] from comment #8)
> The functions that can precisely commute are cos, sin, floor, ceil, round
> (the former is not specialized in Ion at the moment, see also bug 930477).

Well, that's not exactly a proof of commutativity :)  I'd be happy if you could find one, but I think I already asked the Intel FP expert and she said that it didn't strictly commute (but that it was also highly dependent on the C stdlib, since they all produce slightly different results).

> But precision seems more important here, so that's
> clearly a question of tradeoff. What do you think?

Well, it's not so much precision as having identical behavior in the different execution modes.  If we do bug 967709 we'll get the major perf wins associated with sin/cos (others?).

Comment 10

[Jason Orendorff [:jorendorff]]

Here is a proof that the optimization is sound for floor, ceil, and round.

First consider floor. The optimization is sound iff fround(floor(x)) === floorf(x) for all x in float32. What are the cases?

- This equation holds for float32's infinities, NaN, by inspection.

- Assuming floorf(-0) is -0, it holds for -0.

- All x >= 2^23 or <= -2^23 in float32 are integers, and for those, fround(floor(x)) = floorf(x) = x.

- For the rest, including the underflow range, floorf(x) gives a mathematically exact value: the greatest
  integer <= x.  This number is in float32: all the integers in this range (-2^23 to 2^23) are.
  Of course floor(x) returns the same mathematical value, and fround converts the result to float32,
  with no rounding, because the value is in float32.

The same proof goes for ceil and round.

(It is central to the proof is that the fround is a no-op here, and that these functions are not approximations, so sin and cos would require a rather different kind of proof.)

Comment 11

[Benjamin Bouvier [:bbouvier] (inactive)]

I'll take care of taking out all the math functions that don't commute properly.

Comment 12

[Vincent Lefevre]

The problem is here double rounding. Take a function f. You want, for any single-precision x (24-bit significand):

  rnd_24(f(x)) = rnd_24(rnd_53(f(x))).

What happens here is that instead of rounding directly to 24 bits (single precision), you first round it to 53 bits by using the double version of the function, then round it to 24 bits. Instead of having a single rounding, you have two. In directed rounding modes, one can show that one obtains the same result as with a single rounding. But for rounding to nearest, this is not always true. The problem can occur when rnd_53(f(x)) rounds to a midpoint for single precision: if the two successive roundings are done in the same direction, you'll get a different result from rnd_24(f(x)).

So, the problem occurs when f(x) can be very close to a midpoint for single precision... which is what we call a hard-to-round case of f (previously called "worst case"). If p is the final precision (here p = 24), from probabilistic arguments, the distance between the hardest to round case and a midpoint should be around 2^(exponent-2p-small_constant), meaning that in general, you need an intermediate precision of at least 2p+small_constant (the constant depends on the function) to guarantee that you get the same result. Here, p = 24, and with double precision, 53 = 2 * 24 + 5. This means that for some functions, using the double precision will always work, but for other functions, you may get differences for a few arguments[*] (the constant 5 here may be too small). One needs an exhaustive test to find out (not sure whether this has already been done). For single precision, one can do a naive test, e.g. using MPFR. This should be quite fast; so, no need for the methods I've designed for double precision and more.

[*] The failing cases should occur with a probability around 1/2^29, so that random testing as mentioned in comment #8 won't be conclusive in practice.

Note: this assumes that the double version is correctly rounded and you want correctly rounded results. Otherwise I think that you can consider that:
* double rounding is not a problem here: the accuracy of the result will almost not change due to double rounding;
* or on the contrary, double rounding may yield platform-dependent inconsistencies (failures will be platform-dependent), and you regard this as a major problem, so that you should not allow double rounding (except in particular cases such as floor, ceil and round). But IMHO, user programs should be written in a way not to be affected by such inconsistencies.

Comment 13

[Jason Orendorff [:jorendorff]]

(In reply to Vincent Lefevre from comment #12)
> The problem is here double rounding. Take a function f. You want, for any
> single-precision x (24-bit significand):
> 
>   rnd_24(f(x)) = rnd_24(rnd_53(f(x))).
>
> [...] if the two successive roundings are done in the same
> direction, you'll get a different result from rnd_24(f(x)).

Very nicely expressed.

Luke told me earlier today that he has been in touch with Intel and there is a proof that the optimization is safe for Math.sqrt() -- that is, fround(rnd_53(√x))) == rnd_24(√x) for all x in float32. I may have misunderstood him. Anyway, now I am skeptical! The proof would have to depend on some very particular property of square roots and the float32 values.

If Math.sin and Math.fround both always rounded toward zero, the optimization would be OK. But the bias would be distateful, and anyway Math.fround must use roundTiesToEven per spec.

> [...] But IMHO, user programs
> should be written in a way not to be affected by such inconsistencies.

I cannot agree with you there; but thank you very much for your lucid post!

Comment 14

[Vincent Lefevre]

(In reply to Jason Orendorff [:jorendorff] from comment #13)
> Luke told me earlier today that he has been in touch with Intel and there is
> a proof that the optimization is safe for Math.sqrt() -- that is,
> fround(rnd_53(√x))) == rnd_24(√x) for all x in float32. I may have
> misunderstood him. Anyway, now I am skeptical! The proof would have to
> depend on some very particular property of square roots and the float32
> values.

The proof for square root, reciprocal and division is elementary. You can find a study in T. Lang and J.-M. Muller's paper here:

  http://hal.inria.fr/inria-00072593/

which also deals with other algebraic functions (but for the other functions, the bounds are not reached, at least for most precisions).

So, I think that the only particular cases are the round-to-integer functions (floor, ceil, round), sqrt (operations like reciprocal and division are not expressed as functions, so that I think they are out of scope here), and the absolute value (I suppose that it is provided). Not sure because I don't know the set of functions. The good point about sqrt is that it should be correctly rounded on all platforms in practice.

For the other functions, there are choices to be made, as I've explained.

> If Math.sin and Math.fround both always rounded toward zero, the
> optimization would be OK. But the bias would be distateful, and anyway
> Math.fround must use roundTiesToEven per spec.

If support for directed rounding modes is added in some future JavaScript specification, you have to remember that the optimization is OK in these modes...

> > [...] But IMHO, user programs
> > should be written in a way not to be affected by such inconsistencies.
> 
> I cannot agree with you there; but thank you very much for your lucid post!

Well, perhaps JavaScript is (currently) not affected, but with advanced features such as parallelization, code generation, compile-time calculations, replacement of x*y+z by fma (allowed by some languages), replacement of sin(x) and cos(x) [both in the code] by sincos (which may give a different result), etc., inconsistencies tend to occur more and more often if behavior isn't fully specified (which can also prevent some optimizations). There are some examples here: https://www.vinc17.net/research/slides/sieste2010.pdf (slides 10 to 17). Trying to avoid inconsistencies is just a safety measure.

There's also the related question whether a server expects to get consistent results between several clients, which may have different implementations of math functions...

Comment 15

[Benjamin Bouvier [:bbouvier] (inactive)]

Attachment: Adios logf et al.

Comment 16

[Benjamin Bouvier [:bbouvier] (inactive)]

Attachment: Adios logf et al.

Better with the initial test case

Comment 17

[Sean Stangl [:sstangl]]

Comment on attachment 8379073 [details] [diff] [review]
Adios logf et al.

Review of attachment 8379073 [details] [diff] [review]:
-----------------------------------------------------------------

::: js/src/jit/AsmJS.cpp
@@ +4140,3 @@
>        case AsmJSMathBuiltin_floor:  arity = 1; doubleCallee = AsmJSImm_FloorD; floatCallee = AsmJSImm_FloorF;  break;
> +      case AsmJSMathBuiltin_sin:    arity = 1; doubleCallee = AsmJSImm_SinD;   floatCallee = AsmJSImm_Invalid; break;
> +      case AsmJSMathBuiltin_cos:    arity = 1; doubleCallee = AsmJSImm_CosD;   floatCallee = AsmJSImm_Invalid; break;

Is it necessary to rip out sin() and cos()? As Luke said in Comment 7, these are hot. Is this patch anticipating the new sin() implementation?

Comment 18

[Luke Wagner [:luke]]

There is a very promising one in bug 967709.

Comment 19

[Luke Wagner [:luke]]

Oops, I misread your question.  But, "yes".

Comment 20

[Gary Kwong [:gkw] [:nth10sd] (NOT official MoCo now)]

Is this testcase also likely related?

function testMathyFunction(f, inputs) {
    results = [];
    for (var j = 0; j < inputs.length; ++j) {
        for (var k = 0; k < inputs.length; ++k) {
            try {
                results.push(f(inputs[j], inputs[k]))
            } catch (e) {}
        }
    }
    print(uneval(results))
}
f = (function(x, y) {
    return ((Math.log(Math.log1p(~x))) / y)
})
testMathyFunction(f, [Number.MAX_VALUE, 0x80000000])


$ ./js-dbg-32-dm-ts-linux-22650589a724 --fuzzing-safe --ion-parallel-compile=off w7131-cj-in.js
[NaN, NaN, 1.706339210189184e-308, 1.428404023825881e-9]

$ ./js-dbg-32-dm-ts-linux-22650589a724 --fuzzing-safe --ion-parallel-compile=off --ion-eager w7131-cj-in.js
[NaN, NaN, 1.7063392101891837e-308, 1.428404023825881e-9]

Tested on 32-bit deterministic threadsafe debug shell on Linux, on rev 22650589a724.

Comment 21

[Benjamin Bouvier [:bbouvier] (inactive)]

(In reply to Gary Kwong [:gkw] [:nth10sd] from comment #20)
> Is this testcase also likely related?
> 
No, this one is due to the div operation, as in bug 948321

Comment 22

[Benjamin Bouvier [:bbouvier] (inactive)]

https://hg.mozilla.org/integration/mozilla-inbound/rev/5d9d0a9f3e1f

Comment 23

[Ryan VanderMeulen [:RyanVM]]

https://hg.mozilla.org/mozilla-central/rev/5d9d0a9f3e1f

Comment 24

[Benjamin Bouvier [:bbouvier] (inactive)]

For the record, I made up a program [1] that iterates over the entire float32 space and check for the commutativity property of these functions. I think my method is correct (iterating over the uint32 range and read the bytes as a float -- "pun") but I would love to hear comments or suggestions. The conclusion seems to be that only sqrt, round, floor and ceil can indeed commute, at least on my specific configuration (gcc 4.8, ubuntu 64 bits, libc's math lib, etc.).

The program can either show counter-examples for functions which don't commute [2] or count the number of times the property doesn't hold [3]. Results are quite surprising.

[1] https://github.com/bnjbvr/exhaustive32
[2] https://github.com/bnjbvr/exhaustive32/blob/master/check-property.txt
[3] https://github.com/bnjbvr/exhaustive32/blob/master/count.txt

Comment 25

[Benjamin Bouvier [:bbouvier] (inactive)]

Attachment: uplift patch

carrying r+ from sstangl.

https://tbpl.mozilla.org/?tree=Try&rev=35c3a55b88f6
If I understand correctly, the orange is due to a misconfiguration of the tree and I see that the B2G KK emulators builds aren't activated in the beta TBPL builds (plus it seems that the reds are due to harness errors, not jit-tests failures).

I built the ARM simulator shell locally and all tests pass with --tbpl.

Comment 26

[Benjamin Bouvier [:bbouvier] (inactive)]

Comment on attachment 8401860 [details] [diff] [review]
uplift patch

[Approval Request Comment]
Bug caused by (feature/regressing bug #): bugs 918613 and 904918
User impact if declined: precision errors with math functions
Testing completed (on m-c, etc.): m-i, m-c, locally on a clone of m-beta.
Risk to taking this patch (and alternatives if risky): very low
String or IDL/UUID changes made by this patch: n/a

Comment 27

[Ryan VanderMeulen [:RyanVM]]

https://hg.mozilla.org/releases/mozilla-beta/rev/7e9b33204db9