expFloat32x4 function

Float32x4 expFloat32x4(Float32x4 entry)

SIMD Optimized Exponential function.

Author: Graciliano Monteiro Passos

Implementation

Float32x4 expFloat32x4(Float32x4 entry) {
  entry = entry.clamp(_float32x4__87, _float32x4_87);

  final x = entry.x;
  final y = entry.y;
  final z = entry.z;
  final w = entry.w;

  var xInt = _expXToInt(x);
  var yInt = _expXToInt(y);
  var zInt = _expXToInt(z);
  var wInt = _expXToInt(w);

  final xIdx = expIntTableMaxIndex + xInt;
  final yIdx = expIntTableMaxIndex + yInt;
  final zIdx = expIntTableMaxIndex + zInt;
  final wIdx = expIntTableMaxIndex + wInt;

  // Lookup exp(floor(x)).
  // intPartA will have the upper 22 bits, intPartB will have the lower
  // 52 bits.
  final intPartA = Float32x4(expIntTableA[xIdx], expIntTableA[yIdx],
      expIntTableA[zIdx], expIntTableA[wIdx]);
  final intPartB = Float32x4(expIntTableB[xIdx], expIntTableB[yIdx],
      expIntTableB[zIdx], expIntTableB[wIdx]);

  final entryInt = Float32x4(
      xInt.toDouble(), yInt.toDouble(), zInt.toDouble(), wInt.toDouble());

  // Get the fractional part of x, find the greatest multiple of 2^-10 less than
  // x and look up the exp function of it.
  // fracPartA will have the upper 22 bits, fracPartB the lower 52 bits.
  final frac = ((entry - entryInt) * _float32x4_1024);

  final fracIntX = frac.x.toInt();
  final fracIntY = frac.y.toInt();
  final fracIntZ = frac.z.toInt();
  final fracIntW = frac.w.toInt();

  final fracPartA = Float32x4(expFracTableA[fracIntX], expFracTableA[fracIntY],
      expFracTableA[fracIntZ], expFracTableA[fracIntW]);
  final fracPartB = Float32x4(expFracTableB[fracIntX], expFracTableB[fracIntY],
      expFracTableB[fracIntZ], expFracTableB[fracIntW]);

  // epsilon is the difference in x from the nearest multiple of 2^-10.  It
  // has a value in the range 0 <= epsilon < 2^-10.
  // Do the subtraction from x as the last step to avoid possible
  // loss of precision.
  final epsilon = entry -
      (entryInt +
          Float32x4(
                fracIntX.toDouble(),
                fracIntY.toDouble(),
                fracIntZ.toDouble(),
                fracIntW.toDouble(),
              ) /
              _float32x4_1024);

  // Compute z = exp(epsilon) - 1.0 via a minimax polynomial.  z has
  // full double precision (52 bits).  Since z < 2^-10, we will have
  // 62 bits of precision when combined with the contant 1.  This will be
  // used in the last addition below to get proper rounding.

  // Remez generated polynomial.  Converges on the interval [0, 2^-10], error
  // is less than 0.5 ULP
  var o = _float32x4_0_04168701738764507;
  o = o * epsilon + _float32x4_0_1666666505023083;
  o = o * epsilon + _float32x4_0_5000000000042687;
  o = o * epsilon + _float32x4_1_0;
  o = o * epsilon + _float32x4__3_940510424527919E_20;

  // Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
  // expansion.
  // tempA is exact since intPartA and intPartB only have 22 bits each.
  // tempB will have 52 bits of precision.
  final tempA = intPartA * fracPartA;
  final tempB =
      intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;

  // Compute the result.  (1+z)(tempA+tempB).  Order of operations is
  // important.  For accuracy add by increasing size.  tempA is exact and
  // much larger than the others.  If there are extra bits specified from the
  // pow() function, use them.
  final tempC = tempB + tempA;
  final result = tempC * o + tempB + tempA;

  return result;
}