#include "double_extras.h"
#include "d_vec.h"
double
_d_vec_dot_heuristic(const double *vec1, const double *vec2, slong len2,
double *err)
{
double psum = 0, nsum = 0, p, n, d, t;
int pexp, nexp;
slong i;
for (i = 0; i < len2; i++)
{
t = vec1[i] * vec2[i];
if (t >= 0)
psum += t;
else
nsum += t;
}
nsum = -nsum;
if (err != NULL)
{
p = frexp(psum, &pexp);
n = frexp(nsum, &nexp);
p = d_mul_2exp(1.0, pexp - D_BITS);
n = d_mul_2exp(1.0, nexp - D_BITS);
d = FLINT_MAX(p, n);
*err = 2 * len2 * d;
}
return psum - nsum;
}