P1673: `vector_sum_of_squares`: Rescaling algorithm may not be correct

@mhoemmen
I was looking at this and it seems to me this is wrong?

  Scalar scale = init.scaling_factor;
  Scalar ssq = init.scaled_sum_of_squares;
  for (extents<>::size_type i = 0; i < x.extent(0); ++i) {
    if (abs(x(i)) != 0.0) {
      const auto absxi = abs(x(i));
      const auto quotient = scale / absxi;
      if (scale < absxi) {
          ssq = Scalar(1.0) + ssq * quotient * quotient;
          scale = absxi;
      }
      else {
        ssq = ssq + quotient * quotient;
      }
    }
  }

• why does scale change along the way?

• also, shouldn't scale / absxi; instead be absxi/scale? Looking at the spec, section 16.9.2.7:

  Constraints: For all i in the domain of v, and for f and ssq of type T, the expression ssq = ssq + (abs(x[i]) / f)*(abs(x[i]) / f) is well formed.
  

  It says abs(x[i])/f should be well-formed, but not f/abs(x[i]).

• also, if things are correct, in theory one should have that:

  scaling_factor * scaling_factor * scaled_sum_of_squares = sum of squares of abs(x[i]) plus init.scaling_factor * 
  init.scaling_factor * init.scaled_sum_of_squares.
  

  see dlassq. But that does not work for the impl above.

So i tried something on my end:

{ 
  using std::abs;
  using std::max;

  T scale = init.scaling_factor;
  for (std::size_t i = 0; i < x.extent(0); ++i) {
    scale = max(scale, abs(x(i)));
  }

  T ssq = (init.scaling_factor*init.scaling_factor*init.scaled_sum_of_squares)/(scale*scale);
  T s= {};
  for (std::size_t i = 0; i < x.extent(0); ++i) {
    const auto absxi = abs(x(i));
    const auto quotient = absxi/scale;
    ssq = ssq + quotient * quotient;
    s += absxi*absxi;
  }

  std::experimental::linalg::sum_of_squares_result<T> result;
  result.scaled_sum_of_squares = ssq;
  result.scaling_factor = scale;

  // verify that things are consistent according to definition
  const auto lhs = scale*scale*ssq;
  const auto rhs = s+init.scaling_factor*init.scaling_factor*init.scaled_sum_of_squares;
  std::cout << lhs << " " << rhs << '\n';
  EXPECT_NEAR(lhs, rhs, 1e-9);
}

and this works, in the sense that it satisfies:

  scaling_factor * scaling_factor * scaled_sum_of_squares = sum of squares of abs(x[i]) plus init.scaling_factor * 
  init.scaling_factor * init.scaled_sum_of_squares.

Start with include/experimental/__p1673_bits/blas1_vector_sum_of_squares.hpp and the vector_sum_of_squares rescaling loop. Compare its quotient and scale updates with P1673 section 16.9.2.7 and the linked LAPACK dlassq reference. Done means the implementation satisfies the stated scaling-factor invariant, including the initial values.