@Markiyan Hirnyk The statement is true for Maple 14, 15, and 16. If you are basing your counterclaim on results from Maple 13 or an even older major release then I would say that the moral obligation is on you to state that up front.

(For example, Maple 13 is three major releases back, and it would be bizarre to think that the unqualified statement should be that relating to some much older release as opposed to all the later releases.)

We should not forget that the rest of the claim is also true: using the solve command to get float approximations of all complex roots of a higher degree univariate polynomial (which does not factor explicitly) is much slower than using fsolve, in general. It is often worse by a multiple of n the degree of the polynomial, as solve simply calls fsolve n times to do the same work. This is demonstrated in the link given. This part is also true in Maple 13, 12, etc.

This accuracy comparison below is done in Maple 14.01 (Windows, 64bit), but the results are similar in 15.01 and 16.00. For these particular computations done under Digits=12 the results from fsolve are shown to have a greater number of accurate decimal digits, and produce smaller absolute forward error with respect to resubstitution in the original polynomial.

restart:
Digits:=12:

fsol:=[fsolve(z^5-5*z^4-5*z-5, z, complex)]:

forget(evalf):
ssol:=[solve(z^5-5*z^4-5*z-5.0, z)]:

forget(evalf):
Digits:=100:
highsol:=[fsolve(z^5-5*z^4-5*z-5, z, complex)]:

max(map(abs,sort(highsol)-sort(ssol))):
evalf[4](%);
                                         -12
                                 2.729 10   

max(map(abs,sort(highsol)-sort(fsol))):
evalf[4](%);

                                         -16
                                 8.203 10   

[seq(eval(z^5-5*z^4-5*z-5, z=x), x in fsol)]:
evalf[4](sort(%,(a,b)->abs(a)>abs(b)));

 [        -13           -15           -15             -15           -15    
 [5.476 10   , -2.189 10    - 1.555 10    I, -2.189 10    + 1.555 10    I, 

           -16           -16            -16           -16  ]
   8.853 10    + 2.778 10    I, 8.853 10    - 2.778 10    I]

[seq(eval(z^5-5*z^4-5*z-5, z=x), x in ssol)]:
evalf[4](sort(%,(a,b)->abs(a)>abs(b)));

 [         -10           -11           -11             -11           -11    
 [-7.234 10   , -4.471 10    + 5.918 10    I, -4.471 10    - 5.918 10    I, 

            -12           -12             -12           -12  ]
   -1.190 10    - 5.466 10    I, -1.190 10    + 5.466 10    I]

If you are going to make multiple posts like this (which you have just done) then please could you either include in each (both title and body) some keyword that informs readers as to the topics covered?

Otherwise, please couldn't you just roll them together in a single post?

What you've done so far is just submit a bunch of separate posts linking to other sites, plus zip files, each with no indication of the subject of content. That's a lot like spam.

acer

I didn't mention `numoccur` because it brings the extra need to match the nature of the zeros. If the Matrix has exact 0 then numoccur needs to be passed exact 0 as the second argument. And similarly for floating point 0.0.

So using `numoccur` can bring an extra burden of knowing what type of zeroes are to be tested, and what kind of data is in the Matrix. Of course, if you know the type of data for sure, or if you wish to distinguish between 0 and 0.0, then it's not an issue.

> M:=Matrix([[0.0, 0.0],[3.0,0.5]]);
                                    [0.     0. ]
                               M := [          ]
                                    [3.0    0.5]

> numboccur(M, 0);                  
                                       0

> numboccur(M, 0.0);                
                                       2

> 4-add(`if`(x<>0,1,0),x=M);        
                                       2

> 4-rtable_num_elems(M,'NonZero');  
                                       2

> M:=Matrix([[0, 0],[3,1/2]]);      
                                     [0     0 ]
                                M := [        ]
                                     [3    1/2]

> numboccur(M, 0);                
                                       2

> numboccur(M, 0.0);              
                                       0

> 4-add(`if`(x<>0,1,0),x=M);      
                                       2

> 4-rtable_num_elems(M,'NonZero');
                                       2

And then the mixed case gets more awkward still.

> M:=Matrix([[0, 0.0],[3,5]]);
                                     [0    0.]
                                M := [       ]
                                     [3    5 ]

> numboccur(M, 0);            
                                       1

> numboccur(M, 0.0);
                                       1

> 4-rtable_num_elems(M,'NonZero'); 
                                       2

> 4-add(`if`(x<>0,1,0),x=M);       
                                       2

acer

I didn't mention `numoccur` because it brings the extra need to match the nature of the zeros. If the Matrix has exact 0 then numoccur needs to be passed exact 0 as the second argument. And similarly for floating point 0.0.

So using `numoccur` can bring an extra burden of knowing what type of zeroes are to be tested, and what kind of data is in the Matrix. Of course, if you know the type of data for sure, or if you wish to distinguish between 0 and 0.0, then it's not an issue.

> M:=Matrix([[0.0, 0.0],[3.0,0.5]]);
                                    [0.     0. ]
                               M := [          ]
                                    [3.0    0.5]

> numboccur(M, 0);                  
                                       0

> numboccur(M, 0.0);                
                                       2

> 4-add(`if`(x<>0,1,0),x=M);        
                                       2

> 4-rtable_num_elems(M,'NonZero');  
                                       2

> M:=Matrix([[0, 0],[3,1/2]]);      
                                     [0     0 ]
                                M := [        ]
                                     [3    1/2]

> numboccur(M, 0);                
                                       2

> numboccur(M, 0.0);              
                                       0

> 4-add(`if`(x<>0,1,0),x=M);      
                                       2

> 4-rtable_num_elems(M,'NonZero');
                                       2

And then the mixed case gets more awkward still.

> M:=Matrix([[0, 0.0],[3,5]]);
                                     [0    0.]
                                M := [       ]
                                     [3    5 ]

> numboccur(M, 0);            
                                       1

> numboccur(M, 0.0);
                                       1

> 4-rtable_num_elems(M,'NonZero'); 
                                       2

> 4-add(`if`(x<>0,1,0),x=M);       
                                       2

acer

Another ungraceful way to do this particular example posed in the Question, at present, is to substitute a single name for Pi*a, then simplify, and then backsubstitute,

restart:
expr := -2*Pi*sin(Pi*a)/(-1+cos(2*Pi*a)):

subs(g=a*Pi,simplify(algsubs(a*Pi=g,expr)));

                              Pi    
                           ---------
                           sin(Pi a)

In principle this workaround should not be necessary as the Pythagorean identity, upon which this particular simplification can hinge, should hold regardless of such a scaling.

acer

Another ungraceful way to do this particular example posed in the Question, at present, is to substitute a single name for Pi*a, then simplify, and then backsubstitute,

restart:
expr := -2*Pi*sin(Pi*a)/(-1+cos(2*Pi*a)):

subs(g=a*Pi,simplify(algsubs(a*Pi=g,expr)));

                              Pi    
                           ---------
                           sin(Pi a)

In principle this workaround should not be necessary as the Pythagorean identity, upon which this particular simplification can hinge, should hold regardless of such a scaling.

acer

This Question leads to the following case which gets overlooked by the simplify command.

restart:

simplify( -1+cos(a)^2 ); # ok

                                   2
                            -sin(a) 

simplify( -1+cos(a*b)^2 ); # ?? missed ??

                                      2
                         -1 + cos(a b) 

simplify( -1+sin(a*b)^2 ); # ok

                                    2
                           -cos(a b) 

If that middle (problem) case above were fixed, so that the Pythagoras identity were recognized and used, then I could envision the submitter's example succeeding more directly. With that bug fixed then just `simplify, or `simplify` after `expand`, might get the desired result.

I will submit a bug report on that middle case above.

It worked fine for me, on 64bit Maple 15.00 or 15.01 each running on 64bit Linux ubuntu 10.04.

Please, let me ask one question though: are you using 2D Math input? Is it possible that you have an extra space between the int and the (...) bracketed piece? I'd just like to eliminate that possibility first, that it might be a case of inadvertant implicit multiplication.

acer

I don't see how this particular type-check that you've assembled in the Criterium [sic] procedure is useful for the described task. How does it distinguish between floats which are purely real and those which are not?

acer

I don't see how this particular type-check that you've assembled in the Criterium [sic] procedure is useful for the described task. How does it distinguish between floats which are purely real and those which are not?

acer

@Markiyan Hirnyk I have changed the Alias(MM,[18]), which did not appear in my version of pvrbik's procedure but was an additional aside, to Alias(MM,[op([1,1],MM)]) in my comment above so that it is more general. One could pick off the first dimension of an mx1 Matrix in several other ways, as well.

I will mention that you too hard-coded the size 18 into your own version above, as opposed to picking off the dimension programmatically.

Was it the combined local declaration and assignment that concerned you? I haven't checked the documentation to see whether it is properly mentioned. Since pvrbik used it before me, in his version, then I tried it out. Recall that it is his code, not mine. I reckoned that if I could figure it out, based on pvrbik's code and accept it, then other people might too. The echoed printing of a procedure below gives some indication as to how it is parsed,

f:=proc(x::list)
   local a,b,
   c:=x[1],
   d:=[c];
return d[1];
end proc;

                             f := proc(x::list)
                             local a, b, c, d;
                               c := x[1];
                               d := [c];
                               return d[1];
                             end proc;

@Markiyan Hirnyk I have changed the Alias(MM,[18]), which did not appear in my version of pvrbik's procedure but was an additional aside, to Alias(MM,[op([1,1],MM)]) in my comment above so that it is more general. One could pick off the first dimension of an mx1 Matrix in several other ways, as well.

I will mention that you too hard-coded the size 18 into your own version above, as opposed to picking off the dimension programmatically.

Was it the combined local declaration and assignment that concerned you? I haven't checked the documentation to see whether it is properly mentioned. Since pvrbik used it before me, in his version, then I tried it out. Recall that it is his code, not mine. I reckoned that if I could figure it out, based on pvrbik's code and accept it, then other people might too. The echoed printing of a procedure below gives some indication as to how it is parsed,

f:=proc(x::list)
   local a,b,
   c:=x[1],
   d:=[c];
return d[1];
end proc;

                             f := proc(x::list)
                             local a, b, c, d;
                               c := x[1];
                               d := [c];
                               return d[1];
                             end proc;

@Markiyan Hirnyk Just run the code in my previous comment in this answer thread; it runs in Maple 13 if properly pasted in a 1D Maple notation code. (It will not work, if pasted as 2D Math input.)

I don't understand what you are saying about parameters; it reads as nonsense to me. You may be having an off week, but please don't take it out on me.

@Markiyan Hirnyk Just run the code in my previous comment in this answer thread; it runs in Maple 13 if properly pasted in a 1D Maple notation code. (It will not work, if pasted as 2D Math input.)

I don't understand what you are saying about parameters; it reads as nonsense to me. You may be having an off week, but please don't take it out on me.

@Markiyan Hirnyk It's extremely easy to fix it for Maple 13. Just swictch around the R and M declaration/assignments, and have the R be done in terms of M rather than MM.

After such an edit, pvbik's code could become,

M := Matrix(18, 1, [2, 3, 3, 5, 7, 8, 12, 5, 9, -3, 4.1, 7, 7, 7, -3, 9, 3, 8]):

MatLabSort := proc( MM :: Matrix ) :: list( Matrix, list );
local  i, Pos,
    SortedMM := table(),
    M := convert( MM, list ),
    R := convert(Statistics[Rank](M),list);
    
    for i from 1 to nops(R) do
        SortedMM[ R[i] ] := M[i];
        Pos[ R[i] ] := i;
    end do;
    
    return [ convert(convert(SortedMM,list),Vector), convert(Pos,list) ];
    
end proc;

MatLabSort( M );

Or one could use Statistics[Rank](ArrayTools:-Alias(MM,[op([1,1],MM)])), to avoid that inefficiency of the list conversion.