--Copyright The Numerical Algorithms Group Limited 1996.

-- this is here strictly for documentation purposes, nothing is executed
--
--	Here are some problems that Maple and Mathematica cannot
--solve, but SymbMath can do.
--	The following examples came from news on the sci.math.symbolic
--newsgroup in 1991, and were run in Maple V, Mathematica 2.0, or
--SymbMath 2.1.
--
--************************ Example 1 ******************************
--	Maple:
--> int(exp(-a * x^2), x=0..infinity);
--
--                             infinity
--                                /
--                               |               2
--                               |      exp(- a x ) dx
--                               |
--                              /
--                              0
--
--# unevaluated.
--# Declare 'a' non-negative:
--
--> signum(a) := 1;
--
--                                 signum(a) := 1
--
--# The same integral is now evaluated fully:
--
--> int(exp(-a * x^2), x=0..infinity);
--
--                                         1/2
--                                       Pi
--                                   1/2 -----
--                                         1/2
--                                        a
--
--	SymbMath :
--	Input:
--inte(exp(-a*x^2), x from 0 to inf)
--assume(sqrt(a) > 0)
--inte(exp(-a*x^2), x from 0 to inf)
--	Output:
--1/2*a^(-0.5)*sqrt(pi)*erf(inf*sgn(sqrt(a)))
--assumed
--1/2*a^(-0.5)*sqrt(pi)
--
--************************** Example 2 *********************************
--	Maple:
--# Despite the fact that 'n' is declared non-negative...
--> signum(n) := 1;
--                                 signum(n) := 1
--
--# ...this simple proper definite integral (that any freshman calculus
--# student can evaluate!) is left unevaluated:
--> int(x^n, x=0..1);
--
--                                     1
--                                     /
--                                    |   n
--                                    |  x  dx
--                                    |
--                                   /
--                                   0
--
--
--	SymbMath :
--	Input:
--assume(n > -1)
--inte(x^n, x from 0 to 1)
--	Output:
--assumed
--1/(1 + n)
--
--************************** Example 3 *****************************
--	Maple:
--> int(x^n, x=eps..1);
--
--                                          (n + 1)
--                                 1     eps
--                               ----- - ----------
--                               n + 1      n + 1
--
--# ...but the one-sided limit...
--
--> limit(", eps=0, right);
--
--                                               (n + 1)
--                                      1     eps
--                          limit     ----- - ----------
--                          eps -> 0+ n + 1      n + 1
--
--# ...remains unevaluated...
--
--> simplify(");
--
--                                               (n + 1)
--                                      - 1 + eps
--                          limit     - ----------------
--                          eps -> 0+         n + 1
--
--
--# ...no matter what we do!
--
--> eval(");
--
--                                               (n + 1)
--                                      - 1 + eps
--                          limit     - ----------------
--                          eps -> 0+         n + 1
--
--	SymbMath:
--	Input:
--assume(n > -1)
--inte(x^n, x from eps to 1)
--subs(eps=0 to last)
--	Output:
--assumed
--1/(1 + n) - eps^(1 + n)/(1 + n)
--1/(1 + n)
--
--
--*************************** Example 4 ***************************
--	Maple:
--> 0^n;
--            0
--
--# Maple flags the error only when 'n' is replaced by the constant 0:
--
--> 0^0;
--Error, 0^0 is undefined
--
--	SymbMath:
--	Input:
--assume(n > 0)
--0^n
--0^-n
--0^0
--	Output:
--assumed
--0
--discont
--undefined
--
--**************************** Example 5 ******************************
--	Maple:
--> int(x^k, x);
--
--                                     (k + 1)
--                                    x
--                                    --------
--                                      k + 1
--
--	SymbMath:
--	Input:
--inte(x^k*d(x))
--subs(k=-1 to last)
--	Output:
--constant + x^(1 + k)/(1 + k)
--discont
--
--**************************** Example 6 *****************************
--	Maple:
--# The following limit is left unevaluated:
--> limit(x^k/exp(x), x=infinity);
--
--                                               k
--                                              x
--                              limit         ------
--                              x -> infinity exp(x)
--
--
--# We might ask if Maple knows this result for specific (constant) values
--# of the symbolic parameter 'k'.  The answer is a QUALIFIED "yes".
--# Maple knows the result for 'k' equal to 10^8...
--
--> limit(x^(10^8)/exp(x), x=infinity);
--
--                                       0
--
--
--# ...and Maple knows the result for 'k' equal to 10^9...
--
--> limit(x^(10^9)/exp(x), x=infinity);
--
--                                       0
--
--# ...BUT, Maple seems to FORGET the result when 'k' equals 10^10...
--
--> limit(x^(10^10)/exp(x), x=infinity);
--
--                                          10000000000
--                                         x
--                           limit         ------------
--                           x -> infinity    exp(x)
--
--	SymbMath:
--	Input:
--lim(x=inf, x^k/exp(x))
--lim(x=inf, x^(10^10)/exp(x))
--lim(x=inf, x^(10^10000)/exp(x))
--	Output:
--0
--0
--0
--
--****************************** Example 7 ****************************
--	Maple:
--> int(x^m * exp (-b * x), x=0..infinity);
--
--                            infinity
--                               /
--                              |       m
--                              |      x  exp(- b x) dx
--                              |
--                             /
--                             0
--
--
--# As expected, the integral remains unevaluated.  Now declare the signs
--# of the symbolic parameters:
--
--> signum(b) := 1;
--
--                                 signum(b) := 1
--
--> signum(m) := 1;
--
--                                 signum(m) := 1
--
--
--# Upon attempting to compute the integral a second time...
--
--> int(x^m * exp (-b * x), x=0..infinity);
--
--                            infinity
--                               /
--                              |       m
--                              |      x  exp(- b x) dx
--                              |
--                             /
--                             0
--
--# ...THE INTEGRAL CONTINUES TO REMAIN UNEVALUATED, despite the fact that
--# the signs of the parameters 'b' and 'm' were declared PRIOR to the second
--# attempt at computation.
--
--	SymbMath:
--	Input:
--inte(x^n*exp(-a*x), x from 0 to inf)
--	Output:
--inte(x^n*exp(-a*x), x, 0, inf)
--
--
--************************ Example 8 *********************************
--	Mathematica:
--In[1]:= Integrate[1/x,{x,-1,1}]
--Out[1]= -Log[-1]
--
--	Maple:
--has the same problem.
--
--	SymbMath:
--	Input:
--inte(1/x, x from -1 to 1)
--inte(1/x, x from -1 to 2)
--	Output:
--0
--ln(2)
--
--
--*************************** Example 9 *******************************
--	Maple:
--has a problem for int(tan(x), x=0..pi).
--
--	SymbMath:
--	Input:
--inte(tan(x), x from 0 to pi)
--	Output:
--0
--
--**************************** Example 10 ****************************
--	Mathematica:
--cannot evaluate integral of sgn(x).
--
--	Maple:
--# cannot evaluate inegral of signum(x) by int(). Help with a procedure:
--# Load the procedure into Maple:
--
--> read `pvint.txt`;
--
--pvint := proc(f,x,a,b,s)
--         local i1,i2,eps;
--             signum(eps) := 1;
--             i1 := int(f,x = a .. s-eps);
--             i2 := int(f,x = s+eps .. b);
--             simplify(i1+i2);
--             limit(",eps = 0,right)
--         end
--
--
--> pvint(signum(x), x, -1, 1, 0);
--
--                          - eps                   1
--                            /                     /
--                           |                     |
--               limit       |    signum(x) dx +   |  signum(x) dx
--               eps -> 0+   |                     |
--                          /                     /
--                          -1                   eps
--
--# Maple refuses to evaluate this P.V. integral:
--
--> simplify(");
--
--                          - eps                   1
--                            /                     /
--                           |                     |
--               limit       |    signum(x) dx +   |  signum(x) dx
--               eps -> 0+   |                     |
--                          /                     /
--                          -1                   eps
--
--> eval(");
--
--                          - eps                   1
--                            /                     /
--                           |                     |
--               limit       |    signum(x) dx +   |  signum(x) dx
--               eps -> 0+   |                     |
--                          /                     /
--                          -1                   eps
--
--	SymbMath:
--	Input:
--inte(sgn(x), x from -1 to 1)
--inte(sgn(x), x from -1 to 2)
--	Output:
--0
--1
--
--************************* Example 11 ********************************
--Implicit diff. gives 1+y'[x](1+1/y[x])==0; y'[x]==-y[x]/(y[x]+1).
--
--	Mathematica:
-- given this eq. as input to DSolve says that built-in
-- procedure can't solve it.
--
--	 MACSYMA:
--(c1) depends(y,x)$
--
--(c2) ode2(diff(y,x) = -y/(y+1),y,x);
--
--(d2)                   - log(y) - y = x + %c
--
--(c3) method;
--
--(d3)                         separable
--
--
--	SymbMath:
--solve the differential equation by integration inte() or by dsolve().
--	Input:
--d(y)/d(x)*(1+1/y) === -1
--inte(last*d(x))
--Expand=On
--dsolve(d(y)/d(x) === -y/(y+1), y)
--	Output:
--(1 + 1/y)*d(y)/d(x) === -1
--y + ln(y*sgn(y)) === constant - x
--Expand = On
---y - ln(y*sgn(y)) === constant + x
--
--
--************************* Example 12 ********************************
--	Mathematica:
--                 y'[x] = y[x]^(1/2)
--                  y[0] = 0
--
--DSolve could not handle it. (It rarely solves anything!!), and the
--RungeKutta package only gave me the solution
--
--                    y[x]=0
--
--Obviously, there is another solution viz.
--
--                   y[x] = (x/2)^2
--
--	SymbMath:
--	Input:
--dsolve(d(y)/d(x) === sqrt(y), y)
--(last/2)^2
--	Output:
--2*sqrt(y) === constant + x
--y === 1/4*(constant + x)^2
--
--************************* Example 13 ********************************
--	Maple:	
--> sqrt(x*x);
--      x
--
--	Mathematica:
-- Sqrt[a^2] evaluates to Sqrt[a^2].
--
--	SymbMath:
--	Input:
--sqrt(x^2)
--assume(a > 0)
--sqrt(a^2)
--assume(b <0 )
--sqrt(b^2)
--	Output:
--x*sgn(x)
--assumed
--a
--assumed
---b
--
--********************** Example 14 **********************************
--	Maple and Mathematica cannot find the integrals of abs(x).
--
--	SymbMath:
--	Input:
--inte(abs(x), x from -1 to 1)
--inte(abs(x)^5*d(x))
--	Output:
--1
--constant + 1/6*abs(x)^6*sgn(x)
--
-----------------------------------------------------------------------
--
--	The following problems are taken from Swokowski's Calculus
--book.  They cause Mathematica to fail because of a singularity in
--the interior of the interval of integration:
--
--        Section 10.4 Problems 3, 12, 15, 16, 23, 29.
--
--The comments "INTEGRAL IS DIVERGENT" and "Principal Value" come from
--Macsyma.  Mma gives no indication that anything is amiss.
--
--************************ Problem3 *********************************
--	Mathematica:
--In[4]:= Integrate[1/x^2,{x,-3,1}]
--
--          4
--Out[4]= -(-)                            (* INTEGRAL IS DIVERGENT *)
--          3
--
--	SymbMath:
--	Input:
--inte(1/x^2, x from -3 to 1)
--	Output:
--inf
--
--************************** Problem12 ***************************
--	Mathemtica:
--In[6]:= Integrate[x^(-4/3),{x,-1,1}]
--
--Out[6]= -6                              (* INTEGRAL IS DIVERGENT *)
--
--	SymbMath:
--	Input:
--inte(x^(-4/3), x from -1 to 1)
--        Output:
--inf
--
--**************************** Problem15  ************************
--	Mathematica:
--In[7]:= Integrate[1/x,{x,-1,2}]
--
--Out[7]= -Log[-1] + Log[2]
--
--	Maple:
--has the same problem.
--	
--	Macsyma:
--	(c7) integrate(1/x,x,-1,2);
--	Principal Value
--	(d7)                                log(2)
--
--	SymbMath:
--	Input:
--inte(1/x, x from -1 to 2)
--	Output:
--ln(2)
--
--
--********************** Problem16 ************************************
--	Mathematica:
--In[8]:= Integrate[1/(x^2-x-2),{x,0,4}]
--
--        -Log[-2]   Log[2]   Log[5]
--Out[8]= -------- + ------ - ------
--           3         3        3
--
--	Macsyma:
--	(c8) integrate(1/(x^2-x-2),x,0,4);
--	Principal Value
--	                                     log(5)
--	(d8)                               - ------
--	                                       3
--
--
--	SymbMath:
--	Input:
--inte(1/(x^2-x-2), x from 0 to 4)
--	Output:
---1/3*ln(2) + 1/3*ln(2/5)
--
--***************************** Problem23 **********************
--	Mathematica:
--In[10]:= Integrate[(1/x^2)Cos[1/x],{x,-1,2}]
--
--                       1
--Out[10]= Sin[-1] - Sin[-]              (* INTEGRAL IS DIVERGENT *)
--		       2
--
--	SymbMath:
--	Input:
--y=1/x^2*cos(1/x)
--inte(y, x from -1 to 0-zero) + inte(y, x from 0+zero to 2)
--	Output:
--y = x^(-2)*cos(1/x)
--sin(-1) - sin(1/2) + 2*sin(inf)
--
--************************** Problem29 ********************************
--	Mathematica:
--In[12]:= Integrate[1/(x-4)^2,{x,0,Infinity}]
--
--           1
--Out[12]= -(-)                          (* INTEGRAL IS DIVERGENT *)
--           4
--
--	SymbMath:
--	Input:
--y=1/(x-4)^2
--inte(y, x from 0 to 4-zero) + inte(y, x from 4+zero to inf)
--	Output:
--y = (-4 + x)^(-2)
--inf
