These are chat archives for symengine/symengine

8th
Jun 2015
Isuru Fernando
@isuruf
Jun 08 2015 13:04
@certik, for asin(real_double(2.0)) do we want it to return a ComplexDouble or a nan?
Ondřej Čertík
@certik
Jun 08 2015 14:05
Let me think. I personally do not like nan.
So there are two ways around. Either they are in a complex domain, then the answer is something like 1.57079-i*1.3169, or they are in real domain, then it should raise ValueError: math domain error like in Python.
By default, I would suggest a complex domain. That's the default assumption in SymPy as well as SymEngine, that everything is complex.
Unless told otherwise by a user.
Isuru Fernando
@isuruf
Jun 08 2015 14:13
Okay.
This complicates the code for EvaluateDouble and when we have it for MPFR and MPC it's going to be even more complicated.
Ondřej Čertík
@certik
Jun 08 2015 14:17
So for EvalDouble, there I think the user is telling it to use reals, doesn't it?
So there I would just raise an exception.
Because we want it to be fast, when the user knows it is real.
For eval_complex_double, there I would evaluate in the complex plane.
Isuru Fernando
@isuruf
Jun 08 2015 14:18
Yes. I mean for EvaluateDouble which does the automatic simplification
Ondřej Čertík
@certik
Jun 08 2015 14:19
Where is the code?
src/real_double.cpp?
Isuru Fernando
@isuruf
Jun 08 2015 14:19
yes
I just changed the code, let me send a WIP PR
Ondřej Čertík
@certik
Jun 08 2015 14:20
I still don't understand what the problem is. We just want EvaluateRealDouble to raise an exception (should be simple) and EvaluateComplexDouble to evaluate in complex plane (should be simple as well).
Isuru Fernando
@isuruf
Jun 08 2015 14:21
I thought you wanted it to return complex for asin(real_double(2.0))?
Ondřej Čertík
@certik
Jun 08 2015 14:22
It depends. If asin is the general function in SymEngine, then yes, asin(2) should be complex.
Or, to be more specific, asin(2) should be left as is, and asin(2.0) can be evaluated, as complex.
In other words, we can use EvaluateComplexDouble by default, unless the user tells us, then we should use EvaluateRealDouble.
Isuru Fernando
@isuruf
Jun 08 2015 14:23
That's what's in #466
Ondřej Čertík
@certik
Jun 08 2015 14:24

So the EvaluateComplexDouble contains code like:

+    virtual RCP<const Basic> asin(const Basic &x) const {
+        SYMENGINE_ASSERT(is_a<ComplexDouble>(x))
+        return number(std::asin(static_cast<const ComplexDouble &>(x).i));
+    }

that's not complicated, is it?

That's the code that would be used by default.
Isuru Fernando
@isuruf
Jun 08 2015 14:24
EvaluateRealDouble on the other hand
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)));
        }
    }
Ondřej Čertík
@certik
Jun 08 2015 14:25

If the user specifically requests a real evaluation, then this code gets triggered:


+    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)));
+        }
+    }

Which must test the domain, there is no way around it. You can't use std::asin(2.0) (I am pretty sure, correct?). You need to either use std::asin(std::complex<double>(2)) or you just raise a math domain error.

So how would you simplify this whole thing? I don't think there is any way around this.
Isuru Fernando
@isuruf
Jun 08 2015 14:28

The only way to simplify this whole thing, would be to call this function.

RCP<const Number> eval_double(const Basic &b) {
    double d = eval_real_double(b);
    if (isnan(d)) {
        std::complex<double> comp = eval_complex_double(b);
        return complex_double(comp);
    }
    return real_double(d);
}

which is one of the methods we tried earlier for simplification of RealDoubles, but is slower.

Ondřej Čertík
@certik
Jun 08 2015 14:29
I don't like this to be honest.
Isuru Fernando
@isuruf
Jun 08 2015 14:30
What do you suggest for it?
I mean we'll eventually need it for .n() method right?
Ondřej Čertík
@certik
Jun 08 2015 14:31
I like the following workflow: everything is complex by default, and the eval_complex_double is very simple inside, you just evaluate it in complex domain. Everything works, you only get exceptions for things like log(0). So RCP<const Number> eval_double(const Basic &b) just calls eval_complex_double and that's it. Then, if the user is concerned about speed, and the problem turns out to work in real domain only, only then the user calls eval_real_double(). That function will raise exceptions for things like asin(2).
As to .n(), SymPy uses complex by default as well:
In [1]: asin(2)
Out[1]: asin(2)

In [2]: asin(2).n()
Out[2]: 1.5707963267949 - 1.31695789692482⋅ⅈ
Is there anything you don't like on this approach? The code is actually much simpler in complex domain. It's just that I expect std::asin(std::complex<double>(2)) to be slower than std::asin(2.0), so we want the user to allow the faster option, even if less general (and could raise exceptions).
Isuru Fernando
@isuruf
Jun 08 2015 14:33
How does sympy know when you have to work in the real domain, like
In [6]: sin(2).n()
Out[6]: 0.909297426825682
and not to return 0.909297426825682 + 0.0000000000000*I?
Ondřej Čertík
@certik
Jun 08 2015 14:38
I am not sure, I would think it chops small parts. The code is in sympy/core/evalf.py.
It looks to me that it evaluates in complex plane, and if imaginary part is zero within the requested precision, then it returns a real number.
Isuru Fernando
@isuruf
Jun 08 2015 14:52
Thanks, will do it like that.

I've been writing a class RealMPFR.

class RealMPFR : public Number {
public:
    mpfr_t i;
    mpfr_rnd_t rnd_;
}

Is the following okay?

RCP<const Number> RealMPFR::addreal(const Integer &other) const {
    RealMPFR* t = new RealMPFR(rnd_);
    mpfr_add_z(t->i, i, other.i.get_mpz_t(), rnd_);
    return rcp(t);
}
Ondřej Čertík
@certik
Jun 08 2015 14:56
No, you should never use naked new like that.
Rewrite it as follows:
mpfr_t i2;
mpfr_add_z(i2, i, other.i.get_mpz_t(), rnd_); 
return rcp(new RealMPFR(rnd_, i2));
(Also in SymEngine, the instances are immutable, so you should not be changing them once you construct them.) You can use a move constructor to avoid unnecessary copies of i2. I guess it's not a move constructor, just a constructor, that has rvalue references of i2 (&&).
(It would be a move constructor if it has rvalue of the RealMPFR class.)
Isuru Fernando
@isuruf
Jun 08 2015 15:00
Yeah, I didn't want to change it after construction, but couldn't find a way to avoid copying
Ondřej Čertík
@certik
Jun 08 2015 15:00
Isn't mpfr_t just a pointer anyway? So it is cheap to copy.
I assume you must initialize it first then, don't you?
Isuru Fernando
@isuruf
Jun 08 2015 15:01
Yes.
typedef struct {
  mpfr_prec_t  _mpfr_prec;
  mpfr_sign_t  _mpfr_sign;
  mpfr_exp_t   _mpfr_exp;
  mp_limb_t   *_mpfr_d;
} __mpfr_struct;

typedef __mpfr_struct mpfr_t[1];
Ondřej Čertík
@certik
Jun 08 2015 15:02
I would document, that if you must pass initialized mpfr_t into the RealMPFR constructor, and that the class will take care of deallocating it.
It is identical to mpz, the only difference is that in there, we use mpz_class, that takes care of this initialization and deallocation. For mpfr_t, we have two options. We either do it manually as I did above, or we introduce a simple mpfr_class, that will initialize and deallocate it for us. I am actually leaning towards the mpfr_class. We do not need to overload any fancy operators and do some template trickery like mpz_class does, just look how it is implemented, it is a mess --- and it all for just allowing things like z=a+b+c without an overhead. For mpfr_class, we just allocate/deallocate in constructor/destructor, and then provide access to the managed mpfr_t. Then we use the mpfr C interface as before.
Here is how it would work:
mpfr_class i2;
mpfr_add_z(i2.get_mpfr_t(), i, other.i.get_mpz_t(), rnd_); 
return rcp(new RealMPFR(rnd_, i2));
Then the constructor of RealMPFR only accepts the mpfr_class, just like we handle the mpz_class.
Ondřej Čertík
@certik
Jun 08 2015 15:10
Now, if we want to use it like an int or a double, then we would have to also override the arithmetic operators +, -, * and /. But that's where it becomes complicated, because you need to allocate a new mpfr_t for the result, and so if you do things like z=a+b+c, then you are wasting the intermediate temporaries. For this reason, it's cleaner if we simply call the C interface by hand, and do the most efficient way for the given problem. Then there is no waste and the code is reasonably simple.
The other option is to use one of the many C++ interfaces to MPFR, see http://www.mpfr.org/ for a list. But I think a simple mpfr_class like I described above is a very clean solution, without an additional dependency on another library (that might or might not be synchronized with the latest MPFR).
Ondřej Čertík
@certik
Jun 08 2015 15:19
As to the real/complex thing, I looked at Arb, and it also has two modes. A complex mode, which just evaluates everything: http://fredrikj.net/arb/acb.html#c.acb_atan, and a real mode, which returns an "indeterminate interval" (which is just a fancy name for a NaN, see the top of the file): http://fredrikj.net/arb/arb.html#c.arb_asin. So Arb returns a NaN. I would just raise an exception. There are pros and cons, whether you should return a NaN, or raise an exception. The problem with a NaN is that you do not know where it came from. A nice exception tells you exactly where it came from --- and you can still catch it and act accordingly, if that is part of an algorithm. So perhaps the argument for a NaN is speed (it might be faster than an exception), but I think that NaN is kind of an exception anyway, so I personally think exceptions are better for this, unless we can find some convincing argument for NaN. Let me ask Fredrik.
Isuru Fernando
@isuruf
Jun 08 2015 15:22

Exception handling is not done in eval_mpfr, eval_arb, etc.

void bvisit(const Add &x) {
        mpfr_t t;
        mpfr_init2(t, mpfr_get_prec(result_));

        auto d = x.get_args();
        auto p = d.begin();
        apply(result_, *(*p));
        p++;
        for (; p != d.end();  p++) {
            apply(t, *(*p));
            mpfr_add(result_, result_, t, rnd_);
        }
        mpfr_clear(t);
    }

Here t is cleared only if there is no exception. We should do that as well.

Ondřej Čertík
@certik
Jun 08 2015 15:23
That's why we should use the mpfr_class, as I described above. Then it will get deallocated if there is an exception.

Here is how the code would look like:

void bvisit(const Add &x) {
        mpfr_class t(mpfr_get_prec(result_));

        auto d = x.get_args();
        auto p = d.begin();
        apply(result_, *(*p));
        p++;
        for (; p != d.end();  p++) {
            apply(t.get_mpfr_t(), *(*p));
            mpfr_add(result_, result_, t.get_mpfr_t(), rnd_);
        }
    }

Much simpler, safer (no leaks) and still as performant.

Sumith Kulal
@Sumith1896
Jun 08 2015 15:40
I am stuck at FindBOOST.cmake
I looked at Piranha but will a similar one do?
It doesn't use libfind.
Ondřej Čertík
@certik
Jun 08 2015 15:44
I would just find BOOST using our libfind, just like any other library.
Sumith Kulal
@Sumith1896
Jun 08 2015 15:44
But it gives an error. I had pasted that earlier.
Let me get it.
Ondřej Čertík
@certik
Jun 08 2015 15:44
Thanks.
Isuru Fernando
@isuruf
Jun 08 2015 15:46
That's because boost libraries I have are like this,
/usr/lib/x86_64-linux-gnu/libboost_date_time.so.1.55.0
/usr/lib/x86_64-linux-gnu/libboost_filesystem.so.1.55.0
/usr/lib/x86_64-linux-gnu/libboost_iostreams.so.1.55.0
/usr/lib/x86_64-linux-gnu/libboost_system.so.1.55.0
I think it's the same with @Sumith1896
Ondřej Čertík
@certik
Jun 08 2015 15:47
I think there are multiple boost libraries, so you have to look for boost_date_time instead of just boost.
But, I think there might be a lot more logic involved. We might want to just use FindBOOST from cmake itself.
Sumith Kulal
@Sumith1896
Jun 08 2015 15:48
Okay, for that to work I have to remove my FindBOOST.cmake right?
Ondřej Čertík
@certik
Jun 08 2015 15:50
I see, I 100% agree with what you wrote.
@Sumith1896 try to spend some more time on this, and ask questions. This is a good thing to learn anyway, how to use cmake and interface with libraries, and Boost is pretty commonly used. Just ask when you get stuck, so that we can resolve it eventually.
Sumith Kulal
@Sumith1896
Jun 08 2015 15:53
Agreed, thanks
Ondřej Čertík
@certik
Jun 08 2015 16:08
@isuruf what is your opinion about the mpfr_class? Do you think it's a good idea?
Isuru Fernando
@isuruf
Jun 08 2015 16:09
Yes, I think it saves us a lot of trouble when it comes to exceptions, memory management etc.
Ondřej Čertík
@certik
Jun 08 2015 16:09
Ok. Let me know if you need any help with implementing it.
Isuru Fernando
@isuruf
Jun 08 2015 16:12
Thanks
Copy constructor of mpfr_class would copy the mpfr_t right?
Ondřej Čertík
@certik
Jun 08 2015 16:25
Well, it would need to call the MPFR C functions to make a copy of the mpfr_t.
Isuru Fernando
@isuruf
Jun 08 2015 16:25
Yes
Isuru Fernando
@isuruf
Jun 08 2015 16:57
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 = 0;
        mpfr_swap(mp, other.get_mpfr_t());
    }
    ~mpfr_class() {
        if (mp->_mpfr_d != 0) {
            mpfr_clear(mp);
        }
    }
};
Ondřej Čertík
@certik
Jun 08 2015 19:18
You should also define assignment and move operators (operator=)
Otherwise I think that this looks good.
Ondřej Čertík
@certik
Jun 08 2015 19:34
The only questionable thing is the mp->_mpfr_d = 0; usage. I understand that you have to somehow specify that it is not initialized for the move constructor. I think it's fine. @bluescarni, what do you think?
Francesco Biscani
@bluescarni
Jun 08 2015 22:43
@certik I think it's good... might use nullptr instead just to be more C++11ish, but it's a minor nitpick
The thing is that no matter how you do it, you will have to touch the internals of mpfr_t in order to enable move semantics
Another minor thing is that I usually prefer to just swap the members of mpfr_t manually, in case for some reasons mpfr_swap() is not doing just a plain swap (e.g., running some checks in debug mode about the consistency of a null _mpfr_d in an mpfr_t that has size nonzero )
but probably not a big deal in practice