throw(x);
x;

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

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

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

integrate(exp(a*x),x,0,inf);
if a < 0 then -(1/a) elseif equal(a,0)
    then merror("defint: integral is divergent.")
    elseif a > 0 then merror("defint: integral is divergent.");

integrate(x^k,x);
if equal(k,-1) then log(x) elseif notequal(k, -1) then 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 - 1 elseif k > 0 then 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) elseif a+1 > 0 then 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)
    elseif a+1 > 0 then 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 elseif z-1 > 0 then %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 integrate(1/(a^2*cos(t)^2+1),t,0,%pi)
    elseif 2*sqrt(a^2+1)-a^2-2 > 0 then 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 elseif la > 0 then 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(sin(5*x),x,0,inf) elseif s > 0 then 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)
    elseif a-4 > 0 then 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 elseif k > 0 then 1/(k+1);

/* apparently doesn't provoke asksign anymore -- skip it
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) elseif n-1 > 0 then 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)
    elseif 2*sqrt(a^2+1)-a^2-2 > 0 then 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 gamma_incomplete(0,-log(x))-gamma_incomplete(0,-log(2*x)) elseif equal(x,0)
    then 0 elseif x > 0 then 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
    elseif t > 0 then 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 0
    elseif z > 0 then 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 0 = %pi - 1 elseif c > 0 then 2*sin(c/2)-sin(c) = -c+%pi-1;

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

/* apparently doesn't provoke error; skip this one
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^(-a/t),t)
    elseif p > 0 then 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 ''(-1 + x1 + x2) < 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))
    elseif ''(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));

limit(a*x,x,inf);
if a < 0 then minf elseif equal(a, 0) then 0 elseif a > 0 then inf;

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

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

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

tlimit(s/(s^2+1)/sinh(s*T),s,inf);
if T < 0 then 0 elseif equal(T,0)
    then ?merror("expt: undefined: 0 to a negative exponent.") elseif T > 0 then 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),r = %e^(x/2)*sqrt(2*a-%e^x)]
    elseif equal(a-%e^x,0)
    then [r = -%e^(x/2)*sqrt(2*a-%e^x),r = %e^(x/2)*sqrt(2*a-%e^x)]
    elseif a-%e^x > 0 then [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]
    elseif sqrt(a)-sqrt(x) > 0 then [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),
          y = sqrt(e^2*x^2-x^2-a^2*e^2+a^2)] elseif a-e*x > 0
    then [y = -sqrt(e^2*x^2-x^2-a^2*e^2+a^2),
          y = sqrt(e^2*x^2-x^2-a^2*e^2+a^2)];

solve(b^(a*t) = 1,b);
if featurep(t, integer) then [b = 1] elseif featurep(t, noninteger) then [b = 1];

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
    elseif a > 0 then 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)
    elseif 4*b-a^2 > 0 then 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)
    elseif 4*b-a^2 > 0
    then 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
    elseif a > 0 then 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
    elseif 4*b-a^2 > 0 then 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 sin(a)/(1-cos(a)^2*sin(t)^2)$

/* doesn't cause a question, just skip it
(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));
 */

integrate(1/(x^2-a),x);
if a < 0 then atan(x/sqrt(-a))/sqrt(-a) elseif a > 0
    then log((2*x-2*sqrt(a))/(2*x+2*sqrt(a)))/(2*sqrt(a))$

integrate(1/x,x,a,b);
if b - a < 0 then integrate(1/x,x,a,b)
  elseif equal(b - a, 0) then integrate(1/x,x,a,b)
  elseif b - a > 0 then integrate(1/x,x,a,b);

limit((exp(a*x)-1)/(exp(a*x)+1),x,inf);
if a < 0 then -1 elseif equal (a, 0) then 0 elseif a > 0 then 1;

limit(exp(b*a), b, inf);
if a < 0 then 0 elseif equal (a, 0) then 1 elseif a > 0 then inf;

limit(exp(-b*a), b, inf);
if a < 0 then inf elseif equal (a, 0) then 1 elseif a > 0 then 0;

/* not sure if this one is handled correctly by expand_branches --
 * do we get into a cycle of asksign questions??
 *
(foo : %e^(beta*hbar*omega-beta*(hbar*(2*n+1)*omega+delta))
         /(1/(1-%e^-(beta*hbar*omega))^2 +(%e^-(beta*delta)-1)/(1-%e^-(2*beta*hbar*omega))),
 limit (foo, beta, inf));
0;

limit (ratsimp (foo), beta, inf);
0;
 *
 */

limit ((%e^-(b*d)-1)/(1-%e^-(2*b*h*z)), b, inf);
0;

/* a couple of examples inspired by mailing list 2015-03-04
 * tnx David Scherfgen
 */

integrate (x^n + x^m, x);
if equal(m,-1) then integrate(x^n+1/x,x) elseif notequal(m,-1) then
integrate(x^n+x^m,x);

(load (expand_branches), 0);
0;

expand_branches (integrate (x^n + x^m, x));
if equal(m,-1)
    then (if equal(n,-1) then 2*log(x) elseif notequal(n,-1)
              then log(x)+x^(n+1)/(n+1)) elseif notequal(m,-1)
    then (if equal(n,-1) then log(x)+x^(m+1)/(m+1) elseif notequal(n,-1)
              then x^(n+1)/(n+1)+x^(m+1)/(m+1));

/* expected failure due to assume weakness; see SF bug #2915 */

integrate (1/(a*x^2 + b), x);
if a*b < 0 then log((2*a*x-2*sqrt(-a*b))/(2*a*x+2*sqrt(-a*b)))/(2*sqrt(-a*b))
elseif a*b > 0 then atan(a*x/sqrt(a*b))/sqrt(a*b)$