New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Further fixes to integrate Complex module with CSymPy #248

Merged
merged 29 commits into from Jul 23, 2014

Conversation

Projects
None yet
2 participants
@sushant-hiray
Copy link
Contributor

sushant-hiray commented Jul 12, 2014

  • add
  • mul
  • pow
  • tests
  • expand
@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 12, 2014

+1 so far.

@certik

This comment has been minimized.

Copy link

certik commented on src/mul.cpp in 563e1f7 Jul 13, 2014

You forgot ).

@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 16, 2014

I checked into how SymPy handles complex powers and complex numbers raised to number.

In [1]: from sympy import *

In [2]: a = 3 + 2*I

In [3]: b = 4

In [4]: a**b
Out[4]: (3 + 2*I)**4

In [5]: b**a
Out[5]: 4**(3 + 2*I)

In [6]: (a**b).expand()
Out[6]: -119 + 120*I

In [7]: (b**a).expand()
Out[7]: 64*4**(2*I)

Since by default, the complex powers are not evaluated and also complex numbers are not evaluated, so the current changes I made should be sufficient.
I will add tests to check the same.

src/pow.cpp Outdated
@@ -151,9 +153,15 @@ RCP<const Basic> pow(const RCP<const Basic> &a, const RCP<const Basic> &b)
map_basic_basic surd;
surd[exp_new] = div(integer(r), integer(den));
return rcp(new Mul(frac, std::move(surd)));
} else if (is_a<Complex>(*a)) {
// This might raise some issues.

This comment has been minimized.

@certik

certik Jul 16, 2014

Contributor

What kind of issues?

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 17, 2014

Author Contributor

I had written tests and checked how SymPy handles it so I wasn't sure about it then. I'm just tweaking the part in mul a bit so we can perform some manipulations directly.

r2 = mul(mul(I, i2), x);
r3 = mul(mul(I, i3), x);
r2 = add(r1, r2);
assert(eq(r3, r2));

This comment has been minimized.

@certik

certik Jul 16, 2014

Contributor

This looks good. Also check div, sub and pow. And try the expand() with complex numbers, I bet there are a few bugs there.

sushant-hiray added some commits Jul 17, 2014

Updated dict_add_term_new to handle complex coefficients
Now dict doesn't complex terms with power 1 or -1.
These terms are directly incorporated into the coeff.
Rest additions are dependent functions for the same.
Added tests for complex multiplication and subtraction
Tests have been added to check if the complex is multiplied with
coeff instead of adding to dict when its power is 1.
Added tests for complex div. Fixed some issues.
* Fixed printing for mul in case of complex coef.
* Fixed the else condition for add, mul, div, sub in Rational.
* Added tests for complex division
Updated complex printing
In case the real or the imaginary part is zero, it is not printed
Added tests.
virtual bool is_minus_one() const { return ((this->real_ == -1) && (this->imaginary_ == 0)); }
//! \return `true` if both `real_` and `imaginary_` are non-zero
// This is needed while printing inside `mul`.
inline bool is_non_zero() const { return ((this->real_ != 0) && (this->imaginary_ != 0)); }

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 17, 2014

Author Contributor

This is perhaps not the best possible name. Any suggestions?

This comment has been minimized.

@certik

certik Jul 17, 2014

Contributor

Cannot you just do !is_zero()?

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 17, 2014

Author Contributor

I want both the real and imaginary part to be non-zero. It is a bit different than !is_zero(). Which is why the name is not apt.

This comment has been minimized.

@certik

certik Jul 17, 2014

Contributor

Ah, so the condition is that the number is off the x and y-axes (I mean real and imaginary axes). When would you use such a function?

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 17, 2014

Author Contributor

I'm using it to check during printing. So for instance if coeff_ is complex and both real and imaginary are non-zero, then I enclose the coeff_ inside brackets while printing, in other cases it is printed asis.

This comment has been minimized.

@certik

certik Jul 17, 2014

Contributor

I see. You can call it is_reim_non_zero, or is_off_reim_axes, or something like that. I don't like too many _, so maybe you can figure out a better name.

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 17, 2014

Author Contributor

Sure, I like is_reim_non_zero. I will put this for now, until I perhaps can think for a better name.

This comment has been minimized.

@certik

certik Jul 17, 2014

Contributor

Or perhaps is_reim_zero, similar to is_zero, and then just use the ! to obtain nonzero.

@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 18, 2014

@certik I've updated tests.
I also updated the basic printing for complex powers and complex multiplication.
Can you review this PR once?


//! \return `true` if `-1`
virtual bool is_minus_one() const { return ((this->real_ == -1) && (this->imaginary_ == 0)); }
//! \return `true` if both `real_` and `imaginary_` are zero

This comment has been minimized.

@certik

certik Jul 18, 2014

Contributor

Actually, the correct docstring is:

//! \return `true` if either `real_` or `imaginary_` is zero

// `I` is created only once in complex.cpp and reused
// everywhere (faster than creating it all the time):
extern RCP<const Number> I;

This comment has been minimized.

@certik

certik Jul 18, 2014

Contributor

Not in this PR, but I think we should create files constants.h/cpp and add all these global constants in there.

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 19, 2014

Author Contributor

Sure!

src/mul.cpp Outdated
imulnum(outArg(*coef), rcp_static_cast<const Number>(t));
} else if (rcp_static_cast<const Integer>(exp)->is_minus_one()) {
idivnum(outArg(*coef), rcp_static_cast<const Number>(t));
}

This comment has been minimized.

@certik

certik Jul 18, 2014

Contributor

So if exp is 2, and t is complex, then it will not be added into the dictionary at all, if I am reading this right, because the line insert(d, t, exp) below never gets executed.

So that's a bug.

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 19, 2014

Author Contributor

Right, I missed the else condition inside this.

src/mul.cpp Outdated
imulnum(outArg(*coef), rcp_static_cast<const Number>(t));
} else if (rcp_static_cast<const Integer>(it->second)->is_minus_one()) {
idivnum(outArg(*coef), rcp_static_cast<const Number>(t));
}

This comment has been minimized.

@certik

certik Jul 18, 2014

Contributor

If you are adding it to coef, don't you need to remove it from the dict? That's done using the line d.erase(it);.

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 19, 2014

Author Contributor

Correct, I'll fix this.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 18, 2014

We'll have to wrap this into Python and test carefully before merging. Even without running, I am suspecting there are at least 2 bugs, per my comments above.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 18, 2014

So even just using the current wrappers, you can see some stuff:

In [7]: sqrt(-1)**2
Out[7]: -1

In [8]: sqrt(-1)**3
Out[8]: --1^(1/2)

That should be fixed. I think it's just printing, though I think the sqrt(-1) should simplify to I.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 18, 2014

The printing needs to be fixed. In Pow, but also i1 should just print i.

In [1]: from csympy import *

In [2]: I
Out[2]: i1

In [3]: I**2
Out[3]: i1^2

In [4]: type(I)
Out[4]: csympy.lib.csympy_wrapper.Complex

In [5]: I+1
Out[5]: 1 + i1

In [6]: (I+1)**2
Out[6]: 1 + i1^2

In [7]: ((I+1)**2).expand()
Out[7]: 1 + i1^2
@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 18, 2014

Also, as I suspected, expand() is incorrect, i.e. in sympy:

In [6]: ((x+I)**2).expand()
Out[6]: x**2 + 2*I*x - 1

and csympy:

In [14]: ((x+I)**2).expand()
Out[14]: 1 + x^2 + i2x

The 1 should be -1.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 18, 2014

I think that if Complex is pure imaginary, then we should implement the integer powers of it, e.g. i^2, i^4, (5i)^10.... That should be simple.

Add Python wrappers for complex numbers
Test SymPy <-> CSymPy conversion as well.
@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 18, 2014

@sushant-hiray, I wrote the Python wrappers for you: sushant-hiray#1

So just play with this in Python, and make sure that everything looks good and correct.

src/pow.cpp Outdated
res = mulnum(I, minus_one);
imulnum(outArg(overall_coeff), mulnum(im->pow(*pow_new), res));
} else {
Mul::dict_add_term_new(outArg(overall_coeff), d, exp, i2->second);

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 22, 2014

Author Contributor

The main issue with expand not working perfectly is in this case.

So primarily the issue is: things like 2+3*I as treated as a single complex number, one way to go about fixing this is to break 2+3*I into add(i2, mul(i3, I)) and then expand this.

This comment has been minimized.

@certik

certik Jul 22, 2014

Contributor

I think that add(i2, mul(i3, I)) automatically recombines to a single complex number, as it should.

I would leave the current behavior, more on this below.

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 22, 2014

Author Contributor

Aah, my bad it will recombine again.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 22, 2014

Ooops, I broke master. Fixing it now.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 22, 2014

Master is fixed by #262. Tests should now work in this PR. Sorry about that.

@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 22, 2014

No issues about that 😄

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 22, 2014

As to the example:

((1 + I)*x).expand()

I think that's fine. This can't be expanded, as 1+I is now a single complex number. This, btw, is consistent with Mathematica.
As to

In [15]: ((x + 2 + 3*I)**2).expand()
Out[15]: (2 + 3*I)^2 + 4 + 6*I*x + x^2

Here again, it's quite consistent. The issue is if and when to do the (2 + 3*I)^2 expansion. Mathematica does it automatically. I.e. if you have a complex number raised to an integer (positive or negative), it simplifies this automatically, because the answer of (x+I*y)^n is always x+I*y, so the result is clearly simpler (of course, only if x+I*y is of type Complex).

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 22, 2014

To put it another way, the operation Complex^n is just another Complex, so it is just a matter of quickly finding it. As you said before, it's the equation (9) at http://mathworld.wolfram.com/ComplexNumber.html. We have the function that computes the binomial coefficients fast. And then you just put the even ones to the real part, and odd ones to the imaginary part, so I think this is a pretty fast operation. I tried it in Mathematica, and even things like (3+7*I)^50000 are immediate.

And I think this greatly simplifies a lot of things. This would mean that also things like 1/(2+3*I) will get simplified automatically. I think that's fine.

What do you think?

@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 22, 2014

Yes this makes sense.
Given the fact that we have a function computing binomial coefficients, we should perhaps simplify it.
This would also mean less special cases in the code as you pointed out in the case of 1/(2+3*I)

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 22, 2014

Here is the function that you need: https://github.com/sympy/csympy/blob/master/src/pow.cpp#L222 and here is the loop https://github.com/sympy/csympy/blob/master/src/pow.cpp#L287 that uses it. In our case here, the loop will be very simple --- you just power the two rationals x and y to the power specified by the r (result of the multinomial_coefficients_mpz), add the (-1) factor, and put even terms in the new real part and odd terms into the imaginary part. I think that's it.

@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 22, 2014

Thanks, I will look into it tomorrow.
Also builds fail in test_arit.py:

    e = Integer(2)*x
    assert str(e) == "2x"
    e = 2*x
    assert str(e) == "2x"

After the latest commit, these are printed as 2*x.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 22, 2014

Yes, the print tests need to be updated.

Thanks @sushant-hiray, the integer power is the last thing to fix. Given the complexity of this PR already, do you have any objections in merging this first, and implement the integer power in a new PR?

To merge this PR, all tests need to pass.

src/mul.cpp Outdated
@@ -102,12 +104,20 @@ std::string Mul::__str__() const
std::ostringstream o;
if (eq(coef_, minus_one))
o << "-";
else if (is_a<Rational>(*coef_) &&

This comment has been minimized.

@certik

certik Jul 22, 2014

Contributor

You should always put braces around an if statement if you use any else statements, because it is more readable and leads to less errors (there is no ambiguity if there are braces, but without them there might be).

src/mul.cpp Outdated
else if (is_a<Rational>(*coef_) &&
!(rcp_static_cast<const Rational>(coef_)->is_int()))
o << "(" << *coef_ <<")";
else if (is_a<Complex>(*coef_)) {

This comment has been minimized.

@certik

certik Jul 22, 2014

Contributor

braces around else

@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 22, 2014

Sure, I'll fix those print tests and then we can merge the PR once everything passes.

src/mul.cpp Outdated
if (!(rcp_static_cast<const Complex>(coef_)->is_reim_zero()))
o << "(" << *coef_ <<")";
else
o << *coef_;

This comment has been minimized.

@certik

certik Jul 22, 2014

Contributor

put braces here too

src/pow.cpp Outdated
else if (is_a<Rational>(*b) &&
(rcp_static_cast<const Rational>(b)->i.get_num() == 1) &&
(rcp_static_cast<const Rational>(b)->i.get_den() == 2))
return I;

This comment has been minimized.

@certik

certik Jul 22, 2014

Contributor

put braces around these if statements.

src/pow.cpp Outdated
else if (eq(res, integer(2)))
res = minus_one;
else
res = mulnum(I, minus_one);

This comment has been minimized.

@certik

certik Jul 22, 2014

Contributor

put braces around these if statements.

src/pow.cpp Outdated
imulnum(outArg(overall_coeff), i2->second);
else if (exp->is_minus_one())
idivnum(outArg(overall_coeff), i2->second);
else if (rcp_static_cast<const Complex>(i2->second)->is_re_zero() &&

This comment has been minimized.

@certik

certik Jul 22, 2014

Contributor

braces around if statements

src/pow.cpp Outdated
res = I;
else if (eq(res, integer(2)))
res = minus_one;
else

This comment has been minimized.

@certik

certik Jul 22, 2014

Contributor

braces around if statements

Fixed failing tests. PR cleanup
Fixed braces issues around if-else blocks
@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 23, 2014

@certik I've cleaned the PR according to your comments and fixed the tests.
Travis tests now pass here

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 23, 2014

The PR is +1, but we need to run benchmarks before merging. I'll try to get to it soon.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 23, 2014

So this is master (8ffd374):

ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand1
Expanding: (y + x + z + w)^60
104ms
number of terms: 39711
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand1
Expanding: (y + x + z + w)^60
104ms
number of terms: 39711
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand1
Expanding: (y + x + z + w)^60
103ms
number of terms: 39711
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand1
Expanding: (y + x + z + w)^60
107ms
number of terms: 39711
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand2
Expanding: ((y + x + z + w)^15 + w)*(y + x + z + w)^15
1181ms
number of terms: 6272
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand2
Expanding: ((y + x + z + w)^15 + w)*(y + x + z + w)^15
1170ms
number of terms: 6272
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand2
Expanding: ((y + x + z + w)^15 + w)*(y + x + z + w)^15
1168ms
number of terms: 6272
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand2
Expanding: ((y + x + z + w)^15 + w)*(y + x + z + w)^15
1196ms
number of terms: 6272
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand3
Expanding: (z^x + x^y + y^x)^100
31ms
number of terms: 5151
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand3
Expanding: (z^x + x^y + y^x)^100
33ms
number of terms: 5151
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand3
Expanding: (z^x + x^y + y^x)^100
32ms
number of terms: 5151
ondrej@kittiwake:~/repos/csympy/benchmarks(master)$ ./expand3
Expanding: (z^x + x^y + y^x)^100
31ms
number of terms: 5151

And this is this PR (1e4cfbf merged with master 8ffd374):

ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand1
Expanding: (y + x + z + w)^60
108ms
number of terms: 39711
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand1
Expanding: (y + x + z + w)^60
109ms
number of terms: 39711
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand1
Expanding: (y + x + z + w)^60
111ms
number of terms: 39711
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand1
Expanding: (y + x + z + w)^60
112ms
number of terms: 39711
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand2
Expanding: ((y + x + z + w)^15 + w)*(y + x + z + w)^15
1204ms
number of terms: 6272
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand2
Expanding: ((y + x + z + w)^15 + w)*(y + x + z + w)^15
1214ms
number of terms: 6272
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand2
Expanding: ((y + x + z + w)^15 + w)*(y + x + z + w)^15
1170ms
number of terms: 6272
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand2
Expanding: ((y + x + z + w)^15 + w)*(y + x + z + w)^15
1186ms
number of terms: 6272
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand3
Expanding: (z^x + x^y + y^x)^100
31ms
number of terms: 5151
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand3
Expanding: (z^x + x^y + y^x)^100
30ms
number of terms: 5151
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand3
Expanding: (z^x + x^y + y^x)^100
33ms
number of terms: 5151
ondrej@kittiwake:~/repos/csympy/benchmarks(bench)$ ./expand3
Expanding: (z^x + x^y + y^x)^100
33ms
number of terms: 5151

So I think that there is no significant slowdown.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 23, 2014

I ran lots of expand2 benchmarks and the best value for master that I could get is 1146ms. The best value for this PR is 1165ms. So it's a bit slower.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 23, 2014

I think the slow down must be caused by using more of the is_a functions. Those should be sped up along the lines of #127.

I think this PR is fine to merge. Thanks!

certik added a commit that referenced this pull request Jul 23, 2014

Merge pull request #248 from sushant-hiray/complex-fix
Further fixes to integrate Complex module with CSymPy

@certik certik merged commit 31bf2a6 into symengine:master Jul 23, 2014

2 checks passed

continuous-integration/travis-ci The Travis CI build passed
Details
default All builds succeeded on Shippable
Details
@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 23, 2014

@sushant-hiray I created a new issue #263 for the integer powers.

@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 24, 2014

@certik Thanks for merging this.

@sushant-hiray sushant-hiray deleted the sushant-hiray:complex-fix branch Jul 24, 2014

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment