These are chat archives for symengine/symengine

1st
Jul 2014
Ondřej Čertík
@certik
Jul 01 2014 16:38
@thilinarmtb ping me when you update the PR
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 16:38
@certik: Sure, I'll do that.
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 16:44
@certik : If we are to make the exponent positive in an expression like x*2^(-5/3) we need to find the ceiling or floor of -5/3. Can this be done with gmp, finding ceiling or floor of a Rational? There is a ceil function in the mpf class, but I couldn't find anything else.
Ondřej Čertík
@certik
Jul 01 2014 16:45
Why cannot you just check if it is negative?
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 16:48
Yes, I can but x*2^(-5/3) should be represented as (1/4)*x*2^(1/3) right? Not as (1/2^(5/3))*x because latter causes negative non integer exponents in coef
To do that simplification, I think we have to find the ceiling or floor of the negative exponent, or am I wrong?
or the quotient of the rational
Ondřej Čertík
@certik
Jul 01 2014 16:57
I need to think about this a bit
essentially, the goal of the transformation is to always normalize 1/sqrt(n) to sqrt(n)/n
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 16:58
Okay, finding ceiling, floor or quotient will be little bit expensive right?
Ondřej Čertík
@certik
Jul 01 2014 16:58
we just need to figure out the proper transformation for a general 1/n^(a/b)
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 16:58
Yes, that's right.
In the general case, is there any other work around without finding the quotient?
Ondřej Čertík
@certik
Jul 01 2014 16:59
I am looking into it
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:00
If we try to do it manually, as we do in normal arithmetic simplification, what we do inside our brain is that, we multiply 1/n^(a/b) by a power n^t/n^t so that a/b + t is an integer.
Ondřej Čertík
@certik
Jul 01 2014 17:02
Let me see
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:03
Okay
Ondřej Čertík
@certik
Jul 01 2014 17:03
So for a=1, b=2, we get t=1/2
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:03
yes
Ondřej Čertík
@certik
Jul 01 2014 17:03
Obviously, in all this n>0
For n<0, one needs to extract the imaginary unit "i" first, then follow the above rules.
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:04
it need not be true always, I guess. For example we can have something like (-3)^7/3
Ondřej Čertík
@certik
Jul 01 2014 17:05
for (-3)^7/3, you have to extract the imaginary unit appropriately
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:05
if n is negative then b can't be even
you just need to be careful about the principal branch in the complex plane.
anyway, for now we should raise exception if n < 0, as I think we already do
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:07
Of course, I was talking in a real context, in general yes.
Ondřej Čertík
@certik
Jul 01 2014 17:07
(sympy and csympy must work for complex numbers.)
Ok, I agree that we need a/b+t = integer
but there are multiple solutions obviously
I think we probably want "t" to be between 0 and 1?
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:09
Sure, but is there a imaginary unit in (-3)^(1/3) ? it's -1*3^(1/3) isn't it?
Ondřej Čertík
@certik
Jul 01 2014 17:10
I think there is imaginary unit in it: http://www.wolframalpha.com/input/?i=%28-3%29^%281%2F3%29
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:10
Yes, we need to do this as cheaply as possible
Yes, I guess there is.
Ondřej Čertík
@certik
Jul 01 2014 17:12
i.e. (-1)^(1/3) = exp(i*pi)^(1/3) = exp(i*pi/3) = 1/2 + i*sqrt(3)/2
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:14
Hmm, yes. thing is that x^3 + 1 = 0 has one real and two complex solutions, once again I was only thinking about the reals, my bad.
Ondřej Čertík
@certik
Jul 01 2014 17:19
It depends how you define the principal branch
the solutions to x^3+1=0 are of the form x=exp(i*pi + 2*pi*i*n)^(1/3)=exp(i*pi/3 + 2/3*pi*i*n) where n is any integer
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:21
Yes, I get it now, thanks !!
Ondřej Čertík
@certik
Jul 01 2014 17:21
so if you set n=1, you get the real solution
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:21
Yes.
So sage does the simplification as we expect: http://www.sagenb.org/home/pub/5048
I'll try to go through their code and see what they do.
Ondřej Čertík
@certik
Jul 01 2014 17:23
( Public worksheets are currently disabled. )
Btw, you should use cloud.sagemath.com
you can install csympy there and we can even collaborate on it.
If you send me your login there, I can add you to my project
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:23
Ohh, Sorry for that. I will use it, thanks.
Sure
Ondřej Čertík
@certik
Jul 01 2014 17:25
No, you have a good point about the branch cuts --- I can't find the reference for it, I just know that if you do x^(a/b), then you need to write "x" as exp(i*pi*phi), where 0 <= phi < 2*pi (this selects the branch cut) and then you just power this as usual to a/b, or to any complex number, and it will work
As far as Sage, the simplification is done in GiNaC, by this code in ginac/power.cpp:
            // ^(c1,n/m) -> *(c1^q,c1^(n/m-q)), 0<(n/m-q)<1, q integer
            if (basis_is_crational && exponent_is_crational
                && num_exponent->is_real()
                && !num_exponent->is_integer()) {
                const numeric n = num_exponent->numer();
                const numeric m = num_exponent->denom();
                numeric r;
                numeric q = iquo(n, m, r);
                if (r.is_negative()) {
                    r += m;
                    --q;
                }
                if (q.is_zero()) {  // the exponent was in the allowed range 0<(n/m)<1
                    if (num_basis->is_rational() && !num_basis->is_integer()) {
                        // try it for numerator and denominator separately, in order to
                        // partially simplify things like (5/8)^(1/3) -> 1/2*5^(1/3)
                        const numeric bnum = num_basis->numer();
                        const numeric bden = num_basis->denom();
                        const numeric res_bnum = bnum.power(*num_exponent);
                        const numeric res_bden = bden.power(*num_exponent);
                        if (res_bnum.is_integer())
                            return (new mul(power(bden,-*num_exponent),res_bnum))->setflag(status_flags::dynallocated | status_flags::evaluated);
                        if (res_bden.is_integer())
                            return (new mul(power(bnum,*num_exponent),res_bden.inverse()))->setflag(status_flags::dynallocated | status_flags::evaluated);
                    }
                    return this->hold();
                } else {
                    // assemble resulting product, but allowing for a re-evaluation,
                    // because otherwise we'll end up with something like
                    //    (7/8)^(4/3)  ->  7/8*(1/2*7^(1/3))
                    // instead of 7/16*7^(1/3).
                    ex prod = power(*num_basis,r.div(m));
                    return prod*power(*num_basis,q);
                }
            }
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:31
Thanks Ondrej, I'll look in to the code.
Here is my login for cloud.sagemath.org: thilinarmtb@gmail.com
            if (basis_is_crational && exponent_is_crational
                && num_exponent->is_real()
                && !num_exponent->is_integer()) {
                const numeric n = num_exponent->numer();
                const numeric m = num_exponent->denom();
                numeric r;
                numeric q = iquo(n, m, r);
                if (r.is_negative()) {
                    r += m;
                    --q;
                }
Here I think they do what we were talking about.
Ondřej Čertík
@certik
Jul 01 2014 17:38
Right. I think the rule is: n^(a/b) = n^t * n^(a/b-t), and you choose integer t such that 0 < a/b-t < 1. Here n > 0.
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:39
yes, t is an integer here right?
Ondřej Čertík
@certik
Jul 01 2014 17:39
Right.
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:39
Okay
So I'll try to update this in the current PR.
Ondřej Čertík
@certik
Jul 01 2014 17:42
I see, so t = floor(a/b) using Python's floor.
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:43
Yes.
Ondřej Čertík
@certik
Jul 01 2014 17:52
The iquo seems to just divide n/m with a remainder r. So n = q*m+r, and they make sure that r is positive (they don't use gmp). In gmp, r is always positive, per their documentation: https://gmplib.org/manual/Integer-Division.html#Integer-Division
So I would use the function mpz_cdiv_qr
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:53
Yes, thanks I'll use it.
Ondřej Čertík
@certik
Jul 01 2014 17:55
that gives you q and r. Then I think t = r/b and a/b-t=q.
So I think the algorithm is:
mpz_cdiv_qr(q, r, a, b)
return n^q * n^(r/b)
This message was deleted
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 17:56
Yes
Yes, if both exp and base are numbers and exp is a negative integer, we do the above
Ondřej Čertík
@certik
Jul 01 2014 17:59
no, if exp is just negative integer, you leave things be
you do the above if exp is any Rational
i.e. a^(3/2) -> a*a^(1/2), I think we want to do this as well
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:00
Oh yes, sorry for the mistake. We can do this for positive exp also.
Yes, exactly.
Ondřej Čertík
@certik
Jul 01 2014 18:02
I got it backwards above, I think it's t=q and a/b-t=r/b. The algorithm I wrote above should be correct.
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:05
Yes, I think the confusion was due to that in my original chat I used t for the fractional part and then we again use it for q. Algorithm is fine because we only need q and r/b.
Ondřej Čertík
@certik
Jul 01 2014 18:07
So I would call mpz_cdiv_qr(q, r, a, b), then check if q == 0 (if so, then we don't do anything, as then r/b=a/b, i.e. it's already normalized). If q != 0, then you return n^q * n^(r/b). I would implement this simplification in a new PR, unrelated to your fixes. We need to update Pow::canonicalize to make sure that number^rational is always normalized. Then write tests for this. Merge it, and then finish the PR, which should be now easy hopefully.
I like this, this is a well defined, cheap simplification, and it would take care of 1/sqrt(2) vs sqrt(2)/2, which is a very common simplification. It will not take care of things like sqrt(12)/sqrt(6) -> sqrt(2), but that is fine. We'll simply create a function powsimp (or some better name), that will do those.
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:11
Yes, l agree. So this sympy/csympy#215 will have to wait till the new PR gets merged, right?
Ondřej Čertík
@certik
Jul 01 2014 18:11
Also things like 5^(3/2)/5 will then automatically get simplified to 5^(1/2).
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:12
Yes, because 5/5 cancels in coef that's great !!!
(y)
Ondřej Čertík
@certik
Jul 01 2014 18:13
So this takes care of 90% of all the cases, and for the rest one needs to call some function, that will do integer factorization.
As to #215, you can actually finish it too, but don't do any of this simplification there (those are two separate issues). As noted, you will have to modify one of the sec tests (or which one it was).
Sushant Hiray
@sushant-hiray
Jul 01 2014 18:14
So @thilinarmtb you can just comment out the trignometric cases which fail, I'll send a followup PR to fix those
In this way #215 should get merged
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:15
Ahh, yes I'll modify it. Could you elaborate on the other 10% where we need integer factorization?
@sushant-hiray : Without altering the sec test? Do you mean as it is now?
Ondřej Čertík
@certik
Jul 01 2014 18:16
I.e. if #215 gets merged first, we need to modify one of the tests, to handle the 2/sqrt(2) vs sqrt(2) issue. Then, after a new PR that handles these rational exponents normalization (canonicalization), we can (should) change that test back
@thilinarmtb, this test needs to be modified https://github.com/sympy/csympy/pull/215#issuecomment-47274387
It's because the changes in #215 stop simplifying 2/sqrt(2) to sqrt(2).
Sushant Hiray
@sushant-hiray
Jul 01 2014 18:18
so for now @thilinarmtb you can just change 2/sqrt(2) to sqrt(2) in that case
Ondřej Čertík
@certik
Jul 01 2014 18:18
And then the new PR enables this simplification again, and many others.
@sushant-hiray exactly
Sushant Hiray
@sushant-hiray
Jul 01 2014 18:19
yes, also once the complex module is shipped in, we can start checking for negative rational case
Ondřej Čertík
@certik
Jul 01 2014 18:19
exactly
and negative integers as well
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:19
Yes, I'll modify it. I think that's the first test which fails, Not sure I remember correctly, but I'll have to change one or two other test cases as well.
Ondřej Čertík
@certik
Jul 01 2014 18:19
If it's just a few tests, then that's what I would do.
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:19
Great !!!
Yes, not an issue, I'll fix it.
Ondřej Čertík
@certik
Jul 01 2014 18:20
Also, don't forget to incorporate the expand() patch, that I posted in the PR (https://github.com/sympy/csympy/pull/215#issuecomment-47289910)
and tests for it
Sushant Hiray
@sushant-hiray
Jul 01 2014 18:21
so I had purposefully put those tests: as in my hashmap i had sqrt(2), so instead of just checking for sqrt(2) again, i checked for 2/sqrt(2). There might be a couple more tests which could fail because of the same thing
Ondřej Čertík
@certik
Jul 01 2014 18:22
Right. It's just a temporary thing, as we'll enable this again soon.
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:23
Sure Ondrej
Ondřej Čertík
@certik
Jul 01 2014 18:23
I am strongly suspecting that we might not need the old dict_add_term anymore, the new one should be used everywhere
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:24
Hmm, I'll try replacing every dict_add_term with dict_add_term_new
Ondřej Čertík
@certik
Jul 01 2014 18:25
We need to think about each case and for each case construct a test case, that fails in current master, just like I did for expand()
And we also need to think about each case if some code that calls it can't be simplified, because we don't need the "for loop" anymore.
As to cloud.sage, actually, you can add me (ondrej.certik@gmail.com) to your project
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:26
Okay, I'll add you.
Ondřej Čertík
@certik
Jul 01 2014 18:32
Thanks @thilinarmtb for working on this. Both PRs will be very important fixes.
Thilina Rathnayake
@thilinarmtb
Jul 01 2014 18:33
Sure @certik , It's my pleasure.