Float operations gives unexpected result in vectors

[CorentinB78]

Description

Operator /= on a vector of double do not always give the same result as on an array or standard division.

Steps to Reproduce

#include <iostream>
#include <triqs/arrays.hpp>

int main() {
  // just defining the particular double I had problem with
  uint64_t n = 0x3ff40b512eb8870b;
  double x;
  memcpy(&x, &n, sizeof(x));
  std::cout << "x = " << x << std::endl;

  double c = x / x;
  std::cout << "1 - x/x = " << 1 - c << std::endl;

  triqs::arrays::array<double, 1> arr(1);
  triqs::arrays::vector<double> vec(1);

  arr(0) = x;
  arr /= arr(0);
  std::cout << "1 - x/x = " << 1 - arr(0) << std::endl;

  vec(0) = x;
  vec /= vec(0);
  std::cout << "1 - x/x = " << 1 - vec(0) << std::endl;

  vec(0) = x;
  std::cout << "1 - x/x = " << 1 - vec(0) / vec(0) << std::endl;

  return 0;
}

Expected output:

x = 1.25276
1 - x/x = 0
1 - x/x = 0
1 - x/x = 0
1 - x/x = 0

Actual output:

x = 1.25276
1 - x/x = 0
1 - x/x = 0
1 - x/x = 1.11022e-16
1 - x/x = 0

Versions

Triqs 3.0.0 59c511fbfc39837570d419769a27455c051718dc

OS: on my FI workstation centos-release-7-7.1908.0.el7.centos.x86_64

[Wentzell]

Hey Corentin,

Thank you for posting the issue! This is indeed an inconsistency in the implementation of array and vector. I checked and while the array class will just apply /= using foreach, vector will currently calculate the inverse of the rhs and use blas::scal to multiple all elements with this inverse.

These two operations (x/y and x * (1/y)) are not identical, c.f. https://godbolt.org/z/Yc41nf

@parcollet Was there a reason you opted for the two different solutions in vector and array?

[Wentzell]

@CorentinB78 I have made sure for both the unstable and the 3.0.x branch that we are using the array solution for the operator /= instead.