throw(x);
x;

asksign(x);
if x < 0 then neg elseif equal(x,0) then zero else pos;

asksign(log(x));
if log(x) < 0 then neg elseif equal(log(x),0) then zero else pos;

asksign(abs(x));
if equal(x,0) then zero else pos;

integrate(exp(a*x),x,0,inf);
if a < 0 then -1/a elseif equal(a,0) then 0 else 0;

integrate(x^k,x);
if equal(k+1,0) then log(x) else x^(k+1)/(k+1);

integrate(x^k,x,1,b);
if k < 0 then integrate(x^k,x,1,b) elseif equal(k,0)
    then b^(k+1)/(k+1)-1/(k+1) else b^(k+1)/(k+1)-1/(k+1);

integrate(x^a/(1+x),x,0,inf);
if a+1 < 0 then integrate(x^a/(x+1),x,0,inf) elseif equal(a+1,0)
    then 'integrate(x^a/(x+1),x,0,inf) else integrate(x^a/(x+1),x,0,inf);

integrate(x^a/(1+x)^(5/2),x,0,inf);
if a+1 < 0 then 'integrate(x^a/(x+1)^(5/2),x,0,inf) elseif equal(a+1,0)
    then integrate(x^a/(x+1)^(5/2),x,0,inf)
    else integrate(x^a/(x+1)^(5/2),x,0,inf);

integrate(sqrt(1-s^2)/(z-s),s,-1,1);
if z-1 < 0 then integrate(sqrt(1-s^2)/(z-s),s,-1,1) elseif equal(z-1,0)
    then 2*asin(1/z)*z else %pi*z-%pi*sqrt(z^2-1);

(xmax(a,b):=(abs(a-b)+b+a)/2,integrate(xmax(a-x,0),x,0,b));
((b-a)*abs(b-a)-b^2+2*a*b+a*abs(a))/4;

integrate(1/(1+a^2*cos(t)^2),t,0,%pi);
if 2*sqrt(a^2+1)-a^2-2 < 0 then integrate(1/(a^2*cos(t)^2+1),t,0,%pi)
    elseif equal(2*sqrt(a^2+1)-a^2-2,0)
    then ?merror("atan2: atan2(0,0) is undefined.")
    else integrate(1/(a^2*cos(t)^2+1),t,0,%pi);

integrate(exp(-la*t)*la,t,0,inf);
if la < 0 then ?merror("defint: integral is divergent.") elseif equal(la,0)
    then 0 else 1;

integrate(sin(5*x)*exp(-s*x),x,0,inf);
if s < 0 then 'integrate(%e^-(s*x)*sin(5*x),x,0,inf) elseif equal(s,0)
    then 'integrate(%e^-(s*x)*sin(5*x),x,0,inf) else 5/(s^2+25);

integrate(1/x,x,a,4);
if a-4 < 0 then integrate(1/x,x,a,4) elseif equal(a-4,0) then log(4)-log(a)
    else log(4)-log(a);

integrate(x^k,x,0,1);
if k < 0 then integrate(x^k,x,0,1) elseif equal(k,0) then 1/(k+1) else 1/(k+1);

integrate(1/sqrt(i+n),i,1,n);
if n-1 < 0 then integrate(1/sqrt(n+i),i,1,n) elseif equal(n-1,0)
    then 2^(3/2)*sqrt(n)-2*sqrt(n+1) else 2^(3/2)*sqrt(n)-2*sqrt(n+1);

integrate(1/(a^2*sin(x)^2+1),x,0,3*%pi);
if 2*sqrt(a^2+1)-a^2-2 < 0 then integrate(1/(a^2*sin(x)^2+1),x,0,3*%pi)
    elseif equal(2*sqrt(a^2+1)-a^2-2,0)
    then integrate(1/(a^2*sin(x)^2+1),x,0,3*%pi)
    else integrate(1/(a^2*sin(x)^2+1),x,0,3*%pi);

integrate(integrate(p2^2,x1,(-4*%pi)/kp,4*%pi/kp),x2,5*%pi/kp/2,7*%pi/kp/2);
8*%pi^2*p2^2/kp^2;

integrate(1/log(t),t,x,2*x);
if x < 0
    then -log(log(2*x))+expintegral_ei(log(2*x))+log(log(x))
                       -expintegral_ei(log(x))-log(log(-x)^2+%pi^2)/2
                       -%i*atan(%pi/log(-x))+log(log(-2*x)^2+%pi^2)/2
                       +%i*atan(%pi/log(-2*x)) elseif equal(x,0)
    then ?merror("log: encountered log(0).") else integrate(1/log(t),t,x,2*x);

integrate((6*sin(8*x)+6*cos(9*x))*(j-x),x,-%pi,%pi);
(81*j+81*%pi-8)/108-(81*j-81*%pi-8)/108;

integrate(cos(n*x)*sin(n*x),x,0,t);
if t < 0 then 1/(2*n)-cos(n*t)^2/(2*n) elseif equal(t,0) then 0
    else 1/(2*n)-cos(n*t)^2/(2*n);

integrate(1/(1+sin(x)^2),x,0,z);
if z < 0 then atan(sqrt(2)*sin(z)/cos(z))/sqrt(2) elseif equal(z,0)
    then atan(sqrt(2)*sin(z)/cos(z))/sqrt(2)
    else integrate(1/(sin(x)^2+1),x,0,z);

integrate(cos(x)-cos(x-c),x,0,c/2) = -integrate(cos(x-c),x,c,c+%pi/2)-c+%pi;
if c < 0 then 2*sin(c/2)-sin(c) = -c+%pi-1 elseif equal(c,0)
    then 2*sin(c/2)-sin(c) = -c+%pi-1 else 2*sin(c/2)-sin(c) = -c+%pi-1;

integrate(cot(x),x,0,%pi/2);
?merror("defint: integral is divergent.");

integrate(r,z,-sqrt(20-r^2),sqrt(20-r^2));
?merror("defint: lower limit of integration must be real; found ~M",
        -sqrt(20-r^2));

integrate((2+x+x^2)/(c+x+x^2),c);
(x^2+x+2)*log(x^2+x+c);

specint(sqrt(t)*%e^(-a/t-p*t),t);
if p < 0 then 'specint(sqrt(t)*%e^(-p*t-a/t),t) elseif equal(p,0)
    then 'specint(sqrt(t)*%e^(-p*t-a/t),t)
    else specint(sqrt(t)*%e^(-p*t-a/t),t);

(f(x):=x^2/4-2*log((1+exp(x))/2)+x,
 exp((-x^2)/4)*ratsimp(taylor(exp(f(x))*(mu+x*kb)^gamma,x,0,4)),
 integrate(%%,x,minf,inf));
(24*sqrt(%pi)
   *(8*kb^4*mu^gamma*gamma^4-48*kb^4*mu^gamma*gamma^3+88*kb^4*mu^gamma*gamma^2
                            -48*kb^4*mu^gamma*gamma+2*mu^(gamma+4))
 +4*sqrt(%pi)*(96*kb^2*mu^(gamma+2)*gamma^2-96*kb^2*mu^(gamma+2)*gamma)
 +384*sqrt(%pi)*mu^(gamma+4))
 /(192*mu^4);

(aa:1+3*kbt*x/(2*mu)+3*kbt^2*x^2/(8*mu^2)-kbt^3*x^3/(16*mu^3)
     +3*kbt^4*x^4/(128*mu^4),integrate(exp(x)*aa/(1+exp(x))^2,x,minf,inf));
(640*mu^4+80*%pi^2*kbt^2*mu^2+7*%pi^4*kbt^4)/(640*mu^4);

(tmpp:(x1/2+x2/2)*sin(%pi*(x1/2+x2/2-(1-x1-x2-x3)/2+(-x3)/2)/2),
 integrate(integrate(integrate(tmpp,x3,0,-x2-x1+1),x2,0,1-x1),x1,0,1));
if x2+x1-1 < 0
    then (tmpp:(x2/2+x1/2)*sin(%pi*(-x3/2-(-x3-x2-x1+1)/2+x2/2+x1/2)/2),
          integrate(integrate(integrate(tmpp,x3,0,-x2-x1+1),x2,0,1-x1),x1,0,
                    1)) elseif equal(x2+x1-1,0)
    then (tmpp:(x2/2+x1/2)*sin(%pi*(-x3/2-(-x3-x2-x1+1)/2+x2/2+x1/2)/2),
          integrate(integrate(integrate(tmpp,x3,0,-x2-x1+1),x2,0,1-x1),x1,0,
                    1))
    else (tmpp:(x2/2+x1/2)*sin(%pi*(-x3/2-(-x3-x2-x1+1)/2+x2/2+x1/2)/2),
          integrate(integrate(integrate(tmpp,x3,0,-x2-x1+1),x2,0,1-x1),x1,0,
                    1));

limit(a*x,x,inf);
infinity;

limit(x^a,x,inf);
if a < 0 then 0 elseif a > 0 then inf;

limit(x^a,x,minf);
if a < 0 then 0 elseif a > 0 then minf;

declare(a,integer);
done;
limit(x^a,x,minf);
if a < 0 then 0 elseif a > 0 then minf;

tlimit(s/(s^2+1)/sinh(s*T),s,inf);
if T < 0 then 0 elseif equal(T,0)
    then 'limit(s/(s^2*sinh(s*T)+sinh(s*T)),s,inf) else 0;

solve(x = log(sqrt(a^2-r^2)+a),r);
if a-%e^x < 0 then r = -%e^(x/2)*sqrt(2*a-%e^x)
    elseif r = %e^(x/2)*sqrt(2*a-%e^x) then equal(a-%e^x,0)
    elseif r = -%e^(x/2)*sqrt(2*a-%e^x) then r = %e^(x/2)*sqrt(2*a-%e^x)
    else sqrt(a^2-r^2) = %e^x-a;

solve(sqrt(y)+sqrt(x) = sqrt(a),y);
if sqrt(a)-sqrt(x) < 0 then sqrt(y) = sqrt(a)-sqrt(x)
    elseif equal(sqrt(a)-sqrt(x),0) then y = x-2*sqrt(a)*sqrt(x)+a
    else y = x-2*sqrt(a)*sqrt(x)+a;

solve(a^2*e^2+2*a*e*x+x^2+y^2 = (2*a-sqrt(a^2*e^2-2*a*e*x+x^2+y^2))^2,y);
if a-e*x < 0 then sqrt(y^2+x^2-2*a*e*x+a^2*e^2) = a-e*x elseif equal(a-e*x,0)
    then y = -sqrt(e^2*x^2-x^2-a^2*e^2+a^2)
    elseif y = sqrt(e^2*x^2-x^2-a^2*e^2+a^2) then true
    elseif y = -sqrt(e^2*x^2-x^2-a^2*e^2+a^2)
    then y = sqrt(e^2*x^2-x^2-a^2*e^2+a^2);

solve(b^(a*t) = 1,b);
?merror("~:M: improper argument: ~M",declare,a*t);

ode2(a*x+'diff(x,t,2) = b*sin(c*t),x,t);
if a < 0
    then x = -b*sin(c*t)/(c^2-a)+%k1*%e^(%i*sqrt(a)*t)+%k2*%e^-(%i*sqrt(a)*t)
    elseif equal(a,0) then x = -b*sin(c*t)/c^2+%k2*t+%k1
    else x = -b*sin(c*t)/(c^2-a)+%k1*sin(sqrt(a)*t)+%k2*cos(sqrt(a)*t);

ode2(b*x+a*'diff(x,t)+'diff(x,t,2) = 0,x,t);
if 4*b-a^2 < 0
    then x = %k1*%e^((sqrt(a^2-4*b)-a)*t/2)+%k2*%e^((-sqrt(a^2-4*b)-a)*t/2)
    elseif equal(4*b-a^2,0) then x = (%k2*t+%k1)*%e^-(a*t/2)
    else x = %e^-(a*t/2)*(%k1*sin(sqrt(4*b-a^2)*t/2)
                         +%k2*cos(sqrt(4*b-a^2)*t/2));

ode2(b*x+a*'diff(x,t)+'diff(x,t,2) = d*sin(e*t),x,t);
if 4*b-a^2 < 0
    then x = -((d*e^2-b*d)*sin(e*t)+a*d*e*cos(e*t))/(e^4+(a^2-2*b)*e^2+b^2)
           +%k1*%e^((sqrt(a^2-4*b)-a)*t/2)+%k2*%e^((-sqrt(a^2-4*b)-a)*t/2)
    elseif equal(4*b-a^2,0)
    then x = (%k2*t+%k1)*%e^-(a*t/2)-((16*d*e^2-4*a^2*d)*sin(e*t)
                                     +16*a*d*e*cos(e*t))
                                     /(16*e^4+8*a^2*e^2+a^4)
    else x = %e^-(a*t/2)*(%k1*sin(sqrt(4*b-a^2)*t/2)
                         +%k2*cos(sqrt(4*b-a^2)*t/2))
           -((d*e^2-b*d)*sin(e*t)+a*d*e*cos(e*t))/(e^4+(a^2-2*b)*e^2+b^2);

ode2(a*y+'diff(y,x,2) = c*x,y,x);
if a < 0 then y = %k1*%e^(%i*sqrt(a)*x)+%k2*%e^-(%i*sqrt(a)*x)+c*x/a
    elseif equal(a,0) then y = c*x^3/6+%k2*x+%k1
    else y = %k1*sin(sqrt(a)*x)+%k2*cos(sqrt(a)*x)+c*x/a;

ode2(b*y+a*'diff(y,x)+'diff(y,x,2) = c*x,y,x);
if 4*b-a^2 < 0
    then y = %k1*%e^((sqrt(a^2-4*b)-a)*x/2)+%k2*%e^((-sqrt(a^2-4*b)-a)*x/2)
                                           +(b*c*x-a*c)/b^2
    elseif equal(4*b-a^2,0)
    then y = (%k2*x+%k1)*%e^-(a*x/2)+(4*a*c*x-16*c)/a^3
    else y = %e^-(a*x/2)*(%k1*sin(sqrt(4*b-a^2)*x/2)
                         +%k2*cos(sqrt(4*b-a^2)*x/2))
           +(b*c*x-a*c)/b^2;

block([gcd:spmod],expr:sin(a)/(1-cos(a)^2*sin(t)^2),
      assume(cos(a) > 0,sin(a) > 0),integrate(expr,t,0,2*%pi));
if cos(a)^2-1.0 < 0 then -%pi*sin(a)^2/(cos(a)^2-1) elseif cos(a)^2-1.0 > 0
    then block([gcd:spmod],expr:sin(a)/(1-cos(a)^2*sin(t)^2),
               assume(cos(a) > 0,sin(a) > 0),integrate(expr,t,0,2*%pi));

(load(simplify_sum),sum(r*binom(n,r)*binom(m,k-r)/k/binom(m+n,k),r,0,k),
 simplify_sum(intosum(%%)));
?merror(?\~M\ non\-rational\ term\ ratio\ to\ nusum,
        binom(m,-r+k-1)*binom(n,r+1)*(r+1)/(binom(m,k-r)*binom(n,r)*r));
