CGAL/cgal

`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()].

  1. For my application, it is very important that I preserve monotonicity when converting from Quotient<Gmpzf> to double. 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?
  2. 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.