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

Expand integral powers of complex nos #264

Merged
merged 4 commits into from Jul 26, 2014

Conversation

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

sushant-hiray commented Jul 25, 2014

The pow_number needs further checks for overflow.

Updated pow_expand to expand integral powers of complex nos
* Updated print for `add` in case of Complex numbers
* Added `pow_number` function
* Added a test for checking running time
@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 25, 2014

@certik Currently pow_expand checks if the base is of form Add and exp is of form Integer and then only performs exp. If I update this to also include the case when base is Complex, we won't be able to perform dynamic_cast into Add.
The dynamic cast was used to find the size of the dict while printing in the tests.

Should I update this then?
In the tests instead of doing a dynamic cast, we can do a static cast after checking that the final result is an Add.

@sushant-hiray sushant-hiray changed the title Updated pow_expand to expand integral powers of complex nos Expand integral powers of complex nos Jul 25, 2014

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 25, 2014

Pow_number looks good. This needs to be called in the "pow" function. The
pow_expand needs to be updated too, but I don't like that we are currently
repeating code there. We should implement the complex simplification only
once and then just call it. Where is the dynamic cast used?

Sent from my mobile phone.
On Jul 25, 2014 4:47 AM, "Sushant Hiray" notifications@github.com wrote:

@certik https://github.com/certik Currently pow_expand checks if the
base is of form Add and exp is of form Integer and then only performs
exp. If I update this to also include the case when base is Complex, we
won't be able to perform dynamic_cast into Add.
The dynamic cast was used to find the size of the dict while printing in
the tests.

Should I update this then?
In the tests instead of doing a dynamic cast, we can do a static cast
after checking that the final result is an Add.


Reply to this email directly or view it on GitHub
#264 (comment).

@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 25, 2014

I had initially called it in pow function, but then things like (2 + 3*I)^k where k is an integer are expanded by default.
This will perhaps deviate from what SymPy does by default, as SymPy doesn't expand unless called with .expand().
I guess I updated pow_expand. What else do you wish to add in this function?

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 25, 2014

@sushant-hiray we discussed it here:
#248 (comment)
#248 (comment)
where I wrote

And I think this /=expanding integer powers automatically/ 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.

to which you wrote:

Yes this makes sense.

However, if you have some doubts now, we can discuss it further. I am convinced we should expand integer powers (both positive and negative) automatically, because that's what Python does, as well as Mathematica does. SymPy doesn't do it, that maybe stems from the fact that I in SymPy is handled like a symbol with special properties. In CSymPy, we use a better design (in my opinion) that we have a complex number (that have some special properties to multiply), but symbols all act the same way and numbers also all act in the same way.

In our design, the rule of thumb is that if you have any arithmetic operation (+, -, *, /, ^) on Number (be it integer, rational or complex) and the result of the operation is a Number, then you do it automatically (and it will also be fast, because it turns out all these operations are very fast). Sometimes, the result is not a Number, for example 2^(1/2). Then of course we leave it be.

p = mulnum(p, p);
}
return r;
}

This comment has been minimized.

@certik

certik Jul 25, 2014

Contributor

This looks good. This assumes that n>0. Perhaps you can change it to unsigned long, then it's clear. I think it works also for n=0, i.e. it returns 1, so I think that's fine. Then we need a similar function, that does it for negative n --- essentially by calling this pow_number and then dividing by the result (or perhaps first dividing then powering, which might be faster, as you are dividing with a smaller number).

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 25, 2014

Author Contributor

Dividing and powering is a better approach. I'll add the check for the positive and negative power

This comment has been minimized.

@certik

certik Jul 26, 2014

Contributor

Great. I still think the argument n should be unsigned. But otherwise I think we can merge it, this can be fixed in a subsequent PR.

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 26, 2014

Author Contributor

Aah, right I will fix that and then we can merge this.

I wasn't sure if it is safe to convert from signed to unsigned int.
But it is safe so no issues as such.
http://stackoverflow.com/questions/50605/signed-to-unsigned-conversion-in-c-is-it-always-safe

Sushant Hiray,
Senior Undergrad CSE,
IIT Bombay

On Sat, Jul 26, 2014 at 6:41 PM, Ondřej Čertík notifications@github.com
wrote:

In src/pow.cpp:

@@ -322,6 +322,21 @@ void multinomial_coefficients_mpz(
int m, int n, map_vec_mpz &r)
}
}

+RCP pow_number(const RCP &x, long n)
+{

  • RCP r, p;
  • long mask = 1;
  • r = one;
  • p = x;
  • while (mask > 0 && n >= mask) {
  •    if (n & mask)
    
  •        r = mulnum(r, p);
    
  •    mask = mask << 1;
    
  •    p = mulnum(p, p);
    
  • }
  • return r;
    +}

Great. I still think the argument n should be unsigned. But otherwise I
think we can merge it, this can be fixed in a subsequent PR.


Reply to this email directly or view it on GitHub
https://github.com/sympy/csympy/pull/264/files#r15433775.

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 26, 2014

Author Contributor

Oh well you merged this!
We can fix this later on then.

src/pow.cpp Outdated
} else {
res = mulnum(I, minus_one);
}
imulnum(outArg(overall_coeff), mulnum(im->pow(*pow_new), res));

This comment has been minimized.

@certik

certik Jul 25, 2014

Contributor

I don't like that we are repeating this logic here, because we already have it in pow. So we should refactor this into some function and call it from both places.

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 25, 2014

Author Contributor

Yes, I will do this in the upcoming commits.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 25, 2014

If you disagree, let's discuss it. But if you agree, then I would first implement this in pow and test it. Only then we should figure out how to make sure that expand() does the right thing.

sushant-hiray added some commits Jul 25, 2014

Removed redundant code blocks and replaced with a function
Added a new function, pow_complex which contains the code
corresponding to complex^integer
@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 25, 2014

I've handled the case for negative powers.
Also I added a new function pow_complex and replaced the redundant code-blocks for complex^integer by making a function call.

src/pow.cpp Outdated

void pow_complex(const Ptr<RCP<const Number>> &self,
const RCP<const Complex> &base_,
const RCP<const Integer> &exp_)

This comment has been minimized.

@certik

certik Jul 25, 2014

Contributor

The base_ is passed into functions that expect RCP, but exp_ isn't, so you should just declare it as
const Integer &exp_. Essentially, the rule of thumb is: don't use RCP, unless you have to (as in general, things are faster if you don't use RCP).

This comment has been minimized.

@sushant-hiray

sushant-hiray Jul 26, 2014

Author Contributor

Oh sure, I'll keep note of that next time onwards.

@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 25, 2014

Also, I tried this benchmark:

In [2]: %time e = (Integer(2)/3 + I/5)**1000000
CPU times: user 12 s, sys: 28 ms, total: 12 s
Wall time: 12 s

and compared to Sage:

sage: %time e = (2/3+I/5)^1000000
CPU times: user 1.59 s, sys: 19.8 ms, total: 1.61 s
Wall time: 1.61 s

we are a lot slower. Not sure what the problem could be, but after we merge this PR, we should have a look.

Updated function definition for `pow_complex`
Now `pow_complex` takes `exp_` as `const Integer` as opposed to
`RCP<const Integer>` in the former version
@sushant-hiray

This comment has been minimized.

Copy link
Contributor Author

sushant-hiray commented Jul 26, 2014

The difference in computing time for rational and integral (real and imaginary) parts is quite large.

Here are some benchmarks:

In [2]: %time e = (Integer(2)/3 + I/5)**1000000
CPU times: user 1min 4s, sys: 280 ms, total: 1min 4s
Wall time: 1min 4s

In [3]: %time e = (3 + I)**1000000
CPU times: user 276 ms, sys: 0 ns, total: 276 ms
Wall time: 278 ms

In [4]: %time e = (3 + 2*I)**1000000
CPU times: user 296 ms, sys: 0 ns, total: 296 ms
Wall time: 295 ms
@certik

This comment has been minimized.

Copy link
Contributor

certik commented Jul 26, 2014

Yes, the rational version is always slower. I think it looks good for now, so I am merging it. Speed can be improved later.

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

Merge pull request #264 from sushant-hiray/complex-pow
Expand integral powers of complex nos

@certik certik merged commit b759c50 into symengine:master Jul 26, 2014

1 of 2 checks passed

continuous-integration/travis-ci The Travis CI build is in progress
Details
default All builds succeeded on Shippable
Details

@sushant-hiray sushant-hiray deleted the sushant-hiray:complex-pow branch Jul 26, 2014

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