Use the divide command to explicitly check whether the polynomial is divisible by (Px-Py) (over the rational numbers). If you are looking for factorizations over more complicated fields than the rationals, that is also possible, but you need to specify the field extension(s) (if they're not implied by the coefficients). Let me know.

Update: Your polynomial contains terms that have neither a p1 nor a p2. So p1-p2 couldn't possibly be a factor.

You can use a userinfo statement or a print statement to print any infomation about an ongoing computation. That includes plots. If you put a plot inside a userinfo, it needs to be inside print also, like

userinfo(1, mytest, print(plot(...)));

To "activate" the above statement, before you run the procedure do

infolevel[mytest]:= 1:

It is supported: I've uploaded animated GIF files to MaplePrimes many times. I do not use the image->include tool. I always upload files using the green uparrow (Upload File), which is the last thing on the second row of the toolbar. The animation will automatically "play" in the editor and in the post.

I can't understand what you want, and you can't understand what I've given you. Please get someone to help you with your English before posting.

You wrote:

The Maxwell-Boltzmann distribution sampling times are 2002. But you just list five sequences. I don't know why.

Five was just a convenient number that fit on one line on the screen. The way that you are using the word "sampling" is not standard mathematical English. The entire population of possible distributions of 9 quanta among 6 particles has size 2002. But that is neither a sample nor a time.

I am very confuced about the sampling. How can you set sampling times are 2002?

I have no idea why you would want a sample of size 2002. The idea seems crazy to me. But you can change the five to any number that you want. That should be obvious!

'R()' $ 2002;

Now, if you want the whole population, that's another matter. Now (in the worksheet posted below) I have included an option to return the whole population: ouput= population.

You said the computation is exact---not based on random sampling. But latter you said that does either the exact computation or provides a random sampler depending on an input option.

If you say output= exact then the computation is exact! If you say output= sampler then the output is a random sample! Why can't you understand that? Please get someone to help you read the code. The computation where I compared the numbers with those on the website were the exact computations.

What's the difference from random sampling and random sampler?

I am using "sampler" to mean a procedure that generates a random sample. When you call my procedure with the option output= sampler then the object returned is itself a procedure---a procedure which draws one item from the population each time it is called.

Download QuantumDistribution.mw

@H-R It is usually easier to specify a unique characterization of the component that you want to extract than it is to fix its operand number. For your situation, could I characterize the part that you want to extract as "the factor in the sum that has a sinh in it"? Would that be general enough to cover every situation where you would use this and yet precise enough to uniquely specify the part that you want?

For the case at hand, the following extracts the subexpression that you want:

select(hastype, op(1, indets(BC4, specfunc(anything, sum))[]), sinh(algebraic));

This assumes that BC4 has only one sum.

m:= 1:  N:= 2:
Z:= Matrix(m+1, m+1):
Id:= Matrix([1$m+1], scan= band[0]):
A:= Matrix([[Z$N], [Id$N-1]], scan= band[0,1]);

See ?scan (unfortunately there are no examples on this help page).

Here's another solution completely in polar coordinates. Unlike Kitonum's solution (which I don't mean to disparage), I derive the polar form of the necessary formulae for slope and tangent line from first principles.

restart:
R:= theta-> 3 - 3*cos(theta):
theta0:= 3*Pi/4:
polar:= (r,theta)-> {x = r(theta)*cos(theta), y(x) = r(theta)*sin(theta)}:
dy_dx:= PDEtools:-dchange(polar(r,theta), diff(y(x),x)):
line:= PDEtools:-dchange(polar(r,theta), y(x) - eval(y(x), polar(R,theta0)) =
     m*(x - eval(x, polar(R,theta0)))):
line_polar:= solve(line, r(theta)):
m:= eval(dy_dx, [r= R, theta= theta0]):
P1:= plots:-polarplot(R(theta), theta= 0..2*Pi):
P2:= plots:-polarplot(line_polar, theta= Pi/4..Pi):
plots:-display([P1,P2]);

Here's a procedure for it:

Pbar:= proc(P::Matrix(square), k::posint)
local
     m:= op([1,1], P),
     km:= k*m,
     M:= Matrix(km, km, (i,j)-> `if`(i < j, 1/m, 0)),
     n
;
     for n to k do M[(n-1)*m+1..n*m, (n-1)*m+1..m*n]:= P end do;
     M/k
end proc:

Example of use:

Pbar(< 1,2 ; 3,4 >, 3);

One way is by declaring the parameter as type evaln:

mysum:= (data::evaln(list))-> printf("The sum of %a is %3.3f.", data, `+`(eval(data)[])):

listA:= [1,2,3]:  listB:= [4,5,6]:
mysum(listA);
The sum of listA is 6.000.

Note that the parameter itself refers to the assigned name and that the contents of the name are accessed by using eval.

You are missing a left square bracket before b1, and Gamma should be GAMMA.

If that doesn't solve your problem, then I need to know what version of Maple you are using.

When you receive a "lost connection to kernel" error, you need to

1. save your worksheet
2. close your worksheet
3. open it again

Okay, here is one procedure that does either the exact computation or provides a random sampler depending on an input option.

THe result for Bose-Einstein energy level 4 given on that website must be wrong. If you multiply the energy levels times their average values and add them up, you must get 9, right? If you do that for the numbers in the Bose-Einstein column on that website, you will see that you don't get 9. My results come from exact computation on the entire population of possible states.

Download QuantumDistribution.mw

The symbol has no mathematical meaning in Maple, but for typing in mathematical text mode, the ordinary vertical bar works fine |. On my keyboard, it's the shifted version of the backslash key.

P1:= PlanePlot(-3*x+ 2*y + 3*z = 0, [x,y,z], normaloptions= [shape= harpoon], showbasis):
P2:= PlanePlot(
     2*y+3*z = 0, [x,y,z], normaloptions= [shape= harpoon], showbasis,
     planeoptions= [color= red]
):
plots:-display([P1,P2], caption= "");

Here are two ways to accomplish the task with rtable_scanblock. This first one is in the spirit of the examples given at ?rtable_scanblock :

ind:= rtable_scanblock(A, [], (val,ind,res)-> `if`(val >= 3, [ind,res[]], res), []):

That's pretty good timewise and has the benefit of being a one-liner. Timewise it's about midway between Preben's Matrix solution and Pagan's seq solution. But it suffers a bit from the res being constantly appended and passed back and forth. The following is 3 to 4 times faster on the 2000x2000 Matrix that I used for testing. Instead of using res for the results, I store the results in a table.

Ind:= table():
rtable_scanblock(A, [], (val,ind)-> if val >= 3 then :-Ind[ind]:= 0 fi):
ind:= [indices(Ind, nolist)]:

Note that the value 0 that I assign is insignificant; it could be any value. All that matters are the indices of the table that are assigned.

Update: The following variation of Joe Riel's dynamic Array solution (given in a subsequent Answer) fares about the same timewise as my table-based solution:

Ind:= Vector():  k:= 0:
rtable_scanblock(A, [], (val,ind)-> if val >= 3 then :-k:= k+1; :-Ind(k):= ind fi):

The commands diff and convert are not thread safe. See ?index,threadsafe .