exp() with exponent below -88 or so yields nan

[kfjahnke]

Vc version / revision | Operating System | Compiler & Version | Compiler Flags | Assembler & Version | CPU
1.4 Linux g++ 9.3.0, clang++ 10-0-0-4 Haswell i5
Vc built with clang++ 10-0-0-4, fresh checkout. Can supply more info if needed. Can you reproduce this?

Testcase

#include <iostream>
#include "Vc/Vc"

int main ( int argc , char * argv[] )
{
  Vc::SimdArray < float , 8 > e { -87.0f , -87.5f , -88.0f , -88.1f , -88.2f , -88.3f , -88.4f , -88.5f } ;
  auto r = exp ( e ) ;
  std::cout << "exp ( " << e << " ) = " << r << std::endl ;
}

Actual Results

exp ( <-87 -87.5 -88 -88.1 | -88.2 -88.3 -88.4 -88.5> ) = [1.64581e-38, 8.20976e-39, 3.54261e-40, -nan, -nan, -nan, -nan, -nan]

Expected Results

not nan, but 0 if result can't be represented in sp float.

[amyspark]

Can confirm with MSVC 19.29.30038.1, but it only happens when Vc_IMPL is not Scalar (whether AVX, SSE, or any in between doesn't matter).

[amyspark]

Here's a test case porting Geant4's approximation.

#include <Vc/Vc>
#include <iostream>

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);
  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};

  auto r = Vc::exp(e);
  std::cout << "exp ( " << e << " ) = " << r << std::endl;
  auto r2 = geant_exp(e);
  std::cout << "geant_exp ( " << e << " ) = " << r2 << std::endl;
}

When I try with Vc's single precision ldexp approximation, it fails as expected:

exp ( <-87 -87.5 -88 -88.1 | -88.2 -88.3 -88.4 -88.5> ) = [1.64581e-38, 8.20976e-39, 3.54261e-40, -nan, -nan, -nan, -nan, -nan]
geant_exp ( <-87 -87.5 -88 -88.1 | -88.2 -88.3 -88.4 -88.5> ) = <1.64581e-38 8.20976e-39 3.54261e-40 -nan | -nan -nan -nan -nan>

When I force Vc to use the double-precision approximation:

exp ( <-87 -87.5 -88 -88.1 | -88.2 -88.3 -88.4 -88.5> ) = [1.64581e-38, 8.20976e-39, 3.54261e-40, -nan, -nan, -nan, -nan, -nan]
geant_exp ( <-87 -87.5 -88 -88.1 | -88.2 -88.3 -88.4 -88.5> ) = <1.64581e-38 9.98235e-39 6.0546e-39 5.47844e-39 | 4.9571e-39 4.48535e-39 4.05851e-39 3.6723e-39>

The problem seems to be that ldexp's precision is simply not enough for this function. It would seem that when the original was written, the boundaries here

Vc/Vc/common/exponential.h

Lines 40 to 41 in e4b581f

	constexpr float MAXLOGF = 88.72283905206835f;
	constexpr float MINLOGF = -103.278929903431851103f; /* log(2^-149) */

weren't checked against float type precision.

EDIT: after further testing the input values to ldexp on the reporter's range,

  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);

exp ( <-87 -87.5 -88 -88.1 | -88.2 -88.3 -88.4 -88.5> ) = [1.64581e-38, 8.20976e-39, 3.54261e-40, -nan, -nan, -nan, -nan, -nan]
<1.4001 0.849204 1.03014 0.932108 | 0.843408 0.763142 1.38104 1.24962>
<-126 -126 -127 -127 | -127 -127 -128 -128>
<1.17549e-38 1.17549e-38 5.87747e-39 5.87747e-39 | 5.87747e-39 5.87747e-39 2.93874e-39 2.93874e-39>
geant_exp ( <-87 -87.5 -88 -88.1 | -88.2 -88.3 -88.4 -88.5> ) = <1.64581e-38 8.20976e-39 3.54261e-40 -nan | -nan -nan -nan -nan>

Vc hits just below subnormal precision.