|
|
1
|
|
Editor:
Time: 2007/11/18 18:34:41 GMT-8
|
|
Note: update
|
changed:
-
\documentclass{article}
\def\code#1{{\tt #1}}
\def\url#1{{\it #1}}
\def\sage#1{{\bf #1}}
\def\note#1{#1}
\usepackage{axiom}
\begin{document}
\title{A Sage Interface to Axiom}
\author{
AUTHORS OF THIS MODULE:
- William Stein (2006-10): Initial version based on Maxima interface
- Bill Page (2006-10): Axiom version
}
\abstract{
Interface to Axiom
Axiom is a free GPL-compatible (modified BSD license) general purpose
computer algebra system whose development started in 1973 at IBM. It
contains symbolic manipulation algorithms, as well as implementations
of special functions, including elliptic functions and generalized
hypergeometric functions. Moreover, Axiom has implementations of many
functions relating to the invariant theory of the symmetric group $S_n$.
For many links to Axiom documentation see
\url{http://wiki.axiom-developer.org}.
EXAMPLES:
We evaluate a very simple expression in axiom.
\begin{verbatim}
sage: axiom('3 * 5')
15
Type: PositiveInteger
\end{verbatim}
We factor $x^5 - y^5$ in Axiom in several different ways.
The first way yields a Axiom object.
\begin{verbatim}
sage: F = axiom.factor('x^5 - y^5')
sage: type(F)
<class 'axiom.interfaces.axiom.AxiomElement'>
\end{verbatim}
Note that Axiom objects are normally displayed using ``ASCII art'';
to see a normal linear representation of any Axiom object x,
use \code{str(x)}. no emph
\begin{verbatim}
sage: F
4 3 2 2 3 4
- (y - x) (y + x y + x y + x y + x )
\end{verbatim}
You can always use \code{x.str()} to obtain the linear representation
of an object, even without changing the display2d flag. This can
be useful for moving axiom data to other systems.
\begin{verbatim}
sage: F.str()
'-(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4)'
sage: axiom.display2d(False)
sage: F
-(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4)
\end{verbatim}
The \code{axiom.eval} command evaluates an expression in axiom
and returns the result as a string.
\begin{verbatim}
sage: print axiom.eval('factor(x^5 - y^5)')
-(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4)
\end{varbatim}
We can create the polynomial $f$ as a Axiom polynomial, then call
the factor method on it. Notice that the notation \code{f.factor()}
is consistent with how the rest of \sage works.
\begin{verbatim}
sage: f = axiom('x^5 - y^5')
sage: f^2
(x^5 - y^5)^2
sage: f.factor()
-(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4)
\end{verbatim}
Control-C interruption works well with the axiom interface,
because of the excellent implementation of axiom. For example,
try the following sum but with a much bigger range, and hit
control-C.
\begin{verbatim}
sage: axiom('sum(1/x^2, x, 1, 10)')
1968329/1270080
\end{verbatim}
\subsection{Tutorial}
We follow the tutorial at
\url{http://wiki.axiom-developer.org/AxiomTutorial}.
\begin{verbatim}
sage: axiom('1/100 + 1/101')
201/10100
sage: a = axiom('(1 + sqrt(2))^5'); a
(sqrt(2) + 1)^5
sage: a.expand()
29*sqrt(2) + 41
sage: a = axiom('(1 + sqrt(2))^5')
sage: float(a)
82.012193308819747
sage: a.numer()
82.01219330881975
sage: axiom.eval('fpprec : 100')
'100'
sage: a.bfloat()
8.20121933088197564152489730020812442785204843859314941221237124017312418754011041266612384955016056b1
sage: axiom('100!')
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
sage: f = axiom('(x + 3*y + x^2*y)^3')
sage: f.expand()
x^6*y^3 + 9*x^4*y^3 + 27*x^2*y^3 + 27*y^3 + 3*x^5*y^2 + 18*x^3*y^2 + 27*x*y^2 + 3*x^4*y + 9*x^2*y + x^3
sage: f.subst('x=5/z')
(5/z + 25*y/z^2 + 3*y)^3
sage: g = f.subst('x=5/z')
sage: h = g.ratsimp(); h
(27*y^3*z^6 + 135*y^2*z^5 + (675*y^3 + 225*y)*z^4 + (2250*y^2 + 125)*z^3 + (5625*y^3 + 1875*y)*z^2
+ 9375*y^2*z + 15625*y^3)/z^6
sage: h.factor()
(3*y*z^2 + 5*z + 25*y)^3/z^6
sage: eqn = axiom(['a+b*c=1', 'b-a*c=0', 'a+b=5'])
sage: s = eqn.solve('[a,b,c]'); s
[[a = (25*sqrt(79)*%i + 25)/(6*sqrt(79)*%i - 34),b = (5*sqrt(79)*%i + 5)/(sqrt(79)*%i + 11),
c = (sqrt(79)*%i + 1)/10],
[a = (25*sqrt(79)*%i - 25)/(6*sqrt(79)*%i + 34),b = (5*sqrt(79)*%i - 5)/(sqrt(79)*%i - 11),
c = - (sqrt(79)*%i - 1)/10]]
\end{varbatim}
Here is an example of solving an algebraic equation:
\begin{verbatim}
sage: axiom('x^2+y^2=1').solve('y')
[y = - sqrt(1 - x^2),y = sqrt(1 - x^2)]
sage: axiom('x^2 + y^2 = (x^2 - y^2)/sqrt(x^2 + y^2)').solve('y')
[y = - sqrt(( - y^2 - x^2)*sqrt(y^2 + x^2) + x^2),y = sqrt(( - y^2 - x^2)*sqrt(y^2 + x^2) + x^2)]
\end{verbatim}
You can even nicely typeset the solution in latex:
\begin{verbatim}
sage: latex(s)
\left[ \left[ a=\frac{25 \sqrt{79} i+25}{6 \sqrt{79} i-34} , b= \frac{5 \sqrt{79} i+5}{\sqrt{79} i+11} , c=\frac{\sqrt{79} i+1}{10} \right] , \left[ a=\frac{25 \sqrt{79} i-25}{6 \sqrt{79} i+34} , b= \frac{5 \sqrt{79} i-5}{\sqrt{79} i-11} , c=-\frac{\sqrt{79} i-1}{10} \right] \right]
\end{verbatim}
To have the above appear onscreen via \code{xdvi}, type \code{view(s)}.
(TODO: For OS X should create pdf output and use preview instead?)
\begin{verbatim}
sage: e = axiom('sin(u + v) * cos(u)^3'); e
cos(u)^3*sin(v + u)
sage: f = e.trigexpand(); f
cos(u)^3*(cos(u)*sin(v) + sin(u)*cos(v))
sage: f.trigreduce()
(sin(v + 4*u) + sin(v - 2*u))/8 + (3*sin(v + 2*u) + 3*sin(v))/8
sage: w = axiom('3 + k*%i')
sage: f = w^2 + axiom('%e')^w
sage: f.realpart()
%e^3*cos(k) - k^2 + 9
sage: f = axiom('x^3 * %e^(k*x) * sin(w*x)'); f
x^3*%e^(k*x)*sin(w*x)
sage: f.diff('x')
k*x^3*%e^(k*x)*sin(w*x) + 3*x^2*%e^(k*x)*sin(w*x) + w*x^3*%e^(k*x)*cos(w*x)
sage: f.integrate('x')
(((k*w^6 + 3*k^3*w^4 + 3*k^5*w^2 + k^7)*x^3 + (3*w^6 + 3*k^2*w^4 - 3*k^4*w^2 - 3*k^6)*x^2 + ( - 18*k*w^4 - 12*k^3*w^2 + 6*k^5)*x - 6*w^4 + 36*k^2*w^2 - 6*k^4)*%e^(k*x)*sin(w*x) + (( - w^7 - 3*k^2*w^5 - 3*k^4*w^3 - k^6*w)*x^3 + (6*k*w^5 + 12*k^3*w^3 + 6*k^5*w)*x^2 + (6*w^5 - 12*k^2*w^3 - 18*k^4*w)*x - 24*k*w^3 + 24*k^3*w)*%e^(k*x)*cos(w*x))/(w^8 + 4*k^2*w^6 + 6*k^4*w^4 + 4*k^6*w^2 + k^8)
sage: f = axiom('1/x^2')
sage: f.integrate('x', 1, 'inf')
1
sage: g = axiom('f/sinh(k*x)^4')
sage: g.taylor('x', 0, 3)
f/(k^4*x^4) - 2*f/(3*k^2*x^2) + 11*f/45 - 62*k^2*f*x^2/945
sage: axiom.taylor('asin(x)','x',0, 10)
x + x^3/6 + 3*x^5/40 + 5*x^7/112 + 35*x^9/1152
\end{verbatim}
\subsection{Examples involving matrices}
We illustrate computing with the matrix whose $i,j$ entry
is $i/j$, for $i,j=1,\ldots,4$.
\begin{verbatim}
sage: f = axiom.eval('f[i,j] := i/j')
sage: A = axiom('genmatrix(f,4,4)'); A
matrix([1,1/2,1/3,1/4],[2,1,2/3,1/2],[3,3/2,1,3/4],[4,2,4/3,1])
sage: A.determinant()
0
sage: A.echelon()
matrix([1,1/2,1/3,1/4],[0,0,0,0],[0,0,0,0],[0,0,0,0])
sage: A.eigenvalues()
[[0,4],[3,1]]
sage: A.eigenvectors()
[[[0,4],[3,1]],[1,0,0, - 4],[0,1,0, - 2],[0,0,1, - 4/3],[1,2,3,4]]
\end{verbatim}
We can also compute the echelon form in \sage:
\begin{verbatim}
sage: B = matrix(QQ, A)
sage: B.echelon_form()
[ 1 1/2 1/3 1/4]
[ 0 0 0 0]
[ 0 0 0 0]
[ 0 0 0 0]
sage: B.charpoly().factor()
(x - 4) * x^3
\end{verbatim}
\subsection{Laplace Transforms}
We illustrate Laplace transforms:
\begin{verbatim}
sage: _ = axiom.eval("f(t) := t*sin(t)")
sage: axiom("laplace(f(t),t,s)")
2*s/(s^2 + 1)^2
sage: axiom("laplace(delta(t-3),t,s)") #Dirac delta function
%e^-(3*s)
sage: _ = axiom.eval("f(t) == exp(t)*sin(t)")
sage: axiom("laplace(f(t),t,s)")
1/(s^2 - 2*s + 2)
sage: _ = axiom.eval("f(t) := t^5*exp(t)*sin(t)")
sage: axiom("laplace(f(t),t,s)")
360*(2*s - 2)/(s^2 - 2*s + 2)^4 - 480*(2*s - 2)^3/(s^2 - 2*s + 2)^5 + 120*(2*s - 2)^5/(s^2 - 2*s + 2)^6
sage: axiom("laplace(f(t),t,s)")
3 5
360 (2 s - 2) 480 (2 s - 2) 120 (2 s - 2)
--------------- - --------------- + ---------------
2 4 2 5 2 6
(s - 2 s + 2) (s - 2 s + 2) (s - 2 s + 2)
sage: axiom("laplace(diff(x(t),t),t,s)")
s*laplace(x(t),t,s) - x(0)
sage: axiom("laplace(diff(x(t),t,2),t,s)")
-at('diff(x(t),t,1),t = 0) + s^2*laplace(x(t),t,s) - x(0)*s
sage.: axiom("laplace(diff(x(t),t,2),t,s)")
!
d ! 2
- -- (x(t))! + s laplace(x(t), t, s) - x(0) s
dt !
!t = 0
\end{verbatim}
Even better, use \code{view(axiom("laplace(diff(x(t),t,2),t,s)"))} to see
a typeset version.
\subsection{Continued Fractions}
A continued fraction $a + 1/(b + 1/(c + \cdots))$ is
represented in axiom by the list $[a, b, c, \ldots]$.
\begin{verbatim}
sage: axiom("cf((1 + sqrt(5))/2)")
[1,1,1,1,2]
sage: axiom("cf ((1 + sqrt(341))/2)")
[9,1,2,1,2,1,17,1,2,1,2,1,17,1,2,1,2,1,17,2]
\end{verbatim}
\subsection{Special examples}
In this section we illustrate calculations that would be awkward to do
(as far as I know) in non-symbolic computer algebra systems like MAGMA
or GAP.
We compute the gcd of $2x^{n+4} - x^{n+2}$ and $4x^{n+1} + 3x^n$
for arbitrary $n$.
\begin{verbatim}
sage: f = axiom('2*x^(n+4) - x^(n+2)')
sage: g = axiom('4*x^(n+1) + 3*x^n')
sage: f.gcd(g)
x^n
\end{verbatim}
You can plot 3d graphs (via gnuplot):
\begin{verbatim}
sage.: axiom('plot3d(x^2-y^2, [x,-2,2], [y,-2,2], [grid,12,12])')
[displays a 3 dimensional graph]
\end{verbatim}
You can formally evaluate sums (note the \code{nusum} command):
\begin{verbatim}
sage: S = axiom('nusum(exp(1+2*i/n),i,1,n)')
sage.: S
2/n + 3 2/n + 1
%e %e
----------------------- - -----------------------
1/n 1/n 1/n 1/n
(%e - 1) (%e + 1) (%e - 1) (%e + 1)
\end{verbatim}
We formally compute the limit as $n\to\infty$ of $2S/n$ as follows:
\begin{verbatim}
sage: T = S*axiom('2/n')
sage: T.tlimit('n','inf')
%e^3 - %e
\end{verbatim}
\subsection{Miscellaneous}
Obtaining digits of $\pi$:
\begin{verbatim}
sage: axiom.eval('fpprec : 100')
'100'
sage: axiom(pi).bfloat()
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068b0
\end{verbatim}
Defining functions in axiom:
\begin{verbatim}
sage: axiom.eval('fun(a) == a^2')
'fun[a] := a^2'
sage: axiom('fun(10)')
100
\end{verbatim}
\subsection{Interactivity}
Axiom has a non-interactive mode that is initiated via the command
")read file.input".
\subsection{Latex Output}
The latex output of Axiom is not perfect. E.g.,
\begin{verbatim}
sage: axiom.eval('tex(sin(u) + sinh(v^2))')
'$$\\sinhv^2 + \\sinu$$false'
\end{verbatim}
Notice the lack of space after the sin macro, which is a latex syntax
error. In \sage this is automatically fixed via a substition for
trig functions, which may have potentially bad side effects:
\begin{verbatim}
sage: latex(axiom('sin(u) + sinh(v^2)'))
\sinh v^2+\sin u
\end{verbatim}
It would be nice if somebody would fix this problem. One way would
be to improve Axiom by making the fix to Axiom and giving this back
to the Axiom people.
Here's another example:
\begin{verbatim}
sage: g = axiom('exp(3*%i*x)/(6*%i) + exp(%i*x)/(2*%i) + c')
sage: latex(g)
-\frac{i e^{3 i x}}{6}-\frac{i e^{i x}}{2}+c
\end{verbatim}
\subsection{Long Input}
The Axiom interface reads in even very long input (using files) in a
robust manner, as long as you are creating a new object.
\note{Using \code{axiom.eval} for long input
is much less robust, and is not recommended.}
\begin{verbatim}
sage: t = '"%s"'%10^10000 # ten thousand character string.
sage: a = axiom(t)
\end{verbatim}
}
\maketitle
<<axiom.py>>=
r"""
Interface to Axiom
Axiom is a free GPL-compatible (modified BSD license) general purpose
computer algebra system whose development started in 1973 at IBM. It
contains symbolic manipulation algorithms, as well as implementations
of special functions, including elliptic functions and generalized
hypergeometric functions. Moreover, Axiom has implementations of many
functions relating to the invariant theory of the symmetric group $S_n$.
For many links to Axiom documentation see
\url{http://wiki.axiom-developer.org}.
AUTHORS OF THIS MODULE:
- William Stein (2006-10): Initial version based on Maxima interface
- Bill Page (2006-10): Axiom version
If the string "error" (case insensitive) occurs in the output of
anything from axiom, a RuntimeError exception is raised.
EXAMPLES:
We evaluate a very simple expression in axiom.
sage: axiom('3 * 5')
15
Type: PositiveInteger
We factor $x^5 - y^5$ in Axiom in several different ways.
The first way yields a Axiom object.
sage: F = axiom.factor('x^5 - y^5')
sage: type(F)
<class 'axiom.interfaces.axiom.AxiomElement'>
Note that Axiom objects are normally displayed using ``ASCII art'';
to see a normal linear representation of any Axiom object x,
use \code{str(x)}.
sage: F
4 3 2 2 3 4
- (y - x) (y + x y + x y + x y + x )
You can always use \code{x.str()} to obtain the linear representation
of an object, even without changing the display2d flag. This can
be useful for moving axiom data to other systems.
sage: F.str()
'-(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4)'
sage: axiom.display2d(False)
sage: F
-(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4)
The \code{axiom.eval} command evaluates an expression in axiom
and returns the result as a string.
sage: print axiom.eval('factor(x^5 - y^5)')
-(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4)
We can create the polynomial $f$ as a Axiom polynomial, then call
the factor method on it. Notice that the notation \code{f.factor()}
is consistent with how the rest of \sage works.
sage: f = axiom('x^5 - y^5')
sage: f^2
(x^5 - y^5)^2
sage: f.factor()
-(y - x)*(y^4 + x*y^3 + x^2*y^2 + x^3*y + x^4)
Control-C interruption works well with the axiom interface,
because of the excellent implementation of axiom. For example,
try the following sum but with a much bigger range, and hit
control-C.
sage: axiom('sum(1/x^2, x, 1, 10)')
1968329/1270080
\subsection{Tutorial}
We follow the tutorial at
\url{http://wiki.axiom-developer.org/AxiomTutorial}.
sage: axiom('1/100 + 1/101')
201/10100
sage: a = axiom('(1 + sqrt(2))^5'); a
(sqrt(2) + 1)^5
sage: a.expand()
29*sqrt(2) + 41
sage: a = axiom('(1 + sqrt(2))^5')
sage: float(a)
82.012193308819747
sage: a.numer()
82.01219330881975
sage: axiom.eval('fpprec : 100')
'100'
sage: a.bfloat()
8.20121933088197564152489730020812442785204843859314941221237124017312418754011041266612384955016056b1
sage: axiom('100!')
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
sage: f = axiom('(x + 3*y + x^2*y)^3')
sage: f.expand()
x^6*y^3 + 9*x^4*y^3 + 27*x^2*y^3 + 27*y^3 + 3*x^5*y^2 + 18*x^3*y^2 + 27*x*y^2 + 3*x^4*y + 9*x^2*y + x^3
sage: f.subst('x=5/z')
(5/z + 25*y/z^2 + 3*y)^3
sage: g = f.subst('x=5/z')
sage: h = g.ratsimp(); h
(27*y^3*z^6 + 135*y^2*z^5 + (675*y^3 + 225*y)*z^4 + (2250*y^2 + 125)*z^3 + (5625*y^3 + 1875*y)*z^2 + 9375*y^2*z + 15625*y^3)/z^6
sage: h.factor()
(3*y*z^2 + 5*z + 25*y)^3/z^6
sage: eqn = axiom(['a+b*c=1', 'b-a*c=0', 'a+b=5'])
sage: s = eqn.solve('[a,b,c]'); s
[[a = (25*sqrt(79)*%i + 25)/(6*sqrt(79)*%i - 34),b = (5*sqrt(79)*%i + 5)/(sqrt(79)*%i + 11),c = (sqrt(79)*%i + 1)/10],[a = (25*sqrt(79)*%i - 25)/(6*sqrt(79)*%i + 34),b = (5*sqrt(79)*%i - 5)/(sqrt(79)*%i - 11),c = - (sqrt(79)*%i - 1)/10]]
Here is an example of solving an algebraic equation:
sage: axiom('x^2+y^2=1').solve('y')
[y = - sqrt(1 - x^2),y = sqrt(1 - x^2)]
sage: axiom('x^2 + y^2 = (x^2 - y^2)/sqrt(x^2 + y^2)').solve('y')
[y = - sqrt(( - y^2 - x^2)*sqrt(y^2 + x^2) + x^2),y = sqrt(( - y^2 - x^2)*sqrt(y^2 + x^2) + x^2)]
You can even nicely typeset the solution in latex:
sage: latex(s)
\left[ \left[ a=\frac{25 \sqrt{79} i+25}{6 \sqrt{79} i-34} , b= \frac{5 \sqrt{79} i+5}{\sqrt{79} i+11} , c=\frac{\sqrt{79} i+1}{10} \right] , \left[ a=\frac{25 \sqrt{79} i-25}{6 \sqrt{79} i+34} , b= \frac{5 \sqrt{79} i-5}{\sqrt{79} i-11} , c=-\frac{\sqrt{79} i-1}{10} \right] \right]
To have the above appear onscreen via \code{xdvi}, type \code{view(s)}.
(TODO: For OS X should create pdf output and use preview instead?)
sage: e = axiom('sin(u + v) * cos(u)^3'); e
cos(u)^3*sin(v + u)
sage: f = e.trigexpand(); f
cos(u)^3*(cos(u)*sin(v) + sin(u)*cos(v))
sage: f.trigreduce()
(sin(v + 4*u) + sin(v - 2*u))/8 + (3*sin(v + 2*u) + 3*sin(v))/8
sage: w = axiom('3 + k*%i')
sage: f = w^2 + axiom('%e')^w
sage: f.realpart()
%e^3*cos(k) - k^2 + 9
sage: f = axiom('x^3 * %e^(k*x) * sin(w*x)'); f
x^3*%e^(k*x)*sin(w*x)
sage: f.diff('x')
k*x^3*%e^(k*x)*sin(w*x) + 3*x^2*%e^(k*x)*sin(w*x) + w*x^3*%e^(k*x)*cos(w*x)
sage: f.integrate('x')
(((k*w^6 + 3*k^3*w^4 + 3*k^5*w^2 + k^7)*x^3 + (3*w^6 + 3*k^2*w^4 - 3*k^4*w^2 - 3*k^6)*x^2 + ( - 18*k*w^4 - 12*k^3*w^2 + 6*k^5)*x - 6*w^4 + 36*k^2*w^2 - 6*k^4)*%e^(k*x)*sin(w*x) + (( - w^7 - 3*k^2*w^5 - 3*k^4*w^3 - k^6*w)*x^3 + (6*k*w^5 + 12*k^3*w^3 + 6*k^5*w)*x^2 + (6*w^5 - 12*k^2*w^3 - 18*k^4*w)*x - 24*k*w^3 + 24*k^3*w)*%e^(k*x)*cos(w*x))/(w^8 + 4*k^2*w^6 + 6*k^4*w^4 + 4*k^6*w^2 + k^8)
sage: f = axiom('1/x^2')
sage: f.integrate('x', 1, 'inf')
1
sage: g = axiom('f/sinh(k*x)^4')
sage: g.taylor('x', 0, 3)
f/(k^4*x^4) - 2*f/(3*k^2*x^2) + 11*f/45 - 62*k^2*f*x^2/945
sage: axiom.taylor('asin(x)','x',0, 10)
x + x^3/6 + 3*x^5/40 + 5*x^7/112 + 35*x^9/1152
\subsection{Examples involving matrices}
We illustrate computing with the matrix whose $i,j$ entry
is $i/j$, for $i,j=1,\ldots,4$.
sage: f = axiom.eval('f[i,j] := i/j')
sage: A = axiom('genmatrix(f,4,4)'); A
matrix([1,1/2,1/3,1/4],[2,1,2/3,1/2],[3,3/2,1,3/4],[4,2,4/3,1])
sage: A.determinant()
0
sage: A.echelon()
matrix([1,1/2,1/3,1/4],[0,0,0,0],[0,0,0,0],[0,0,0,0])
sage: A.eigenvalues()
[[0,4],[3,1]]
sage: A.eigenvectors()
[[[0,4],[3,1]],[1,0,0, - 4],[0,1,0, - 2],[0,0,1, - 4/3],[1,2,3,4]]
We can also compute the echelon form in \sage:
sage: B = matrix(QQ, A)
sage: B.echelon_form()
[ 1 1/2 1/3 1/4]
[ 0 0 0 0]
[ 0 0 0 0]
[ 0 0 0 0]
sage: B.charpoly().factor()
(x - 4) * x^3
\subsection{Laplace Transforms}
We illustrate Laplace transforms:
sage: _ = axiom.eval("f(t) := t*sin(t)")
sage: axiom("laplace(f(t),t,s)")
2*s/(s^2 + 1)^2
sage: axiom("laplace(delta(t-3),t,s)") #Dirac delta function
%e^-(3*s)
sage: _ = axiom.eval("f(t) == exp(t)*sin(t)")
sage: axiom("laplace(f(t),t,s)")
1/(s^2 - 2*s + 2)
sage: _ = axiom.eval("f(t) := t^5*exp(t)*sin(t)")
sage: axiom("laplace(f(t),t,s)")
360*(2*s - 2)/(s^2 - 2*s + 2)^4 - 480*(2*s - 2)^3/(s^2 - 2*s + 2)^5 + 120*(2*s - 2)^5/(s^2 - 2*s + 2)^6
sage: axiom("laplace(f(t),t,s)")
3 5
360 (2 s - 2) 480 (2 s - 2) 120 (2 s - 2)
--------------- - --------------- + ---------------
2 4 2 5 2 6
(s - 2 s + 2) (s - 2 s + 2) (s - 2 s + 2)
sage: axiom("laplace(diff(x(t),t),t,s)")
s*laplace(x(t),t,s) - x(0)
sage: axiom("laplace(diff(x(t),t,2),t,s)")
-at('diff(x(t),t,1),t = 0) + s^2*laplace(x(t),t,s) - x(0)*s
sage.: axiom("laplace(diff(x(t),t,2),t,s)")
!
d ! 2
- -- (x(t))! + s laplace(x(t), t, s) - x(0) s
dt !
!t = 0
Even better, use \code{view(axiom("laplace(diff(x(t),t,2),t,s)"))} to see
a typeset version.
\subsection{Continued Fractions}
A continued fraction $a + 1/(b + 1/(c + \cdots))$ is
represented in axiom by the list $[a, b, c, \ldots]$.
sage: axiom("cf((1 + sqrt(5))/2)")
[1,1,1,1,2]
sage: axiom("cf ((1 + sqrt(341))/2)")
[9,1,2,1,2,1,17,1,2,1,2,1,17,1,2,1,2,1,17,2]
\subsection{Special examples}
In this section we illustrate calculations that would be awkward to do
(as far as I know) in non-symbolic computer algebra systems like MAGMA
or GAP.
We compute the gcd of $2x^{n+4} - x^{n+2}$ and $4x^{n+1} + 3x^n$
for arbitrary $n$.
sage: f = axiom('2*x^(n+4) - x^(n+2)')
sage: g = axiom('4*x^(n+1) + 3*x^n')
sage: f.gcd(g)
x^n
You can plot 3d graphs (via gnuplot):
sage.: axiom('plot3d(x^2-y^2, [x,-2,2], [y,-2,2], [grid,12,12])')
[displays a 3 dimensional graph]
You can formally evaluate sums (note the \code{nusum} command):
sage: S = axiom('nusum(exp(1+2*i/n),i,1,n)')
sage.: S
2/n + 3 2/n + 1
%e %e
----------------------- - -----------------------
1/n 1/n 1/n 1/n
(%e - 1) (%e + 1) (%e - 1) (%e + 1)
We formally compute the limit as $n\to\infty$ of $2S/n$ as follows:
sage: T = S*axiom('2/n')
sage: T.tlimit('n','inf')
%e^3 - %e
\subsection{Miscellaneous}
Obtaining digits of $\pi$:
sage: axiom.eval('fpprec : 100')
'100'
sage: axiom(pi).bfloat()
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068b0
Defining functions in axiom:
sage: axiom.eval('fun(a) == a^2')
'fun[a] := a^2'
sage: axiom('fun(10)')
100
\subsection{Interactivity}
Axiom has a non-interactive mode that is initiated via the command
")read file.input".
\subsection{Latex Output}
The latex output of Axiom is not perfect. E.g.,
sage: axiom.eval('tex(sin(u) + sinh(v^2))')
'$$\\sinhv^2 + \\sinu$$false'
Notice the lack of space after the sin macro, which is a latex syntax
error. In \sage this is automatically fixed via a substition for
trig functions, which may have potentially bad side effects:
sage: latex(axiom('sin(u) + sinh(v^2)'))
\sinh v^2+\sin u
It would be nice if somebody would fix this problem. One way would
be to improve Axiom by making the fix to Axiom and giving this back
to the Axiom people.
Here's another example:
sage: g = axiom('exp(3*%i*x)/(6*%i) + exp(%i*x)/(2*%i) + c')
sage: latex(g)
-\frac{i e^{3 i x}}{6}-\frac{i e^{i x}}{2}+c
\subsection{Long Input}
The Axiom interface reads in even very long input (using files) in a
robust manner, as long as you are creating a new object.
\note{Using \code{axiom.eval} for long input
is much less robust, and is not recommended.}
sage: t = '"%s"'%10^10000 # ten thousand character string.
sage: a = axiom(t)
"""
#*****************************************************************************
# Copyright (C) 2005 William Stein <wstein@gmail.com>
#
# Distributed under the terms of the GNU General Public License (GPL)
#
# This code is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# The full text of the GPL is available at:
#
# http://www.gnu.org/licenses/
#*****************************************************************************
import os, re
from expect import Expect, ExpectElement, FunctionElement, ExpectFunction, tmp
from sage.misc.misc import verbose
from sage.misc.multireplace import multiple_replace
cnt = 0
seq = 0
from sage.misc.all import pager, verbose, DOT_SAGE, SAGE_ROOT
COMMANDS_CACHE = '%s/axiom_commandlist_cache.sobj'%DOT_SAGE
import sage.server.support
# The Axiom commands ")what thing det" ")show Matrix" and ")display
# op det" commands, gives a list of all identifiers that begin in
# a certain way. This could maybe be useful somehow... (?) Also
# axiom has a lot a lot of ways for getting documentation from the
# system -- this could also be useful.
class Axiom(Expect):
"""
Interface to the Axiom interpreter.
"""
def __call__(self, x):
import sage.rings.all
if sage.rings.all.is_Infinity(x):
return Expect.__call__(self, 'inf')
else:
return Expect.__call__(self, x)
def __init__(self, script_subdirectory=None, logfile=None, server=None):
"""
Create an instance of the Axiom interpreter.
"""
eval_using_file_cutoff = 200
self.__eval_using_file_cutoff = eval_using_file_cutoff
Expect.__init__(self,
name = 'axiom',
prompt = '\([0-9]+\) -> ',
command = "axiom -nox",
maxread = 1, # CRUCIAL to use less buffering for Axiom!
script_subdirectory = script_subdirectory,
restart_on_ctrlc = False,
verbose_start = False,
init_code = [],
logfile = logfile,
eval_using_file_cutoff=eval_using_file_cutoff)
def __getattr__(self, attrname):
if attrname[:1] == "_":
raise AttributeError
return AxiomExpectFunction(self, attrname)
def _start(self):
# For some reason sending a single input line at startup avoids
# lots of weird timing issues when doing doctests.
Expect._start(self)
#out = self._eval_line(')set output algebra off', reformat=False)
#out = self._eval_line(')set output tex on', reformat="False")
out = self._eval_line(')set message autoload off', reformat=False)
self._expect.expect(self._prompt)
#out = self._expect.before
#print "out 0 = '%s'"%out
#self(1)
def _eval_line_using_file(self, line, tmp):
F = open(tmp, 'w')
F.write(line)
F.close()
if self._expect is None:
self._start()
# For some reason this trivial comp
# keeps certain random freezes from occuring. Do not remove this.
# The space before the \n is also important.
self._expect.sendline(')read "%s"\n'%tmp)
self._expect.expect(self._prompt)
return ''
def __reduce__(self):
return reduce_load_Axiom, tuple([])
def _quit_string(self):
return ')lisp (quit)'
def _eval_line(self, line, reformat=True, allow_use_file=False,
wait_for_prompt=True):
if not wait_for_prompt:
return Expect._eval_line(self, line)
line = line.rstrip().rstrip(';')
if line == '':
return ''
global seq
seq += 1
#start = SAGE_START + str(seq)
#end = SAGE_END + str(seq)
#line = '%s;\n%s; %s;'%(start, line, end)
if self._expect is None:
self._start()
if allow_use_file and self.__eval_using_file_cutoff and \
len(line) > self.__eval_using_file_cutoff:
return self._eval_line_using_file(line, tmp)
try:
E = self._expect
# debug
# print "in = '%s'"%line
E.sendline(line)
self._expect.expect(self._prompt)
out = self._expect.before
# debug
# print "out = '%s'"%out
except KeyboardInterrupt:
self._keyboard_interrupt()
if 'Incorrect syntax:' in out:
raise RuntimeError, out
if not reformat:
return out
if 'error' in out:
return out
#out = out.lstrip()
i = out.find('\n')
out = out[i+1:]
outs = out.split("\n")
i = 0
outline = ''
for line in outs:
line = line.rstrip()
# print "'%s'"%line
if line[:4] == ' (':
i = line.find('(')
i += line[i:].find(')')
if line[i+1:] == "":
i = 0
outs = outs[1:]
break;
out = "\n".join(line[i+1:] for line in outs[1:])
return out
###########################################
# Interactive help
###########################################
def help(self, s):
if sage.server.support.EMBEDDED_MODE:
os.system('asq "describe(%s); "< /dev/null'%s)
else:
os.system('asq "describe(%s);"'%s)
def example(self, s):
if sage.server.support.EMBEDDED_MODE:
os.system('axiom -ht -nogr "example(%s);" < /dev/null'%s)
else:
os.system('axiom -ht "example(%s);"'%s)
describe = help
def demo(self, s):
if sage.server.support.EMBEDDED_MODE:
os.system('axiom -ht "demo(%s);" < /dev/null'%s)
else:
os.system('axiom -ht "demo(%s);"'%s)
def completions(self, s):
"""
Return all commands that complete the command starting with the
string s. This is like typing s[tab] in the maple interpreter.
"""
s = self.eval('apropos(%s)'%s).replace('\\ - ','-')
return [x for x in s[1:-1].split(',') if x[0] != '?']
def _commands(self):
"""
Return list of all commands defined in Axiom.
"""
try:
return self.__commands
except AttributeError:
self.__commands = sum([self.completions(chr(97+n)) for n in range(26)], [])
return self.__commands
def trait_names(self, verbose=True, use_disk_cache=True):
try:
return self.__trait_names
except AttributeError:
import sage.misc.persist
if use_disk_cache:
try:
self.__trait_names = sage.misc.persist.load(COMMANDS_CACHE)
return self.__trait_names
except IOError:
pass
if verbose:
print "\nBuilding Axiom command completion list (this takes"
print "a few seconds only the first time you do it)."
print "To force rebuild later, delete %s."%COMMANDS_CACHE
v = self._commands()
self.__trait_names = v
sage.misc.persist.save(v, COMMANDS_CACHE)
return v
def _object_class(self):
return AxiomElement
def _true_symbol(self):
return 'true'
def _false_symbol(self):
return 'false'
def function(self, args, defn, repr=None, latex=None):
"""
Return the Axiom function with given arguments
and definition.
INPUT:
args -- a string with variable names separated by commas
defn -- a string (or Axiom expression) that defines
a function of the arguments in Axiom.
repr -- an optional string; if given, this is how the function will print.
EXAMPLES:
sage: f = axiom.function('x', 'sin(x)')
sage: f(3.2)
-0.058374143427580086
sage: f = axiom.function('x,y', 'sin(x)+cos(y)')
sage: f(2,3.5)
sin(2) - 0.9364566872907963
sage: f
sin(x)+cos(y)
sage: g = f.integrate('z'); g
(cos(y) + sin(x))*z
sage: g(1,2,3)
3*(cos(2) + sin(1))
The function definition can be an axiom object:
sage: an_expr = axiom('sin(x)*gamma(x)')
sage: t = axiom.function('x', an_expr)
sage: t
gamma(x)*sin(x)
sage: t(2)
sin(2)
sage: float(t(2))
0.90929742682568171
sage: loads(t.dumps())
gamma(x)*sin(x)
"""
name = self._next_var_name()
defn = str(defn)
args = str(args)
axiom.eval('%s(%s) := %s'%(name, args, defn))
if repr is None:
repr = defn
f = AxiomFunction(self, name, repr, args, latex)
return f
def set(self, var, value):
"""
Set the variable var to the given value.
"""
cmd = '%s := %s'%(var, value)
out = self._eval_line(cmd, reformat=False)
#out = self._eval_line(cmd, reformat=False, allow_use_file=True)
if out.find("error") != -1:
raise TypeError, "Error executing code in Axiom\nCODE:\n\t%s\nAxiom ERROR:\n\t%s"%(cmd, out)
def get(self, var):
"""
Get the string value of the variable var.
"""
s = self._eval_line('%s'%var)
return s
#def clear(self, var):
# """
# Clear the variable named var.
# """
# if self._expect is None:
# return
# self._expect.sendline('kill(%s);'%var)
# self._expect.expect(self._prompt)
def console(self):
axiom_console()
def plot2d(self, *args):
r"""
Plot a 2d graph using Axiom / Viewman.
axiom.plot2d(f, '[var, min, max]', options)
INPUT:
f -- a string representing a function (such as f="sin(x)")
[var, xmin, xmax]
options -- an optional string representing plot2d options in gnuplot format
EXAMPLES:
sage.: axiom.plot2d('sin(x)','[x,-5,5]')
sage.: opts = '[gnuplot_term, ps], [gnuplot_out_file, "sin-plot.eps"]'
sage.: axiom.plot2d('sin(x)','[x,-5,5]',opts)
The eps file is saved in the current directory.
"""
self('plot2d(%s)'%(','.join([str(x) for x in args])))
def plot2d_parametric(self, r, var, trange, nticks=50, options=None):
r"""
Plots r = [x(t), y(t)] for t = tmin...tmax using gnuplot with options
INPUT:
r -- a string representing a function (such as r="[x(t),y(t)]")
var -- a string representing the variable (such as var = "t")
trange -- [tmin, tmax] are numbers with tmin<tmax
nticks -- int (default: 50)
options -- an optional string representing plot2d options in gnuplot format
EXAMPLES:
sage.: axiom.plot2d_parametric(["sin(t)","cos(t)"], "t",[-3.1,3.1])
sage.: opts = '[gnuplot_preamble, "set nokey"], [gnuplot_term, ps], [gnuplot_out_file, "circle-plot.eps"]'
sage.: axiom.plot2d_parametric(["sin(t)","cos(t)"], "t", [-3.1,3.1], options=opts)
The eps file is saved to the current working directory.
Here is another fun plot:
sage.: axiom.plot2d_parametric(["sin(5*t)","cos(11*t)"], "t", [0,2*pi()], nticks=400)
"""
tmin = trange[0]
tmax = trange[1]
cmd = "plot2d([parametric, %s, %s, [%s, %s, %s], [nticks, %s]]"%( \
r[0], r[1], var, tmin, tmax, nticks)
if options is None:
cmd += ")"
else:
cmd += ", %s)"%options
self(cmd)
def plot3d(self, *args):
r"""
Plot a 3d graph using Axiom / Viewman.
axiom.plot3d(f, '[x, xmin, xmax]', '[y, ymin, ymax]', '[grid, nx, ny]', options)
INPUT:
f -- a string representing a function (such as f="sin(x)")
[var, min, max]
EXAMPLES:
sage.: axiom.plot3d('1 + x^3 - y^2', '[x,-2,2]', '[y,-2,2]', '[grid,12,12]')
sage.: axiom.plot3d('sin(x)*cos(y)', '[x,-2,2]', '[y,-2,2]', '[grid,30,30]')
sage.: opts = '[gnuplot_term, ps], [gnuplot_out_file, "sin-plot.eps"]'
sage.: axiom.plot3d('sin(x+y)', '[x,-5,5]', '[y,-1,1]', opts)
The eps file is saved in the current working directory.
"""
self('plot3d(%s)'%(','.join([str(x) for x in args])))
def plot3d_parametric(self, r, vars, urange, vrange, options=None):
r"""
Plot a 3d parametric graph with r=(x,y,z), x = x(u,v), y = y(u,v), z = z(u,v),
for u = umin...umax, v = vmin...vmax using gnuplot with options.
INPUT:
x, y, z -- a string representing a function (such as x="u^2+v^2", ...)
vars is a list or two strings representing variables (such as vars = ["u","v"])
urange -- [umin, umax]
vrange -- [vmin, vmax] are lists of numbers with
umin < umax, vmin < vmax
options -- optional string representing plot2d options in gnuplot format
OUTPUT:
displays a plot on screen or saves to a file
EXAMPLES:
sage.: axiom.plot3d_parametric(["v*sin(u)","v*cos(u)","v"], ["u","v"],[-3.2,3.2],[0,3])
sage.: opts = '[gnuplot_term, ps], [gnuplot_out_file, "sin-cos-plot.eps"]'
sage.: axiom.plot3d_parametric(["v*sin(u)","v*cos(u)","v"], ["u","v"],[-3.2,3.2],[0,3],opts)
The eps file is saved in the current working directory.
Here is a torus:
sage.: _ = axiom.eval("expr_1: cos(y)*(10.0+6*cos(x)); expr_2: sin(y)*(10.0+6*cos(x)); expr_3: -6*sin(x);") # optional
sage.: axiom.plot3d_parametric(["expr_1","expr_2","expr_3"], ["x","y"],[0,6],[0,6])
Here is a Mobius strip:
sage.: x = "cos(u)*(3 + v*cos(u/2))"
sage.: y = "sin(u)*(3 + v*cos(u/2))"
sage.: z = "v*sin(u/2)"
sage.: axiom.plot3d_parametric([x,y,z],["u","v"],[-3.1,3.2],[-1/10,1/10])
"""
umin = urange[0]
umax = urange[1]
vmin = vrange[0]
vmax = vrange[1]
cmd = 'plot3d([%s, %s, %s], [%s, %s, %s], [%s, %s, %s]'%(
r[0], r[1], r[2], vars[0], umin, umax, vars[1], vmin, vmax)
if options is None:
cmd += ')'
else:
cmd += ', %s)'%options
axiom(cmd)
def de_solve(axiom, de, vars, ics=None):
"""
Solves a 1st or 2nd order ordinary differential equation (ODE)
in two variables, possibly with initial conditions.
INPUT:
de -- a string representing the ODE
vars -- a list of strings representing the two variables.
ics -- a triple of numbers [a,b1,b2] representing
y(a)=b1, y'(a)=b2
EXAMPLES:
sage.: axiom.de_solve('diff(y,x,2) + 3*x = y', ['x','y'], [1,1,1])
y = 3*x - 2*%e^(x - 1)
sage.: axiom.de_solve('diff(y,x,2) + 3*x = y', ['x','y'])
y = %k1*%e^x + %k2*%e^ - x + 3*x
sage.: axiom.de_solve('diff(y,x) + 3*x = y', ['x','y'])
y = (%c - 3*( - x - 1)*%e^ - x)*%e^x
sage.: axiom.de_solve('diff(y,x) + 3*x = y', ['x','y'],[1,1])
y = - %e^ - 1*(5*%e^x - 3*%e*x - 3*%e)
"""
if not isinstance(vars, str):
str_vars = '%s, %s'%(vars[1], vars[0])
else:
str_vars = vars
axiom.eval('depends(%s)'%str_vars)
m = axiom(de)
a = 'ode2(%s, %s)'%(m.name(), str_vars)
if ics != None:
if len(ics) == 3:
cmd = "ic2("+a+",%s=%s,%s=%s,diff(%s,%s)=%s);"%(vars[0],ics[0], vars[1],ics[1], vars[1], vars[0], ics[2])
return axiom(cmd)
if len(ics) == 2:
return axiom("ic1("+a+",%s=%s,%s=%s);"%(vars[0],ics[0], vars[1],ics[1]))
return axiom(a+";")
def de_solve_laplace(self, de, vars, ics=None):
"""
Solves an ordinary differential equation (ODE) using Laplace transforms.
INPUT:
de -- a string representing the ODE
(e.g., de = "diff(f(x),x,2)=diff(f(x),x)+sin(x)")
vars -- a list of strings representing the variables
(e.g., vars = ["x","f"])
ics -- a list of numbers representing initial conditions,
with symbols allowed which are represented by strings
(eg, f(0)=1, f'(0)=2 is ics = [0,1,2])
EXAMPLES:
sage.: axiom.clear('x'); axiom.clear('f')
sage.: axiom.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)", ["x","f"], [0,1,2])
f(x) = x*%e^x + %e^x
sage.: axiom.clear('x'); axiom.clear('f')
sage.: f = axiom.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)", ["x","f"])
sage.: f
!
x d ! x x
f(x) = x %e (-- (f(x))! ) - f(0) x %e + f(0) %e
dx !
!x = 0
\note{The second equation sets the values of $f(0)$ and
$f'(0)$ in axiom, so subsequent ODEs involving these
variables will have these initial conditions automatically
imposed.}
"""
if not (ics is None):
d = len(ics)
for i in range(0,d-1):
ic = 'atvalue(diff(%s(%s), %s, %s), %s = %s, %s)'%(
vars[1], vars[0], vars[0], i, vars[0], ics[0], ics[1+i])
axiom.eval(ic)
return axiom('desolve(%s, %s(%s))'%(de, vars[1], vars[0]))
def solve_linear(self, eqns,vars):
"""
Wraps axiom's linsolve.
INPUT:
eqns is a list of m strings, each rperesenting a linear question
in m <= n variables
vars is a list of n strings, each representing a variable
EXAMPLES:
sage: eqns = ["x + z = y","2*a*x - y = 2*a^2","y - 2*z = 2"]
sage: vars = ["x","y","z"]
sage: axiom.solve_linear(eqns, vars)
[x = a + 1,y = 2*a,z = a - 1]
"""
eqs = "["
for i in range(len(eqns)):
if i<len(eqns)-1:
eqs = eqs + eqns[i]+","
if i==len(eqns)-1:
eqs = eqs + eqns[i]+"]"
vrs = "["
for i in range(len(vars)):
if i<len(vars)-1:
vrs = vrs + vars[i]+","
if i==len(vars)-1:
vrs = vrs + vars[i]+"]"
return self('linsolve(%s, %s)'%(eqs, vrs))
def unit_quadratic_integer(self, n):
r"""
Finds a unit of the ring of integers of the quadratic number
field $\Q(\sqrt{n})$, $n>1$, using the qunit axiom command.
EXAMPLE:
sage: u = axiom.unit_quadratic_integer(101)
sage: u.parent()
Number Field in a with defining polynomial x^2 - 101
sage: u
a + 10
sage: u = axiom.unit_quadratic_integer(13)
sage: u
5*a + 18
sage: u.parent()
Number Field in a with defining polynomial x^2 - 13
"""
from sage.rings.all import QuadraticField, Integer
# Take square-free part so sqrt(n) doesn't get simplified further by axiom
# (The original version of this function would yield wrong answers if
# n is not squarefree.)
n = Integer(n).square_free_part()
if n < 1:
raise ValueError, "n (=%s) must be >= 1"%n
s = str(self('qunit(%s)'%n)).lower()
r = re.compile('sqrt\(.*\)')
s = r.sub('a', s)
a = QuadraticField(n, 'a').gen()
return eval(s)
def plot_list(self, ptsx, ptsy, options=None):
r"""
Plots a curve determined by a sequence of points.
INPUT:
ptsx -- [x1,...,xn], where the xi and yi are real,
ptsy -- [y1,...,yn]
options -- a string representing axiom plot2d options.
The points are (x1,y1), (x2,y2), etc.
EXAMPLES:
sage.: zeta_ptsx = [ (pari(1/2 + i*I/10).zeta().real()).precision(1) for i in range (70,150)]
sage.: zeta_ptsy = [ (pari(1/2 + i*I/10).zeta().imag()).precision(1) for i in range (70,150)]
sage.: axiom.plot_list(zeta_ptsx, zeta_ptsy)
sage.: opts='[gnuplot_preamble, "set nokey"], [gnuplot_term, ps], [gnuplot_out_file, "zeta.eps"]'
sage.: axiom.plot_list(zeta_ptsx, zeta_ptsy, opts)
"""
cmd = 'plot2d([discrete,%s, %s]'%(ptsx, ptsy)
if options is None:
cmd += ')'
else:
cmd += ', %s)'%options
self(cmd)
def plot_multilist(self, pts_list, options=None):
r"""
Plots a list of list of points pts_list=[pts1,pts2,...,ptsn],
where each ptsi is of the form [[x1,y1],...,[xn,yn]]
x's must be integers and y's reals
options is a string representing axiom plot2d options.
EXAMPLES:
sage.: xx = [ i/10.0 for i in range (-10,10)]
sage.: yy = [ i/10.0 for i in range (-10,10)]
sage.: x0 = [ 0 for i in range (-10,10)]
sage.: y0 = [ 0 for i in range (-10,10)]
sage.: zeta_ptsx1 = [ (pari(1/2+i*I/10).zeta().real()).precision(1) for i in range (10)]
sage.: zeta_ptsy1 = [ (pari(1/2+i*I/10).zeta().imag()).precision(1) for i in range (10)]
sage.: axiom.plot_multilist([[zeta_ptsx1,zeta_ptsy1],[xx,y0],[x0,yy]])
sage.: zeta_ptsx1 = [ (pari(1/2+i*I/10).zeta().real()).precision(1) for i in range (10,150)]
sage.: zeta_ptsy1 = [ (pari(1/2+i*I/10).zeta().imag()).precision(1) for i in range (10,150)]
sage.: axiom.plot_multilist([[zeta_ptsx1,zeta_ptsy1],[xx,y0],[x0,yy]])
sage.: opts='[gnuplot_preamble, "set nokey"]'
sage.: axiom.plot_multilist([[zeta_ptsx1,zeta_ptsy1],[xx,y0],[x0,yy]],opts)
"""
n = len(pts_list)
cmd = '['
for i in range(n):
if i < n-1:
cmd = cmd+'[discrete,'+str(pts_list[i][0])+','+str(pts_list[i][1])+'],'
if i==n-1:
cmd = cmd+'[discrete,'+str(pts_list[i][0])+','+str(pts_list[i][1])+']]'
# debug
# print cmd
if options is None:
self('plot2d('+cmd+')')
else:
self('plot2d('+cmd+','+options+')')
class AxiomElement(ExpectElement):
def __call__(self, x):
self._check_valid()
P = self.parent()
return P('%s[%s]'%(self.name(), x))
def _cmp_(self, other):
"""
EXAMPLES:
sage: a = axiom(1); b = axiom(2)
sage: a == b
False
sage: a < b
True
sage: a > b
False
sage: b < a
False
sage: b > a
True
We can also compare more complicated object such as functions:
sage: f = axiom('sin(x)'); g = axiom('cos(x)')
sage: -f == g.diff('x')
True
"""
# thanks to David Joyner for telling me about using "is".
P = self.parent()
if P.eval("is (%s < %s)"%(self.name(), other.name())) == P._true_symbol():
return -1
elif P.eval("is (%s > %s)"%(self.name(), other.name())) == P._true_symbol():
return 1
elif P.eval("is (%s = %s)"%(self.name(), other.name())) == P._true_symbol():
return 0
else:
return -1 # everything is supposed to be comparable in Python, so we define
# the comparison thus when no comparable in interfaced system.
def numer(self):
P = self.parent()
return P('numeric(%s)'%self._name)
def real(self):
return self.realpart()
def imag(self):
return self.imagpart()
def str(self):
self._check_valid()
P = self.parent()
return P.get('%s::InputForm'%self._name)
def __repr__(self):
self._check_valid()
P = self.parent()
return P.get(self._name)
def diff(self, var='x', n=1):
"""
Return the n-th derivative of self.
INPUT:
var -- variable (default: 'x')
n -- integer (default: 1)
OUTPUT:
n-th derivative of self with respect to the variable var
EXAMPLES:
sage: f = axiom('x^2')
sage: f.diff()
2*x
sage: f.diff('x')
2*x
sage: f.diff('x', 2)
2
sage: axiom('sin(x^2)').diff('x',4)
16*x^4*sin(x^2) - 12*sin(x^2) - 48*x^2*cos(x^2)
sage: f = axiom('x^2 + 17*y^2')
sage: f.diff('x')
2*x
sage: f.diff('y')
34*y
"""
return ExpectElement.__getattr__(self, 'differentiate')(var, n)
derivative = diff
def nintegral(self, var='x', a=0, b=1,
desired_relative_error='1e-8',
maximum_num_subintervals=200):
r"""
Return a numerical approximation to the integral
of self from a to b.
INPUT:
var -- variable to integrate with respect to
a -- lower endpoint of integration
b -- upper endpoint of integration
desired_relative_error -- (default: '1e-8') the desired
relative error
maximum_num_subintervals -- (default: 200) axiom number
of subintervals
OUTPUT:
-- approximation to the integral
-- estimated absolute error of the approximation
-- the number of integrand evaluations
-- an error code:
0 -- no problems were encountered
1 -- too many subintervals were done
2 -- excessive roundoff error
3 -- extremely bad integrand behavior
4 -- failed to converge
5 -- integral is probably divergent or slowly convergent
6 -- the input is invalid
EXAMPLES:
sage: axiom('exp(-sqrt(x))').nintegral('x',0,1)
(0.5284822353142306, 4.1633141378838445E-11, 231, 0)
Note that GP also does numerical integration, and can do
so to very high precision very quickly:
sage: gp('intnum(x=0,1,exp(-sqrt(x)))')
0.5284822353142307136179049194 # 32-bit
0.52848223531423071361790491935415653021 # 64-bit
sage: _ = gp.set_precision(80)
sage: gp('intnum(x=0,1,exp(-sqrt(x)))')
0.52848223531423071361790491935415653021675547587292866196865279321015401702040079
"""
from sage.rings.all import Integer
v = self.quad_qags(var, a, b, desired_relative_error,
maximum_num_subintervals)
return v[0], v[1], Integer(v[2]), Integer(v[3])
def integral(self, var='x', min=None, max=None):
r"""
Return the integral of self with respect to the variable x.
INPUT:
var -- variable
min -- default: None
max -- default: None
Returns the definite integral if xmin is not None,
otherwise returns an indefinite integral.
EXAMPLES:
sage: axiom('x^2+1').integral()
x^3/3 + x
sage: axiom('x^2+ 1 + y^2').integral('y')
y^3/3 + x^2*y + y
sage: axiom('x / (x^2+1)').integral()
log(x^2 + 1)/2
sage: axiom('1/(x^2+1)').integral()
atan(x)
sage.: axiom('1/(x^2+1)').integral('x', 0, infinity)
%pi/2
sage: axiom('x/(x^2+1)').integral('x', -1, 1)
0
sage: f = axiom('exp(x^2)').integral('x',0,1); f
-sqrt(%pi)*%i*erf(%i)/2
sage: f.numer() # I wonder how to get a real number (~1.463)??
-0.8862269254527579*%i*erf(%i)
"""
I = ExpectElement.__getattr__(self, 'integrate')
if min is None:
return I(var)
else:
if max is None:
raise ValueError, "neither or both of min/max must be specified."
return I(var, min, max)
integrate = integral
def __float__(self):
return float(str(self.numer()))
def __len__(self):
"""
Return the length of a list.
EXAMPLES:
sage: v = axiom('[x^i for i in 0..5]')
sage: len(v)
6
"""
self._check_valid()
return int(self.parent().eval('#(%s)'%self.name()))
def __getattr__(self, attrname):
if attrname[:1] == "_":
raise AttributeError
return AxiomFunctionElement(self, attrname)
def __getitem__(self, n):
r"""
Return the n-th element of this list.
\note{Lists are 0-based when accessed via the \sage interface,
not 1-based as they are in the Axiom interpreter.}
EXAMPLES:
sage: v = axiom('[i*x^i for i in 0..5]'); v
[0,x,2*x^2,3*x^3,4*x^4,5*x^5]
sage: v[3]
3*x^3
sage: v[0]
0
sage: v[10]
Traceback (most recent call last):
...
IndexError: n = (10) must be between 0 and 5
"""
n = int(n)
if n < 0 or n >= len(self):
raise IndexError, "n = (%s) must be between %s and %s"%(n, 0, len(self)-1)
return ExpectElement.__getitem__(self, n+1)
def subst(self, val):
P = self.parent()
return P('subst(%s, %s)'%(self.name(), val))
def comma(self, args):
self._check_valid()
P = self.parent()
return P('%s, %s'%(self.name(), args))
def _latex_(self):
self._check_valid()
P = self.parent()
s = axiom._eval_line('outputAsTex(%s)'%self.name(), reformat=False)
if not '$$' in s:
raise RuntimeError, "Error texing axiom object."
i = s.find('$$')
j = s.rfind('$$')
s = s[i+2:j]
s = multiple_replace({'\r\n':' ',
'\\%':'',
'\\arcsin ':'\\sin^{-1} ',
'\\arccos ':'\\cos^{-1} ',
'\\arctan ':'\\tan^{-1} '}, s)
return s
def trait_names(self):
return self.parent().trait_names()
def _matrix_(self, R):
r"""
If self is a Axiom matrix, return the corresponding \sage
matrix over the \sage ring $R$.
This may or may not work depending in how complicated the
entries of self are! It only works if the entries of self
can be coerced as strings to produce meaningful elements
of $R$.
EXAMPLES:
sage: _ = axiom.eval("f[i,j] := i/j")
sage: A = axiom('genmatrix(f,4,4)'); A
matrix([1,1/2,1/3,1/4],[2,1,2/3,1/2],[3,3/2,1,3/4],[4,2,4/3,1])
sage: A._matrix_(QQ)
[ 1 1/2 1/3 1/4]
[ 2 1 2/3 1/2]
[ 3 3/2 1 3/4]
[ 4 2 4/3 1]
You can also use the \code{matrix} command (which is defined
in \code{sage.misc.functional}):
sage: matrix(QQ, A)
[ 1 1/2 1/3 1/4]
[ 2 1 2/3 1/2]
[ 3 3/2 1 3/4]
[ 4 2 4/3 1]
"""
from sage.matrix.all import MatrixSpace
self._check_valid()
P = self.parent()
nrows = int(P.eval('length(%s)'%self.name()))
if nrows == 0:
return MatrixSpace(R, 0, 0)(0)
ncols = int(P.eval('length(%s[1])'%self.name()))
M = MatrixSpace(R, nrows, ncols)
s = self.str().replace('matrix','').replace(',',"','").\
replace("]','[","','").replace('([',"['").replace('])',"']")
s = eval(s)
return M([R(x) for x in s])
def partial_fraction_decomposition(self, var='x'):
"""
Return the partial fraction decomposition of self with respect to
the variable var.
EXAMPLES:
sage: f = axiom('1/((1+x)*(x-1))')
sage: f.partial_fraction_decomposition('x')
1 1
--------- - ---------
2 (x - 1) 2 (x + 1)
"""
return self.partfrac(var)
class AxiomFunctionElement(FunctionElement):
def _sage_doc_(self):
return self._obj.parent().help(self._name)
class AxiomExpectFunction(ExpectFunction):
def _sage_doc_(self):
M = self._parent
return M.help(self._name)
class AxiomFunction(AxiomElement):
def __init__(self, parent, name, defn, args, latex):
AxiomElement.__init__(self, parent, name, is_name=True)
self.__defn = defn
self.__args = args
self.__latex = latex
def __call__(self, *x):
self._check_valid()
P = self.parent()
if len(x) == 1:
x = '(%s)'%x
return P('%s%s'%(self.name(), x))
def __repr__(self):
return self.__defn
def _latex_(self):
if self.__latex is None:
return '\\mbox{\\rm %s}'%self.__defn
else:
return self.__latex
def integrate(self, var):
return self.integral(var)
def integral(self, var):
self._check_valid()
P = self.parent()
f = P('integrate(%s(%s), %s)'%(self.name(), self.__args, var))
if var in self.__args.split(','):
args = self.__args
else:
args = self.__args + ',' + var
return P.function(args, str(f))
def is_AxiomElement(x):
return isinstance(x, AxiomElement)
# An instance
axiom = Axiom(script_subdirectory=None)
def reduce_load_Axiom():
return axiom
import os
def axiom_console():
os.system('axiom')
def __doctest_cleanup():
import sage.interfaces.quit
sage.interfaces.quit.expect_quitall()
@
\end{document}
