[GiNaC-list] Floating to Rational

[Vladimir V. Kisil]

>>>>> On Wed, 19 Jul 2023 10:42:35 +0100, "Vladimir V. Kisil" <V.Kisil at leeds.ac.uk> said:
    VVK>     Yes, see the code below as an example.

    It is easy to see that rationalize function from CLN produces huge
  fractions due to floating point errors. I am including a second method
  based on the continued fraction approximations. Its behaviour can be
  adjusted through parameters tolerance and cf_size.
-- 
Vladimir V. Kisil                 http://www1.maths.leeds.ac.uk/~kisilv/
  Book:      Geometry of Mobius Maps       https://doi.org/10.1142/p835
  Soft:      Geometry of cycles         http://moebinv.sourceforge.net/
  Jupyter notebooks:        https://github.com/vvkisil?tab=repositories

#include <iostream>
#include <ginac/ginac.h>
#include <math.h>   
using namespace std;
using namespace GiNaC;

int tolerance = 1E6;
int cf_size = 15;

numeric rationalize(numeric x) {
	if (ex_to<numeric>(x).is_rational())
		return x;
	x = ex_to<numeric>(evalf(x));
	int sign = (x.is_positive() ? +1 : -1);
	double xf= abs(x.to_double());
    lst coeff = lst{};
	for (int i = 0; i < cf_size; ++i ) {
        if ( xf > tolerance )
            break;
        else {
			int base = floor(xf);
            coeff.append(numeric(base));
			 if (xf - base == 0)
				 break;
            xf = 1/(xf-base);
		}
	}
	
	lst::const_reverse_iterator i = coeff.rbegin();

	numeric res = ex_to<numeric>(*i);
	++i;
	while ( i != coeff.rend()) {
		res =  ex_to<numeric>(*i) + 1/res;
		++i;
	}
	return sign*res;
}
	
struct map_rationalize_ex : public map_function {
         ex operator()(const ex &e)
         {
			 if (is_a<numeric>(e))
				 return rationalize(ex_to<numeric>(e).real()) + I * rationalize(ex_to<numeric>(e).imag());
			 else if (e.nops() > 0)
				 return e.map(*this);
			 else
				 return e;
		 }
};

int main()
{
	symbol x("x"), y("y");

	ex e = 0.5*sin((-.6+I*.3)*x +exp((.25-I*.8)*y+.3333-I*numeric(1,7)))+pow(.4*x, .6);

	map_rationalize_ex rationalize_ex;

	cout << e << endl << "is rationalised to" << endl << rationalize_ex(e) << endl;
	
	return 0;
}