M:=8:  K:=6:  K1:=K-2*K2:

Sum(M!/(M-K1-K2)! * K!/(K1! * K2)! * 1/M^K * 1/2^K2, K2=0..K/2);

value(%);

I do not understand why you need something to simplify. You can easily to define  this matrix as done below and calculate its determinant with any precision for any  k  and  alpha :

Example:

L:=[J,Y,`I`,K]:

A:=Matrix(4,(i,j)->`if`(i=1 or i=2,cat(Bessel,L[j])(i+1,2*k*sqrt(alpha)/(alpha-1)),cat(Bessel,L[j])(i+1,2*k/(alpha-1)))):

F:=unapply(Matrix(4,{seq(seq(`if`(j=3 and (i=2 or i=4),(i,j)=-A[i,j], (i,j)=A[i,j]),j=1..4),i=1..4)}), k,alpha):

F(k,alpha);

evalf[20](LinearAlgebra[Determinant](F(2,3)));

 Edit. The code and the example was corrected (at first I missed minus signs)

I corrected your syntax as I understood it:

f := x->piecewise(abs(x)<1,1-x,0);

U := (x,t)->2/3*f(x+t)+1/3*f(x-2*t)+1/(3*Pi)*sin(Pi*(x+t))+1/(3*Pi)*sin(Pi*(x-2*t));

plot3d(U(x,t), x=0..3, t=0..3, axes=normal, numpoints=10000);

The original problem has an obvious meaning, if one or both parameters  H1  and  H2  are 0 or negative, but the proposed solution does not work for these cases. Here is the solution in the form of a procedure that works for all cases. The procedure is very simple, it is based on direct calculation of the volume by a single integration over the cross sectional area (each a cross-section - this is the usual circle segment).

restart;

V:=proc(R,H1,H2)

local r, H;

if is(H1^2+H2^2>R^2) then error "Should be H1^2+H2^2<=R^2" fi;

r:=int((R^2-x^2)*arccos(H1/sqrt(R^2-x^2))-H1*sqrt(R^2-x^2-H1^2), x=H2..sqrt(R^2-H1^2));

if is(H1>=0) then return simplify(r) else

H:=R-sqrt(R^2-H1^2);

simplify(r+Pi*H^2*(R-H/3)) fi;

end proc:

Examples of use:

V(1,1/2,1/2);

 evalf(%);

V(1,0,1/2);

V(1,1/2,0); 

V(1,-1/2,1/2);

 evalf(%);

V(1,1/2,-1/2);

 evalf(%); 

V(1,0,0); 

V(1,-1/2,-1/2);

 evalf(%);

V(1,-1/sqrt(2),-1/sqrt(2));

 evalf(%);

Here are solutions the same examples by Rouben's formula:

Vwedge :=(R,a,b)->-( (1/3)*Pi*R^3-(1/6)*Pi*(2*R+a)*(R-a)^2-(1/6)*Pi*(2*R+b)*(R-b)^2-(1/3)*a*(3*R^2-a^2)*arcsin(b/sqrt(R^2-a^2))-(1/3)*b*(3*R^2-b^2)*arcsin(a/sqrt(R^2-b^2))-(2/3)*a*b*sqrt(R^2-a^2-b^2)+(1/3)*R^3*arctan((R^2+R*b-a^2)/(sqrt(R^2-a^2-b^2)*a))-(1/3)*R^3*arctan((R^2-R*b-a^2)/(sqrt(R^2-a^2-b^2)*a))):

Vwedge(1,1/2,1/2);

 evalf(%);

Vwedge(1,0,1/2);

Vwedge(1,1/2,0);

Vwedge(1,-1/2,1/2);

 evalf(%);

Vwedge(1,1/2,-1/2);

 evalf(%);

Vwedge(1,0,0);

Vwedge(1,-1/2,-1/2);

 evalf(%);

Vwedge(1,-1/sqrt(2),-1/sqrt(2));

Edit.  The procedure  V  worked incorrectly if  H1<0 . So the procedure and all the examples was corrected. Now all right. The number of examples was increased. 

plot3d([sin(u)*(7+cos(u/3-2*v)+2*cos(u/3+v)), cos(u)*(7+cos(u/3-2*v)+2*cos(u/3+v)),

sin(u/3-2*v)+2*sin(u/3+v)], u=-Pi..Pi, v=-Pi..Pi,  axes=normal, orientation=[30,30], lightmodel=light4);

Edited: the extra commas was removed. 

To the points of not too often flashed for every point allotted 10 frames:

L := [seq([i, i^2]$10, i = 1 .. 4)];

Frames := seq(plots[display](plots[pointplot](L[2 .. i], style = point, symbol = solidcircle, symbolsize = 12, color = red)), i = 1 .. nops(L)):

plots[display](Frames, insequence = true, view = [0 .. 4.5, 0 .. 16]);

This will be true if we assume that all constants are positive:

temp := (1/2)*(a*r*t-b*p)*sqrt(p+v*sqrt(a*r))*sqrt(b)/(sqrt(a*b*r*t)*sqrt(a*r*t*(p+v*sqrt(a*r))))+(1/2)*(a*r*v*t-v*b*p+p*t*sqrt(a*r)-b*v^2*sqrt(a*r))*sqrt(b)/(sqrt(a*b*r*t*(p+v*sqrt(a*r)))*sqrt(t*(p+v*sqrt(a*r))))+(1/2)*b*(p+v*sqrt(a*r))/(a*r*t);

normal(expand(rationalize(simplify(temp))))  assuming positive;

                                               1

You can decide in short, if you already know the answer:

is(temp = 1)  assuming positive;

                                             true

I think that in order to achieve complete clarity in finding all real roots of the equation  a=1-(1+x)*exp(-x)  for any real  a  is useful to make the simplest study of the function defined of the right-hand side of the equation:

f:=x->1-(1+x)*exp(-x):

g:=diff(f(x), x);

map(solve, [g=0, g<0, g>0]);

minimize(f(x), x=-infinity..0, location),  maximize(f(x), x=-infinity..0, location);

minimize(f(x), x=0..infinity, location),  maximize(f(x), x=0..infinity, location);

plot([f(x), 1], x=-1.5..5, thickness=[2,1], linestyle=[1,2], scaling=constrained);

We can see (and the plot confirms) that the function  f  is strictly decreasing in the domain  RealRange(-infinity,0) and takes all the values from  RealRange(0,infinity)  and  f  is strictly increasing in the domain  RealRange(0,infinity) and takes all the values from  RealRange(0, Open(1)) . 

On the basis of this information written procedure  SS , which finds all the real roots of the equation with a given accuracy (the numbers of digits):

restart;

SS:=proc(a::realcons, n::posint:=10)

local eq, b;

Digits:=n;

b:=evalf(a);

eq:=a=1-(1+x)*exp(-x);

if b<0 then return `No real solutions` elif

b=0. then return x=0 elif

b>0 and b<1 then return x[1]=fsolve(eq, x=-infinity..0), x[2]=fsolve(eq, x=0..infinity) else return x=fsolve(eq, x=-infinity..0) fi;

end proc:

Examples of use:

SS(-2);  SS(0);  SS(1/Pi);  SS(1);  SS(3, 20);

This system can be solved instantly (and we obtain all solutions) by  solve  command:

restart;

eq[1]:= 0.223569*c_1+2.35589*c_2*c_1^2+0.002356*c_1*c_2^2:

eq[2]:= 1.277899*c_1*c_3-2.350023*c_2*c_3^2+7.5856*c_3*c_2^2:

eq[3]:= 3.225989*c_1^2-2.35589*c_3*c_1^2-7.28356*c_3*c_2^3:

solve({eq[1], eq[2], eq[3]}, {c_1, c_2, c_3});

Your example can be easily generalized to any natural number, if you use the binary number system. In it  any number  S such that  1<= S <2^n  can be recorded as the sequence of no more than  n  of  "0"  and  "1" .

The first procedure for a given number  S  returns the list of the number of dollars in each bag that your condition is true. The second procedure returns the corresponding partition for any number.

MinBags:=proc(S)

[seq(2^n, n=0..ilog[2](S))];

end proc:

Parts:=proc(s)

local L;

L:=convert(s, base, 2);

convert(s,symbol)=convert([seq(`if`(L[i]<>0,convert(L[i]*2^(i-1),symbol),NULL), i=1..nops(L))], `+`);

end proc:

Examples.

Your example:

MinBags(15);

for n from 1 to 15 do

Parts(n);

od;

Another example:

MinBags(137);

for n from 100 to 110 do

Parts(n);

od;

You can easily solve the problem by converting the output into the list, and reversing it:

A:=<1,1,1,1; 1,2,3,4; 4,3,2,1>;

LinearAlgebra[NullSpace](A);

ListTools[Reverse](convert(%, list));

L:={1,2,3,4,5,6,7,8}:

P:=combinat[choose](L, 3):

k:=0:

for p in P do

a:=p[1]; b:=p[2]; c:=p[3];

if a+b>c and a+c>b and b+c>a then k:=k+1 fi;

od:

k;

                        22

Addition:  faster code for large sets.

Example for  L:={$ 1..500}:

restart;

ts:=time():

k:=0:

for a from 1 to 498 do

for b from a+1 to 499 do

for c from b+1 to 500 do

if a+b>c and a+c>b and b+c>a then k:=k+1 fi;

od: od: od:

k;

time()-ts;

                                 10323125

                                    23.875

The first code fails in this example.

solves the problem for  specific parameters values.

Example:

Optimization[QPSolve](2*x+5*y+3*x^2+3*x*y+2*y^2-x*z+4*y*z+2*z^2, {x+y+z=1, 2*x-3*y+z=5});

         [-4.18080357142857, [x = 0.410714285714286, y = -0.897321428571428, z = 1.48660714285714]]

ben maths  You wrote about the mission to give a better code. Here is the shorter and faster code: 

Simp38:=proc(f,x0,xn,n)

local  h:=(xn-x0)/n,  x:=i->x0+i*h;

3*h/8*(f(x0)+f(xn)+3*add(f(x(i)),i=1..n-1,3)+3*add(f(x(i)),i=2..n-1,3)+2*add(f(x(i)),i=3..n-1,3));

end proc:

Test:

ts:=time():

simp38(x->x^5-5.15*x^2+8.55*x-4.05045,1,5,3000);

int(x^5-5.15*x^2+8.55*x-4.05045,x=1..5);  

time()-ts;

ts:=time():

Simp38(x->x^5-5.15*x^2+8.55*x-4.05045,1,5,3000);

int(x^5-5.15*x^2+8.55*x-4.05045,x=1..5);

time()-ts; 

plots[inequal](not (x>4 and y>4), x=-10..10, y=-10..10, color=yellow);