\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/input easter.input}
\author{The Axiom Team}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{License}
<<license>>=
--Copyright The Numerical Algorithms Group Limited 1994.
@
<<*>>=
<<license>>
)cl all
-- ----------[ A x i o m ]----------
-- ---------- Initialization ----------
)set bre resume
)set messages autoload off
)set messages time off
)set quit unprotected
)set streams calculate 7
-- ---------- Numbers ----------
-- Let's begin by playing with numbers: infinite precision integers
factorial(50)
factor(%)
-- Infinite precision rational numbers
1/2 + 1/3 + 1/4 + 1/5 + 1/6 + 1/7 + 1/8 + 1/9 + 1/10
-- Arbitrary precision floating point numbers
digits(50)
-- This number is nearly an integer
exp(sqrt(163.)*%pi)
digits(20)
-- Special functions
besselJ(2, 1 + %i)
-- Complete decimal expansion of a rational number
decimal(1/7)
-- Continued fractions
continuedFraction(3.1415926535)
-- Simplify an expression with nested square roots
sqrt(2*sqrt(3) + 4)
simplify(%)
-- Try a more complicated example (from the Putnam exam)
sqrt(14 + 3*sqrt(3 + 2*sqrt(5 - 12*sqrt(3 - 2*sqrt(2)))))
simplify(%)
-- Cardinal numbers
2*Aleph(0) - 3
-- ---------- Statistics ----------
-- ---------- Algebra ----------
-- Numbers are nice, but symbols allow for variability---try some high school
-- algebra: rational simplification
(x**2 - 4)/(x**2 + 4*x + 4)
-- This example requires more sophistication
(%e**x - 1)/(%e**(x/2) + 1)
normalize(%)
-- Expand and factor polynomials
(x + 1)**20
D(%, x)
factor(%)
x**100 - 1
factor(%)
-- Factor polynomials over finite fields and field extensions
p:= x**4 - 3*x**2 + 1
factor(p)
phi:= rootOf(phi**2 - phi - 1)
factor(p, [phi])
factor(p :: Polynomial(PrimeField(5)))
expand(%)
-- Partial fraction decomposition
(x**2 + 2*x + 3)/(x**3 + 4*x**2 + 5*x + 2)
padicFraction(
   partialFraction(numerator(%) :: UnivariatePolynomial(x, Fraction Integer),
                   factor(denominator(%) :: Polynomial Integer) ::
                      Factored UnivariatePolynomial(x, Fraction Integer)))
-- ---------- Inequalities ----------
-- ---------- Trigonometry ----------
-- Trigonometric manipulations---these are typically difficult for students
r:= cos(3*x)/cos(x)
real(complexNormalize(%))
real(normalize(simplify(complexNormalize(r))))
-- Use rewrite rules
sincosAngles:= rule _
  (cos((n | integer?(n)) * x) == _
      cos((n - 1)*x) * cos(x) - sin((n - 1)*x) * sin(x); _
   sin((n | integer?(n)) * x) == _
      sin((n - 1)*x) * cos(x) + cos((n - 1)*x) * sin(x) )
sincosAngles r
r:= 'r
-- ---------- Determining Zero Equivalence ----------
-- The following expressions are all equal to zero
sqrt(997) - (997**3)**(1/6)
sqrt(999983) - (999983**3)**(1/6)
(2**(1/3) + 4**(1/3))**3 - 6*(2**(1/3) + 4**(1/3)) - 6
simplify(%)
-- This expression is zero for x, y > 0 and n not equal to zero
x**(1/n)*y**(1/n) - (x*y)**(1/n)
normalize(%)
-- See Joel Moses, ``Algebraic Simplification: A Guide for the Perplexed'',
-- CACM, Volume 14, Number 8, August 1971
expr:= log(tan(1/2*x + %pi/4)) - asinh(tan(x))
complexNormalize(%)
-- Use a roundabout method---show that expr is a constant equal to zero
D(expr, x)
simplify(real(complexNormalize(expand(simplify(%)))))
normalize(eval(expr, x = 0))
log((2*sqrt(r) + 1)/sqrt(4*r + 4*sqrt(r) + 1))
simplify(%)
(4*r + 4*sqrt(r) + 1)**(sqrt(r)/(2*sqrt(r) + 1)) _
   * (2*sqrt(r) + 1)**(1/(2*sqrt(r) + 1)) - 2*sqrt(r) - 1
normalize(%)
-- ---------- The Complex Domain ----------
-- Complex functions---separate into their real and imaginary parts
rectform(z) == real(z) + %i*imag(z)
rectform(log(3 + 4*%i))
simplify(rectform(tan(x + %i*y)))
-- Check for branch abuse.  See David R. Stoutemyer, ``Crimes and Misdemeanors
-- in the Computer Algebra Trade'', Notices of the AMS, Volume 38, Number 7,
-- September 1991.  This first expression can simplify to sqrt(x y)/sqrt(x),
-- but no further in general (consider what happens when x, y = -1).
sqrt(x*y*abs(z)**2) / (sqrt(x)*abs(z))
-- If z = -1, sqrt(1/z) is not equal to 1/sqrt(z)
sqrt(1/z) - 1/sqrt(z)
-- If z = 3 pi i, log(exp(z)) is not equal to z
log(%e**z)
normalize(%)
-- The principal value of this expression is (10 - 4 pi) i
log(%e**(10*%i))
normalize(%)
-- If z = pi, arctan(tan(z)) is not equal to z
atan(tan(z))
-- If z = 2 pi i, sqrt(exp(z)) is not equal to exp(z/2)
sqrt(%e**z) - %e**(z/2)
-- ---------- Equations ----------
-- Manipulate an equation using a natural syntax
(x = 0)/2 + 1
-- Solve various nonlinear equations---this cubic polynomial has all real roots
radicalSolve(3*x**3 - 18*x**2 + 33*x - 19 = 0, x)
map(e +-> lhs(e) = rectform(rhs(e)), %)
-- Some simple seeming problems can have messy answers
eqn:= x**4 + x**3 + x**2 + x + 1 = 0
radicalSolve(eqn, x)
-- Check one of the answers
eval(eqn, %.1)
%e**(2*x) + 2*%e**x + 1 = z
solve(%, x)
-- This equation is already factored and so *should* be easy to solve
(x + 1) * (sin(x)**2 + 1)**2 * cos(3*x)**3 = 0
solve(%, x)
-- The following equations have an infinite number of solutions (let n be an
-- arbitrary integer): z = 0 [+ n 2 pi i]
solve(%e**z = 1, z)
-- x = pi/4 [+ n pi]
solve(sin(x) = cos(x), x)
solve(tan(x) = 1, x)
-- x = 0, 0 [+ n pi, + n 2 pi]
solve(sin(x) = tan(x), x)
-- This equation has no solutions
solve(sqrt(x**2 + 1) = x - 2, x)
-- Solve a system of linear equations
eq1:=   x +   y +   z =  6
eq2:= 2*x +   y + 2*z = 10
eq3:=   x + 3*y +   z = 10
-- Note that the solution is parametric
solve([eq1, eq2, eq3], [x, y, z])
-- Solve a system of nonlinear equations
eq1:= x**2*y + 3*y*z - 4 = 0
eq2:= -3*x**2*z + 2*y**2 + 1 = 0
eq3:= 2*y*z**2 - z**2 - 1 = 0
-- Solving this by hand would be a nightmare
solve([eq1, eq2, eq3], [x, y, z])
-- ---------- Matrix Algebra ----------
m:= matrix([[a, b], [1, a*b]])
-- Invert the matrix
minv:= inverse(m)
m * minv
-- Define a Vandermonde matrix (useful for doing polynomial interpolations)
matrix([[1,    1,    1,    1   ], _
        [w,    x,    y,    z   ], _
        [w**2, x**2, y**2, z**2], _
        [w**3, x**3, y**3, z**3]])
determinant(%)
-- The following formula implies a general result
factor(%)
-- Compute the eigenvalues of a matrix from its characteristic polynomial
m:= matrix([[ 5, -3, -7], _
            [-2,  1,  2], _
            [ 2, -3, -4]])
characteristicPolynomial(m, lambda)
solve(% = 0, lambda)
m:= 'm
-- ---------- Tensors ----------
-- ---------- Sums and Products ----------
-- Sums: finite and infinite
summation(k**3, k = 1..n)
sum(k**3, k = 1..n)
limit(sum(1/k**2 + 1/k**3, k = 1..n), n = %plusInfinity)
-- Products
product(k, k = 1..n)
-- ---------- Calculus ----------
-- Limits---start with a famous example
limit((1 + 1/n)**n, n = %plusInfinity)
limit((1 - cos(x))/x**2, x = 0)
-- Apply the chain rule---this is important for PDEs and many other
-- applications
y:= operator('y)
x:= operator('x)
D(y(x(t)), t, 2)
)clear properties x y
-- ---------- Indefinite Integrals ----------
1/(x**3 + 2)
-- This would be very difficult to do by hand
integrate(%, x)
D(%, x)
-- This example involves several symbolic parameters
integrate(1/(a + b*cos(x)), x)
map(simplify, map(f +-> D(f, x), %))
-- Calculus on a non-smooth (but well defined) function
D(abs(x), x)
integrate(abs(x), x)
-- Calculus on a piecewise defined function
a(x) == if x < 0 then -x else x
D(a(x), x)
integrate(a(x), x)
)clear properties a
-- The following two integrals should be equivalent.  The correct solution is
-- [(1 + x)^(3/2) + (1 - x)^(3/2)] / 3
integrate(x/(sqrt(1 + x) + sqrt(1 - x)), x)
integrate((sqrt(1 + x) - sqrt(1 - x))/2, x)
-- ---------- Definite Integrals ----------
-- The following two functions have a pole at zero
integrate(1/x, x = -1..1)
integrate(1/x**2, x = -1..1)
-- Different branches of the square root need to be chosen in the intervals
-- [0, 1] and [1, 2].  The correct results are 4/3, [4 - sqrt(8)]/3,
-- [8 - sqrt(8)]/3, respectively.
integrate(sqrt(x + 1/x - 2), x = 0..1)
integrate(sqrt(x + 1/x - 2), x = 0..1, "noPole")
integrate(sqrt(x + 1/x - 2), x = 1..2)
integrate(sqrt(x + 1/x - 2), x = 1..2, "noPole")
integrate(sqrt(x + 1/x - 2), x = 0..2)
integrate(sqrt(x + 1/x - 2), x = 0..2, "noPole")
-- Contour integrals
integrate(cos(x)/(x**2 + a**2), x = %minusInfinity..%plusInfinity)
integrate(cos(x)/(x**2 + a**2), x = %minusInfinity..%plusInfinity, "noPole")
-- Integrand with a branch point
integrate(t**(a - 1)/(1 + t), t = 0..%plusInfinity)
integrate(t**(a - 1)/(1 + t), t = 0..%plusInfinity, "noPole")
-- Multiple integrals: volume of a tetrahedron
integrate(integrate(integrate(1, z = 0..c*(1 - x/a - y/b)), _
                    y = 0..b*(1 - x/a)), _
          x = 0..a)
-- ---------- Series ----------
-- Taylor series---this first example comes from special relativity
1/sqrt(1 - (v/c)**2)
series(%, v = 0)
1/%**2
tsin:= series(sin(x), x = 0)
tcos:= series(cos(x), x = 0)
-- Note that additional terms will be computed as needed
tsin/tcos
series(tan(x), x = 0)
-- Look at the Taylor series around x = 1
)set streams calculate 1
log(x)**a*exp(-b*x)
series(%, x = 1)
)set streams calculate 7
-- Compare the Taylor series of two different formulations of a function
taylor(log(sinh(z)) + log(cosh(z + w)), z = 0)
% - taylor(log(sinh(z) * cosh(z + w)), z = 0)
-- Power series (compute the general formula)
log(sin(x)/x)
series(%, x = 0)
exp(-x)*sin(x)
series(%, x = 0)
-- Derive an explicit Taylor series solution of y as a function of x from the
-- following implicit relation
y:= operator('y)
x = sin(y(x)) + cos(y(x))
seriesSolve(%, y, x = 1, 0)
)clear properties y
-- Pade (rational function) approximation
pade(1, 1, taylor(exp(-x), x = 0))
-- ---------- Transforms ----------
-- Laplace and inverse Laplace transforms
laplace(cos((w - 1)*t), t, s)
inverseLaplace(%, s, t)
-- ---------- Difference and Differential Equations ----------
-- Second order linear recurrence equation
r:= operator('r)
r(n + 2) - 2 * r(n + 1) + r(n) = 2
[%, r(0) = 1, r(1) = m]
)clear properties r
-- Second order ODE with initial conditions---solve first using Laplace
-- transforms
f:= operator('f)
ode:= D(f(t), t, 2) + 4*f(t) = sin(2*t)
map(e +-> laplace(e, t, s), %)
-- Now, solve the ODE directly
solve(ode, f, t = 0, [0, 0])
-- First order linear ODE
y:= operator('y)
x**2 * D(y(x), x) + 3*x*y(x) = sin(x)/x
solve(%, y, x)
-- Nonlinear ODE
D(y(x), x, 2) + y(x)*D(y(x), x)**3 = 0
solve(%, y, x)
-- A simple parametric ODE
D(y(x, a), x) = a*y(x, a)
solve(%, y, x)
-- ODE with boundary conditions.  This problem has nontrivial solutions
-- y(x) = A sin([pi/2 + n pi] x) for n an arbitrary integer.
solve(D(y(x), x, 2) + k**2*y(x) = 0, y, x)
-- bc(%, x = 0, y = 0, x = 1, D(y(x), x) = 0)
-- System of two linear, constant coefficient ODEs
x:= operator('x)
system:= [D(x(t), t) = x(t) - y(t), D(y(t), t) = x(t) + y(t)]
-- Check the answer
-- Triangular system of two ODEs
system:= [D(x(t), t) = x(t) * (1 + cos(t)/(2 + sin(t))), _
          D(y(t), t) = x(t) - y(t)]
-- Try solving this system one equation at a time
solve(system.1, x, t)
isTimes(subst(%.basis.1, cos(t) = sqrt(1 - sin(t)**2)))
reduce(*, cons(subst(
   factors(factor(subst(%.1**2, sin(t) = u) :: Polynomial Integer)).1.factor,
                     u = sin(t)),
               rest(%)))
x(t) = C1 * %
solve(map(e +-> subst(e, %), system.2), y, t)
y(t) = simplify(%.particular) + C2 * %.basis.1
)clear properties x y
-- ---------- Operators ----------
-- Linear differential operator
DD:= operator("D") :: Operator(Expression Integer)
evaluate(DD, e +-> D(e, x))$Operator(Expression Integer)
L:= (DD - 1) * (DD + 2)
g:= operator('g)
L(f(x))
subst(L(subst(g(y), y = x)), x = y)
subst(L(subst(A * sin(z**2), z = x)), x = z)
-- Truncated Taylor series operator
T:= (f, xx, a) +-> subst((DD**0)(f(x)), x = a)/factorial(0) * (xx - a)**0 + _
                   subst((DD**1)(f(x)), x = a)/factorial(1) * (xx - a)**1 + _
                   subst((DD**2)(f(x)), x = a)/factorial(2) * (xx - a)**2
T(f, x, a)
T(g, y, b)
Sin:= operator("sin") :: Operator(Expression Integer)
evaluate(Sin, x +-> sin(x))$Operator(Expression Integer)
T(Sin, z, c)
-- ---------- Programming ----------
-- Write a simple program to compute Legendre polynomials
p(n, x) == 1/(2**n*factorial(n)) * D((x**2 - 1)**n, x, n)
for i in 0..4 repeat {  output("");    output(concat(["p(", string(i), ", x) = "]));    output(p(i, x))}
eval(p(4, x), x = 1)
-- Now, perform the same computation using a recursive definition
pp(0, x) == 1
pp(1, x) == x
pp(n, x) == ((2*n - 1)*x*pp(n - 1, x) - (n - 1)*pp(n - 2, x))/n
for i in 0..4 repeat {   output("");   output(concat(["pp(", string(i), ", x) = "]));   output(pp(i, x))}
)clear properties p pp
-- ---------- Translation ----------
-- Horner's rule---this is important for numerical algorithms
a:= operator('a)
sum(a(i)*x**i, i = 1..5)
p:= factor(%)
-- Convert the result into FORTRAN syntax
)set fortran ints2floats off
outputAsFortran('p = p)
-- ---------- Boolean Logic ----------
-- Simplify logical expressions
true and false
x or (not x)
x or y or (x and y)
)set bre nobr
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}