GSoC 2015 : Week 3

This week, I continued my work from the previous two weeks to complete the floating point support for SymEngine.

RealDouble was the wrapper for double used in SymEngine and it used to return NaN for numerical evaluations that resulted in a complex number. Checks were added to ensure that if the result was complex, a ComplexDouble was returned. So now, `asin` for RealDouble looks like this

    virtual RCP<const Basic> asin(const Basic &x) const {
        SYMENGINE_ASSERT(is_a<RealDouble>(x))
        double d = static_cast<const RealDouble &>(x).i;
        if (d <= 1.0 && d >= -1.0) {
            return number(std::asin(d));
        } else {
            return number(std::asin(std::complex<double>(d)));
        }
    }

Here, when the result is a double, a RealMPFR is returned, while when it is complex, RealDouble is converted to std::complex<double>, evaluated and then a ComplexDouble is returned.

Implementing RealMPFR took more time than I thought, because of several design decisions that were needed to be taken.

Basic classes in SymEngine are immutable classes. So, when constructing a RealMPFR class, the mpfr_t value has to be passed into the RealMPFR's constructor. Which means when working with RealMPFR, we would first construct a mpfr_t, initialize it, set a value to the mpfr_t and then construct the RealMPFR. Therefore when exceptions are raised and caught, memory leaks happen, because the mpfr_t is not managed. This is also in EvalMPFR class where if we had an expression like 3 + x that we want to evaluate numerically, but when trying to evaluate x numerically, an exception is raised and a memory leak happens. Solution was to implement a managed mpfr_class that acts like mpz_class. 

To avoid adding a new dependency, I implemented a simple mpfr_class that would manage the mpfr_t inside it. 

class mpfr_class {
private:
    mpfr_t mp;
public:
    mpfr_ptr get_mpfr_t() { return mp; }
    mpfr_srcptr get_mpfr_t() const { return mp; }

    mpfr_class(mpfr_t m) {
        mpfr_init2(mp, mpfr_get_prec(m));
        mpfr_set(mp, m, MPFR_RNDN);
    }
    mpfr_class(mpfr_prec_t prec = 53) {
        mpfr_init2(mp, prec);
    }
mpfr_class(const mpfr_class& other) {
        mpfr_init2(mp, mpfr_get_prec(other.get_mpfr_t()));
        mpfr_set(mp, other.get_mpfr_t(), MPFR_RNDN);
    }
    mpfr_class(mpfr_class&& other) {
        mp->_mpfr_d = nullptr;
        mpfr_swap(mp, other.get_mpfr_t());
    }
    mpfr_class& operator=(const mpfr_class& other) {
        mpfr_set_prec(mp, mpfr_get_prec(other.get_mpfr_t()));
        mpfr_set(mp, other.get_mpfr_t(), MPFR_RNDN);
        return *this;
    }
    mpfr_class& operator=(mpfr_class&& other) {
        mpfr_swap(mp, other.get_mpfr_t());
        return *this;
    }
    ~mpfr_class() {
        if (mp->_mpfr_d != nullptr) {
            mpfr_clear(mp);
        }
    }
};

To add the move constructors and destructors, it was needed to know if a mpfr_class has been initialized or not. This was achieved by setting _mpfr_d pointer to null when move constructor is called to avoid adding a data member to mpfr_class. Thanks to +Ondřej Čertík  and +Francesco Biscani  for the ideas and comments. 

Similarly an mpc_class was introduced to manage a mpc_t object. Another decision taken was to use a default rounding mode for rounding when operations on RealMPFR's and ComplexMPC's are done.

When adding two RealMPFR's of different precisions, it was decided to promote the RealMPFR with low precision to higher precision. An example of this is given below.

RCP<const Number> RealMPFR::addreal(const RealMPFR &other) const {
    mpfr_class t(std::max(get_prec(), other.get_prec()));
    mpfr_add(t.get_mpfr_t(), i.get_mpfr_t(), other.i.get_mpfr_t(), MPFR_RNDN);
    return rcp(new RealMPFR(std::move(t)));
}

To add SymEngine into Sage, I had to add CMake into Sage as well, since CMake is a dependency of SymEngine. Problem with including CMake was that it fails to build with OS X 10.10 due to a bug in a header installed in OS X 10.10 (/usr/include/dispatch/object.h). This header does not give a fault when compiling with clang but with gcc which is included in Sage this gives an error. This seemed to be a problem in building bundled curl in cmake, but realized that it was not so afterwards. After building curl and linking with it in cmake, I got the same error on OS X. I got access to OS X 10.10 only this week, but will try next week as time permits.

For next week, I am going to implement the wrappers for SymEngine for Sage. (i.e. conversions to and from Sage).