@Matt C Anderson  
Your method for constructing the set of desired objects using  b := `union`(b, {a}) construction is extremely inefficient. This isn't particularly evident in your example of finding twin primes, as they are quite rare, but in other cases, the difference compared to other methods would be very significant.

The examples below compare the construction of the set of all odd primes less than 1,000,000 using your method and another method:

restart;
primes:=proc(N) # First way (your method)
local a,b;
b:={};
for a by 2 to N do
if isprime(a)  then b:=`union`(b,{a})
end if;
end do;
return b;
end proc:

primes1:=proc(N) # Second way
local a,b,k;
k:=0:
for a by 2 to N do
if isprime(a)  then k:=k+1; b[k]:=a 
end if;
end do;
return convert(b,set);
end proc:

CodeTools:-Usage(primes(1000000)):

#  memory used=23.16GiB, alloc change=28.49MiB, cpu time=91.33s, real time=96.09s, gc time=3.86s
 

CodeTools:-Usage(primes1(1000000)):

#  memory used=217.08MiB, alloc change=40.00MiB, cpu time=2.28s, real time=2.38s, gc time=171.88ms

Thanks everyone for the helpful comments. I hope this long-standing bug is finally fixed in future Maple releases.

@C_R  If you need to simplify a subexpression (in your example, this is the expression in the numerator under the integral sign), then you can the simplifying procedure apply to the corresponding operand.

restart;
expr := Int(sqrt(-2*alpha^2*x - 2*alpha^2 + 4)/sqrt(-alpha^2 + 1)/(-alpha^2*x - alpha^2 + 2), alpha = 0 .. z);

applyop(e->content(e)*primpart(e), [1,1], expr);

If you open the brackets in the expression  dio1 , you will notice that all of their coefficients for  dio  and  dio1  are the same, except for one.

restart; 
dio := 3*x^2*y^2*z^2-36*x^2*y^2*z-30*x^2*y*z^2-24*x*y^2*z^2+126*x^2*y^2+432*x^2*y*z+93*x^2*z^2+360*x*y^2*z+312*x*y*z^2+66*y^2*z^2-1692*x^2*y-1476*x^2*z-144*x*y^2-5040*x*y*z-1104*x*z^2-1080*y^2*z-948*y*z^2+6120*x^2+21096*x*y+18576*x*z+4554*y^2+16056*y*z+3540*z^2-79848*x-69300*y-61272*z+268452;
dio1:=add(((x-i)*(y-i-1)*(z-i-2))^2, i=1..7,3);
dio2:=expand(dio1);
is(dio-dio1=0);
coeffs(dio);
coeffs(dio2);

This example is fairly easy to solve by hand without Maple. First, we remove the irrationality from the denominator, then simplify the expression under the external root and represent it as a square.

Below are the main step-by-step transformations:

There is no error in Maple 2018.2. Maple simply returns this integral unevaluated:

@janhardo Yes, you're right. The  ~  command isn't needed here.

@sand15  Thank you.

@imparter 
 

restart;
with(plots): with(plottools):
h:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))):
K1:=((4/h(z)^4)-(sin(alpha)/F)-h(z)^2+Nb*h(z)^4):
T1:=Int(K1,z=0..1):
L1:=0.2: F:=10:
d1:=0.2:
T2:=evalf(T1):
T3:=unapply(T2,alpha,Nb,delta2):
T:=textplot([[0.037,4.18,"Color Representation",color=black],[0.04,4.14,"Nb=0.02",color=blue],[0.04,4.1,"Nb=0.03",color=red],[0.04,4.06,"Nb=0.04",color=green]], font=[times,bold,16]):
Lines:=line([0.033,4.14],[0.035,4.14],color=blue,thickness=10),line([0.033,4.1],[0.035,4.1],color=red,thickness=10),line([0.033,4.06],[0.035,4.06],color=green,thickness=10):                     
T11:=textplot([[0.067,4.18,"Line Style"],[0.07,4.14,typeset(alpha=Pi,"/",6)],[0.07,4.1,typeset(alpha=Pi,"/",4),line=longdash],[0.07,4.06,typeset(alpha=Pi,"/",3),line=dash]], font=[times,bold,16]):
Lines1:=line([0.03+0.03,4.14],[0.035+0.03,4.14],linestyle=solid,color=blue),line([0.03+0.03,4.1],[0.035+0.03,4.1],linestyle=longdash,color=red),line([0.03+0.03,4.06],[0.035+0.03,4.06],linestyle=spacedash,color=green):
P:=plot([seq(seq(T3(alpha, Nb, delta2), Nb in [0.1, 0.2,0.3]), alpha in [Pi/6, Pi/4, Pi/3])], delta2 = 0.02.. 0.1,  titlefont = ["ROMAN", 15], labels = ["δ1", "T3"], labeldirections = ["horizontal", "vertical"], labelfont = ["SYMBOL", 11], linestyle = [solid, solid,solid,longdash, longdash,longdash, spacedash,spacedash, spacedash], symbol = [BOX, CROSS, CIRCLE, BOX, CROSS, CIRCLE, BOX, CROSS, CIRCLE], color = [blue,red, green, blue,red, green,blue,red, green], axes = boxed,size = [1000, 1000]):
plots:-display(P,T,T11,Lines,Lines1);

@janhardo  Yes,  minimize  works symbolically, but  the Optimization package - only numerically.

@janhardo  

f:=419*x^2 + 116*x*y - 426*x*z + 78*y^2 - 142*y*z + 133*z^2 - 1604*x - 682*y + 1086*z + 2306:
minimize(f, location);

@Alfred_F  I think Maple does this exactly symbolically, for example

is(sqrt(13)+sqrt(468)=sqrt(637));

                                 true

But

is(sqrt(13)+sqrt(468)=sqrt(637)+10^(-100));

                                 false

@dharr  But if we use the  is  command to check for identity, we don't get the desired result:

is(sqrt(x+2*sqrt(x-1))=sqrt(x-1)+1);
is(sqrt(x+2*sqrt(x-1))=sqrt(x-1)+1) assuming x>=1 ;

When trying to calculate this sum symbolically, we get an incorrect result (returns  0 ):

restart;
sum(sum(`if`(igcd(m,n)=1, 1/(m^2*n^2), 0), m=1..infinity), n=1..infinity);

                                   0

@nm  Maple probably thinks that  B  is no simpler than  A . 
It's also worth noting that  expand  is a standard command in Maple for breaking an algebraic fraction (with a sum in the numerator) into the sum of several fractions:

expand((a+b+c)/d);
expand((a+b+c)/(d+e));