exp: adjust boundaries for single-precision floating point

[amyspark]

👋

In this MR, I adjust the boundaries of exp to account for the precision allowed by float.

I've used the following snippet of code to determine the exact values. These are, in each case, the last value of boundary.

Open for snippet

#include <Vc/Vc>
#include <iostream>
#include <limits>
#include <cmath>

const auto LOG2E = 1.4426950408889634073599;
const auto PX1exp = 1.26177193074810590878E-4;
const auto PX2exp = 3.02994407707441961300E-2;
const auto PX3exp = 9.99999999999999999910E-1;
const auto QX1exp = 3.00198505138664455042E-6;
const auto QX2exp = 2.52448340349684104192E-3;
const auto QX3exp = 2.27265548208155028766E-1;
const auto QX4exp = 2.00000000000000000009E0;
const auto EXP_LIMIT = 708;

template <typename A, int B>
Vc::SimdArray<A, B> geant_exp(Vc::SimdArray<A, B> initial_x) {
  Vc::SimdArray<A, B> x(initial_x);
  Vc::SimdArray<A, B> px(Vc::floor(LOG2E * x + 0.5));
  const auto n (Vc::simd_cast<Vc::SimdArray<int, B>>(px));
  x -= px * 6.93145751953125E-1;
  x -= px * 1.42860682030941723212E-6;
  const auto xx = x * x;

  // px = x * P(x**2).
  px = PX1exp;
  px *= xx;
  px += PX2exp;
  px *= xx;
  px += PX3exp;
  px *= x;

  // Evaluate Q(x**2).
  Vc::SimdArray<A, B> qx(QX1exp);
  qx *= xx;
  qx += QX2exp;
  qx *= xx;
  qx += QX3exp;
  qx *= xx;
  qx += QX4exp;

  // e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
  x = px / (qx - px);
  x = 1.0 + 2.0 * x;

  // Build 2^n in double.
  //x *= G4ExpConsts::uint642dp((((uint64_t)n) + 1023) << 52);

  // std::cout << x << std::endl;
  // std::cout << n << std::endl;
  // const auto a = Vc::ldexp(Vc::SimdArray<double, B>(1), n);
  // std::cout << a << std::endl;
  x = Vc::ldexp(x, n);
  // x = Vc::ldexp(Vc::simd_cast<Vc::SimdArray<double, B>>(x), n);
  // for (auto i = 0; i != 8; i++) {
  //   x[i] *= std::powf(2, n[i]);
  // }

  x(initial_x > EXP_LIMIT) = std::numeric_limits<A>::infinity();
  x(initial_x < -EXP_LIMIT) = 0;

  return x;
}

int main() {
  {
    Vc::SimdArray<float, 8> e{-87.0f, -87.5f, -88.0f, -88.1f,
                              -88.2f, -88.3f, -88.4f, -88.5f};

    std::cout << "exp ( " << e << " ) = " << Vc::exp(e) << std::endl;
    std::cout << "geant_exp ( " << e << " ) = " << geant_exp(e) << std::endl;
    std::cout << "geant_exp !double! ( " << e << " ) = " << geant_exp(Vc::simd_cast<Vc::SimdArray<double, 8>>(e)) << std::endl;
  }

  {
    bool limit = false;

    Vc::SimdArray<float, 4> e{-87.0f, -87.75f, -88.5f, 1.0f};

    while (!limit) {
      auto r = Vc::exp(e);
      std::cout.precision(std::numeric_limits<double>::max_digits10);

      std::cout << "exp ( " << e << " ) = " << r << " " << e[1] << std::endl;

      if (std::isfinite(r[0]) && std::isnan(r[1])) {
        e[0] = e[0];
        e[2] = e[1];
        const auto midpoint = e[0] + (e[2] - e[0]) / 2.f;
        if (e[0] == midpoint) {
          limit = true;
        }
        else {
          e[1] = midpoint;
        }
      } else if (std::isfinite(r[1]) && std::isnan(r[2])) {
        e[0] = e[1];
        e[2] = e[2];
        const auto midpoint = e[0] + (e[2] - e[0]) / 2.f;
        if (midpoint == e[2]) {
          limit = true;
        }
        else {
          e[1] = midpoint;
        }
      } else if (std::isfinite(r[0]) && std::isfinite(r[2]) ||
                 std::isnan(r[0]) && std::isnan(r[2])) {
        assert(!"Oops!");
      }
    }

    {
      bool limit = false;
      float min = -128.f;
      float max = -125.f;

      while (!limit) {
        const auto boundary = min + (max - min) / 2.0f;
        const auto res = std::logf(std::powf(2.f, boundary));

        std::cout << e[1] << " => " << min << " " << boundary << " " << max
                  << " => " << res << std::endl;

        if (res <= e[1]) {
          if (min == boundary || res == e[1]) {
            limit = true;
          }
          else {
            min = boundary;
          }
        } else {
          if (max == boundary || res == e[1]) {
            limit = true;
          }
          else {
          max = boundary;
          }
        }
      }
    }
  }

  {
    bool limit = false;

    Vc::SimdArray<float, 4> e{80.f, 90.f, 100.f, 1.0f};

    while (!limit) {
      auto r = Vc::exp(e);
      std::cout.precision(std::numeric_limits<double>::max_digits10);

      std::cout << "exp ( " << e << " ) = " << r << " " << e[1] << std::endl;

      if (std::isfinite(r[0]) && !std::isfinite(r[1])) {
        e[0] = e[0];
        e[2] = e[1];
        const auto midpoint = e[0] + (e[2] - e[0]) / 2.f;
        if (e[0] == midpoint) {
          limit = true;
        }
        else {
          e[1] = midpoint;
        }
      } else if (std::isfinite(r[1]) && !std::isfinite(r[2])) {
        e[0] = e[1];
        e[2] = e[2];
        const auto midpoint = e[0] + (e[2] - e[0]) / 2.f;
        if (midpoint == e[2]) {
          limit = true;
        } else {
          e[1] = midpoint;
        }
      } else if (std::isfinite(r[0]) && std::isfinite(r[2]) ||
                 !std::isfinite(e[0]) && !std::isfinite(e[2])) {
        assert(!"Oops!");
      }
    }

    {
      bool limit = false;
      float min = 127.f;
      float max = 128.f;

      while (!limit) {
        const auto boundary = min + (max - min) / 2.0f;
        const auto res = std::logf(std::powf(2.f, boundary));

        std::cout << e[1] << " => " << min << " "<<  boundary << " "<< max << " => " << res << std::endl;

        if (std::isfinite(res)) {
          if (min == boundary || res == e[1]) {
            limit =true;
          }
          else {
            min = boundary;
          }
        } else {
          if (max == boundary || res == e[1]) {
            limit = true;
          }
          else {
            max = boundary;
          }
        }
      }
    }
  }
}

Please let me know if I can improve it!

Fixes #273

[bernhardmgruber]

Hi again! Thank you for the great work! I will ask my colleague @amadio to review this, since he is also working on Geant4, to which you referred in the linked issue. This might take a few weeks though, since he is currently on vacation. You might try to play with the CI for now to get it compiling. Thank you again!

[amadio]

I think the change is ok, but the commit needs to be fixed before merging, as there's a comment terminator inside the comment that needs to be removed.

[amyspark]

I'm getting CI failures but I can't seem to find the actual error.

[amadio]

You can see the errors on CDash: http://cdash.cern.ch/index.php?project=Vc&date=2021-08-12, though I also only see warnings there for your builds (check with short commit hash 408...).

[bernhardmgruber]

@amyspark please rebase on 1.4 to get the CI working.

[bernhardmgruber]

LGTM. I would like @amadio's feedback as well though.