Protect against potential floating-point overflows when determining b…

…eta.

Otherwise, the numerator of the midpoint calculation in the binary
search could overflow, even if the value of the midpoint itself was
fine. Admittedly, this should not be possible as the operands should be
well below the float limits, but better safe than sorry.

include/qdtsne/gaussian.hpp

@@ -117,11 +117,11 @@ void compute_gaussian_perplexity(NeighborList<Index_, Float_>& neighbors, Float_
                             if (max_beta == max_value) {
                                 beta *= static_cast<Float_>(2);
                             } else {
-                                beta = (beta + max_beta) / static_cast<Float_>(2);
+                                beta += (max_beta - beta) / static_cast<Float_>(2); // i.e., midpoint that avoids problems with potential overflow.
                             }
                         } else {
                             max_beta = beta;
-                            beta = (beta + min_beta) / static_cast<Float_>(2);
+                            beta += (min_beta - beta) / static_cast<Float_>(2); // i.e., midpoint that avoids problems with potential overflow.
                         }
                     }