SSE2 Optimization for ParamEvent::SetTarget

Source/modules/webaudio/AudioParamTimeline.cpp

 /*
  * Copyright (C) 2011 Google Inc. All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions
  * are met:
  *
  * 1.  Redistributions of source code must retain the above copyright
  *     notice, this list of conditions and the following disclaimer.
  * 2.  Redistributions in binary form must reproduce the above copyright
@@ skipping 331 matching lines @@
 
             case ParamEvent::SetTarget:
                 {
                     currentTime = fillToTime;
 
                     // Exponential approach to target value with given time constant.
                     float target = event.value();
                     float timeConstant = event.timeConstant();
                     float discreteTimeConstant = static_cast<float>(AudioUtilities::discreteTimeConstantForSampleRate(timeConstant, controlRate));
 
+#if CPU(X86) || CPU(X86_64)
+                    // Necessary for unrolling..

[Raymond Toy, 2015/08/13 15:42:56]

As in the other panner node CL, can you include a brief comment on how you
derived these constants?

[adrian.belgun, 2015/08/14 07:49:21]

Done.

+                    const float c0 = discreteTimeConstant;
+                    const float c1 = 2 * c0 - c0 * c0;
+                    const float c2 = 3 * c0 - 3 * c0 * c0 + c0 * c0 * c0;
+                    const float c3 = 4 * c0 - 6 * c0 * c0 + 4 * c0 * c0 * c0 - c0 * c0 * c0 * c0;

[Raymond Toy, 2015/08/13 15:42:56]

Would it help accuracy if the constants were computed using Horner's rule like
so:

c1 = c0*(2-c0)
c2 = c0*((c0-3)*c0+3)
c3 = c0*(c0*((4-c0)*c0-6)+4)

If that's not enough, would it help if the computations were done using double
precision (but leaving the c<n> values as floats)?

[adrian.belgun, 2015/08/14 07:49:21]

Experimented also with doing the computations with doubles and also with
Horner's rule. Got the same result.
Will keep the expanded version for clarity.

[Raymond Toy, 2015/08/14 16:00:03]

You mean you computed the c's using doubles and also Horner's rule and still got
the same result?

Horner's rule has fewer operations (not really important here), and also
generally has much better roundoff errors.  I prefer Horner's rule here for
these reasons.

[adrian.belgun, 2015/08/17 07:58:14]

Yes, applied Horner's rule and for testing to make sure that the computations
were made on doubles made target, timeConstant,  discreteTimeConstant and c*
doubles. Still a 3.9e-5 error was needed to pass the test.

The issue is caused by the line operation 
  value += delta * c3;
If we exchange this with 
  value += delta * c0;
  delta = target - value;
  value += delta * c0;
  delta = target - value;
  value += delta * c0;
  delta = target - value;
  value += delta * c0;
The error needed to pass the test got reduced to 2.84e-5 (but still higher due
to SSE ops), but the performance gain got completely canceled.
Changed as per request.

[Raymond Toy, 2015/08/17 15:54:50]

Interesting. Thanks for doing the tests.

I think this change in threshold is acceptable.  Note also that the test
computes the expected result using the simple (non-SSE2) equation, so there
might be some round-off there as well.  The results might also be better than
indicated if we used a relative error instead of absolute.  I think the curve
starts at a gain of 100 so absolute errors might be larger than the relative
error.

+
+                    float delta;
+                    __m128 vC = _mm_set_ps(c2, c1, c0, 0);
+                    __m128 vDelta, vValue, vResult;
+
+                    // Process in batches of 4.
+                    unsigned fillToFrameTrunc = writeIndex + ((fillToFrame - writeIndex) / 4) * 4;
+                    for (; writeIndex < fillToFrameTrunc; writeIndex += 4) {
+                        delta = target - value;
+                        vDelta = _mm_set_ps1(delta);
+                        vValue = _mm_set_ps1(value);
+
+                        vResult = _mm_add_ps(vValue, _mm_mul_ps(vDelta, vC));
+                        _mm_storeu_ps(values + writeIndex, vResult);
+
+                        value += delta * c3;
+                    }
+#endif
+                    // Serially process remaining values
                     for (; writeIndex < fillToFrame; ++writeIndex) {
                         values[writeIndex] = value;
                         value += (target - value) * discreteTimeConstant;
                     }
 
                     break;
                 }
 
             case ParamEvent::SetValueCurve:
                 {
@@ skipping 96 matching lines @@
     // to the end of the values buffer.
     for (; writeIndex < numberOfValues; ++writeIndex)
         values[writeIndex] = value;
 
     return value;
 }
 
 } // namespace blink
 
 #endif // ENABLE(WEB_AUDIO)