`to_double(Quotient<Gmpzf>)` is not monotonic
Opened this issue · 2 comments
Issue Details
The function CGAL::to_double<Quotient<Gmpzf>> is not monotonic (even weakly), i.e., it can be true that a > b && to_double(a) < to_double(b) even when both a and b are mathematically in the interval [0, numeric_limits<double>::max()].
- For my application, it is very important that I preserve monotonicity when converting from
Quotient<Gmpzf>todouble. This conversion is not in a hot loop so performance is a secondary concern. Is there a recommended way to do this with the existing API? - I realise that there is nowhere in the documentation that provides any sort of monotonicity guarantee for
to_double, but IMHO this is surprising behaviour and it should either be corrected (if possible with not too much effort) or documented explicitly.
Source Code
#include <CGAL/Gmpzf.h>
#include <CGAL/Quotient.h>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <sstream>
auto construct(std::stringstream& base, const int& e) -> CGAL::Gmpzf {
CGAL::Gmpz m;
base >> m;
CGAL::Gmpzf result_base(m);
CGAL::Gmpzf result_exp(pow(2, e));
return result_base * result_exp;
}
auto main() -> int {
using CGAL::Gmpzf, CGAL::Quotient;
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10);
// Numerator of a
auto num_a_base = std::stringstream(
"495465331884540240762104420639278860096116250934219950844827973437034498625"
); //*2^-379
auto num_a_exp = -379;
auto num_a = construct(num_a_base, num_a_exp);
std::cout << "num_a = ";
CGAL::print(std::cout, num_a);
std::cout << '\n';
// Denominator of a
auto denom_a_base = std::stringstream("2245694908428994193174821578352066558595985337089");
auto denom_a_exp = -299;
auto denom_a = construct(denom_a_base, denom_a_exp);
std::cout << "denom_a = ";
CGAL::print(std::cout, denom_a);
std::cout << '\n';
// Numerator of b
auto num_b_base = std::stringstream("3447327498334006902041169");
auto num_b_exp = -215;
auto num_b = construct(num_b_base, num_b_exp);
std::cout << "num_b = ";
CGAL::print(std::cout, num_b);
std::cout << '\n';
// Denominator of b
auto denom_b_base = std::stringstream("1");
auto denom_b_exp = -141;
auto denom_b = construct(denom_b_base, denom_b_exp);
std::cout << "denom_b = ";
CGAL::print(std::cout, denom_b);
std::cout << '\n';
auto a = Quotient<Gmpzf>(num_a, denom_a), b = Quotient<Gmpzf>(num_b, denom_b);
auto a_d = CGAL::to_double(a), b_d = CGAL::to_double(b);
std::cout << "a = " << a << ", b = " << b << '\n';
std::cout << "double(a) = " << a_d << ", double(b) = " << b_d << '\n';
std::cout << "a > b: " << (a > b ? "True" : "False") << '\n';
std::cout << "a_d >= b_d: " << (a_d >= b_d ? "True" : "False") << '\n';
}Output
num_a = 495465331884540240762104420639278860096116250934219950844827973437034498625*2^-379
denom_a = 2245694908428994193174821578352066558595985337089*2^-299
num_b = 3447327498334006902041169*2^-215
denom_b = 1*2^-141
a = 4.0238790231336247e-40/2.2048652400042302e-42, b = 6.5468663604502424e-41/3.5873240686715317e-43
double(a) = 182.49999819154047, double(b) = 182.4999981915405
a > b: True
a_d >= b_d: False
Environment
- Operating system: Ubuntu 24.04
- Compiler: g++ 13.3.0
- Release or debug mode: Both release and debug
- Specific flags used (if any): -frounding-math
- CGAL version: 6.0.1
- Boost version: 1.88.0
- Other libraries versions if used (Eigen, TBB, etc.): Eigen 3.4.0
Ideally to_double would return the nearest double, and then it would naturally be monotonic. The current implementations for Quotient and Quotient<Gmpzf> are not very good. If speed is not an issue, I think I would convert the Quotient to a rational (Gmpq, mpq_class or mpq_rational), by converting the numerator and denominator to rationals and dividing them, and then apply to_double to that.
By the way, is there a particular reason why you are using Gmpzf? Mpzf has similar functionality and is faster.
Thank you for the information. Is there a specific reason that Mpzf is faster than Gmpzf? I was able to resolve my problem for the moment by constructing an mpq_t from the mantissa of the numerator and denominator, and then using mul_2exp on the difference of exponents. It does not look like Mpzf currently provides access to the mantissa and exponent separately.