@tomleslie 

You have removed a hell of a thorn from my side!

@acer 

I like your second solution.
I vote up

@uzairmukadam 

"Is it because we need the denominator out before starting the elementary operations?"
Yes.
In the present case all the denominators are equal and thus if SA is the SmithForm  of matrix A then the SmithForm SB of this matrix A with all these terms multiplied by the "something" is proportional to "something" times SA.
If this "something" is a polynomial whose leading term as a coefficient equal to C, than SB =(1/C) *~ SA :

restart;
with(LinearAlgebra):
A := Matrix([[1,2*x,2*x^2+2*x],[1,6*x,6*x^2+6*x],[1,3,x]]): # help pages' example
S := factor~(SmithForm(A));

p := randpoly(x, degree=3, coeffs = rand(-5..5));
B  := p*~A:
SB := factor~(SmithForm(B, x)):
simplify(SB /~ p);

Things are more complicated if not all the denominatores of the elements of A are equal.
Let N the matrix defined by N[i,j] = numer(A[i,j])
In this case the structure of Smith forms of A and N are the same (for instance zeros are at the same places), but the previous proportionnality no longer holds.
See for Instance SmithForm.mw

@acer 

There are indeed two safer ways to use the aliases and I am very eager to receive this kind of advices.

About your final question "Why do you use concatenated names instead of indexed names?"
I found this (IMO) a safe way to avoid conflicts, for instance:

T := 1:  # here t represents a variable
:
:
plot( ..., t=0..T); # here I use the same T defined above
:
:
# fine length of the tube named T and of elements 1 and 2

L__T := 28 :
L[1] := 10: L[2]:=20

Another situation where I found using "__" useful is when the "main name" is protected; for instance:

# Temperatures of elements 1, 2, ...
 Temperature[1] := 10
Error, attempting to assign to array/table `Temperature` which is protected.  Try declaring `local Temperature`; see ?protect for details.

Temperature__1 := 10;


# accelerations

gamma[r] := 10;
Error, attempting to assign to array/table `gamma` which is protected.  Try declaring `local gamma`; see ?protect for details.

gamma__r := 10;

As soon as I have decided to work with variable names like L__T instead of L[..], it seemed "natural"  to use the concatenation symbol || to construct them programmatically. 

@acer 

Thanks for this explanation.

About my "slightly suspicious objective":
I study a many degreees of freedom mechanical system. Each of these dof (a mass) is conventionnally refered to by a bame (not a number). 
For instance one can have names [A, B, C].

TO enhance readibility of the differential system whivh governs the dynamics of these masses, I thought to avoid representing the dependent variable.
Here is an example

Download purpose.mw

@Carl Love 

Thank you for this technical detail about your method

Regarding my Answer: I need to have the same bounds for the dimensions to be transposed, even if transposition of rectangular matrices  is well defined (which makes your answer all the more clever)

restart
A := Array(1..3, 1..2, 1..2, (i, j, k) -> i+2*j+3*k):
TA := Array(1..3, 1..2, 1..2, (i, j, k) -> A[j, i, k]):

Error, Array index out of range
      TA[() .. (), () .. (), 1], TA[() .. (), () .. (), 2]

@Carl Love 

I vote up

@vv 

I'm afraid you're right.
Unless one of a, b or c is 0, in such case solutions are easy to find, I've spent a lot of time trying to find a trick.

restart:
sys := [x*y*z + y*z + y = a, x*y*z + x*z + z = b, x*y*z + x*y + x = c]:
print~(sys):
                      x y z + y z + y = a
                      x y z + x z + z = b
                      x y z + x y + x = c

# an eample of a particular case

solve(eval(sys, a=0), [x, y, z])
       [[                 2                           ]  
       [[      c        -b  + b c - b + c          b  ]  
       [[x = - -, y = - -----------------, z = - -----], 
       [[      b                c                b - c]  

                                  ]
         [                    b  ]]
         [x = c, y = 0, z = -----]]
         [                  c + 1]]

# only the second solution can lead to (x, y, z) integers (provided c+1 divides b)
# but not the first one:
#   as b divides c, let's rewrite it while setting c=K*b (K non null integer)
#   to see that z cannot be an integer

simplify(eval(%[1], c=K*b))
           [              K b + K - b - 1        1  ]
           [x = -K, y = - ---------------, z = -----]
           [                     K             K - 1]

One can also prove some facts such that

• a=b=c implies x=y=z, but it remains to prove x is an integer (unless a=b=c=0 which gives x=y=z=0)
• a=b implies x=c-a, y=z=0  (the other solutions can't be real)

I'd put some hope in observing that a+b+c involves the symmetric functions of the root of a 3rd degree equation

add(sys);

       3 x y z + x y + x z + y z + x + y + z = a + b + c

# let sigma__p the sum of the products of p roots of a 3rd degree equation

rel := a+b+c = 3*sigma__1+sigma__2-sigma__3

          a + b + c = 3 sigma__1 + sigma__2 - sigma__3

# thus x, y, and z must be the integer roots of function F below

F := (U-x)*(U-y)*(U-z);
F := collect(expand(F), U)
                    (U - x) (U - y) (U - z)
        3                 2                              
       U  + (-x - y - z) U  + (x y + x z + y z) U - x y z

... but it seemed to be a dead end.

As I know you are (were) involved in IMO, have you not ever heard of this problem?

@tomleslie 

At first glance I have had the same reaction as yours.
It was so obvious that I made myself to wonder if the question of the OP was not this one "How to choose a triple (a, b, c) of integers such that x, y and z would be integers too?"

Or stated differently: given a triple (a, b, c) of integers, and three functions f, g and h of (x, y, z), how can I check if the solution of the system {f(x,y,z)=a, g(x,y,z)=b, h(x,y,z)=c} wrt (x,y,z) is a integer?

Thiinvolve solving third degree equations in integer numbers, which is far from simple.
Unless there is some trick (the problem looks like some asked in mathematics problem books).

@Jaime_mc2 

MY MISTAKE !!!

I've been lured by the way you did the computations.
I hadn't realized that 
In fact if H(n, y) is the nth Hermite polynomial relatively to the measure mu(y)dy, you have written 
HN(n, y) = sqrt(H(n, y)*mu(y)) which makes < HN(n, y), HN(m, y)>=0 wrt the measure dy.
So the HN(n, y)  are the Hermite-Gauss functions... which I hadn't realized first sight

Apologies.

I've just edit my answer

The code you give reminds me of classical problems in probability:

• u7 is what is called the (type ) Gram-Charlier expansion of the density of probability of some RV X.
  In this case the coefficients a0...a6 are determined from the moments of X.
• Sometimes X obey an algebraic/ODE/PDE relation (on aODE in your case ?) and a method to build the pdf/cdf of X is to use "chaos polynomials", for instance Hermite's.
  Here again we obtaine something close to what you wrote.

Any connection with this ?

@Mariusz Iwaniuk 

How do you deduce from your first table that the formula

sum(m!*x^(m - n)*LaguerreL(m, n - m, x)^2/n!, n = 0 .. infinity) = exp(z)

given by the OP is right?

@acer 

Excellent, I vote up!

@acer 

Thank you acer for this very clear answer

@acer 

Hi, 
In my posted remark I used fsolve, but my first guess was to use RootFinding:-NextZero instead.
It seems that the latter is a little bit faster than the former.
Unfortunately RootFinding:-NextZero doesn't seeem capable to to find the two highest positive zeros of ChebyshevT(N, x) for N >= 25 (Maple 2015).
Do you think that  fsolve is more versatile than RootFinding:-NextZero or that some tuning if this latter is possible.

(just by curiosity)