If the question is as I suppose (but that guy seems to be too lazy
to answer the implicite question), then I think there is not extremum
except at the boundaries.

E0:=Int(q*exp((1-b)^2* z^2/lambda )/((1+b^2*q^2+b^2*z^2)*(q^2+z^2)),
  z=-infinity..infinity);

E:=eval(E0, lambda=-alpha^2); # since 0 < lambda ===> infinite
E1:=PDETools[dchange](z = x/(1-b)*alpha, %, [x]) assuming 0 < alpha:

That depends on sign( b-1 ), but changes only the sign (use 'Flip').

S:=[b = 1/(s^(1/2)*beta), q = p^(1/2)*beta, alpha = beta-1/(s^(1/2))];

  E1 assuming b-1 < 0;
  eval(%, S);
  PS:=simplify(%) ;

                       infinity
                      /                    2   1/2
                     |             s exp(-x ) p
              PS :=  |          --------------------- dx
                     |                2    2
                    /           (p + x ) (x  + p + s)
                      -infinity

  plot3d(PS, s=0..2,p=20..1, axes=boxed);

This plot suggests, that the value is decreasing with increasing p
(but I am too lazy too for trying a proof), which should give that 
there is no finite minimum.

I think he wants to minimize an expectation value w.r.t. that q
and that the integral (up to factors) can be written as 

  Int( exp(-x^2) * f(x,q), x = -infinity .. infinity), 
  f(x,q)=q/polynom(x,q)

where that polynomial in x has 4 distinct complex roots.

Maple does not give the integral in closed form (even after using
(real) partial fractions), i.e. the special function(s) involved
are possibly not of known form. And even if it would, it is not
clear that if would lead to some usefull result for q_min.

And 0 = diff(theIntegral,q) looks not that nice to solve ...

May be VariationalCalculus can be used. But I guess, one needs
sound Math, not only Maple, for a usefull description of the
desired aim.

I always forget the 'distributed' option for collect ...

One can also use coeftayl(p,c0=0,1) or coeff(p,c0,1)

Have not looked into ListTools, but the following seems to be what you are looking for:

  g:=[seq('(e[i],f[i])',i=1..8)]; 

        g := [3, 5, 7, 9, 6, 3, 6, 7, 5, 4, 4, 6, 3, 9, 2, 6]

This is not a polynomial, but a Taylor series - make it a polynomial:

  p:=convert(%,polynom);

Then you can use pattern matching

  pattern:= F*c0+G*c1+H*c2+K*c3;
  match(p = pattern , {c0,c1,c2,c3}, 's'); 
  s;

The polynomial function you are looking for then are obvious and
can be done like this:

  L:=convert(s,list);

  rhs(L[1]): collect(%,x);
  f:= unapply( %, x );

                        4           6  2  2             3  4          2   8
    f := x ->  1 - 2/3 x  a - 1/45 x  a  R  + (-1/2520 a  R  + 8/315 a ) x

The same for the rest ...

Thank you both, I see.

Ok, I can not solve it, the question is from groups.google.de/group/sci.math.symbolic/browse_frm/thread/456afbdf28007cbf/
where the notations are a bit different from Manzoni's post (but his final task is just the problem for me)

you should not mix lower case and upper case for sums:

  S:=Sum((A*x[i]+B)^2, i = 1 .. n); expand(%); value(%); 
  #subs(sum=Sum, %); # what I prefer here ...

                /  n        \         /  n       \
                |-----      |         |-----     |
              2 | \        2|         | \        |    2
             A  |  )   x[i] | + 2 A B |  )   x[i]| + B  n
                | /         |         | /        |
                |-----      |         |-----     |
                \i = 1      /         \i = 1     /

  interface(version);

      Classic Worksheet Interface, Maple 11.02, Windows, Nov 10 2007 Build ID 330022

  assume(u::real, 0 <v, 0 < lambda); # your assumptions
  (u+i*v / u+i*v+lambda)^4;          # your expression
  subs(i=I,%):                       # the correct imag I = sqrt(-1)
  Im(%);                             # does not work
  evalc(%);                          # but that's the better way to use it
  simplify(%,size); lprint(%);       # let it look nicer

           4*(u^2+(lambda-v)*u-v)*v*(1+u)*(u^2+(lambda+v)*u+v)*(u+lambda)/u^3

  Sum((A*x[i]+B)^2, i = 1 .. n) = A^2*(sum(x[i]^2, i = 1 .. n))+2*A*B*(sum(x[i], i = 1 .. n))+n*B^2;
  subs(sum=Sum, %);
  expand(%); 
  value(%); 
  is(%);

                                 true

where I prefer to test  lhs(%) - rhs(%) to be zero

first: if you have a question you should break it up into the essential part - I will not read 85 pages

second: essentially you ask how to handle an array as inpu, yes? You should try it with some simple code, say 2 or 3 lines

third: you have only 3 or 4 dimensions for your CS - in that case I simply would not invest much time and use 3 or 4 parameters cs1, ... , cs4 and if I really would need that as an array, then I would write a simple wrapper (with array notation) to call the function (with parameter notation)

for that Y3:= subs(CS(1)=cs1, CS(2)=cs2, CS(3)=cs3, CS(4)=cs4, Y3) should work and one can proceed as suggested

finally you should aware, that the resulting Fortran code may cause numerical problems - at least it will be not quite clear that it is 'stable'

so perhaps it is reasonable to break it up before (for better control) and not to pump all into 1 function - but that is on you, you have to know why you want it your way

at least you need intensive tests

i would try to evaluate, which gives a LerchPhi and to convert to hypergeometrics, which is a 2F1 and there are papers by Nico Temme for large parameters - may be they cover it, otherwise i would look into his references for other roads

rationalize(A); Int(%,x=0..infinity); value(%);

                                   1/2
                               -2 2

For Numerics the modulo function is not quite handy :-(

The following should do, what you expect using it as plot('p'(x),x=-6..4):

  p:=proc(x)
    local r,m;
    m:=63;
    r:=irem(floor(16*x),m);
    `if`(0<=r,r,r+m);
  end proc;

Some tests:

  p(21.2);
  p(0);
  p(4);
  p(4-1/16);

                                  24
                                  0
                                  1
                                  0

  4-1/16; %*16; `mod`(%,63);

                                  0

of course "convert(%,Si)" works - and use Pi, not pi ...

a multiplication sign is missing (trapped by the interface?):

g*(w-c)+(1/4)*(a-g*w+g*r)^2/(b*(p-k))-(a-g*w+2*g*r)*(a-g*w+g*r)/(2*b*(p-k));

..................................................^^^......................