isequal is relatively slow, about 2X slower than '=='. For a 100-element vector the absolute difference is still small, about 400 microseconds, but I think it is worth testing using isempty before using all. I made that change and another small tweak or two and checked in the changeset here (https://hg.savannah.gnu.org/hgweb/octave/rev/e02f500a0980).
Marking as fixed and closing report.
|
@Markus: thanks.
A patch is included.
(file #45248)
|
Would this do what you need?
|
If you take
you get true, which is not correct. How do I check whether sort(x(imag(x)>0)) is the same vector of sort(conj(x(imag(x)<0))))?
|
I will do asap.
|
Okay, then I think this is fine to work up in to a true changeset and commit it.
|
As far as I know, eig (that is, the undelying LAPACK routines) gives exact complex conjugates for real matrices. This is confirmed by experiments like
Therefore, I would not introduce the tolerance.
|
That seems to work. Would the call to eig ever return values which are not exactly complex conjugates such that we need to use a tolerance on the equality test? In other words,
|
What about the following?
|
Yes, but is it on the right track? Using your example from comment #2 I find that the largest imaginary value is 10 eps which is quite small. At least there are no horrendous errors like 256i as in the first example.
It seems that sort followed by a test to replace small values with 0 would work. The algorithm in poly.m is
which has a multiplication and a subtraction per round of the algorithm. To me that says worst case of 2 eps per cycle or 2 x n x eps for the algorithm. It's not perfect, but for eps one could use the largest value in v.
Of course, this is starting to look like a fair amount of work so if detecting the conjugate pairs is easier then we could do that as you suggested in comment #0.
|
@Rik: your idea does not work for the following example
i.e., elements of the vector with the same absolute value.
|
The algorithm in poly.m uses a for loop and repeated multiplications. It could be something as simple as, with the inversion of matrices, using balance() to put the largest factors first, etc.
In this example, if I use sort on the input x then the conjugates are paired and the result of poly is entirely real without spurious imaginary terms.
See the session below
Should we just sort the input vector 'v' after the input validation in poly.m?
|
If you try the following
you see that N has spurious imaginary parts (quite large, in absolute value). This should not happen, since the roots x are complex conjugate pairs. We should first detect this fact (I thought to imag (sum (sort (x))) == 0, it is probably possible to do without computations, but only with sorting and comparison). Then, either we use the formula
+verbatim
(y-x(1))*(y-conj(x(1)) = y^2-2*real(x(1))*y+abs(x(1))^2
-verbatim-
for each complex conjugate pair or we simply remove the spurious imaginary parts after the computation of the coefficients.
|