Fix accuracy in angle computation

[mcara]

This PR is an alternative to #223 and it builds upon it. While I kept significant simplifications in the vector formulas (dot product of two cross products) proposed by @perrygreenfield in #223, I also reduced the number of calls to the sqrt function to normalize the cross products entering the dot product.

This was done because it turns out that the main source of "erratic values" in the 2nd double in the quad comes from the sqrt function. I have tried to keep the same algorithm and increase accuracy of computations but it still did not solve the issue completely: it only reduced the occurrence of these erratic values but did not eliminate them.

On the other hand, a complete replacement of the algorithm in the sqrt function with the Heron's method resulted in a much more stable algorithm (albeit it is likely slower than the original multiply-add algorithm). I kept the checks that @perrygreenfield introduced for the case when the 1st double is +/- 1 but do not use fabs as values smaller than 1 in magnitude could still be legitimate values.

Also, the formula for computing conjugate angles was slightly improved in the accuracy.

I have arrived at the conclusion that the "erratic values" are not due to bugs in the qd library but rather a limitation of the accuracy of quad numbers: the original paper by the authors of the qd library says that at least 212 bits should be significant. This means that quad numbers have about 64 significant digits, in agreement with observed magnitudes of the erratic values in the second double (<1e-65). Instead on relying on a hard-coded value, this PR uses epsilon() value provided by the quad library (which is hardcoded though).

Closes #222

[mcara]

Also, there are many other functions that compute unnecessary cross products and those computations should be reviewed in a separate PR.

[pllim]

Add that test? 😉

[mcara]

@pllim Working on it right now. Patience, my friend.

[mcara]

@pllim See second commit

[codecov]

Codecov Report

Base: 65.09% // Head: 65.09% // No change to project coverage 👍

Coverage data is based on head (5979731) compared to base (1aca4c3).
Patch has no changes to coverable lines.

❗ Current head 5979731 differs from pull request most recent head 3306362. Consider uploading reports for the commit 3306362 to get more accurate results

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #224   +/-   ##
=======================================
  Coverage   65.09%   65.09%           
=======================================
  Files           8        8           
  Lines        1252     1252           
=======================================
  Hits          815      815           
  Misses        437      437           

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[pllim]

I am not worthy to review this but I appreciate you adding a test. 😸

[perrygreenfield]

I also wish there was a way to partition changes made due to PEP8 changes and actual code changes. The PEP8 changes make it harder to see the code changes. Two separate PRs?

[perrygreenfield]

I don't understand this part where the product of the dot products of each vector is computed. In principle each dot product should be positive (or zero), but never negative, and yet there is a test for a negative value. If there is a way for the dot product to be negative through a bug, then both might be negative. And what does this have to do with normalizing with regard to the dot of A and B? I am missing something.

[mcara]

And what does this have to do with normalizing with regard to the dot of A and B? I am missing something.

It has nothing to do with normalizing with regard to the dot of A and B. Yes, you are missing the comment above the function:

normalized_dot_qd returns dot product of normalized input vectors.

The original code was doing this:

Anormalized = A / SQRT(A.A)
Bnormalized = B / SQRT(B.B)
AB = Anormalized.Bnormalized

Modified code does this:

AB = A.B / SQRT((A.A)*(B.B))

This reduces the number of calls to the SQRT function - the one that has the tendency t introduce those random bits in the end.

[perrygreenfield]

So I was missing something.

[mcara]

There are no PEP8 changes here. Unfortunately some lines had trailing white spaces and my editor strips those spaces on Save.

[perrygreenfield]

Looks Stupendous To Me

[perrygreenfield]

So I was missing something.