axiom-commit Mailing List for Axiom Computer Algebra System
Brought to you by:
daly,
starseeker
Menu
▾
▴
axiom-commit — Journal Commits to Axiom SVN
You can subscribe to this list here.
| 2006 |
Jan
|
Feb
|
Mar
|
Apr
|
May
(12) |
Jun
(30) |
Jul
(5) |
Aug
(44) |
Sep
(54) |
Oct
(83) |
Nov
(138) |
Dec
(77) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2007 |
Jan
(8) |
Feb
(19) |
Mar
(36) |
Apr
(57) |
May
(62) |
Jun
(95) |
Jul
(51) |
Aug
(79) |
Sep
(39) |
Oct
(25) |
Nov
(8) |
Dec
(33) |
| 2008 |
Jan
(35) |
Feb
(37) |
Mar
(41) |
Apr
(52) |
May
(42) |
Jun
(55) |
Jul
(42) |
Aug
(47) |
Sep
(45) |
Oct
(24) |
Nov
(53) |
Dec
(102) |
| 2009 |
Jan
(39) |
Feb
(58) |
Mar
(48) |
Apr
(77) |
May
(269) |
Jun
(111) |
Jul
(21) |
Aug
(102) |
Sep
(54) |
Oct
(73) |
Nov
(42) |
Dec
(26) |
| 2010 |
Jan
(45) |
Feb
(35) |
Mar
(15) |
Apr
(25) |
May
(80) |
Jun
(62) |
Jul
(28) |
Aug
(27) |
Sep
(31) |
Oct
(40) |
Nov
(9) |
Dec
(21) |
| 2011 |
Jan
(21) |
Feb
(16) |
Mar
(12) |
Apr
(9) |
May
(23) |
Jun
(2) |
Jul
(19) |
Aug
(18) |
Sep
(20) |
Oct
(5) |
Nov
(19) |
Dec
(19) |
| 2012 |
Jan
(20) |
Feb
(21) |
Mar
(12) |
Apr
(19) |
May
(14) |
Jun
(3) |
Jul
|
Aug
(1) |
Sep
|
Oct
(1) |
Nov
|
Dec
|
| 2013 |
Jan
(12) |
Feb
(16) |
Mar
(65) |
Apr
(3) |
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| 2016 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
(1) |
| 2020 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
(1) |
Nov
|
Dec
|
| 2022 |
Jan
|
Feb
|
Mar
|
Apr
|
May
(1) |
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| S | M | T | W | T | F | S |
|---|---|---|---|---|---|---|
|
|
|
|
|
|
|
1
(6) |
|
2
|
3
(2) |
4
|
5
(2) |
6
|
7
|
8
(4) |
|
9
(4) |
10
|
11
(1) |
12
|
13
|
14
(8) |
15
(2) |
|
16
(4) |
17
|
18
|
19
|
20
|
21
(2) |
22
|
|
23
|
24
|
25
|
26
|
27
|
28
(2) |
29
(2) |
|
30
|
|
|
|
|
|
|
|
From: <gi...@ax...> - 2007-09-29 15:59:01
|
changelog | 2 +
src/input/Makefile.pamphlet | 6 +-
src/input/pfaffian.input.pamphlet | 483 -------------------------------------
3 files changed, 5 insertions(+), 486 deletions(-)
New commits:
commit 929c79aa86c9a3e375ab255b7378211c082a1df8
Author: Tim Daly <root@localhost.localdomain>
Date: Sat Sep 22 22:25:01 2007 -0400
20070929 tpd src/input/Makefile pfaffian regression test removed
20070929 tpd src/input/pfaffian.input regression test removed
|
|
From: <da...@us...> - 2007-09-29 15:58:56
|
Revision: 752
http://axiom.svn.sourceforge.net/axiom/?rev=752&view=rev
Author: daly
Date: 2007-09-29 08:58:56 -0700 (Sat, 29 Sep 2007)
Log Message:
-----------
20070929 tpd src/input/Makefile pfaffian regression test removed
20070929 tpd src/input/pfaffian.input regression test removed
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/input/Makefile.pamphlet
Removed Paths:
-------------
trunk/axiom/src/input/pfaffian.input.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-28 03:13:09 UTC (rev 751)
+++ trunk/axiom/changelog 2007-09-29 15:58:56 UTC (rev 752)
@@ -1,3 +1,5 @@
+20070929 tpd src/input/Makefile pfaffian regression test removed
+20070929 tpd src/input/pfaffian.input regression test removed
20070927 tpd src/input/Makefile pfaffian regression test added
20070927 mxr src/input/pfaffian.input regression test added
20070920 tpd src/input/Makefile add bug101.input regression test
Modified: trunk/axiom/src/input/Makefile.pamphlet
===================================================================
--- trunk/axiom/src/input/Makefile.pamphlet 2007-09-28 03:13:09 UTC (rev 751)
+++ trunk/axiom/src/input/Makefile.pamphlet 2007-09-29 15:58:56 UTC (rev 752)
@@ -340,7 +340,7 @@
padic.regress parabola.regress pascal1.regress pascal.regress \
patch51.regress page.regress \
patmatch.regress pat.regress perman.regress perm.regress \
- pfaffian.regress pfr1.regress pfr.regress pmint.regress \
+ pfr1.regress pfr.regress pmint.regress \
poly1.regress polycoer.regress poly.regress psgenfcn.regress \
quat1.regress quat.regress r20abugs.regress r20bugs.regress \
r21bugsbig.regress r21bugs.regress radff.regress radix.regress \
@@ -593,7 +593,7 @@
${OUT}/pascal.input \
${OUT}/patch51.input \
${OUT}/patmatch.input ${OUT}/perman.input \
- ${OUT}/perm.input ${OUT}/pfaffian.input \
+ ${OUT}/perm.input \
${OUT}/pfr.input ${OUT}/pfr1.input \
${OUT}/pinch.input ${OUT}/plotfile.input ${OUT}/pollevel.input \
${OUT}/pmint.input ${OUT}/polycoer.input \
@@ -870,7 +870,7 @@
${DOC}/pat.input.dvi ${DOC}/patch51.input.dvi \
${DOC}/patmatch.input.dvi \
${DOC}/pdecomp0.as.dvi ${DOC}/perman.input.dvi \
- ${DOC}/perm.input.dvi ${DOC}/pfaffian.input.dvi \
+ ${DOC}/perm.input.dvi \
${DOC}/pfr1.input.dvi \
${DOC}/pfr.input.dvi ${DOC}/pinch.input.dvi \
${DOC}/plotfile.input.dvi ${DOC}/plotlist.input.dvi \
Deleted: trunk/axiom/src/input/pfaffian.input.pamphlet
===================================================================
--- trunk/axiom/src/input/pfaffian.input.pamphlet 2007-09-28 03:13:09 UTC (rev 751)
+++ trunk/axiom/src/input/pfaffian.input.pamphlet 2007-09-29 15:58:56 UTC (rev 752)
@@ -1,483 +0,0 @@
-\documentclass{article}
-\usepackage{amsfonts}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/input pfaffian}
-\author{Timothy Daly, Gunter Rote, Martin Rubey}
-\maketitle
-\begin{abstract}
-In mathematics, the determinant of a skew-symmetric matrix can always
-be written as the square of a polynomial in the matrix entries. This
-polynomial is called the Pfaffian of the matrix. The Pfaffian is
-nonvanishing only for $2n \times 2n$ skew-symmetric matrices, in which
-case it is a polynomial of degree $n$.
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{Examples}
-$$
-{\rm Pfaffian}\left[
-\begin{array}{cc}
-0 & a\\
--a & 0
-\end{array}
-\right] = a
-$$
-
-$$
-{\rm Pfaffian}\left[
-\begin{array}{cccc}
-0 & a & b & c\\
--a & 0 & d & e\\
--b & -d & 0 & f\\
--c & -e & -f & 0
-\end{array}
-\right] = af -be + dc
-$$
-
-$$
-{\rm Pfaffian}\left[
-\begin{array}{cccc}
-\begin{array}{cc}
-0 & \lambda_1\\
--\lambda_1 & 0
-\end{array} & 0 & \cdots & 0\\
-0 &
-\begin{array}{cc}
-0 & \lambda_2\\
--\lambda_2 & 0
-\end{array} & & 0\\
-\vdots & & \ddots & \vdots\\
-0 & 0 & \cdots &
-\begin{array}{cc}
-0 & \lambda_n\\
--\lambda_n & 0
-\end{array}
-\end{array}
-\right] = \lambda_1\lambda_2\cdots\lambda_n
-$$
-\section{Formal definition}
-
-Let $A = \{a_{i,j}\}$ be a $2n \times 2n$ skew-symmetric matrix.
-The Pfaffian of A is defined by the equation
-
-$$Pf(A) = \frac{1}{s^n n!}\sum_{\sigma\in{}S_{2n}}sgn(\sigma)
-\prod_{i=1}^n{a_{\sigma(2i-1),\sigma(2i)}}$$
-
-where $S_{2n}$ is the symmetric group and $sgn(\sigma)$
-is the signature of $\sigma$.
-
-One can make use of the skew-symmetry of $A$ to avoid summing over all
-possible permutations. Let $\Pi$ be the set of all partitions of
-$\{1, 2, \ldots, 2n\}$ into pairs without regard to order.
-There are $(2n - 1)!!$ such partitions. An element
-$\alpha \in \Pi$, can be written as
-
-$$\alpha=\{(i_1,j_1),(i_2,j_2),\cdots,(i_n,j_n)\}$$
-
-with $i_k < j_k$ and $i_1 < i_2 < \cdots < i_n$. Let
-
-$$
-\pi=
-\left[
-\begin{array}{cccccc}
-1 & 2 & 3 & 4 & \cdots & 2n \\
-i_1 & j_1 & i_2 & j_2 & \cdots & j_{n}
-\end{array}
-\right]
-$$
-
-be a corresponding permutation. This depends only on the partition $\alpha$
-and not on the particular choice of $\Pi$. Given a partition $\alpha$ as above
-define
-
-$$A_\alpha =sgn(\pi)a_{i_1,j_1}a_{i_2,j_2}\cdots a_{i_n,j_n}$$
-
-The Pfaffian of $A$ is then given by
-
-$$Pf(A)=\sum_{\alpha\in\Pi} A_\alpha$$
-
-The Pfaffian of a $n\times n$ skew-symmetric matrix for n odd is
-defined to be zero.
-
-\subsection{Alternative definition}
-
-One can associate to any skew-symmetric $2n \times 2n$
-matrix $A=\{a_{ij}\}$ a bivector
-
-$$\omega=\sum_{i<j} a_{ij}~e^i\wedge e^j$$
-
-where $\{e^1, e^2, \ldots, e^{2n}\}$ is the standard basis of
-$\mathbb{R}^{2n}$. The Pfaffian is then defined by the equation
-
-$$\frac{1}{n!}\omega^n = \mbox{Pf}(A)~e^1\wedge e^2\wedge\cdots\wedge e^{2n}$$
-
-here $\omega^n$ denotes the wedge product of $n$ copies of
-$\omega^n$ with itself.
-
-\subsection{Derivation from Determinant}
-
-The Pfaffian can be derived from the determinant for a skew-symmetric
-matrix $A$ as follows. Using Laplace's formula we can write the
-determinant as
-
-$$\det(A)=(-1)^{p+1}a_{p1}A_{p1} +
-(-1)^{p+2}a_{p2}A_{p2}+\cdots+(-1)^{n+p}a_{pn}A_{pn}$$
-
-where $A_{pi}$ is the $p,i$ minor matrix of the matrix $A$. We further
-use Laplace's formula to note that
-
-$$\det(A[A_{ij}]) = \vert A \vert^n$$
-
-since this determinant is that of an $n \times n$ matrix whose only
-non-zero elements are the diagonals (each with value det(A)) and
-$[A_{ij}]$ is a matrix whose $i,j$th component is the corresponding $i,j$
-minor matrix. In this way, following a proof by Parameswaran, we can
-write the determinant as,
-
-$$\det(A)\equiv\Delta_n=
-\left|
-\begin{array}{cccc}
-a_{11}&a_{12}&\cdots&a_{1n}\\
-a_{21}&a_{22}&\cdots&a_{2n}\\
-\vdots&&&\vdots\\
-a_{n1}&a_{n2}&\cdots&a_{nn}
-\end{array}
-\right|
-$$
-
-The minor of
-$$
-\left|
-\begin{array}{cc}
-a_{11}&a_{12}\\
-a_{21}&a_{22}
-\end{array}
-\right|
-$$
-would be $\Delta_{n - 2}$. With this notation
-
-$$\left|
-\begin{array}{cccc}
-1&0&\cdots&0\\
-0&1&\cdots&0\\
-a_{31}&a_{32}&\cdots&a_{3n}\\
-\cdots&\cdots&\cdots&\cdots\\
-a_{n1}&a_{n2}&\cdots&a_{nn}
-\end{array}
-\right|
-\times
-\left|
-\begin{array}{cccc}
-A_{11}&A_{21}&\cdots&A_{n1}\\
-A_{12}&A_{22}&\cdots&A_{n2}\\
-A_{13}&A_{23}&\cdots&A_{n3}\\
-\cdots&\cdots&\cdots&\cdots\\
-A_{1n}&A_{2n}&\cdots&A_{nn}
-\end{array}
-\right| =
-\left|
-\begin{array}{cc}
-A_{11}&A_{21}\\
-A_{12}&a_{22}
-\end{array}
-\right|
-\Delta_n^{n-2}
-$$
-
-Thus
-
-$$\Delta_{n-2}\Delta_n^{n-1}=
-\left|
-\begin{array}{cc}
-A_{11}&A_{21}\\
-A_{12}&a_{22}
-\end{array}
-\right|
-\Delta_n^{n-2}
-$$
-
-Of course, it was only arbitrarily that we chose to remove the first
-two rows, and more generically we can write
-
-$$A_{rr}A_{ss} - A_{rs}A_{sr} = \Delta_{rs,rs}\Delta_n$$
-
-where $\Delta_{rs,rs}$ is the determinant of the original matrix with
-the rows $r$ and $s$, as well as the columns $r$ and $s$ removed. The
-equation above simplifies in the skew-symmetric case to
-
-$$A_{rs} = \sqrt{\Delta_{rs,rs}\Delta_n}$$
-
-We now plug this back into the original formula for the determinant,
-
-$$\Delta_n=
-(-1)^{p+1}a_{p1}\sqrt{\Delta_{p1,p1}\Delta_n} +
-(-1)^{p+2}a_{p2}\sqrt{\Delta_{p2,p2}\Delta_n} +
-\cdots +
-(-1)^{n+p}a_{pn}\sqrt{\Delta_{pn,pn}\Delta_n}
-$$
-
-or with slight manipulation,
-
-$$\sqrt{\Delta_n}=(-1)^{p+1}
-\left( \ a_{p1}\sqrt{\Delta_{p1,p1}} -
-a_{p2}\sqrt{\Delta_{p2,p2}} +
-\cdots +
-(-1)^{n-1}a_{pn}\sqrt{\Delta_{pn,pn}} \
-\right)
-$$
-
-The determinant is thus the square of the right hand side, and so we
-identify the right hand side as the Pfaffian.
-
-\section{Identities}
-
-For a $2n \times 2n$ skew-symmetric matrix $A$ and an arbitrary
-$2n \times 2n$ matrix B,
-\begin{itemize}
-\item $Pf(A)^2 = \det(A)$
-\item $Pf(BAB^T) = \det(B)Pf(A)$
-\item $Pf(\lambda{}A) = \lambda^nPf(A)$
-\item $Pf(A^T) = ( - 1)^nPf(A)$
-\item For a block-diagonal matrix
-$$A_1\oplus A_2=
-\left[
-\begin{array}{cc}
-A_1 & 0 \\
-0 & A_2
-\end{array}
-\right]
-$$
-$$Pf(A1\oplus A2) = Pf(A_11)Pf(A_2)$$
-\item For an arbitrary $n \times n$ matrix $M$:
-$$\mbox{Pf}
-\left[
-\begin{array}{cc}
-0 & M \\
--M^T & 0
-\end{array}
-\right] =
-(-1)^{n(n-1)/2}\det M
-$$
-\end{itemize}
-\section{Applications}
-
-The Pfaffian is an invariant polynomial of a skew-symmetric matrix
-(note that it is not invariant under a general change of basis but
-rather under a proper orthogonal transformation). As such, it is
-important in the theory of characteristic classes. In particular, it
-can be used to define the Euler class of a Riemannian manifold which
-is used in the generalized Gauss-Bonnet theorem.
-
-The number of perfect matchings in a planar graph turns out to be the
-absolute value of a Pfaffian, hence is polynomial time
-computable. This is surprising given that for a general graph, the
-problem is very difficult (so called \#P-complete). This result is used
-to calculate the partition function of Ising models of spin glasses in
-physics, respectively of Markov random fields in machine learning
-(Globerson and Jaakkola, 2007), where the underlying graph is
-planar. Recently it is also used to derive efficient algorithms for
-some otherwise seemingly intractable problems, including the efficient
-simulation of certain types of restricted quantum computation.
-
-The calculation of the number of possible ways to tile a standard
-chessboard or 8-by-8 checkerboard with 32 dominoes is a simple example
-of a problem which may be solved through the use of the Pfaffian
-technique. There are 12,988,816 possible ways to tile a chessboard in
-this manner. Specifically, 12988816 is the number of possible ways to
-cover an 8-by-8 square with 32 1-by-2 rectangles. 12988816 is a square
-number: $12988816 = 3604^2$). Note that 12988816 can be written in the
-form: $2\times 1802^2 + 2\times 1802^2$,
-where all the numbers have a digital root of 2.
-
-More generally, the number of ways to cover a $2n \times 2n$ square with
-$2n^2$ dominoes (as calculated independently by Temperley and M.E. Fisher and
-Kasteleyn in 1961) is given by
-
-$$
-\prod_{j=1}^N
-\prod_{k=1}^N
-\left ( 4\cos^2 \frac{\pi j}{2n + 1} +
-4\cos^2 \frac{\pi k}{2n + 1} \right )
-$$
-
-This technique can be applied in many mathematics-related subjects,
-for example, in the classical, 2-dimensional computation of the
-dimer-dimer correlator function in quantum mechanics.
-
-\section{History}
-
-The term Pfaffian was introduced by Arthur Cayley, who used the term
-in 1852: ``The permutants of this class (from their connection with the
-researches of Pfaff on differential equations) I shall term
-Pfaffians.'' The term honors German mathematician Johann Friedrich
-Pfaff.
-
-\section{Axiom code}
-I have hacked together an algorithm to compute a Pfaffian, using an algorithm
-of Gunter Rote. Currently it's only an .input script, but if it's useful for
-somebody else than myself, we could make it a little more professional.
-
-Martin
-
-<<*>>=
-)spool pfaffian.output
-)set message test on
-)set message auto off
-)clear all
-
---S 1 of 12
-B0 n == matrix [[(if i=j+1 and odd? j then -1 else _
- if i=j-1 and odd? i then 1 else 0) _
- for j in 1..n] for i in 1..n]
---R
---R Type: Void
---E 1
-
---S 2 of 12
-PfChar(lambda, A) ==
- n := nrows A
- (n = 2) => lambda^2 + A.(1,2)
- M := subMatrix(A, 3, n, 3, n)
- r := subMatrix(A, 1, 1, 3, n)
- s := subMatrix(A, 3, n, 2, 2)
-
- p := PfChar(lambda, M)
- d := degree(p, lambda)
-
- B := B0(n-2)
- C := r*B
- g := [(C*s).(1,1), A.(1,2), 1]
- if d >= 4 then
- B := M*B
- for i in 4..d by 2 repeat
- C := C*B
- g := cons((C*s).(1,1), g)
- g := reverse! g
-
- res := 0
- for i in 0..d by 2 for j in 2..d+2 repeat
- c := coefficient(p, lambda, i)
- for e in first(g, j) for k in 2..-d by -2 repeat
- res := res + c * e * lambda^(k+i)
-
- res
---R
---R Type: Void
---E 2
-
---S 3 of 12
-pfaffian A == eval(PfChar(l, A), l=0)
---R
---R Type: Void
---E 3
-
---S 4 of 12
-m:Matrix(Integer):=[[0,15],[-15,0]]
---R
---R + 0 15+
---R (4) | |
---R +- 15 0 +
---R Type: Matrix Integer
---E 4
-
---S 5 of 12
-pfaffian m
---R
---R (5) 15
---R Type: Polynomial Integer
---E 5
-
---S 6 of 12
-m1:Matrix(Polynomial(Integer)):=[[0,a,b,c],[-a,0,d,e],[-b,-d,0,f],[-c,-e,-f,0]]
---R
---R
---R + 0 a b c+
---R | |
---R |- a 0 d e|
---R (6) | |
---R |- b - d 0 f|
---R | |
---R +- c - e - f 0+
---R Type: Matrix Polynomial Integer
---E 6
-
---S 7 of 12
-pfaffian m1
---R
---R Compiling function B0 with type PositiveInteger -> Matrix Integer
---R
---R (7) a f - b e + c d
---R Type: Polynomial Integer
---E 7
-
---S 8 of 12
-(a,b,c,d,e,f):=(3,5,7,11,13,17)
---R
---R (8) 17
---R Type: PositiveInteger
---E 8
-
---S 9 of 12
-m1
---R
---R
---R + 0 a b c+
---R | |
---R |- a 0 d e|
---R (9) | |
---R |- b - d 0 f|
---R | |
---R +- c - e - f 0+
---R Type: Matrix Polynomial Integer
---E 9
-
---S 10 of 12
-a*f-b*e+d*c
---R
---R
---R (10) 63
---R Type: PositiveInteger
---E 10
-
---S 11 of 12
-n:=pfaffian m1
---R
---R (11) a f - b e + c d
---R Type: Polynomial Integer
---E 9
-
---S 12 of 12
-eval(n,['a,'b,'c,'d,'e,'f]::List(Symbol),[a,b,c,d,e,f])
---R
---R (12) 63
---R Type: Polynomial Integer
---E 12
-)spool
-)lisp (bye)
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} {\bf ``http://en.wikipedia.org/wiki/Pfaffian''}
-\bibitem{2} ``The statistics of dimers on a lattice, Part I'', Physica,
-vol.27, 1961, pp.1209-25, P.W. Kasteleyn.
-\bibitem{3} Propp, James (2004),
-``Lambda-determinants and domino-tilings'', arXiv:math.CO/0406301.
-\bibitem{4} Globerson, Amir and Tommi Jaakkola (2007),
-``Approximate inference using planar graph decomposition'',
-Advances in Neural Information Processing Systems 19, MIT Press.
-\bibitem{5} ``The Games and Puzzles Journal'',
-No.14, 1996, pp.204-5, Robin J. Chapman, University of Exeter
-\bibitem{6} ``Domino Tilings and Products of Fibonacci and Pell numbers'',
-Journal of Integer Sequences, Vol.5, 2002, Article 02.1.2,
-James A. Sellers, The Pennsylvania State University
-\bibitem{7} ``The Penguin Dictionary of Curious and Interesting Numbers'',
-revised ed., 1997, ISBN 0-14-026149-4, David Wells, p.182.
-\bibitem{8} ``A Treatise on the Theory of Determinants'',
-1882, Macmillan and Co., Thomas Muir, Online
-\bibitem{9} ``Skew-Symmetric Determinants'',
-The American Mathematical Monthly, vol. 61, no.2., 1954, p.116,
-S. Parameswaran Online-Subscription
-\end{thebibliography}
-\end{document}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <gi...@ax...> - 2007-09-28 03:13:16
|
changelog | 2 +
src/input/Makefile.pamphlet | 12 +-
src/input/pfaffian.input.pamphlet | 483 +++++++++++++++++++++++++++++++++++++
3 files changed, 492 insertions(+), 5 deletions(-)
New commits:
commit 26f1835ed63b5046b731579766e1705d3bb1bb97
Author: Tim Daly <root@localhost.localdomain>
Date: Sat Sep 22 02:32:08 2007 -0400
20070927 mxr add pfaffian code
|
|
From: <da...@us...> - 2007-09-28 03:13:08
|
Revision: 751
http://axiom.svn.sourceforge.net/axiom/?rev=751&view=rev
Author: daly
Date: 2007-09-27 20:13:09 -0700 (Thu, 27 Sep 2007)
Log Message:
-----------
20070927 tpd src/input/Makefile pfaffian regression test added
20070927 mxr src/input/pfaffian.input regression test added
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/input/Makefile.pamphlet
Added Paths:
-----------
trunk/axiom/src/input/pfaffian.input.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-21 15:55:49 UTC (rev 750)
+++ trunk/axiom/changelog 2007-09-28 03:13:09 UTC (rev 751)
@@ -1,3 +1,5 @@
+20070927 tpd src/input/Makefile pfaffian regression test added
+20070927 mxr src/input/pfaffian.input regression test added
20070920 tpd src/input/Makefile add bug101.input regression test
20070920 tpd src/input/bug101.input test laplace(log(z),z,w) bug 101
20070920 wxh src/algebra/laplace.spad fix laplace(log(z),z,w) bug 101
Modified: trunk/axiom/src/input/Makefile.pamphlet
===================================================================
--- trunk/axiom/src/input/Makefile.pamphlet 2007-09-21 15:55:49 UTC (rev 750)
+++ trunk/axiom/src/input/Makefile.pamphlet 2007-09-28 03:13:09 UTC (rev 751)
@@ -338,9 +338,9 @@
oct.regress ode.regress odpol.regress op1.regress \
opalg.regress operator.regress op.regress ovar.regress \
padic.regress parabola.regress pascal1.regress pascal.regress \
- patch51.regress \
+ patch51.regress page.regress \
patmatch.regress pat.regress perman.regress perm.regress \
- pfr1.regress pfr.regress page.regress pmint.regress \
+ pfaffian.regress pfr1.regress pfr.regress pmint.regress \
poly1.regress polycoer.regress poly.regress psgenfcn.regress \
quat1.regress quat.regress r20abugs.regress r20bugs.regress \
r21bugsbig.regress r21bugs.regress radff.regress radix.regress \
@@ -592,8 +592,9 @@
${OUT}/pascal1.input \
${OUT}/pascal.input \
${OUT}/patch51.input \
- ${OUT}/patmatch.input ${OUT}/perman.input \
- ${OUT}/perm.input ${OUT}/pfr.input ${OUT}/pfr1.input \
+ ${OUT}/patmatch.input ${OUT}/perman.input \
+ ${OUT}/perm.input ${OUT}/pfaffian.input \
+ ${OUT}/pfr.input ${OUT}/pfr1.input \
${OUT}/pinch.input ${OUT}/plotfile.input ${OUT}/pollevel.input \
${OUT}/pmint.input ${OUT}/polycoer.input \
${OUT}/poly1.input ${OUT}/psgenfcn.input \
@@ -869,7 +870,8 @@
${DOC}/pat.input.dvi ${DOC}/patch51.input.dvi \
${DOC}/patmatch.input.dvi \
${DOC}/pdecomp0.as.dvi ${DOC}/perman.input.dvi \
- ${DOC}/perm.input.dvi ${DOC}/pfr1.input.dvi \
+ ${DOC}/perm.input.dvi ${DOC}/pfaffian.input.dvi \
+ ${DOC}/pfr1.input.dvi \
${DOC}/pfr.input.dvi ${DOC}/pinch.input.dvi \
${DOC}/plotfile.input.dvi ${DOC}/plotlist.input.dvi \
${DOC}/pollevel.input.dvi ${DOC}/poly1.input.dvi \
Added: trunk/axiom/src/input/pfaffian.input.pamphlet
===================================================================
--- trunk/axiom/src/input/pfaffian.input.pamphlet (rev 0)
+++ trunk/axiom/src/input/pfaffian.input.pamphlet 2007-09-28 03:13:09 UTC (rev 751)
@@ -0,0 +1,483 @@
+\documentclass{article}
+\usepackage{amsfonts}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/input pfaffian}
+\author{Timothy Daly, Gunter Rote, Martin Rubey}
+\maketitle
+\begin{abstract}
+In mathematics, the determinant of a skew-symmetric matrix can always
+be written as the square of a polynomial in the matrix entries. This
+polynomial is called the Pfaffian of the matrix. The Pfaffian is
+nonvanishing only for $2n \times 2n$ skew-symmetric matrices, in which
+case it is a polynomial of degree $n$.
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{Examples}
+$$
+{\rm Pfaffian}\left[
+\begin{array}{cc}
+0 & a\\
+-a & 0
+\end{array}
+\right] = a
+$$
+
+$$
+{\rm Pfaffian}\left[
+\begin{array}{cccc}
+0 & a & b & c\\
+-a & 0 & d & e\\
+-b & -d & 0 & f\\
+-c & -e & -f & 0
+\end{array}
+\right] = af -be + dc
+$$
+
+$$
+{\rm Pfaffian}\left[
+\begin{array}{cccc}
+\begin{array}{cc}
+0 & \lambda_1\\
+-\lambda_1 & 0
+\end{array} & 0 & \cdots & 0\\
+0 &
+\begin{array}{cc}
+0 & \lambda_2\\
+-\lambda_2 & 0
+\end{array} & & 0\\
+\vdots & & \ddots & \vdots\\
+0 & 0 & \cdots &
+\begin{array}{cc}
+0 & \lambda_n\\
+-\lambda_n & 0
+\end{array}
+\end{array}
+\right] = \lambda_1\lambda_2\cdots\lambda_n
+$$
+\section{Formal definition}
+
+Let $A = \{a_{i,j}\}$ be a $2n \times 2n$ skew-symmetric matrix.
+The Pfaffian of A is defined by the equation
+
+$$Pf(A) = \frac{1}{s^n n!}\sum_{\sigma\in{}S_{2n}}sgn(\sigma)
+\prod_{i=1}^n{a_{\sigma(2i-1),\sigma(2i)}}$$
+
+where $S_{2n}$ is the symmetric group and $sgn(\sigma)$
+is the signature of $\sigma$.
+
+One can make use of the skew-symmetry of $A$ to avoid summing over all
+possible permutations. Let $\Pi$ be the set of all partitions of
+$\{1, 2, \ldots, 2n\}$ into pairs without regard to order.
+There are $(2n - 1)!!$ such partitions. An element
+$\alpha \in \Pi$, can be written as
+
+$$\alpha=\{(i_1,j_1),(i_2,j_2),\cdots,(i_n,j_n)\}$$
+
+with $i_k < j_k$ and $i_1 < i_2 < \cdots < i_n$. Let
+
+$$
+\pi=
+\left[
+\begin{array}{cccccc}
+1 & 2 & 3 & 4 & \cdots & 2n \\
+i_1 & j_1 & i_2 & j_2 & \cdots & j_{n}
+\end{array}
+\right]
+$$
+
+be a corresponding permutation. This depends only on the partition $\alpha$
+and not on the particular choice of $\Pi$. Given a partition $\alpha$ as above
+define
+
+$$A_\alpha =sgn(\pi)a_{i_1,j_1}a_{i_2,j_2}\cdots a_{i_n,j_n}$$
+
+The Pfaffian of $A$ is then given by
+
+$$Pf(A)=\sum_{\alpha\in\Pi} A_\alpha$$
+
+The Pfaffian of a $n\times n$ skew-symmetric matrix for n odd is
+defined to be zero.
+
+\subsection{Alternative definition}
+
+One can associate to any skew-symmetric $2n \times 2n$
+matrix $A=\{a_{ij}\}$ a bivector
+
+$$\omega=\sum_{i<j} a_{ij}~e^i\wedge e^j$$
+
+where $\{e^1, e^2, \ldots, e^{2n}\}$ is the standard basis of
+$\mathbb{R}^{2n}$. The Pfaffian is then defined by the equation
+
+$$\frac{1}{n!}\omega^n = \mbox{Pf}(A)~e^1\wedge e^2\wedge\cdots\wedge e^{2n}$$
+
+here $\omega^n$ denotes the wedge product of $n$ copies of
+$\omega^n$ with itself.
+
+\subsection{Derivation from Determinant}
+
+The Pfaffian can be derived from the determinant for a skew-symmetric
+matrix $A$ as follows. Using Laplace's formula we can write the
+determinant as
+
+$$\det(A)=(-1)^{p+1}a_{p1}A_{p1} +
+(-1)^{p+2}a_{p2}A_{p2}+\cdots+(-1)^{n+p}a_{pn}A_{pn}$$
+
+where $A_{pi}$ is the $p,i$ minor matrix of the matrix $A$. We further
+use Laplace's formula to note that
+
+$$\det(A[A_{ij}]) = \vert A \vert^n$$
+
+since this determinant is that of an $n \times n$ matrix whose only
+non-zero elements are the diagonals (each with value det(A)) and
+$[A_{ij}]$ is a matrix whose $i,j$th component is the corresponding $i,j$
+minor matrix. In this way, following a proof by Parameswaran, we can
+write the determinant as,
+
+$$\det(A)\equiv\Delta_n=
+\left|
+\begin{array}{cccc}
+a_{11}&a_{12}&\cdots&a_{1n}\\
+a_{21}&a_{22}&\cdots&a_{2n}\\
+\vdots&&&\vdots\\
+a_{n1}&a_{n2}&\cdots&a_{nn}
+\end{array}
+\right|
+$$
+
+The minor of
+$$
+\left|
+\begin{array}{cc}
+a_{11}&a_{12}\\
+a_{21}&a_{22}
+\end{array}
+\right|
+$$
+would be $\Delta_{n - 2}$. With this notation
+
+$$\left|
+\begin{array}{cccc}
+1&0&\cdots&0\\
+0&1&\cdots&0\\
+a_{31}&a_{32}&\cdots&a_{3n}\\
+\cdots&\cdots&\cdots&\cdots\\
+a_{n1}&a_{n2}&\cdots&a_{nn}
+\end{array}
+\right|
+\times
+\left|
+\begin{array}{cccc}
+A_{11}&A_{21}&\cdots&A_{n1}\\
+A_{12}&A_{22}&\cdots&A_{n2}\\
+A_{13}&A_{23}&\cdots&A_{n3}\\
+\cdots&\cdots&\cdots&\cdots\\
+A_{1n}&A_{2n}&\cdots&A_{nn}
+\end{array}
+\right| =
+\left|
+\begin{array}{cc}
+A_{11}&A_{21}\\
+A_{12}&a_{22}
+\end{array}
+\right|
+\Delta_n^{n-2}
+$$
+
+Thus
+
+$$\Delta_{n-2}\Delta_n^{n-1}=
+\left|
+\begin{array}{cc}
+A_{11}&A_{21}\\
+A_{12}&a_{22}
+\end{array}
+\right|
+\Delta_n^{n-2}
+$$
+
+Of course, it was only arbitrarily that we chose to remove the first
+two rows, and more generically we can write
+
+$$A_{rr}A_{ss} - A_{rs}A_{sr} = \Delta_{rs,rs}\Delta_n$$
+
+where $\Delta_{rs,rs}$ is the determinant of the original matrix with
+the rows $r$ and $s$, as well as the columns $r$ and $s$ removed. The
+equation above simplifies in the skew-symmetric case to
+
+$$A_{rs} = \sqrt{\Delta_{rs,rs}\Delta_n}$$
+
+We now plug this back into the original formula for the determinant,
+
+$$\Delta_n=
+(-1)^{p+1}a_{p1}\sqrt{\Delta_{p1,p1}\Delta_n} +
+(-1)^{p+2}a_{p2}\sqrt{\Delta_{p2,p2}\Delta_n} +
+\cdots +
+(-1)^{n+p}a_{pn}\sqrt{\Delta_{pn,pn}\Delta_n}
+$$
+
+or with slight manipulation,
+
+$$\sqrt{\Delta_n}=(-1)^{p+1}
+\left( \ a_{p1}\sqrt{\Delta_{p1,p1}} -
+a_{p2}\sqrt{\Delta_{p2,p2}} +
+\cdots +
+(-1)^{n-1}a_{pn}\sqrt{\Delta_{pn,pn}} \
+\right)
+$$
+
+The determinant is thus the square of the right hand side, and so we
+identify the right hand side as the Pfaffian.
+
+\section{Identities}
+
+For a $2n \times 2n$ skew-symmetric matrix $A$ and an arbitrary
+$2n \times 2n$ matrix B,
+\begin{itemize}
+\item $Pf(A)^2 = \det(A)$
+\item $Pf(BAB^T) = \det(B)Pf(A)$
+\item $Pf(\lambda{}A) = \lambda^nPf(A)$
+\item $Pf(A^T) = ( - 1)^nPf(A)$
+\item For a block-diagonal matrix
+$$A_1\oplus A_2=
+\left[
+\begin{array}{cc}
+A_1 & 0 \\
+0 & A_2
+\end{array}
+\right]
+$$
+$$Pf(A1\oplus A2) = Pf(A_11)Pf(A_2)$$
+\item For an arbitrary $n \times n$ matrix $M$:
+$$\mbox{Pf}
+\left[
+\begin{array}{cc}
+0 & M \\
+-M^T & 0
+\end{array}
+\right] =
+(-1)^{n(n-1)/2}\det M
+$$
+\end{itemize}
+\section{Applications}
+
+The Pfaffian is an invariant polynomial of a skew-symmetric matrix
+(note that it is not invariant under a general change of basis but
+rather under a proper orthogonal transformation). As such, it is
+important in the theory of characteristic classes. In particular, it
+can be used to define the Euler class of a Riemannian manifold which
+is used in the generalized Gauss-Bonnet theorem.
+
+The number of perfect matchings in a planar graph turns out to be the
+absolute value of a Pfaffian, hence is polynomial time
+computable. This is surprising given that for a general graph, the
+problem is very difficult (so called \#P-complete). This result is used
+to calculate the partition function of Ising models of spin glasses in
+physics, respectively of Markov random fields in machine learning
+(Globerson and Jaakkola, 2007), where the underlying graph is
+planar. Recently it is also used to derive efficient algorithms for
+some otherwise seemingly intractable problems, including the efficient
+simulation of certain types of restricted quantum computation.
+
+The calculation of the number of possible ways to tile a standard
+chessboard or 8-by-8 checkerboard with 32 dominoes is a simple example
+of a problem which may be solved through the use of the Pfaffian
+technique. There are 12,988,816 possible ways to tile a chessboard in
+this manner. Specifically, 12988816 is the number of possible ways to
+cover an 8-by-8 square with 32 1-by-2 rectangles. 12988816 is a square
+number: $12988816 = 3604^2$). Note that 12988816 can be written in the
+form: $2\times 1802^2 + 2\times 1802^2$,
+where all the numbers have a digital root of 2.
+
+More generally, the number of ways to cover a $2n \times 2n$ square with
+$2n^2$ dominoes (as calculated independently by Temperley and M.E. Fisher and
+Kasteleyn in 1961) is given by
+
+$$
+\prod_{j=1}^N
+\prod_{k=1}^N
+\left ( 4\cos^2 \frac{\pi j}{2n + 1} +
+4\cos^2 \frac{\pi k}{2n + 1} \right )
+$$
+
+This technique can be applied in many mathematics-related subjects,
+for example, in the classical, 2-dimensional computation of the
+dimer-dimer correlator function in quantum mechanics.
+
+\section{History}
+
+The term Pfaffian was introduced by Arthur Cayley, who used the term
+in 1852: ``The permutants of this class (from their connection with the
+researches of Pfaff on differential equations) I shall term
+Pfaffians.'' The term honors German mathematician Johann Friedrich
+Pfaff.
+
+\section{Axiom code}
+I have hacked together an algorithm to compute a Pfaffian, using an algorithm
+of Gunter Rote. Currently it's only an .input script, but if it's useful for
+somebody else than myself, we could make it a little more professional.
+
+Martin
+
+<<*>>=
+)spool pfaffian.output
+)set message test on
+)set message auto off
+)clear all
+
+--S 1 of 12
+B0 n == matrix [[(if i=j+1 and odd? j then -1 else _
+ if i=j-1 and odd? i then 1 else 0) _
+ for j in 1..n] for i in 1..n]
+--R
+--R Type: Void
+--E 1
+
+--S 2 of 12
+PfChar(lambda, A) ==
+ n := nrows A
+ (n = 2) => lambda^2 + A.(1,2)
+ M := subMatrix(A, 3, n, 3, n)
+ r := subMatrix(A, 1, 1, 3, n)
+ s := subMatrix(A, 3, n, 2, 2)
+
+ p := PfChar(lambda, M)
+ d := degree(p, lambda)
+
+ B := B0(n-2)
+ C := r*B
+ g := [(C*s).(1,1), A.(1,2), 1]
+ if d >= 4 then
+ B := M*B
+ for i in 4..d by 2 repeat
+ C := C*B
+ g := cons((C*s).(1,1), g)
+ g := reverse! g
+
+ res := 0
+ for i in 0..d by 2 for j in 2..d+2 repeat
+ c := coefficient(p, lambda, i)
+ for e in first(g, j) for k in 2..-d by -2 repeat
+ res := res + c * e * lambda^(k+i)
+
+ res
+--R
+--R Type: Void
+--E 2
+
+--S 3 of 12
+pfaffian A == eval(PfChar(l, A), l=0)
+--R
+--R Type: Void
+--E 3
+
+--S 4 of 12
+m:Matrix(Integer):=[[0,15],[-15,0]]
+--R
+--R + 0 15+
+--R (4) | |
+--R +- 15 0 +
+--R Type: Matrix Integer
+--E 4
+
+--S 5 of 12
+pfaffian m
+--R
+--R (5) 15
+--R Type: Polynomial Integer
+--E 5
+
+--S 6 of 12
+m1:Matrix(Polynomial(Integer)):=[[0,a,b,c],[-a,0,d,e],[-b,-d,0,f],[-c,-e,-f,0]]
+--R
+--R
+--R + 0 a b c+
+--R | |
+--R |- a 0 d e|
+--R (6) | |
+--R |- b - d 0 f|
+--R | |
+--R +- c - e - f 0+
+--R Type: Matrix Polynomial Integer
+--E 6
+
+--S 7 of 12
+pfaffian m1
+--R
+--R Compiling function B0 with type PositiveInteger -> Matrix Integer
+--R
+--R (7) a f - b e + c d
+--R Type: Polynomial Integer
+--E 7
+
+--S 8 of 12
+(a,b,c,d,e,f):=(3,5,7,11,13,17)
+--R
+--R (8) 17
+--R Type: PositiveInteger
+--E 8
+
+--S 9 of 12
+m1
+--R
+--R
+--R + 0 a b c+
+--R | |
+--R |- a 0 d e|
+--R (9) | |
+--R |- b - d 0 f|
+--R | |
+--R +- c - e - f 0+
+--R Type: Matrix Polynomial Integer
+--E 9
+
+--S 10 of 12
+a*f-b*e+d*c
+--R
+--R
+--R (10) 63
+--R Type: PositiveInteger
+--E 10
+
+--S 11 of 12
+n:=pfaffian m1
+--R
+--R (11) a f - b e + c d
+--R Type: Polynomial Integer
+--E 9
+
+--S 12 of 12
+eval(n,['a,'b,'c,'d,'e,'f]::List(Symbol),[a,b,c,d,e,f])
+--R
+--R (12) 63
+--R Type: Polynomial Integer
+--E 12
+)spool
+)lisp (bye)
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} {\bf ``http://en.wikipedia.org/wiki/Pfaffian''}
+\bibitem{2} ``The statistics of dimers on a lattice, Part I'', Physica,
+vol.27, 1961, pp.1209-25, P.W. Kasteleyn.
+\bibitem{3} Propp, James (2004),
+``Lambda-determinants and domino-tilings'', arXiv:math.CO/0406301.
+\bibitem{4} Globerson, Amir and Tommi Jaakkola (2007),
+``Approximate inference using planar graph decomposition'',
+Advances in Neural Information Processing Systems 19, MIT Press.
+\bibitem{5} ``The Games and Puzzles Journal'',
+No.14, 1996, pp.204-5, Robin J. Chapman, University of Exeter
+\bibitem{6} ``Domino Tilings and Products of Fibonacci and Pell numbers'',
+Journal of Integer Sequences, Vol.5, 2002, Article 02.1.2,
+James A. Sellers, The Pennsylvania State University
+\bibitem{7} ``The Penguin Dictionary of Curious and Interesting Numbers'',
+revised ed., 1997, ISBN 0-14-026149-4, David Wells, p.182.
+\bibitem{8} ``A Treatise on the Theory of Determinants'',
+1882, Macmillan and Co., Thomas Muir, Online
+\bibitem{9} ``Skew-Symmetric Determinants'',
+The American Mathematical Monthly, vol. 61, no.2., 1954, p.116,
+S. Parameswaran Online-Subscription
+\end{thebibliography}
+\end{document}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <gi...@ax...> - 2007-09-21 15:55:58
|
changelog | 3 ++
src/algebra/laplace.spad.pamphlet | 50 +++++++++++++++++++++++++-----
src/input/Makefile.pamphlet | 9 ++++--
src/input/bug101.input.pamphlet | 59 +++++++++++++++++++++++++++++++++++++
4 files changed, 109 insertions(+), 12 deletions(-)
New commits:
commit 541a0f407052e4d11418747a7f87d5d456dffc85
Author: Tim Daly <root@localhost.localdomain>
Date: Tue Sep 18 07:39:30 2007 -0400
20070920 tpd src/input/Makefile add bug101.input regression test
20070920 tpd src/input/bug101.input test laplace(log(z),z,w) bug 101
20070920 wxh src/algebra/laplace.spad fix laplace(log(z),z,w) bug 101
|
|
From: <da...@us...> - 2007-09-21 15:55:48
|
Revision: 750
http://axiom.svn.sourceforge.net/axiom/?rev=750&view=rev
Author: daly
Date: 2007-09-21 08:55:49 -0700 (Fri, 21 Sep 2007)
Log Message:
-----------
20070920 tpd src/input/Makefile add bug101.input regression test
20070920 tpd src/input/bug101.input test laplace(log(z),z,w) bug 101
20070920 wxh src/algebra/laplace.spad fix laplace(log(z),z,w) bug 101
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/algebra/laplace.spad.pamphlet
trunk/axiom/src/input/Makefile.pamphlet
Added Paths:
-----------
trunk/axiom/src/input/bug101.input.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-16 22:49:24 UTC (rev 749)
+++ trunk/axiom/changelog 2007-09-21 15:55:49 UTC (rev 750)
@@ -1,3 +1,6 @@
+20070920 tpd src/input/Makefile add bug101.input regression test
+20070920 tpd src/input/bug101.input test laplace(log(z),z,w) bug 101
+20070920 wxh src/algebra/laplace.spad fix laplace(log(z),z,w) bug 101
20070916 tpd src/input/Makefile add bug103.input regression test
20070916 tpd src/input/bug103.input test solve(z=z,z) bug fix
20070916 tpd src/algebra/polycat.spad solve(z=z,z) bug fix
Modified: trunk/axiom/src/algebra/laplace.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/laplace.spad.pamphlet 2007-09-16 22:49:24 UTC (rev 749)
+++ trunk/axiom/src/algebra/laplace.spad.pamphlet 2007-09-21 15:55:49 UTC (rev 750)
@@ -161,41 +161,68 @@
lapkernel(f, t, tt, ss) ==
(k := retractIfCan(f)@Union(K, "failed")) case "failed" => "failed"
- empty?(arg := argument(k::K)) or not empty? rest arg => "failed"
+ empty?(arg := argument(k::K)) => "failed"
+ is?(op := operator k, "%diff"::SE) =>
+ not( #arg = 3) => "failed"
+ not(is?(arg.3, t)) => "failed"
+ fint := eval(arg.1, arg.2, tt)
+ s := name operator (kernels(ss).1)
+ ss * locallaplace(fint, t, tt, s, ss) - eval(fint, tt = 0)
+ not (empty?(rest arg)) => "failed"
member?(t, variables(a := first(arg) / tt)) => "failed"
is?(op := operator k, "Si"::SE) => atan(a / ss) / ss
is?(op, "Ci"::SE) => log((ss**2 + a**2) / a**2) / (2 * ss)
is?(op, "Ei"::SE) => log((ss + a) / a) / ss
+ -- digamma (or Gamma) needs SpecialFunctionCategory
+ -- which we do not have here
+ -- is?(op, "log"::SE) => (digamma(1) - log(a) - log(ss)) / ss
"failed"
+ -- Below we try to apply one of the texbook rules for computing
+ -- Laplace transforms, either reducing problem to simpler cases
+ -- or using one of known base cases
locallaplace(f, t, tt, s, ss) ==
zero? f => 0
-- one? f => inv ss
(f = 1) => inv ss
+
+ -- laplace(f(t)/t,t,s)
+ -- = integrate(laplace(f(t),t,v), v = s..%plusInfinity)
(x := tdenom(f, tt)) case F =>
g := locallaplace(x::F, t, tt, vv := new()$SE, vvv := vv::F)
(x := intlaplace(f, ss, g, vv, vvv)) case F => x::F
oplap(f, tt, ss)
+
+ -- Use linearity
(u := mkPlus f) case List(F) =>
+/[locallaplace(g, t, tt, s, ss) for g in u::List(F)]
(rec := splitConstant(f, t)).const ^= 1 =>
rec.const * locallaplace(rec.nconst, t, tt, s, ss)
+
+ -- laplace(t^n*f(t),t,s) = (-1)^n*D(laplace(f(t),t,s), s, n))
(v := atn(f, t)) case Record(coef:F, deg:PI) =>
vv := v::Record(coef:F, deg:PI)
is?(la := locallaplace(vv.coef, t, tt, s, ss), oplap) => oplap(f,tt,ss)
(-1$Integer)**(vv.deg) * differentiate(la, s, vv.deg)
+
+ -- Complex shift rule
(w := aexp(f, t)) case Record(coef:F, coef1:F, coef0:F) =>
ww := w::Record(coef:F, coef1:F, coef0:F)
exp(ww.coef0) * locallaplace(ww.coef,t,tt,s,ss - ww.coef1)
+
+ -- Try base cases
(x := lapkernel(f, t, tt, ss)) case F => x::F
- -- last chance option: try to use the fact that
- -- laplace(f(t),t,s) = s laplace(g(t),t,s) - g(0) where dg/dt = f(t)
- elem?(int := lfintegrate(f, t)) and (rint := retractIfCan int) case F =>
- fint := rint :: F
- -- to avoid infinite loops, we don't call laplace recursively
- -- if the integral has no new logs and f is an algebraic function
- empty?(logpart int) and algebraic?(f, t) => oplap(fint, tt, ss)
- ss * locallaplace(fint, t, tt, s, ss) - eval(fint, tt = 0)
+
+-- -- The following does not seem to help computing transforms, but
+-- -- quite frequently leads to loops, so I (wh) disabled it for now
+-- -- last chance option: try to use the fact that
+-- -- laplace(f(t),t,s) = s laplace(g(t),t,s) - g(0) where dg/dt = f(t)
+-- elem?(int := lfintegrate(f, t)) and (rint := retractIfCan int) case F =>
+-- fint := rint :: F
+-- -- to avoid infinite loops, we don't call laplace recursively
+-- -- if the integral has no new logs and f is an algebraic function
+-- empty?(logpart int) and algebraic?(f, t) => oplap(fint, tt, ss)
+-- ss * locallaplace(fint, t, tt, s, ss) - eval(fint, tt = 0)
oplap(f, tt, ss)
setProperty(oplap,SPECIALDIFF,dvlap@((List F,SE)->F) pretend None)
@@ -228,6 +255,7 @@
++ Laplace transform of \spad{f(s)}
++ using t as the new variable or "failed" if unable to find
++ a closed form.
+ ++ Handles only rational \spad{f(s)}.
Implementation ==> add
@@ -246,8 +274,12 @@
ilt(expr,var,t) ==
expr = 0 => 0
r := univariate(expr,kernel(var))
+
+ -- Check that r is a rational function such that degree of
+ -- the numarator is lower that degree of denominator
not(numer(r) quo denom(r) = 0) => "failed"
not( freeOf?(numer r,var) and freeOf?(denom r,var)) => "failed"
+
ilt1(r,t::F)
hintpac := TranscendentalHermiteIntegration(F, UP)
Modified: trunk/axiom/src/input/Makefile.pamphlet
===================================================================
--- trunk/axiom/src/input/Makefile.pamphlet 2007-09-16 22:49:24 UTC (rev 749)
+++ trunk/axiom/src/input/Makefile.pamphlet 2007-09-21 15:55:49 UTC (rev 750)
@@ -287,7 +287,8 @@
arrows.regress assign.regress atansqrt.regress \
asec.regress bags.regress bbtree.regress \
binary.regress bop.regress bstree.regress bouquet.regress \
- bug100.regress bug103.regress bug10069.regress \
+ bug100.regress bug101.regress \
+ bug103.regress bug10069.regress \
bugs.regress bug10312.regress bug6357.regress bug9057.regress \
calcprob.regress \
calculus2.regress calculus.regress cardinal.regress card.regress \
@@ -502,7 +503,8 @@
${OUT}/bags.input ${OUT}/bbtree.input ${OUT}/bern.input \
${OUT}/bernpoly.input ${OUT}/binary.input ${OUT}/bop.input \
${OUT}/bouquet.input ${OUT}/bstree.input ${OUT}/bug6357.input \
- ${OUT}/bug9057.input ${OUT}/bug100.input ${OUT}/bug103.input \
+ ${OUT}/bug9057.input ${OUT}/bug100.input ${OUT}/bug101.input \
+ ${OUT}/bug103.input \
${OUT}/bug10069.input ${OUT}/bug10312.input \
${OUT}/calcprob.input ${OUT}/calculus.input \
${OUT}/cardinal.input ${OUT}/card.input ${OUT}/carten.input \
@@ -678,7 +680,8 @@
${DOC}/bernpoly.input.dvi ${DOC}/binary.input.dvi \
${DOC}/bop.input.dvi ${DOC}/bouquet.input.dvi \
${DOC}/bstree.input.dvi ${DOC}/bug10069.input.dvi \
- ${DOC}/bug100.input.dvi ${DOC}/bug103.input.dvi \
+ ${DOC}/bug100.input.dvi ${DOC}/bug101.input.dvi \
+ ${DOC}/bug103.input.dvi \
${DOC}/bug10312.input.dvi ${DOC}/bug6357.input.dvi \
${DOC}/bug9057.input.dvi ${DOC}/bugs.input.dvi \
${DOC}/c02aff.input.dvi ${DOC}/c02agf.input.dvi \
Added: trunk/axiom/src/input/bug101.input.pamphlet
===================================================================
--- trunk/axiom/src/input/bug101.input.pamphlet (rev 0)
+++ trunk/axiom/src/input/bug101.input.pamphlet 2007-09-21 15:55:49 UTC (rev 750)
@@ -0,0 +1,59 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/input bug101.input}
+\author{Timothy Daly}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+@
+The call
+\begin{verbatim}
+laplace(log(z),z,w)
+\end{verbatim}
+causes an infinite loop between lfintegrate and monomialIntegrate.
+
+Waldek modified the code to fail and return unevaluated.
+
+A more reasonable return value, not supported by Axiom, would be:
+\begin{verbatim}
+ -log(w)-gamma
+ -------------
+ w
+\end{verbatim}
+
+<<*>>=
+)spool bug101.output
+)set message test on
+)set message auto off
+)clear all
+
+--S 1 of 2
+laplace(log(z),z,w)
+--R
+--R (1) laplace(log(z),z,w)
+--R Type: Expression Integer
+--E 1
+
+--S 2 of 2
+laplace(log(z),w,z)
+--R
+--R log(z)
+--R (2) ------
+--R z
+--R Type: Expression Integer
+--E 2
+@
+<<*>>=
+)spool
+)lisp (bye)
+
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <gi...@ax...> - 2007-09-16 22:49:32
|
changelog | 4 +
src/algebra/catdef.spad.pamphlet | 1001 +++++------
src/algebra/polycat.spad.pamphlet | 3641 ++++++++++++++++---------------------
src/input/Makefile.pamphlet | 6 +-
src/input/bug103.input.pamphlet | 82 +
5 files changed, 2140 insertions(+), 2594 deletions(-)
New commits:
commit 0d25893fa3d2f6ae66a9c12a00464f88e44bf294
Author: Tim Daly <root@localhost.localdomain>
Date: Fri Sep 14 15:43:32 2007 -0400
20070916 tpd src/input/Makefile add bug103.input regression test
commit 8089767155e0288a92013d6c25ea0f95adb96cea
Author: Tim Daly <root@localhost.localdomain>
Date: Fri Sep 14 15:41:15 2007 -0400
20070916 tpd src/input/Makefile add bug103.input regression test
20070916 tpd src/input/bug103.input test solve(z=z,z) bug fix
20070916 tpd src/algebra/polycat.spad solve(z=z,z) bug fix
20070916 tpd src/algebra/catdef.spad add zero? to exquo
|
|
From: <da...@us...> - 2007-09-16 22:49:21
|
Revision: 749
http://axiom.svn.sourceforge.net/axiom/?rev=749&view=rev
Author: daly
Date: 2007-09-16 15:49:24 -0700 (Sun, 16 Sep 2007)
Log Message:
-----------
20070916 tpd src/input/Makefile add bug103.input regression test
20070916 tpd src/input/bug103.input test solve(z=z,z) bug fix
20070916 tpd src/algebra/polycat.spad solve(z=z,z) bug fix
20070916 tpd src/algebra/catdef.spad add zero? to exquo
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/algebra/catdef.spad.pamphlet
trunk/axiom/src/algebra/polycat.spad.pamphlet
trunk/axiom/src/input/Makefile.pamphlet
Added Paths:
-----------
trunk/axiom/src/input/bug103.input.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-16 03:12:14 UTC (rev 748)
+++ trunk/axiom/changelog 2007-09-16 22:49:24 UTC (rev 749)
@@ -1,3 +1,7 @@
+20070916 tpd src/input/Makefile add bug103.input regression test
+20070916 tpd src/input/bug103.input test solve(z=z,z) bug fix
+20070916 tpd src/algebra/polycat.spad solve(z=z,z) bug fix
+20070916 tpd src/algebra/catdef.spad add zero? to exquo
20070915 tpd merge bug100 branch
20070915 tpd src/input/Makefile add bug100.input regression test
20070915 tpd src/input/bug100.input test integrate((z^a+1)^b,z) infinite loop
Modified: trunk/axiom/src/algebra/catdef.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/catdef.spad.pamphlet 2007-09-16 03:12:14 UTC (rev 748)
+++ trunk/axiom/src/algebra/catdef.spad.pamphlet 2007-09-16 22:49:24 UTC (rev 749)
@@ -1319,6 +1319,7 @@
x quo y == divide(x,y).quotient --divide must be user-supplied
x rem y == divide(x,y).remainder
x exquo y ==
+ zero? x => 0
zero? y => "failed"
qr:=divide(x,y)
zero?(qr.remainder) => qr.quotient
@@ -1544,601 +1545,507 @@
@
\subsubsection{EUCDOM-;sizeLess?;2SB;1}
<<EUCDOM-;sizeLess?;2SB;1>>=
-(DEFUN |EUCDOM-;sizeLess?;2SB;1| (|x| |y| |$|)
- (COND
- ((SPADCALL |y| (QREFELT |$| 8)) (QUOTE NIL))
- ((SPADCALL |x| (QREFELT |$| 8)) (QUOTE T))
- ((QUOTE T)
- (|<| (SPADCALL |x| (QREFELT |$| 10)) (SPADCALL |y| (QREFELT |$| 10))))))
+(DEFUN |EUCDOM-;sizeLess?;2SB;1| (|x| |y| $)
+ (COND
+ ((SPADCALL |y| (QREFELT $ 8)) (QUOTE NIL))
+ ((SPADCALL |x| (QREFELT $ 8)) (QUOTE T))
+ ((QUOTE T)
+ (< (SPADCALL |x| (QREFELT $ 10)) (SPADCALL |y| (QREFELT $ 10))))))
@
\subsubsection{EUCDOM-;quo;3S;2}
<<EUCDOM-;quo;3S;2>>=
-(DEFUN |EUCDOM-;quo;3S;2| (|x| |y| |$|)
- (QCAR (SPADCALL |x| |y| (QREFELT |$| 13))))
+(DEFUN |EUCDOM-;quo;3S;2| (|x| |y| $)
+ (QCAR (SPADCALL |x| |y| (QREFELT $ 13))))
@
\subsubsection{EUCDOM-;rem;3S;3}
<<EUCDOM-;rem;3S;3>>=
-(DEFUN |EUCDOM-;rem;3S;3| (|x| |y| |$|)
- (QCDR (SPADCALL |x| |y| (QREFELT |$| 13))))
+(DEFUN |EUCDOM-;rem;3S;3| (|x| |y| $)
+ (QCDR (SPADCALL |x| |y| (QREFELT $ 13))))
@
\subsubsection{EUCDOM-;exquo;2SU;4}
<<EUCDOM-;exquo;2SU;4>>=
-(DEFUN |EUCDOM-;exquo;2SU;4| (|x| |y| |$|)
- (PROG (|qr|)
- (RETURN
- (SEQ
- (COND
- ((SPADCALL |y| (QREFELT |$| 8)) (CONS 1 "failed"))
- ((QUOTE T)
- (SEQ
- (LETT |qr|
- (SPADCALL |x| |y| (QREFELT |$| 13))
- |EUCDOM-;exquo;2SU;4|)
- (EXIT
- (COND
- ((SPADCALL (QCDR |qr|) (QREFELT |$| 8)) (CONS 0 (QCAR |qr|)))
- ((QUOTE T) (CONS 1 "failed")))))))))))
+(DEFUN |EUCDOM-;exquo;2SU;4| (|x| |y| $)
+ (PROG (|qr|)
+ (RETURN
+ (SEQ
+ (COND
+ ((SPADCALL |x| (QREFELT $ 8)) (CONS 0 (|spadConstant| $ 16)))
+ ((SPADCALL |y| (QREFELT $ 8)) (CONS 1 "failed"))
+ ((QUOTE T)
+ (SEQ
+ (LETT |qr|
+ (SPADCALL |x| |y| (QREFELT $ 13))
+ |EUCDOM-;exquo;2SU;4|)
+ (EXIT
+ (COND
+ ((SPADCALL (QCDR |qr|) (QREFELT $ 8)) (CONS 0 (QCAR |qr|)))
+ ((QUOTE T) (CONS 1 "failed")))))))))))
@
\subsubsection{EUCDOM-;gcd;3S;5}
<<EUCDOM-;gcd;3S;5>>=
-(DEFUN |EUCDOM-;gcd;3S;5| (|x| |y| |$|)
- (PROG (|#G13| |#G14|)
- (RETURN
- (SEQ
- (LETT |x| (SPADCALL |x| (QREFELT |$| 18)) |EUCDOM-;gcd;3S;5|)
- (LETT |y| (SPADCALL |y| (QREFELT |$| 18)) |EUCDOM-;gcd;3S;5|)
- (SEQ G190
- (COND
- ((NULL
- (COND
- ((SPADCALL |y| (QREFELT |$| 8)) (QUOTE NIL))
- ((QUOTE T) (QUOTE T))))
- (GO G191)))
- (SEQ
- (PROGN
- (LETT |#G13| |y| |EUCDOM-;gcd;3S;5|)
- (LETT |#G14| (SPADCALL |x| |y| (QREFELT |$| 19)) |EUCDOM-;gcd;3S;5|)
- (LETT |x| |#G13| |EUCDOM-;gcd;3S;5|)
- (LETT |y| |#G14| |EUCDOM-;gcd;3S;5|))
- (EXIT
- (LETT |y| (SPADCALL |y| (QREFELT |$| 18)) |EUCDOM-;gcd;3S;5|)))
- NIL
- (GO G190)
- G191
- (EXIT NIL))
- (EXIT |x|)))))
+(DEFUN |EUCDOM-;gcd;3S;5| (|x| |y| $)
+ (PROG (|#G13| |#G14|)
+ (RETURN
+ (SEQ
+ (LETT |x| (SPADCALL |x| (QREFELT $ 19)) |EUCDOM-;gcd;3S;5|)
+ (LETT |y| (SPADCALL |y| (QREFELT $ 19)) |EUCDOM-;gcd;3S;5|)
+ (SEQ G190
+ (COND
+ ((NULL
+ (COND
+ ((SPADCALL |y| (QREFELT $ 8)) (QUOTE NIL))
+ ((QUOTE T) (QUOTE T))))
+ (GO G191)))
+ (SEQ
+ (PROGN
+ (LETT |#G13| |y| |EUCDOM-;gcd;3S;5|)
+ (LETT |#G14| (SPADCALL |x| |y| (QREFELT $ 20)) |EUCDOM-;gcd;3S;5|)
+ (LETT |x| |#G13| |EUCDOM-;gcd;3S;5|)
+ (LETT |y| |#G14| |EUCDOM-;gcd;3S;5|))
+ (EXIT
+ (LETT |y| (SPADCALL |y| (QREFELT $ 19)) |EUCDOM-;gcd;3S;5|)))
+ NIL
+ (GO G190)
+ G191
+ (EXIT NIL))
+ (EXIT |x|)))))
@
\subsubsection{EUCDOM-;unitNormalizeIdealElt}
<<EUCDOM-;unitNormalizeIdealElt>>=
-(DEFUN |EUCDOM-;unitNormalizeIdealElt| (|s| |$|)
- (PROG (|#G16| |u| |c| |a|)
- (RETURN
- (SEQ
- (PROGN
- (LETT |#G16| (SPADCALL (QVELT |s| 2) (QREFELT |$| 22)) |EUCDOM-;unitNormalizeIdealElt|)
- (LETT |u| (QVELT |#G16| 0) |EUCDOM-;unitNormalizeIdealElt|)
- (LETT |c| (QVELT |#G16| 1) |EUCDOM-;unitNormalizeIdealElt|)
- (LETT |a| (QVELT |#G16| 2) |EUCDOM-;unitNormalizeIdealElt|)
- |#G16|)
- (EXIT
- (COND
- ((SPADCALL |a| (QREFELT |$| 23)) |s|)
- ((QUOTE T)
- (VECTOR
- (SPADCALL |a| (QVELT |s| 0) (QREFELT |$| 24))
- (SPADCALL |a| (QVELT |s| 1) (QREFELT |$| 24))
- |c|))))))))
+(DEFUN |EUCDOM-;unitNormalizeIdealElt| (|s| $)
+ (PROG (|#G16| |u| |c| |a|)
+ (RETURN
+ (SEQ
+ (PROGN
+ (LETT |#G16|
+ (SPADCALL (QVELT |s| 2) (QREFELT $ 23))
+ |EUCDOM-;unitNormalizeIdealElt|)
+ (LETT |u| (QVELT |#G16| 0) |EUCDOM-;unitNormalizeIdealElt|)
+ (LETT |c| (QVELT |#G16| 1) |EUCDOM-;unitNormalizeIdealElt|)
+ (LETT |a| (QVELT |#G16| 2) |EUCDOM-;unitNormalizeIdealElt|)
+ |#G16|)
+ (EXIT
+ (COND
+ ((SPADCALL |a| (|spadConstant| $ 24) (QREFELT $ 25)) |s|)
+ ((QUOTE T)
+ (VECTOR
+ (SPADCALL |a| (QVELT |s| 0) (QREFELT $ 26))
+ (SPADCALL |a| (QVELT |s| 1) (QREFELT $ 26))
+ |c|))))))))
@
\subsubsection{EUCDOM-;extendedEuclidean;2SR;7}
<<EUCDOM-;extendedEuclidean;2SR;7>>=
-(DEFUN |EUCDOM-;extendedEuclidean;2SR;7| (|x| |y| |$|)
- (PROG (|s3| |s2| |qr| |s1|)
- (RETURN
- (SEQ
- (LETT |s1|
- (|EUCDOM-;unitNormalizeIdealElt|
- (VECTOR (|spadConstant| |$| 25) (|spadConstant| |$| 26) |x|) |$|)
- |EUCDOM-;extendedEuclidean;2SR;7|)
- (LETT |s2|
- (|EUCDOM-;unitNormalizeIdealElt|
- (VECTOR (|spadConstant| |$| 26) (|spadConstant| |$| 25) |y|) |$|)
- |EUCDOM-;extendedEuclidean;2SR;7|)
- (EXIT
- (COND
- ((SPADCALL |y| (QREFELT |$| 8)) |s1|)
- ((SPADCALL |x| (QREFELT |$| 8)) |s2|)
- ((QUOTE T)
- (SEQ
- (SEQ G190
- (COND
- ((NULL
- (COND
- ((SPADCALL (QVELT |s2| 2) (QREFELT |$| 8))
- (QUOTE NIL))
- ((QUOTE T) (QUOTE T))))
- (GO G191)))
- (SEQ
- (LETT |qr|
- (SPADCALL (QVELT |s1| 2) (QVELT |s2| 2) (QREFELT |$| 13))
- |EUCDOM-;extendedEuclidean;2SR;7|)
- (LETT |s3|
- (VECTOR
- (SPADCALL
- (QVELT |s1| 0)
- (SPADCALL
- (QCAR |qr|)
- (QVELT |s2| 0)
- (QREFELT |$| 24))
- (QREFELT |$| 27))
- (SPADCALL
- (QVELT |s1| 1)
- (SPADCALL
- (QCAR |qr|)
- (QVELT |s2| 1)
- (QREFELT |$| 24))
- (QREFELT |$| 27))
- (QCDR |qr|))
- |EUCDOM-;extendedEuclidean;2SR;7|)
- (LETT |s1| |s2| |EUCDOM-;extendedEuclidean;2SR;7|)
- (EXIT
- (LETT |s2|
- (|EUCDOM-;unitNormalizeIdealElt| |s3| |$|)
- |EUCDOM-;extendedEuclidean;2SR;7|)))
- NIL
- (GO G190)
- G191
- (EXIT NIL))
- (COND
- ((NULL (SPADCALL (QVELT |s1| 0) (QREFELT |$| 8)))
- (COND
- ((NULL (SPADCALL (QVELT |s1| 0) |y| (QREFELT |$| 28)))
- (SEQ
- (LETT |qr|
- (SPADCALL (QVELT |s1| 0) |y| (QREFELT |$| 13))
- |EUCDOM-;extendedEuclidean;2SR;7|)
- (QSETVELT |s1| 0 (QCDR |qr|))
- (QSETVELT |s1| 1
- (SPADCALL
- (QVELT |s1| 1)
- (SPADCALL (QCAR |qr|) |x| (QREFELT |$| 24))
- (QREFELT |$| 29)))
- (EXIT
- (LETT |s1|
- (|EUCDOM-;unitNormalizeIdealElt| |s1| |$|)
- |EUCDOM-;extendedEuclidean;2SR;7|)))))))
- (EXIT |s1|)))))))))
+(DEFUN |EUCDOM-;extendedEuclidean;2SR;7| (|x| |y| $)
+ (PROG (|s3| |s2| |qr| |s1|)
+ (RETURN
+ (SEQ
+ (LETT |s1|
+ (|EUCDOM-;unitNormalizeIdealElt|
+ (VECTOR (|spadConstant| $ 24) (|spadConstant| $ 16) |x|)
+ $)
+ |EUCDOM-;extendedEuclidean;2SR;7|)
+ (LETT |s2|
+ (|EUCDOM-;unitNormalizeIdealElt|
+ (VECTOR (|spadConstant| $ 16) (|spadConstant| $ 24) |y|)
+ $)
+ |EUCDOM-;extendedEuclidean;2SR;7|)
+ (EXIT
+ (COND
+ ((SPADCALL |y| (QREFELT $ 8)) |s1|)
+ ((SPADCALL |x| (QREFELT $ 8)) |s2|)
+ ((QUOTE T)
+ (SEQ
+ (SEQ
+ G190
+ (COND
+ ((NULL
+ (COND
+ ((SPADCALL (QVELT |s2| 2) (QREFELT $ 8)) (QUOTE NIL))
+ ((QUOTE T) (QUOTE T))))
+ (GO G191)))
+ (SEQ
+ (LETT |qr|
+ (SPADCALL (QVELT |s1| 2) (QVELT |s2| 2) (QREFELT $ 13))
+ |EUCDOM-;extendedEuclidean;2SR;7|)
+ (LETT |s3|
+ (VECTOR
+ (SPADCALL (QVELT |s1| 0)
+ (SPADCALL (QCAR |qr|) (QVELT |s2| 0) (QREFELT $ 26))
+ (QREFELT $ 27))
+ (SPADCALL (QVELT |s1| 1)
+ (SPADCALL (QCAR |qr|) (QVELT |s2| 1) (QREFELT $ 26))
+ (QREFELT $ 27))
+ (QCDR |qr|))
+ |EUCDOM-;extendedEuclidean;2SR;7|)
+ (LETT |s1| |s2| |EUCDOM-;extendedEuclidean;2SR;7|)
+ (EXIT
+ (LETT |s2|
+ (|EUCDOM-;unitNormalizeIdealElt| |s3| $)
+ |EUCDOM-;extendedEuclidean;2SR;7|)))
+ NIL
+ (GO G190)
+ G191
+ (EXIT NIL))
+ (COND
+ ((NULL (SPADCALL (QVELT |s1| 0) (QREFELT $ 8)))
+ (COND
+ ((NULL (SPADCALL (QVELT |s1| 0) |y| (QREFELT $ 28)))
+ (SEQ
+ (LETT |qr|
+ (SPADCALL (QVELT |s1| 0) |y| (QREFELT $ 13))
+ |EUCDOM-;extendedEuclidean;2SR;7|)
+ (QSETVELT |s1| 0 (QCDR |qr|))
+ (QSETVELT |s1| 1
+ (SPADCALL (QVELT |s1| 1)
+ (SPADCALL (QCAR |qr|) |x| (QREFELT $ 26)) (QREFELT $ 29)))
+ (EXIT
+ (LETT |s1|
+ (|EUCDOM-;unitNormalizeIdealElt| |s1| $)
+ |EUCDOM-;extendedEuclidean;2SR;7|)))))))
+ (EXIT |s1|)))))))))
@
\subsubsection{EUCDOM-;extendedEuclidean;3SU;8}
<<EUCDOM-;extendedEuclidean;3SU;8>>=
-(DEFUN |EUCDOM-;extendedEuclidean;3SU;8| (|x| |y| |z| |$|)
- (PROG (|s| |w| |qr|)
- (RETURN
- (SEQ
- (COND
- ((SPADCALL |z| (QREFELT |$| 8))
- (CONS 0 (CONS (|spadConstant| |$| 26) (|spadConstant| |$| 26))))
- ((QUOTE T)
- (SEQ
- (LETT |s|
- (SPADCALL |x| |y| (QREFELT |$| 32))
- |EUCDOM-;extendedEuclidean;3SU;8|)
- (LETT |w|
- (SPADCALL |z| (QVELT |s| 2) (QREFELT |$| 33))
- |EUCDOM-;extendedEuclidean;3SU;8|)
- (EXIT
- (COND
- ((QEQCAR |w| 1) (CONS 1 "failed"))
- ((SPADCALL |y| (QREFELT |$| 8))
- (CONS 0
- (CONS
- (SPADCALL (QVELT |s| 0) (QCDR |w|) (QREFELT |$| 24))
- (SPADCALL (QVELT |s| 1) (QCDR |w|) (QREFELT |$| 24)))))
- ((QUOTE T)
- (SEQ
- (LETT |qr|
- (SPADCALL
- (SPADCALL (QVELT |s| 0) (QCDR |w|) (QREFELT |$| 24))
- |y|
- (QREFELT |$| 13))
- |EUCDOM-;extendedEuclidean;3SU;8|)
- (EXIT
- (CONS
- 0
- (CONS
- (QCDR |qr|)
- (SPADCALL
- (SPADCALL
- (QVELT |s| 1)
- (QCDR |w|)
- (QREFELT |$| 24))
- (SPADCALL
- (QCAR |qr|)
- |x|
- (QREFELT |$| 24))
- (QREFELT |$| 29))))))))))))))))
+(DEFUN |EUCDOM-;extendedEuclidean;3SU;8| (|x| |y| |z| $)
+ (PROG (|s| |w| |qr|)
+ (RETURN
+ (SEQ
+ (COND
+ ((SPADCALL |z| (QREFELT $ 8))
+ (CONS 0 (CONS (|spadConstant| $ 16) (|spadConstant| $ 16))))
+ ((QUOTE T)
+ (SEQ
+ (LETT |s|
+ (SPADCALL |x| |y| (QREFELT $ 32))
+ |EUCDOM-;extendedEuclidean;3SU;8|)
+ (LETT |w|
+ (SPADCALL |z| (QVELT |s| 2) (QREFELT $ 33))
+ |EUCDOM-;extendedEuclidean;3SU;8|)
+ (EXIT
+ (COND
+ ((QEQCAR |w| 1) (CONS 1 "failed"))
+ ((SPADCALL |y| (QREFELT $ 8))
+ (CONS 0
+ (CONS (SPADCALL (QVELT |s| 0) (QCDR |w|) (QREFELT $ 26))
+ (SPADCALL (QVELT |s| 1) (QCDR |w|) (QREFELT $ 26)))))
+ ((QUOTE T)
+ (SEQ
+ (LETT |qr|
+ (SPADCALL
+ (SPADCALL (QVELT |s| 0) (QCDR |w|) (QREFELT $ 26))
+ |y|
+ (QREFELT $ 13))
+ |EUCDOM-;extendedEuclidean;3SU;8|)
+ (EXIT
+ (CONS 0
+ (CONS (QCDR |qr|)
+ (SPADCALL
+ (SPADCALL (QVELT |s| 1) (QCDR |w|) (QREFELT $ 26))
+ (SPADCALL (QCAR |qr|) |x| (QREFELT $ 26))
+ (QREFELT $ 29))))))))))))))))
@
\subsubsection{EUCDOM-;principalIdeal;LR;9}
<<EUCDOM-;principalIdeal;LR;9>>=
-(DEFUN |EUCDOM-;principalIdeal;LR;9| (|l| |$|)
- (PROG (|uca| |v| |u| #1=#:G83663 |vv| #2=#:G83664)
- (RETURN
- (SEQ
- (COND
- ((SPADCALL |l| NIL (QREFELT |$| 38))
- (|error| "empty list passed to principalIdeal"))
- ((SPADCALL (CDR |l|) NIL (QREFELT |$| 38))
- (SEQ
- (LETT |uca|
- (SPADCALL (|SPADfirst| |l|) (QREFELT |$| 22))
- |EUCDOM-;principalIdeal;LR;9|)
- (EXIT (CONS (LIST (QVELT |uca| 0)) (QVELT |uca| 1)))))
- ((SPADCALL (CDR (CDR |l|)) NIL (QREFELT |$| 38))
- (SEQ
- (LETT |u|
- (SPADCALL
- (|SPADfirst| |l|)
- (SPADCALL |l| (QREFELT |$| 39))
- (QREFELT |$| 32))
- |EUCDOM-;principalIdeal;LR;9|)
- (EXIT
- (CONS (LIST (QVELT |u| 0) (QVELT |u| 1)) (QVELT |u| 2)))))
- ((QUOTE T)
- (SEQ
- (LETT |v|
- (SPADCALL (CDR |l|) (QREFELT |$| 42))
- |EUCDOM-;principalIdeal;LR;9|)
- (LETT |u|
- (SPADCALL (|SPADfirst| |l|) (QCDR |v|) (QREFELT |$| 32))
- |EUCDOM-;principalIdeal;LR;9|)
- (EXIT
- (CONS
- (CONS
- (QVELT |u| 0)
- (PROGN
- (LETT #1# NIL |EUCDOM-;principalIdeal;LR;9|)
- (SEQ
- (LETT |vv| NIL |EUCDOM-;principalIdeal;LR;9|)
- (LETT #2# (QCAR |v|) |EUCDOM-;principalIdeal;LR;9|)
- G190
- (COND
- ((OR
- (ATOM #2#)
- (PROGN
- (LETT |vv|
- (CAR #2#)
- |EUCDOM-;principalIdeal;LR;9|)
- NIL))
- (GO G191)))
- (SEQ
- (EXIT
- (LETT #1#
- (CONS
- (SPADCALL
- (QVELT |u| 1)
- |vv|
- (QREFELT |$| 24))
- #1#)
- |EUCDOM-;principalIdeal;LR;9|)))
- (LETT #2# (CDR #2#) |EUCDOM-;principalIdeal;LR;9|)
- (GO G190)
- G191
- (EXIT (NREVERSE0 #1#)))))
- (QVELT |u| 2))))))))))
+(DEFUN |EUCDOM-;principalIdeal;LR;9| (|l| $)
+ (PROG (|uca| |v| |u| #0=#:G1497 |vv| #1=#:G1498)
+ (RETURN
+ (SEQ
+ (COND
+ ((SPADCALL |l| NIL (QREFELT $ 38))
+ (|error| "empty list passed to principalIdeal"))
+ ((SPADCALL (CDR |l|) NIL (QREFELT $ 38))
+ (SEQ
+ (LETT |uca|
+ (SPADCALL (|SPADfirst| |l|) (QREFELT $ 23))
+ |EUCDOM-;principalIdeal;LR;9|)
+ (EXIT (CONS (LIST (QVELT |uca| 0)) (QVELT |uca| 1)))))
+ ((SPADCALL (CDR (CDR |l|)) NIL (QREFELT $ 38))
+ (SEQ
+ (LETT |u|
+ (SPADCALL (|SPADfirst| |l|)
+ (SPADCALL |l| (QREFELT $ 39)) (QREFELT $ 32))
+ |EUCDOM-;principalIdeal;LR;9|)
+ (EXIT (CONS (LIST (QVELT |u| 0) (QVELT |u| 1)) (QVELT |u| 2)))))
+ ((QUOTE T)
+ (SEQ
+ (LETT |v|
+ (SPADCALL (CDR |l|) (QREFELT $ 42))
+ |EUCDOM-;principalIdeal;LR;9|)
+ (LETT |u|
+ (SPADCALL (|SPADfirst| |l|) (QCDR |v|) (QREFELT $ 32))
+ |EUCDOM-;principalIdeal;LR;9|)
+ (EXIT
+ (CONS
+ (CONS (QVELT |u| 0)
+ (PROGN
+ (LETT #0# NIL |EUCDOM-;principalIdeal;LR;9|)
+ (SEQ
+ (LETT |vv| NIL |EUCDOM-;principalIdeal;LR;9|)
+ (LETT #1# (QCAR |v|) |EUCDOM-;principalIdeal;LR;9|)
+ G190
+ (COND
+ ((OR (ATOM #1#)
+ (PROGN
+ (LETT |vv| (CAR #1#) |EUCDOM-;principalIdeal;LR;9|) NIL))
+ (GO G191)))
+ (SEQ
+ (EXIT
+ (LETT #0#
+ (CONS (SPADCALL (QVELT |u| 1) |vv| (QREFELT $ 26))
+ #0#)
+ |EUCDOM-;principalIdeal;LR;9|)))
+ (LETT #1# (CDR #1#)
+ |EUCDOM-;principalIdeal;LR;9|)
+ (GO G190)
+ G191
+ (EXIT (NREVERSE0 #0#)))))
+ (QVELT |u| 2))))))))))
+
@
\subsubsection{EUCDOM-;expressIdealMember;LSU;10}
<<EUCDOM-;expressIdealMember;LSU;10>>=
-(DEFUN |EUCDOM-;expressIdealMember;LSU;10| (|l| |z| |$|)
- (PROG (#1=#:G83681 #2=#:G83682 |pid| |q| #3=#:G83679 |v| #4=#:G83680)
- (RETURN
- (SEQ
- (COND
- ((SPADCALL |z| (|spadConstant| |$| 26) (QREFELT |$| 44))
- (CONS
- 0
- (PROGN
- (LETT #1# NIL |EUCDOM-;expressIdealMember;LSU;10|)
- (SEQ
- (LETT |v| NIL |EUCDOM-;expressIdealMember;LSU;10|)
- (LETT #2# |l| |EUCDOM-;expressIdealMember;LSU;10|)
- G190
- (COND
- ((OR
- (ATOM #2#)
- (PROGN
- (LETT |v|
- (CAR #2#)
- |EUCDOM-;expressIdealMember;LSU;10|)
- NIL))
- (GO G191)))
- (SEQ
- (EXIT
- (LETT #1#
- (CONS (|spadConstant| |$| 26) #1#)
- |EUCDOM-;expressIdealMember;LSU;10|)))
- (LETT #2# (CDR #2#) |EUCDOM-;expressIdealMember;LSU;10|)
- (GO G190)
- G191
- (EXIT (NREVERSE0 #1#))))))
- ((QUOTE T)
- (SEQ
- (LETT |pid|
- (SPADCALL |l| (QREFELT |$| 42))
- |EUCDOM-;expressIdealMember;LSU;10|)
- (LETT |q|
- (SPADCALL |z| (QCDR |pid|) (QREFELT |$| 33))
- |EUCDOM-;expressIdealMember;LSU;10|)
- (EXIT
- (COND
- ((QEQCAR |q| 1) (CONS 1 "failed"))
- ((QUOTE T)
- (CONS
- 0
- (PROGN
- (LETT #3# NIL |EUCDOM-;expressIdealMember;LSU;10|)
- (SEQ
- (LETT |v| NIL |EUCDOM-;expressIdealMember;LSU;10|)
- (LETT #4# (QCAR |pid|) |EUCDOM-;expressIdealMember;LSU;10|)
- G190
- (COND
- ((OR
- (ATOM #4#)
- (PROGN
- (LETT |v|
- (CAR #4#)
- |EUCDOM-;expressIdealMember;LSU;10|)
- NIL))
- (GO G191)))
- (SEQ
- (EXIT
- (LETT #3#
- (CONS
- (SPADCALL (QCDR |q|) |v| (QREFELT |$| 24))
- #3#)
- |EUCDOM-;expressIdealMember;LSU;10|)))
- (LETT #4#
- (CDR #4#)
- |EUCDOM-;expressIdealMember;LSU;10|)
- (GO G190)
- G191
- (EXIT (NREVERSE0 #3#)))))))))))))))
+(DEFUN |EUCDOM-;expressIdealMember;LSU;10| (|l| |z| $)
+ (PROG (#0=#:G1513 #1=#:G1514 |pid| |q| #2=#:G1515 |v| #3=#:G1516)
+ (RETURN
+ (SEQ
+ (COND
+ ((SPADCALL |z| (|spadConstant| $ 16) (QREFELT $ 25))
+ (CONS 0
+ (PROGN
+ (LETT #0# NIL |EUCDOM-;expressIdealMember;LSU;10|)
+ (SEQ
+ (LETT |v| NIL |EUCDOM-;expressIdealMember;LSU;10|)
+ (LETT #1# |l| |EUCDOM-;expressIdealMember;LSU;10|)
+ G190
+ (COND
+ ((OR (ATOM #1#)
+ (PROGN
+ (LETT |v| (CAR #1#) |EUCDOM-;expressIdealMember;LSU;10|) NIL))
+ (GO G191)))
+ (SEQ
+ (EXIT
+ (LETT #0#
+ (CONS (|spadConstant| $ 16) #0#)
+ |EUCDOM-;expressIdealMember;LSU;10|)))
+ (LETT #1# (CDR #1#) |EUCDOM-;expressIdealMember;LSU;10|)
+ (GO G190)
+ G191
+ (EXIT (NREVERSE0 #0#))))))
+ ((QUOTE T)
+ (SEQ
+ (LETT |pid|
+ (SPADCALL |l| (QREFELT $ 42))
+ |EUCDOM-;expressIdealMember;LSU;10|)
+ (LETT |q|
+ (SPADCALL |z| (QCDR |pid|) (QREFELT $ 33))
+ |EUCDOM-;expressIdealMember;LSU;10|)
+ (EXIT
+ (COND
+ ((QEQCAR |q| 1) (CONS 1 "failed"))
+ ((QUOTE T)
+ (CONS 0
+ (PROGN
+ (LETT #2# NIL |EUCDOM-;expressIdealMember;LSU;10|)
+ (SEQ
+ (LETT |v| NIL |EUCDOM-;expressIdealMember;LSU;10|)
+ (LETT #3# (QCAR |pid|) |EUCDOM-;expressIdealMember;LSU;10|)
+ G190
+ (COND
+ ((OR (ATOM #3#)
+ (PROGN
+ (LETT |v| (CAR #3#) |EUCDOM-;expressIdealMember;LSU;10|)
+ NIL))
+ (GO G191)))
+ (SEQ
+ (EXIT
+ (LETT #2#
+ (CONS (SPADCALL (QCDR |q|) |v| (QREFELT $ 26))
+ #2#)
+ |EUCDOM-;expressIdealMember;LSU;10|)))
+ (LETT #3# (CDR #3#) |EUCDOM-;expressIdealMember;LSU;10|)
+ (GO G190)
+ G191
+ (EXIT (NREVERSE0 #2#)))))))))))))))
@
\subsubsection{EUCDOM-;multiEuclidean;LSU;11}
<<EUCDOM-;multiEuclidean;LSU;11>>=
-(DEFUN |EUCDOM-;multiEuclidean;LSU;11| (|l| |z| |$|)
- (PROG (|n| |l1| |l2| #1=#:G83565 #2=#:G83702 #3=#:G83688 #4=#:G83686
- #5=#:G83687 #6=#:G83566 #7=#:G83701 #8=#:G83691 #9=#:G83689
- #10=#:G83690 |u| |v1| |v2|)
- (RETURN
- (SEQ
- (LETT |n| (LENGTH |l|) |EUCDOM-;multiEuclidean;LSU;11|)
- (EXIT
- (COND
- ((ZEROP |n|) (|error| "empty list passed to multiEuclidean"))
- ((EQL |n| 1) (CONS 0 (LIST |z|)))
- ((QUOTE T)
- (SEQ
- (LETT |l1|
- (SPADCALL |l| (QREFELT |$| 47))
- |EUCDOM-;multiEuclidean;LSU;11|)
- (LETT |l2|
- (SPADCALL |l1| (QUOTIENT2 |n| 2) (QREFELT |$| 49))
- |EUCDOM-;multiEuclidean;LSU;11|)
- (LETT |u|
- (SPADCALL
- (PROGN
- (LETT #5# NIL |EUCDOM-;multiEuclidean;LSU;11|)
- (SEQ
- (LETT #1# NIL |EUCDOM-;multiEuclidean;LSU;11|)
- (LETT #2# |l1| |EUCDOM-;multiEuclidean;LSU;11|)
- G190
- (COND
- ((OR
- (ATOM #2#)
- (PROGN
- (LETT #1#
- (CAR #2#)
- |EUCDOM-;multiEuclidean;LSU;11|)
- NIL))
- (GO G191)))
- (SEQ
- (EXIT
- (PROGN
- (LETT #3# #1# |EUCDOM-;multiEuclidean;LSU;11|)
- (COND
- (#5#
- (LETT #4#
- (SPADCALL #4# #3# (QREFELT |$| 24))
- |EUCDOM-;multiEuclidean;LSU;11|))
- ((QUOTE T)
- (PROGN
- (LETT #4#
- #3#
- |EUCDOM-;multiEuclidean;LSU;11|)
- (LETT #5#
- (QUOTE T)
- |EUCDOM-;multiEuclidean;LSU;11|)))))))
- (LETT #2# (CDR #2#) |EUCDOM-;multiEuclidean;LSU;11|)
- (GO G190)
- G191
- (EXIT NIL))
- (COND (#5# #4#) ((QUOTE T) (|spadConstant| |$| 25))))
- (PROGN
- (LETT #10# NIL |EUCDOM-;multiEuclidean;LSU;11|)
- (SEQ
- (LETT #6# NIL |EUCDOM-;multiEuclidean;LSU;11|)
- (LETT #7# |l2| |EUCDOM-;multiEuclidean;LSU;11|)
- G190
- (COND
- ((OR
- (ATOM #7#)
- (PROGN
- (LETT #6#
- (CAR #7#)
- |EUCDOM-;multiEuclidean;LSU;11|)
- NIL))
- (GO G191)))
- (SEQ
- (EXIT
- (PROGN
- (LETT #8# #6# |EUCDOM-;multiEuclidean;LSU;11|)
- (COND
- (#10#
- (LETT #9#
- (SPADCALL #9# #8# (QREFELT |$| 24))
- |EUCDOM-;multiEuclidean;LSU;11|))
- ((QUOTE T)
- (PROGN
- (LETT #9#
- #8#
- |EUCDOM-;multiEuclidean;LSU;11|)
- (LETT #10#
- (QUOTE T)
- |EUCDOM-;multiEuclidean;LSU;11|)))))))
- (LETT #7# (CDR #7#) |EUCDOM-;multiEuclidean;LSU;11|)
- (GO G190)
- G191
- (EXIT NIL))
- (COND
- (#10# #9#)
- ((QUOTE T) (|spadConstant| |$| 25))))
- |z|
- (QREFELT |$| 50))
- |EUCDOM-;multiEuclidean;LSU;11|)
- (EXIT
- (COND
- ((QEQCAR |u| 1) (CONS 1 "failed"))
- ((QUOTE T)
- (SEQ
- (LETT |v1|
- (SPADCALL |l1| (QCDR (QCDR |u|)) (QREFELT |$| 51))
- |EUCDOM-;multiEuclidean;LSU;11|)
- (EXIT
- (COND
- ((QEQCAR |v1| 1) (CONS 1 "failed"))
- ((QUOTE T)
- (SEQ
- (LETT |v2|
- (SPADCALL
- |l2|
- (QCAR (QCDR |u|))
- (QREFELT |$| 51))
- |EUCDOM-;multiEuclidean;LSU;11|)
- (EXIT
- (COND
- ((QEQCAR |v2| 1) (CONS 1 "failed"))
- ((QUOTE T)
- (CONS
- 0
- (SPADCALL
- (QCDR |v1|)
- (QCDR |v2|)
- (QREFELT |$| 52))))))))))))))))))))))
+(DEFUN |EUCDOM-;multiEuclidean;LSU;11| (|l| |z| $)
+ (PROG (|n| |l1| |l2| #0=#:G1405 #1=#:G1535 #2=#:G1522 #3=#:G1520
+ #4=#:G1521 #5=#:G1406 #6=#:G1536 #7=#:G1525 #8=#:G1523 #9=#:G1524
+ |u| |v1| |v2|)
+ (RETURN
+ (SEQ
+ (LETT |n| (LENGTH |l|) |EUCDOM-;multiEuclidean;LSU;11|)
+ (EXIT
+ (COND
+ ((ZEROP |n|) (|error| "empty list passed to multiEuclidean"))
+ ((EQL |n| 1) (CONS 0 (LIST |z|)))
+ ((QUOTE T)
+ (SEQ
+ (LETT |l1|
+ (SPADCALL |l| (QREFELT $ 46)) |EUCDOM-;multiEuclidean;LSU;11|)
+ (LETT |l2|
+ (SPADCALL |l1| (QUOTIENT2 |n| 2) (QREFELT $ 48))
+ |EUCDOM-;multiEuclidean;LSU;11|)
+ (LETT |u|
+ (SPADCALL
+ (PROGN
+ (LETT #4# NIL |EUCDOM-;multiEuclidean;LSU;11|)
+ (SEQ
+ (LETT #0# NIL |EUCDOM-;multiEuclidean;LSU;11|)
+ (LETT #1# |l1| |EUCDOM-;multiEuclidean;LSU;11|)
+ G190
+ (COND
+ ((OR (ATOM #1#)
+ (PROGN
+ (LETT #0# (CAR #1#) |EUCDOM-;multiEuclidean;LSU;11|)
+ NIL))
+ (GO G191)))
+ (SEQ
+ (EXIT
+ (PROGN
+ (LETT #2# #0# |EUCDOM-;multiEuclidean;LSU;11|)
+ (COND
+ (#4#
+ (LETT #3#
+ (SPADCALL #3# #2# (QREFELT $ 26))
+ |EUCDOM-;multiEuclidean;LSU;11|))
+ ((QUOTE T)
+ (PROGN
+ (LETT #3# #2# |EUCDOM-;multiEuclidean;LSU;11|)
+ (LETT #4# (QUOTE T) |EUCDOM-;multiEuclidean;LSU;11|)))))))
+ (LETT #1# (CDR #1#) |EUCDOM-;multiEuclidean;LSU;11|)
+ (GO G190)
+ G191
+ (EXIT NIL))
+ (COND (#4# #3#) ((QUOTE T) (|spadConstant| $ 24))))
+ (PROGN
+ (LETT #9# NIL |EUCDOM-;multiEuclidean;LSU;11|)
+ (SEQ
+ (LETT #5# NIL |EUCDOM-;multiEuclidean;LSU;11|)
+ (LETT #6# |l2| |EUCDOM-;multiEuclidean;LSU;11|)
+ G190
+ (COND
+ ((OR (ATOM #6#)
+ (PROGN
+ (LETT #5# (CAR #6#) |EUCDOM-;multiEuclidean;LSU;11|)
+ NIL))
+ (GO G191)))
+ (SEQ
+ (EXIT
+ (PROGN
+ (LETT #7# #5# |EUCDOM-;multiEuclidean;LSU;11|)
+ (COND
+ (#9#
+ (LETT #8#
+ (SPADCALL #8# #7# (QREFELT $ 26))
+ |EUCDOM-;multiEuclidean;LSU;11|))
+ ((QUOTE T)
+ (PROGN
+ (LETT #8# #7# |EUCDOM-;multiEuclidean;LSU;11|)
+ (LETT #9# (QUOTE T) |EUCDOM-;multiEuclidean;LSU;11|)))))))
+ (LETT #6# (CDR #6#) |EUCDOM-;multiEuclidean;LSU;11|)
+ (GO G190)
+ G191
+ (EXIT NIL))
+ (COND (#9# #8#) ((QUOTE T) (|spadConstant| $ 24))))
+ |z| (QREFELT $ 49))
+ |EUCDOM-;multiEuclidean;LSU;11|)
+ (EXIT
+ (COND
+ ((QEQCAR |u| 1) (CONS 1 "failed"))
+ ((QUOTE T)
+ (SEQ
+ (LETT |v1|
+ (SPADCALL |l1| (QCDR (QCDR |u|)) (QREFELT $ 50))
+ |EUCDOM-;multiEuclidean;LSU;11|)
+ (EXIT
+ (COND
+ ((QEQCAR |v1| 1) (CONS 1 "failed"))
+ ((QUOTE T)
+ (SEQ
+ (LETT |v2|
+ (SPADCALL |l2| (QCAR (QCDR |u|)) (QREFELT $ 50))
+ |EUCDOM-;multiEuclidean;LSU;11|)
+ (EXIT
+ (COND
+ ((QEQCAR |v2| 1) (CONS 1 "failed"))
+ ((QUOTE T)
+ (CONS 0
+ (SPADCALL
+ (QCDR |v1|)
+ (QCDR |v2|)
+ (QREFELT $ 51))))))))))))))))))))))
@
\subsubsection{EuclideanDomain\&}
<<EuclideanDomainAmp>>=
-(DEFUN |EuclideanDomain&| (|#1|)
- (PROG (|DV$1| |dv$| |$| |pv$|)
- (RETURN
- (PROGN
- (LETT |DV$1| (|devaluate| |#1|) . #1=(|EuclideanDomain&|))
- (LETT |dv$| (LIST (QUOTE |EuclideanDomain&|) |DV$1|) . #1#)
- (LETT |$| (GETREFV 54) . #1#)
- (QSETREFV |$| 0 |dv$|)
- (QSETREFV |$| 3 (LETT |pv$| (|buildPredVector| 0 0 NIL) . #1#))
- (|stuffDomainSlots| |$|)
- (QSETREFV |$| 6 |#1|)
- |$|))))
+(DEFUN |EuclideanDomain&| (|#1|)
+ (PROG (DV$1 |dv$| $ |pv$|)
+ (RETURN
+ (PROGN
+ (LETT DV$1 (|devaluate| |#1|) . #0=(|EuclideanDomain&|))
+ (LETT |dv$| (LIST (QUOTE |EuclideanDomain&|) DV$1) . #0#)
+ (LETT $ (GETREFV 53) . #0#)
+ (QSETREFV $ 0 |dv$|)
+ (QSETREFV $ 3 (LETT |pv$| (|buildPredVector| 0 0 NIL) . #0#))
+ (|stuffDomainSlots| $)
+ (QSETREFV $ 6 |#1|)
+ $))))
@
\subsubsection{EUCDOM-;MAKEPROP}
<<EUCDOM-;MAKEPROP>>=
-(MAKEPROP
- (QUOTE |EuclideanDomain&|)
- (QUOTE |infovec|)
- (LIST
- (QUOTE
- #(NIL NIL NIL NIL NIL NIL
- (|local| |#1|)
- (|Boolean|)
- (0 . |zero?|)
- (|NonNegativeInteger|)
- (5 . |euclideanSize|)
- |EUCDOM-;sizeLess?;2SB;1|
- (|Record| (|:| |quotient| |$|) (|:| |remainder| |$|))
- (10 . |divide|)
- |EUCDOM-;quo;3S;2|
- |EUCDOM-;rem;3S;3|
- (|Union| |$| (QUOTE "failed"))
- |EUCDOM-;exquo;2SU;4|
- (16 . |unitCanonical|)
- (21 . |rem|)
- |EUCDOM-;gcd;3S;5|
- (|Record| (|:| |unit| |$|) (|:| |canonical| |$|) (|:| |associate| |$|))
- (27 . |unitNormal|)
- (32 . |one?|)
- (37 . |*|)
- (43 . |One|)
- (47 . |Zero|)
- (51 . |-|)
- (57 . |sizeLess?|)
- (63 . |+|)
- (|Record| (|:| |coef1| |$|) (|:| |coef2| |$|) (|:| |generator| |$|))
- |EUCDOM-;extendedEuclidean;2SR;7|
- (69 . |extendedEuclidean|)
- (75 . |exquo|)
- (|Record| (|:| |coef1| |$|) (|:| |coef2| |$|))
- (|Union| 34 (QUOTE "failed"))
- |EUCDOM-;extendedEuclidean;3SU;8|
- (|List| 6)
- (81 . |=|)
- (87 . |second|)
- (|Record| (|:| |coef| 41) (|:| |generator| |$|))
- (|List| |$|)
- (92 . |principalIdeal|)
- |EUCDOM-;principalIdeal;LR;9|
- (97 . |=|)
- (|Union| 41 (QUOTE "failed"))
- |EUCDOM-;expressIdealMember;LSU;10|
- (103 . |copy|)
- (|Integer|)
- (108 . |split!|)
- (114 . |extendedEuclidean|)
- (121 . |multiEuclidean|)
- (127 . |concat|)
- |EUCDOM-;multiEuclidean;LSU;11|))
- (QUOTE
- #(|sizeLess?| 133 |rem| 139 |quo| 145 |principalIdeal| 151
- |multiEuclidean| 156 |gcd| 162 |extendedEuclidean| 168 |exquo| 181
- |expressIdealMember| 187))
- (QUOTE NIL)
- (CONS
- (|makeByteWordVec2| 1 (QUOTE NIL))
- (CONS
- (QUOTE #())
- (CONS
- (QUOTE #())
- (|makeByteWordVec2| 53
- (QUOTE
- (1 6 7 0 8 1 6 9 0 10 2 6 12 0 0 13 1 6 0 0 18 2 6 0 0 0 19 1 6
- 21 0 22 1 6 7 0 23 2 6 0 0 0 24 0 6 0 25 0 6 0 26 2 6 0 0 0 27
- 2 6 7 0 0 28 2 6 0 0 0 29 2 6 30 0 0 32 2 6 16 0 0 33 2 37 7 0
- 0 38 1 37 6 0 39 1 6 40 41 42 2 6 7 0 0 44 1 37 0 0 47 2 37 0 0
- 48 49 3 6 35 0 0 0 50 2 6 45 41 0 51 2 37 0 0 0 52 2 0 7 0 0 11
- 2 0 0 0 0 15 2 0 0 0 0 14 1 0 40 41 43 2 0 45 41 0 53 2 0 0 0 0
- 20 3 0 35 0 0 0 36 2 0 30 0 0 31 2 0 16 0 0 17 2 0 45 41 0
- 46))))))
- (QUOTE |lookupComplete|)))
+(MAKEPROP
+ (QUOTE |EuclideanDomain&|)
+ (QUOTE |infovec|)
+ (LIST
+ (QUOTE #(NIL NIL NIL NIL NIL NIL (|local| |#1|) (|Boolean|) (0 . |zero?|)
+ (|NonNegativeInteger|) (5 . |euclideanSize|) |EUCDOM-;sizeLess?;2SB;1|
+ (|Record| (|:| |quotient| $) (|:| |remainder| $)) (10 . |divide|)
+ |EUCDOM-;quo;3S;2| |EUCDOM-;rem;3S;3| (16 . |Zero|)
+ (|Union| $ (QUOTE "failed")) |EUCDOM-;exquo;2SU;4| (20 . |unitCanonical|)
+ (25 . |rem|) |EUCDOM-;gcd;3S;5|
+ (|Record| (|:| |unit| $) (|:| |canonical| $) (|:| |associate| $))
+ (31 . |unitNormal|) (36 . |One|) (40 . =) (46 . *) (52 . -)
+ (58 . |sizeLess?|) (64 . +)
+ (|Record| (|:| |coef1| $) (|:| |coef2| $) (|:| |generator| $))
+ |EUCDOM-;extendedEuclidean;2SR;7|
+ (70 . |extendedEuclidean|) (76 . |exquo|)
+ (|Record| (|:| |coef1| $) (|:| |coef2| $))
+ (|Union| 34 (QUOTE "failed")) |EUCDOM-;extendedEuclidean;3SU;8|
+ (|List| 6) (82 . =) (88 . |second|)
+ (|Record| (|:| |coef| 41) (|:| |generator| $))
+ (|List| $) (93 . |principalIdeal|) |EUCDOM-;principalIdeal;LR;9|
+ (|Union| 41 (QUOTE "failed")) |EUCDOM-;expressIdealMember;LSU;10|
+ (98 . |copy|) (|Integer|) (103 . |split!|) (109 . |extendedEuclidean|)
+ (116 . |multiEuclidean|) (122 . |concat|) |EUCDOM-;multiEuclidean;LSU;11|))
+ (QUOTE
+ #(|sizeLess?| 128 |rem| 134 |quo| 140 |principalIdeal| 146
+ |multiEuclidean| 151 |gcd| 157 |extendedEuclidean| 163
+ |exquo| 176 |expressIdealMember| 182))
+ (QUOTE NIL)
+ (CONS (|makeByteWordVec2| 1 (QUOTE NIL))
+ (CONS (QUOTE #())
+ (CONS (QUOTE #())
+ (|makeByteWordVec2| 52 (QUOTE (1 6 7 0 8 1 6 9 0 10 2 6 12 0 0 13 0
+ 6 0 16 1 6 0 0 19 2 6 0 0 0 20 1 6 22 0 23 0 6 0 24 2 6 7 0 0 25 2 6 0
+ 0 0 26 2 6 0 0 0 27 2 6 7 0 0 28 2 6 0 0 0 29 2 6 30 0 0 32 2 6 17 0 0
+ 33 2 37 7 0 0 38 1 37 6 0 39 1 6 40 41 42 1 37 0 0 46 2 37 0 0 47 48 3
+ 6 35 0 0 0 49 2 6 44 41 0 50 2 37 0 0 0 51 2 0 7 0 0 11 2 0 0 0 0 15 2
+ 0 0 0 0 14 1 0 40 41 43 2 0 44 41 0 52 2 0 0 0 0 21 3 0 35 0 0 0 36 2 0
+ 30 0 0 31 2 0 17 0 0 18 2 0 44 41 0 45))))))
+ (QUOTE |lookupComplete|)))
@
<<EUCDOM-.lsp BOOTSTRAP>>=
Modified: trunk/axiom/src/algebra/polycat.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/polycat.spad.pamphlet 2007-09-16 03:12:14 UTC (rev 748)
+++ trunk/axiom/src/algebra/polycat.spad.pamphlet 2007-09-16 22:49:24 UTC (rev 749)
@@ -568,8 +568,10 @@
unit(s := squareFree p) * */[f.factor for f in factors s]
content(p,v) == content univariate(p,v)
primitivePart p ==
+ zero? p => p
unitNormal((p exquo content p) ::%).canonical
primitivePart(p,v) ==
+ zero? p => p
unitNormal((p exquo content(p,v)) ::%).canonical
if R has OrderedSet then
p:% < q:% ==
@@ -617,171 +619,125 @@
<<POLYCAT.lsp BOOTSTRAP>>=
-(|/VERSIONCHECK| 2)
+(/VERSIONCHECK 2)
(SETQ |PolynomialCategory;CAT| (QUOTE NIL))
(SETQ |PolynomialCategory;AL| (QUOTE NIL))
-(DEFUN |PolynomialCategory| (|&REST| #1=#:G101841 |&AUX| #2=#:G101839)
- (DSETQ #2# #1#)
- (LET (#3=#:G101840)
- (COND
- ((SETQ #3# (|assoc| (|devaluateList| #2#) |PolynomialCategory;AL|))
- (CDR #3#))
- (T
- (SETQ |PolynomialCategory;AL|
- (|cons5|
- (CONS
- (|devaluateList| #2#)
- (SETQ #3# (APPLY (FUNCTION |PolynomialCategory;|) #2#)))
- |PolynomialCategory;AL|))
- #3#))))
+(DEFUN |PolynomialCategory| (&REST #0=#:G1430 &AUX #1=#:G1428)
+ (DSETQ #1# #0#)
+ (LET (#2=#:G1429)
+ (COND
+ ((SETQ #2# (|assoc| (|devaluateList| #1#) |PolynomialCategory;AL|))
+ (CDR #2#))
+ (T
+ (SETQ |PolynomialCategory;AL|
+ (|cons5|
+ (CONS (|devaluateList| #1#)
+ (SETQ #2# (APPLY (FUNCTION |PolynomialCategory;|) #1#)))
+ |PolynomialCategory;AL|))
+ #2#))))
-(DEFUN |PolynomialCategory;| (|t#1| |t#2| |t#3|)
- (PROG (#1=#:G101838)
- (RETURN
- (PROG1
- (LETT #1#
- (|sublisV|
- (PAIR
- (QUOTE (|t#1| |t#2| |t#3|))
- (LIST
- (|devaluate| |t#1|)
- (|devaluate| |t#2|)
- (|devaluate| |t#3|)))
- (COND
- (|PolynomialCategory;CAT|)
- ((QUOTE T)
- (LETT |PolynomialCategory;CAT|
- (|Join|
- (|PartialDifferentialRing| (QUOTE |t#3|))
- (|FiniteAbelianMonoidRing| (QUOTE |t#1|) (QUOTE |t#2|))
- (|Evalable| (QUOTE |$|))
- (|InnerEvalable| (QUOTE |t#3|) (QUOTE |t#1|))
- (|InnerEvalable| (QUOTE |t#3|) (QUOTE |$|))
- (|RetractableTo| (QUOTE |t#3|))
- (|FullyLinearlyExplicitRingOver| (QUOTE |t#1|))
- (|mkCategory|
- (QUOTE |domain|)
- (QUOTE (
- ((|degree| ((|NonNegativeInteger|) |$| |t#3|)) T)
- ((|degree|
- ((|List| (|NonNegativeInteger|))
- |$|
- (|List| |t#3|))) T)
- ((|coefficient|
- (|$| |$| |t#3| (|NonNegativeInteger|))) T)
- ((|coefficient|
- (|$|
- |$|
- (|List| |t#3|)
- (|List| (|NonNegativeInteger|)))) T)
- ((|monomials| ((|List| |$|) |$|)) T)
- ((|univariate|
- ((|SparseUnivariatePolynomial| |$|) |$| |t#3|)) T)
- ((|univariate|
- ((|SparseUnivariatePolynomial| |t#1|) |$|)) T)
- ((|mainVariable| ((|Union| |t#3| "failed") |$|)) T)
- ((|minimumDegree|
- ((|NonNegativeInteger|) |$| |t#3|)) T)
- ((|minimumDegree|
- ((|List| (|NonNegativeInteger|))
- |$|
- (|List| |t#3|))) T)
- ((|monicDivide|
- ((|Record|
- (|:| |quotient| |$|)
- (|:| |remainder| |$|))
- |$|
- |$|
- |t#3|)) T)
- ((|monomial| (|$| |$| |t#3| (|NonNegativeInteger|))) T)
- ((|monomial|
- (|$|
- |$|
- (|List| |t#3|)
- (|List| (|NonNegativeInteger|)))) T)
- ((|multivariate|
- (|$| (|SparseUnivariatePolynomial| |t#1|) |t#3|)) T)
- ((|multivariate|
- (|$| (|SparseUnivariatePolynomial| |$|) |t#3|)) T)
- ((|isPlus| ((|Union| (|List| |$|) "failed") |$|)) T)
- ((|isTimes| ((|Union| (|List| |$|) "failed") |$|)) T)
- ((|isExpt|
- ((|Union|
- (|Record|
- (|:| |var| |t#3|)
- (|:| |exponent| (|NonNegativeInteger|)))
- "failed")
- |$|)) T)
- ((|totalDegree| ((|NonNegativeInteger|) |$|)) T)
- ((|totalDegree|
- ((|NonNegativeInteger|) |$| (|List| |t#3|))) T)
- ((|variables| ((|List| |t#3|) |$|)) T)
- ((|primitiveMonomials| ((|List| |$|) |$|)) T)
- ((|resultant| (|$| |$| |$| |t#3|))
- (|has| |t#1| (|CommutativeRing|)))
- ((|discriminant| (|$| |$| |t#3|))
- (|has| |t#1| (|CommutativeRing|)))
- ((|content| (|$| |$| |t#3|))
- (|has| |t#1| (|GcdDomain|)))
- ((|primitivePart| (|$| |$|))
- (|has| |t#1| (|GcdDomain|)))
- ((|primitivePart| (|$| |$| |t#3|))
- (|has| |t#1| (|GcdDomain|)))
- ((|squareFree| ((|Factored| |$|) |$|))
- (|has| |t#1| (|GcdDomain|)))
- ((|squareFreePart| (|$| |$|)) (
- |has| |t#1| (|GcdDomain|)))))
- (QUOTE (
- ((|OrderedSet|) (|has| |t#1| (|OrderedSet|)))
- ((|ConvertibleTo| (|InputForm|))
- (AND
- (|has| |t#3| (|ConvertibleTo| (|InputForm|)))
- (|has| |t#1| (|ConvertibleTo| (|InputForm|)))))
- ((|ConvertibleTo| (|Pattern| (|Integer|)))
- (AND
- (|has| |t#3|
- (|ConvertibleTo| (|Pattern| (|Integer|))))
- (|has| |t#1|
- (|ConvertibleTo| (|Pattern| (|Integer|))))))
- ((|ConvertibleTo| (|Pattern| (|Float|)))
- (AND
- (|has| |t#3|
- (|ConvertibleTo| (|Pattern| (|Float|))))
- (|has| |t#1|
- (|ConvertibleTo| (|Pattern| (|Float|))))))
- ((|PatternMatchable| (|Integer|))
- (AND
- (|has| |t#3| (|PatternMatchable| (|Integer|)))
- (|has| |t#1| (|PatternMatchable| (|Integer|)))))
- ((|PatternMatchable| (|Float|))
- (AND
- (|has| |t#3| (|PatternMatchable| (|Float|)))
- (|has| |t#1| (|PatternMatchable| (|Float|)))))
- ((|GcdDomain|) (|has| |t#1| (|GcdDomain|)))
- (|canonicalUnitNormal|
- (|has| |t#1| (ATTRIBUTE |canonicalUnitNormal|)))
- ((|PolynomialFactorizationExplicit|)
- (|has| |t#1| (|PolynomialFactorizationExplicit|)))))
- (QUOTE (
- (|Factored| |$|)
- (|List| |$|)
- (|List| |t#3|)
- (|NonNegativeInteger|)
- (|SparseUnivariatePolynomial| |$|)
- (|SparseUnivariatePolynomial| |t#1|)
- (|List| (|NonNegativeInteger|))))
- NIL))
- . #2=(|PolynomialCategory|)))))
- . #2#)
- (SETELT #1# 0
- (LIST
- (QUOTE |PolynomialCategory|)
- (|devaluate| |t#1|)
- (|devaluate| |t#2|)
- (|devaluate| |t#3|)))))))
+(DEFUN |PolynomialCategory;| (|t#1| |t#2| |t#3|)
+ (PROG (#0=#:G1427)
+ (RETURN
+ (PROG1
+ (LETT #0#
+ (|sublisV|
+ (PAIR (QUOTE (|t#1| |t#2| |t#3|)) (LIST (|devaluate| |t#1|) (|devaluate| |t#2|) (|devaluate| |t#3|)))
+ (COND
+ (|PolynomialCategory;CAT|)
+ ((QUOTE T)
+ (LETT |PolynomialCategory;CAT|
+ (|Join|
+ (|PartialDifferentialRing| (QUOTE |t#3|))
+ (|FiniteAbelianMonoidRing| (QUOTE |t#1|) (QUOTE |t#2|))
+ (|Evalable| (QUOTE $))
+ (|InnerEvalable| (QUOTE |t#3|) (QUOTE |t#1|))
+ (|InnerEvalable| (QUOTE |t#3|) (QUOTE $))
+ (|RetractableTo| (QUOTE |t#3|))
+ (|FullyLinearlyExplicitRingOver| (QUOTE |t#1|))
+ (|mkCategory| (QUOTE |domain|)
+ (QUOTE
+ (((|degree| ((|NonNegativeInteger|) $ |t#3|)) T)
+ ((|degree| ((|List| (|NonNegativeInteger|)) $ (|List| |t#3|))) T)
+ ((|coefficient| ($ $ |t#3| (|NonNegativeInteger|))) T)
+ ((|coefficient| ($ $ (|List| |t#3|)
+ (|List| (|NonNegativeInteger|)))) T)
+ ((|monomials| ((|List| $) $)) T)
+ ((|univariate| ((|SparseUnivariatePolynomial| $) $ |t#3|)) T)
+ ((|univariate| ((|SparseUnivariatePolynomial| |t#1|) $)) T)
+ ((|mainVariable| ((|Union| |t#3| "failed") $)) T)
+ ((|minimumDegree| ((|NonNegativeInteger|) $ |t#3|)) T)
+ ((|minimumDegree| ((|List| (|NonNegativeInteger|)) $
+ (|List| |t#3|))) T)
+ ((|monicDivide|
+ ((|Record| (|:| |quotient| $) (|:| |remainder| $)) $ $ |t#3|))
+ T)
+ ((|monomial| ($ $ |t#3| (|NonNegativeInteger|))) T)
+ ((|monomial| ($ $ (|List| |t#3|) (|List| (|NonNegativeInteger|))))
+ T)
+ ((|multivariate| ($ (|SparseUnivariatePolynomial| |t#1|) |t#3|))
+ T)
+ ((|multivariate| ($ (|SparseUnivariatePolynomial| $) |t#3|)) T)
+ ((|isPlus| ((|Union| (|List| $) "failed") $)) T)
+ ((|isTimes| ((|Union| (|List| $) "failed") $)) T)
+ ((|isExpt|
+ ((|Union|
+ (|Record| (|:| |var| |t#3|)
+ (|:| |exponent| (|NonNegativeInteger|)))
+ "failed") $))
+ T)
+ ((|totalDegree| ((|NonNegativeInteger|) $)) T)
+ ((|totalDegree| ((|NonNegativeInteger|) $ (|List| |t#3|))) T)
+ ((|variables| ((|List| |t#3|) $)) T)
+ ((|primitiveMonomials| ((|List| $) $)) T)
+ ((|resultant| ($ $ $ |t#3|)) (|has| |t#1| (|CommutativeRing|)))
+ ((|discriminant| ($ $ |t#3|)) (|has| |t#1| (|CommutativeRing|)))
+ ((|content| ($ $ |t#3|)) (|has| |t#1| (|GcdDomain|)))
+ ((|primitivePart| ($ $)) (|has| |t#1| (|GcdDomain|)))
+ ((|primitivePart| ($ $ |t#3|)) (|has| |t#1| (|GcdDomain|)))
+ ((|squareFree| ((|Factored| $) $)) (|has| |t#1| (|GcdDomain|)))
+ ((|squareFreePart| ($ $)) (|has| |t#1| (|GcdDomain|)))))
+ (QUOTE
+ (((|OrderedSet|) (|has| |t#1| (|OrderedSet|)))
+ ((|ConvertibleTo| (|InputForm|))
+ (AND (|has| |t#3| (|ConvertibleTo| (|InputForm|)))
+ (|has| |t#1| (|ConvertibleTo| (|InputForm|)))))
+ ((|ConvertibleTo| (|Pattern| (|Integer|)))
+ (AND (|has| |t#3| (|ConvertibleTo| (|Pattern| (|Integer|))))
+ (|has| |t#1| (|ConvertibleTo| (|Pattern| (|Integer|))))))
+ ((|ConvertibleTo| (|Pattern| (|Float|)))
+ (AND (|has| |t#3| (|ConvertibleTo| (|Pattern| (|Float|))))
+ (|has| |t#1| (|ConvertibleTo| (|Pattern| (|Float|))))))
+ ((|PatternMatchable| (|Integer|))
+ (AND
+ (|has| |t#3| (|PatternMatchable| (|Integer|)))
+ (|has| |t#1| (|PatternMatchable| (|Integer|)))))
+ ((|PatternMatchable| (|Float|))
+ (AND
+ (|has| |t#3| (|PatternMatchable| (|Float|)))
+ (|has| |t#1| (|PatternMatchable| (|Float|)))))
+ ((|GcdDomain|) (|has| |t#1| (|GcdDomain|)))
+ (|canonicalUnitNormal|
+ (|has| |t#1| (ATTRIBUTE |canonicalUnitNormal|)))
+ ((|PolynomialFactorizationExplicit|)
+ (|has| |t#1| (|PolynomialFactorizationExplicit|)))))
+ (QUOTE
+ ((|Factored| $)
+ (|List| $)
+ (|List| |t#3|)
+ (|NonNegativeInteger|)
+ (|SparseUnivariatePolynomial| $)
+ (|SparseUnivariatePolynomial| |t#1|)
+ (|List| (|NonNegativeInteger|))))
+ NIL))
+ . #1=(|PolynomialCategory|)))))
+ . #1#)
+ (SETELT #0# 0
+ (LIST (QUOTE |PolynomialCategory|)
+ (|devaluate| |t#1|) (|devaluate| |t#2|) (|devaluate| |t#3|)))))))
@
\section{POLYCAT-.lsp BOOTSTRAP}
@@ -797,1924 +753,1521 @@
(|/VERSIONCHECK| 2)
-(DEFUN |POLYCAT-;eval;SLS;1| (|p| |l| |$|)
- (PROG (#1=#:G101870 #2=#:G101860 #3=#:G101868 #4=#:G101869
- |lvar| #5=#:G101866 |e| #6=#:G101867)
- (RETURN
- (SEQ
- (COND
- ((NULL |l|) |p|)
- ((QUOTE T)
- (SEQ
- (SEQ
- (EXIT
- (SEQ
- (LETT |e| NIL |POLYCAT-;eval;SLS;1|)
- (LETT #1# |l| |POLYCAT-;eval;SLS;1|)
- G190
- (COND
- ((OR
- (ATOM #1#)
- (PROGN (LETT |e| (CAR #1#) |POLYCAT-;eval;SLS;1|) NIL))
- (GO G191)))
- (SEQ
- (EXIT
- (COND
- ((QEQCAR
- (SPADCALL
- (SPADCALL |e| (QREFELT |$| 11))
- (QREFELT |$| 13))
- 1)
- (PROGN
- (LETT #2#
- (|error| "cannot find a variable to evaluate")
- |POLYCAT-;eval;SLS;1|)
- (GO #2#))))))
- (LETT #1# (CDR #1#) |POLYCAT-;eval;SLS;1|)
- (GO G190)
- G191
- (EXIT NIL)))
- #2#
- (EXIT #2#))
- (LETT |lvar|
- (PROGN
- (LETT #3# NIL |POLYCAT-;eval;SLS;1|)
- (SEQ
- (LETT |e| NIL |POLYCAT-;eval;SLS;1|)
- (LETT #4# |l| |POLYCAT-;eval;SLS;1|)
- G190
- (COND
- ((OR
- (ATOM #4#)
- (PROGN
- (LETT |e| (CAR #4#) |POLYCAT-;eval;SLS;1|) NIL))
- (GO G191)))
- (SEQ
- (EXIT
- (LETT #3#
- (CONS
- (SPADCALL
- (SPADCALL |e| (QREFELT |$| 11))
- (QREFELT |$| 14))
- #3#)
- |POLYCAT-;eval;SLS;1|)))
- (LETT #4# (CDR #4#) |POLYCAT-;eval;SLS;1|)
- (GO G190)
- G191
- (EXIT (NREVERSE0 #3#))))
+
+(/VERSIONCHECK 2)
+
+(DEFUN |POLYCAT-;eval;SLS;1| (|p| |l| $)
+ (PROG (#0=#:G1444 #1=#:G1438 #2=#:G1445 #3=#:G1446 |lvar| #4=#:G1447
+ |e| #5=#:G1448)
+ (RETURN
+ (SEQ
+ (COND
+ ((NULL |l|) |p|)
+ ((QUOTE T)
+ (SEQ
+ (SEQ
+ (EXIT
+ (SEQ
+ (LETT |e| NIL |POLYCAT-;eval;SLS;1|)
+ (LETT #0# |l| |POLYCAT-;eval;SLS;1|)
+ G190
+ (COND
+ ((OR (ATOM #0#)
+ (PROGN (LETT |e| (CAR #0#) |POLYCAT-;eval;SLS;1|) NIL))
+ (GO G191)))
+ (SEQ
+ (EXIT
+ (COND
+ ((QEQCAR
+ (SPADCALL (SPADCALL |e| (QREFELT $ 11)) (QREFELT $ 13)) 1)
+ (PROGN
+ (LETT #1#
+ (|error| "cannot find a variable to evaluate")
|POLYCAT-;eval;SLS;1|)
- (EXIT
- (SPADCALL
- |p|
- |lvar|
- (PROGN
- (LETT #5# NIL |POLYCAT-;eval;SLS;1|)
- (SEQ
- (LETT |e| NIL |POLYCAT-;eval;SLS;1|)
- (LETT #6# |l| |POLYCAT-;eval;SLS;1|)
- G190
- (COND
- ((OR
- (ATOM #6#)
- (PROGN
- (LETT |e| (CAR #6#) |POLYCAT-;eval;SLS;1|)
- NIL))
- (GO G191)))
- (SEQ
- (EXIT
- (LETT #5#
- (CONS (SPADCALL |e| (QREFELT |$| 15)) #5#)
- |POLYCAT-;eval;SLS;1|)))
- (LET...
[truncated message content] |
|
From: <gi...@ax...> - 2007-09-16 03:12:20
|
changelog | 4 ++
src/algebra/intef.spad.pamphlet | 43 ++++++++++++++++++++-
src/input/Makefile.pamphlet | 6 ++-
src/input/bug100.input.pamphlet | 80 +++++++++++++++++++++++++++++++++++++++
4 files changed, 129 insertions(+), 4 deletions(-)
New commits:
commit 606d07277b3fa49ee06808bccdcf38d1ce27dd83
Author: Tim Daly <root@localhost.localdomain>
Date: Thu Sep 13 21:32:20 2007 -0400
20070915 tpd merge bug100 branch
20070915 tpd src/input/Makefile add bug100.input regression test
20070915 tpd src/input/bug100.input test integrate((z^a+1)^b,z) infinite loop
20070915 wxh src/algebra/intef.spad fix integrate((z^a+1)^b,z) infinite loop
commit 3d342b4465dc00ef5d3b22d8739767747117f17a
Author: Tim Daly <root@localhost.localdomain>
Date: Thu Sep 13 21:30:13 2007 -0400
fix bug 100 integrate((z^a+1)^b,z)
|
|
From: <da...@us...> - 2007-09-16 03:12:14
|
Revision: 748
http://axiom.svn.sourceforge.net/axiom/?rev=748&view=rev
Author: daly
Date: 2007-09-15 20:12:14 -0700 (Sat, 15 Sep 2007)
Log Message:
-----------
20070915 tpd merge bug100 branch
20070915 tpd src/input/Makefile add bug100.input regression test
20070915 tpd src/input/bug100.input test integrate((z^a+1)^b,z) infinite loop
20070915 wxh src/algebra/intef.spad fix integrate((z^a+1)^b,z) infinite loop
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/algebra/intef.spad.pamphlet
trunk/axiom/src/input/Makefile.pamphlet
Added Paths:
-----------
trunk/axiom/src/input/bug100.input.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-15 16:38:22 UTC (rev 747)
+++ trunk/axiom/changelog 2007-09-16 03:12:14 UTC (rev 748)
@@ -1,3 +1,7 @@
+20070915 tpd merge bug100 branch
+20070915 tpd src/input/Makefile add bug100.input regression test
+20070915 tpd src/input/bug100.input test integrate((z^a+1)^b,z) infinite loop
+20070915 wxh src/algebra/intef.spad fix integrate((z^a+1)^b,z) infinite loop
20070915 tpd src/algebra/carten minor edit for regression cleanup
20070914 wxh src/hyper/hyper fix bad bracing of )hd change
20070914 tpd src/algebra/fraction.spad remove double )spool command
Modified: trunk/axiom/src/algebra/intef.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/intef.spad.pamphlet 2007-09-15 16:38:22 UTC (rev 747)
+++ trunk/axiom/src/algebra/intef.spad.pamphlet 2007-09-16 03:12:14 UTC (rev 748)
@@ -227,9 +227,41 @@
algint(f, t, y, differentiate(#1, differentiate(#1, x),
differentiate(t::F, x)::UP))
+@
+Bug \#100 is an infinite loop that eventually kills Axiom
+from the input
+\begin{verbatim}
+ integrate((z^a+1)^b,z)
+\end{verbatim}
+
+Line 2 of this function used to read:
+\begin{verbatim}
+ symbolIfCan(k := kmax(l := union(l, varselect(kernels g, x))))
+\end{verbatim}
+
+The loop occurs when the call to union causes
+\begin{verbatim}
+ a log(z)
+ %e
+\end{verbatim}
+to get added to the list every time. This gives the argument to kmax
+\begin{verbatim}
+ a log(z)
+ arg1= [z,%e ]
+\end{verbatim}
+and the result being
+\begin{verbatim}
+ a log(z)
+ %e
+\end{verbatim}
+We keep coming back to process this term, which ends up
+putting the same term back on the list and we loop.
+Waldek's solution is to remove the union call.
+
+<<package INTEF ElementaryIntegration>>=
lfextendedint(f, x, g) ==
empty?(l := varselect(kernels f, x)) => [x::F * f, 0]
- symbolIfCan(k := kmax(l := union(l, varselect(kernels g, x))))
+ symbolIfCan(k := kmax(l))
case SE =>
map(multivariate(#1, k), extendedint(univariate(f, k),
univariate(g, k)))
@@ -238,9 +270,16 @@
has?(operator k, ALGOP) => alglfextint(f, k, l, x, g)
unkextint(f, x, g)
+@
+This is part of the fix for bug 100. Line 2 of this function used to read:
+\begin{verbatim}
+ symbolIfCan(k := kmax(l := union(l, vark(lu, x)))) case SE =>
+\end{verbatim}
+See the above discussion for why this causes an infinite loop.
+<<package INTEF ElementaryIntegration>>=
lflimitedint(f, x, lu) ==
empty?(l := varselect(kernels f, x)) => [x::F * f, empty()]
- symbolIfCan(k := kmax(l := union(l, vark(lu, x)))) case SE =>
+ symbolIfCan(k := kmax(l)) case SE =>
map(multivariate(#1, k), limitedint(univariate(f, k),
[univariate(u, k) for u in lu]))
is?(k, "exp"::SE) => explimint(f, x, k, lu)
Modified: trunk/axiom/src/input/Makefile.pamphlet
===================================================================
--- trunk/axiom/src/input/Makefile.pamphlet 2007-09-15 16:38:22 UTC (rev 747)
+++ trunk/axiom/src/input/Makefile.pamphlet 2007-09-16 03:12:14 UTC (rev 748)
@@ -287,7 +287,7 @@
arrows.regress assign.regress atansqrt.regress \
asec.regress bags.regress bbtree.regress \
binary.regress bop.regress bstree.regress bouquet.regress \
- bug10069.regress \
+ bug100.regress bug10069.regress \
bugs.regress bug10312.regress bug6357.regress bug9057.regress \
calcprob.regress \
calculus2.regress calculus.regress cardinal.regress card.regress \
@@ -502,7 +502,8 @@
${OUT}/bags.input ${OUT}/bbtree.input ${OUT}/bern.input \
${OUT}/bernpoly.input ${OUT}/binary.input ${OUT}/bop.input \
${OUT}/bouquet.input ${OUT}/bstree.input ${OUT}/bug6357.input \
- ${OUT}/bug9057.input ${OUT}/bug10069.input ${OUT}/bug10312.input \
+ ${OUT}/bug9057.input ${OUT}/bug100.input \
+ ${OUT}/bug10069.input ${OUT}/bug10312.input \
${OUT}/calcprob.input ${OUT}/calculus.input \
${OUT}/cardinal.input ${OUT}/card.input ${OUT}/carten.input \
${OUT}/cclass.input ${OUT}/cdraw.input ${OUT}/char.input \
@@ -677,6 +678,7 @@
${DOC}/bernpoly.input.dvi ${DOC}/binary.input.dvi \
${DOC}/bop.input.dvi ${DOC}/bouquet.input.dvi \
${DOC}/bstree.input.dvi ${DOC}/bug10069.input.dvi \
+ ${DOC}/bug100.input.dvi \
${DOC}/bug10312.input.dvi ${DOC}/bug6357.input.dvi \
${DOC}/bug9057.input.dvi ${DOC}/bugs.input.dvi \
${DOC}/c02aff.input.dvi ${DOC}/c02agf.input.dvi \
Added: trunk/axiom/src/input/bug100.input.pamphlet
===================================================================
--- trunk/axiom/src/input/bug100.input.pamphlet (rev 0)
+++ trunk/axiom/src/input/bug100.input.pamphlet 2007-09-16 03:12:14 UTC (rev 748)
@@ -0,0 +1,80 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/input bug100.input}
+\author{Timothy Daly}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+@
+This code used to cause an infinite loop that eventually killed Axiom.
+The bug exists in the {\tt intef.spad}, in ElementaryIntegration (INTEF),
+in the routines {\tt lfextendedint} and {\tt lflimitedint}.
+
+The input
+\begin{verbatim}
+ integrate((z^a+1)^b,z)
+\end{verbatim}
+triggers the bug because of line 2 of each function.
+Line 2 of these functions used to read:
+\begin{verbatim}
+ symbolIfCan(k := kmax(l := union(l, varselect(kernels g, x))))
+\end{verbatim}
+
+The loop occurs when the call to union causes
+\begin{verbatim}
+ a log(z)
+ %e
+\end{verbatim}
+to get added to the list every time. This gives the argument to kmax
+\begin{verbatim}
+ a log(z)
+ arg1= [z,%e ]
+\end{verbatim}
+and the result being
+\begin{verbatim}
+ a log(z)
+ %e
+\end{verbatim}
+We keep coming back to process this term, which ends up
+putting the same term back on the list and we loop.
+Waldek's solution is to remove the union call.
+
+<<*>>=
+)spool bug100.output
+)set message test on
+)set message auto off
+)clear all
+
+--S 1 of 1
+integrate((z^a+1)^b,z)
+--R
+--R z
+--R ++ a b
+--R (1) | (%I + 1) d%I
+--R ++
+--R Type: Union(Expression Integer,...)
+--E 1
+@
+If we substitute z for %I we get:
+\begin{verbatim}
+ z
+ ++ a b
+ | (z + 1) dz
+ ++
+
+\end{verbatim}
+which is the original input integral.
+<<*>>=
+)spool
+)lisp (bye)
+
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <gi...@ax...> - 2007-09-15 16:38:31
|
changelog | 1 +
src/algebra/carten.spad.pamphlet | 96 +++++++++++++++++++-------------------
2 files changed, 49 insertions(+), 48 deletions(-)
New commits:
commit f47f592674e1a73ee0b07042c483c48c2033a79e
Author: Tim Daly <root@localhost.localdomain>
Date: Thu Sep 13 15:28:11 2007 -0400
20070915 tpd src/algebra/carten minor edit for regression cleanup
|
|
From: <da...@us...> - 2007-09-15 16:38:31
|
Revision: 747
http://axiom.svn.sourceforge.net/axiom/?rev=747&view=rev
Author: daly
Date: 2007-09-15 09:38:22 -0700 (Sat, 15 Sep 2007)
Log Message:
-----------
20070915 tpd src/algebra/carten minor edit for regression cleanup
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/algebra/carten.spad.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-14 18:26:42 UTC (rev 746)
+++ trunk/axiom/changelog 2007-09-15 16:38:22 UTC (rev 747)
@@ -1,3 +1,4 @@
+20070915 tpd src/algebra/carten minor edit for regression cleanup
20070914 wxh src/hyper/hyper fix bad bracing of )hd change
20070914 tpd src/algebra/fraction.spad remove double )spool command
20070914 tpd src/algebra/kl.spad remove double )spool command
Modified: trunk/axiom/src/algebra/carten.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/carten.spad.pamphlet 2007-09-14 18:26:42 UTC (rev 746)
+++ trunk/axiom/src/algebra/carten.spad.pamphlet 2007-09-15 16:38:22 UTC (rev 747)
@@ -108,7 +108,7 @@
)set message test on
)set message auto off
)clear all
---S 1
+--S 1 of 48
CT := CARTEN(i0 := 1, 2, Integer)
--R
--R
@@ -116,7 +116,7 @@
--R Type: Domain
--E 1
---S 2
+--S 2 of 48
t0: CT := 8
--R
--R
@@ -124,7 +124,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 2
---S 3
+--S 3 of 48
rank t0
--R
--R
@@ -132,7 +132,7 @@
--R Type: NonNegativeInteger
--E 3
---S 4
+--S 4 of 48
v: DirectProduct(2, Integer) := directProduct [3,4]
--R
--R
@@ -140,7 +140,7 @@
--R Type: DirectProduct(2,Integer)
--E 4
---S 5
+--S 5 of 48
Tv: CT := v
--R
--R
@@ -148,7 +148,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 5
---S 6
+--S 6 of 48
m: SquareMatrix(2, Integer) := matrix [ [1,2],[4,5] ]
--R
--R
@@ -158,7 +158,7 @@
--R Type: SquareMatrix(2,Integer)
--E 6
---S 7
+--S 7 of 48
Tm: CT := m
--R
--R
@@ -168,7 +168,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 7
---S 8
+--S 8 of 48
n: SquareMatrix(2, Integer) := matrix [ [2,3],[0,1] ]
--R
--R
@@ -178,7 +178,7 @@
--R Type: SquareMatrix(2,Integer)
--E 8
---S 9
+--S 9 of 48
Tn: CT := n
--R
--R
@@ -188,7 +188,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 9
---S 10
+--S 10 of 48
t1: CT := [2, 3]
--R
--R
@@ -196,7 +196,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 10
---S 11
+--S 11 of 48
rank t1
--R
--R
@@ -204,7 +204,7 @@
--R Type: PositiveInteger
--E 11
---S 12
+--S 12 of 48
t2: CT := [t1, t1]
--R
--R
@@ -214,7 +214,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 12
---S 13
+--S 13 of 48
t3: CT := [t2, t2]
--R
--R
@@ -224,7 +224,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 13
---S 14
+--S 14 of 48
tt: CT := [t3, t3]; tt := [tt, tt]
--R
--R
@@ -238,7 +238,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 14
---S 15
+--S 15 of 48
rank tt
--R
--R
@@ -246,7 +246,7 @@
--R Type: PositiveInteger
--E 15
---S 16
+--S 16 of 48
Tmn := product(Tm, Tn)
--R
--R
@@ -260,7 +260,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 16
---S 17
+--S 17 of 48
Tmv := contract(Tm,2,Tv,1)
--R
--R
@@ -268,7 +268,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 17
---S 18
+--S 18 of 48
Tm*Tv
--R
--R
@@ -276,7 +276,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 18
---S 19
+--S 19 of 48
Tmv = m * v
--R
--R
@@ -284,7 +284,7 @@
--R Type: Equation CartesianTensor(1,2,Integer)
--E 19
---S 20
+--S 20 of 48
t0()
--R
--R
@@ -292,7 +292,7 @@
--R Type: PositiveInteger
--E 20
---S 21
+--S 21 of 48
t1(1+1)
--R
--R
@@ -300,7 +300,7 @@
--R Type: PositiveInteger
--E 21
---S 22
+--S 22 of 48
t2(2,1)
--R
--R
@@ -308,7 +308,7 @@
--R Type: PositiveInteger
--E 22
---S 23
+--S 23 of 48
t3(2,1,2)
--R
--R
@@ -316,7 +316,7 @@
--R Type: PositiveInteger
--E 23
---S 24
+--S 24 of 48
Tmn(2,1,2,1)
--R
--R
@@ -324,7 +324,7 @@
--R Type: NonNegativeInteger
--E 24
---S 25
+--S 25 of 48
t0[]
--R
--R
@@ -332,7 +332,7 @@
--R Type: PositiveInteger
--E 25
---S 26
+--S 26 of 48
t1[2]
--R
--R
@@ -340,7 +340,7 @@
--R Type: PositiveInteger
--E 26
---S 27
+--S 27 of 48
t2[2,1]
--R
--R
@@ -348,7 +348,7 @@
--R Type: PositiveInteger
--E 27
---S 28
+--S 28 of 48
t3[2,1,2]
--R
--R
@@ -356,7 +356,7 @@
--R Type: PositiveInteger
--E 28
---S 29
+--S 29 of 48
Tmn[2,1,2,1]
--R
--R
@@ -364,7 +364,7 @@
--R Type: NonNegativeInteger
--E 29
---S 30
+--S 30 of 48
cTmn := contract(Tmn,1,2)
--R
--R
@@ -374,7 +374,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 30
---S 31
+--S 31 of 48
trace(m) * n
--R
--R
@@ -384,7 +384,7 @@
--R Type: SquareMatrix(2,Integer)
--E 31
---S 32
+--S 32 of 48
contract(Tmn,1,2) = trace(m) * n
--R
--R
@@ -394,7 +394,7 @@
--R Type: Equation CartesianTensor(1,2,Integer)
--E 32
---S 33
+--S 33 of 48
contract(Tmn,1,3) = transpose(m) * n
--R
--R
@@ -404,7 +404,7 @@
--R Type: Equation CartesianTensor(1,2,Integer)
--E 33
---S 34
+--S 34 of 48
contract(Tmn,1,4) = transpose(m) * transpose(n)
--R
--R
@@ -414,7 +414,7 @@
--R Type: Equation CartesianTensor(1,2,Integer)
--E 34
---S 35
+--S 35 of 48
contract(Tmn,2,3) = m * n
--R
--R
@@ -424,7 +424,7 @@
--R Type: Equation CartesianTensor(1,2,Integer)
--E 35
---S 36
+--S 36 of 48
contract(Tmn,2,4) = m * transpose(n)
--R
--R
@@ -434,7 +434,7 @@
--R Type: Equation CartesianTensor(1,2,Integer)
--E 36
---S 37
+--S 37 of 48
contract(Tmn,3,4) = trace(n) * m
--R
--R
@@ -444,7 +444,7 @@
--R Type: Equation CartesianTensor(1,2,Integer)
--E 37
---S 38
+--S 38 of 48
tTmn := transpose(Tmn,1,3)
--R
--R
@@ -458,7 +458,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 38
---S 39
+--S 39 of 48
transpose Tmn
--R
--R
@@ -472,7 +472,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 39
---S 40
+--S 40 of 48
transpose Tm = transpose m
--R
--R
@@ -482,7 +482,7 @@
--R Type: Equation CartesianTensor(1,2,Integer)
--E 40
---S 41
+--S 41 of 48
rTmn := reindex(Tmn, [1,4,2,3])
--R
--R
@@ -496,7 +496,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 41
---S 42
+--S 42 of 48
tt := transpose(Tm)*Tn - Tn*transpose(Tm)
--R
--R
@@ -506,7 +506,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 42
---S 43
+--S 43 of 48
Tv*(tt+Tn)
--R
--R
@@ -514,7 +514,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 43
---S 44
+--S 44 of 48
reindex(product(Tn,Tn),[4,3,2,1])+3*Tn*product(Tm,Tm)
--R
--R
@@ -528,7 +528,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 44
---S 45
+--S 45 of 48
delta: CT := kroneckerDelta()
--R
--R
@@ -538,7 +538,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 45
---S 46
+--S 46 of 48
contract(Tmn, 2, delta, 1) = reindex(Tmn, [1,3,4,2])
--R
--R
@@ -552,7 +552,7 @@
--R Type: Equation CartesianTensor(1,2,Integer)
--E 46
---S 47
+--S 47 of 48
epsilon:CT := leviCivitaSymbol()
--R
--R
@@ -562,7 +562,7 @@
--R Type: CartesianTensor(1,2,Integer)
--E 47
---S 48
+--S 48 of 48
contract(epsilon*Tm*epsilon, 1,2) = 2 * determinant m
--R
--R
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <gi...@ax...> - 2007-09-14 18:27:04
|
changelog | 1 +
src/hyper/hyper.pamphlet | 4 ++--
2 files changed, 3 insertions(+), 2 deletions(-)
New commits:
commit 0df0ac99dcf617f9c2ac041a77357d5fb2d38734
Author: Tim Daly <root@localhost.localdomain>
Date: Wed Sep 12 19:16:23 2007 -0400
20070914 wxh src/hyper/hyper fix bad bracing of )hd change
|
|
From: <da...@us...> - 2007-09-14 18:26:39
|
Revision: 746
http://axiom.svn.sourceforge.net/axiom/?rev=746&view=rev
Author: daly
Date: 2007-09-14 11:26:42 -0700 (Fri, 14 Sep 2007)
Log Message:
-----------
20070914 wxh src/hyper/hyper fix bad bracing of )hd change
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/hyper/hyper.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-14 15:04:42 UTC (rev 745)
+++ trunk/axiom/changelog 2007-09-14 18:26:42 UTC (rev 746)
@@ -1,3 +1,4 @@
+20070914 wxh src/hyper/hyper fix bad bracing of )hd change
20070914 tpd src/algebra/fraction.spad remove double )spool command
20070914 tpd src/algebra/kl.spad remove double )spool command
20070914 tpd src/algebra/lindep.spad remove double )spool command
Modified: trunk/axiom/src/hyper/hyper.pamphlet
===================================================================
--- trunk/axiom/src/hyper/hyper.pamphlet 2007-09-14 15:04:42 UTC (rev 745)
+++ trunk/axiom/src/hyper/hyper.pamphlet 2007-09-14 18:26:42 UTC (rev 746)
@@ -948,13 +948,13 @@
fprintf(stderr, "(HyperDoc) Warning: Not connected to AXIOM Server!\n");
MenuServerOpened = 0;
}
- else
+ else {
/* In order to allow hyperdoc restarts from the console we clean up
* the socket on exit */
atexit(&clean_socket);
MenuServerOpened = 1;
+ }
-
/*
* If I have opened the MenuServer socket, then I should also try to open
* the SpadServer socket, so I can send stuff right to SPAD.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <gi...@ax...> - 2007-09-14 15:11:25
|
changelog | 4 ++++
src/algebra/fraction.spad.pamphlet | 2 --
src/algebra/kl.spad.pamphlet | 2 --
src/algebra/lindep.spad.pamphlet | 2 --
src/algebra/radix.spad.pamphlet | 2 --
5 files changed, 4 insertions(+), 8 deletions(-)
New commits:
commit 3aa984ac83988f404ff83bc0f9c21053f30285e6
Author: Tim Daly <root@localhost.localdomain>
Date: Wed Sep 12 16:23:25 2007 -0400
20070914 tpd src/algebra/fraction.spad remove double )spool command
20070914 tpd src/algebra/kl.spad remove double )spool command
20070914 tpd src/algebra/lindep.spad remove double )spool command
20070914 tpd src/algebra/radix.spad remove double )spool command
|
|
From: <da...@us...> - 2007-09-14 15:04:38
|
Revision: 745
http://axiom.svn.sourceforge.net/axiom/?rev=745&view=rev
Author: daly
Date: 2007-09-14 08:04:42 -0700 (Fri, 14 Sep 2007)
Log Message:
-----------
20070914 tpd src/algebra/fraction.spad remove double )spool command
20070914 tpd src/algebra/kl.spad remove double )spool command
20070914 tpd src/algebra/lindep.spad remove double )spool command
20070914 tpd src/algebra/radix.spad remove double )spool command
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/algebra/fraction.spad.pamphlet
trunk/axiom/src/algebra/kl.spad.pamphlet
trunk/axiom/src/algebra/lindep.spad.pamphlet
trunk/axiom/src/algebra/radix.spad.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-14 14:22:14 UTC (rev 744)
+++ trunk/axiom/changelog 2007-09-14 15:04:42 UTC (rev 745)
@@ -1,3 +1,7 @@
+20070914 tpd src/algebra/fraction.spad remove double )spool command
+20070914 tpd src/algebra/kl.spad remove double )spool command
+20070914 tpd src/algebra/lindep.spad remove double )spool command
+20070914 tpd src/algebra/radix.spad remove double )spool command
20070914 jap adapt changes for )hd restart to Axiom sources
20070914 jap Jose Alfredo Portes <doy...@gm...>
20070914 wxh src/sman/bookvol6 enable restart of hyperdoc with )hd
Modified: trunk/axiom/src/algebra/fraction.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/fraction.spad.pamphlet 2007-09-14 14:22:14 UTC (rev 744)
+++ trunk/axiom/src/algebra/fraction.spad.pamphlet 2007-09-14 15:04:42 UTC (rev 745)
@@ -513,8 +513,6 @@
--R Type: Fraction Complex Integer
--E 12
)spool
-
-)spool
)lisp (bye)
@
<<Fraction.help>>=
Modified: trunk/axiom/src/algebra/kl.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/kl.spad.pamphlet 2007-09-14 14:22:14 UTC (rev 744)
+++ trunk/axiom/src/algebra/kl.spad.pamphlet 2007-09-14 15:04:42 UTC (rev 745)
@@ -304,8 +304,6 @@
--R Type: List Expression Integer
--E 19
)spool
-
-)spool
)lisp (bye)
@
<<Kernel.help>>=
Modified: trunk/axiom/src/algebra/lindep.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/lindep.spad.pamphlet 2007-09-14 14:22:14 UTC (rev 744)
+++ trunk/axiom/src/algebra/lindep.spad.pamphlet 2007-09-14 15:04:42 UTC (rev 745)
@@ -164,8 +164,6 @@
--R Type: Union(Vector Fraction Integer,...)
--E 8
)spool
-
-)spool
)lisp (bye)
@
<<IntegerLinearDependence.help>>=
Modified: trunk/axiom/src/algebra/radix.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/radix.spad.pamphlet 2007-09-14 14:22:14 UTC (rev 744)
+++ trunk/axiom/src/algebra/radix.spad.pamphlet 2007-09-14 15:04:42 UTC (rev 745)
@@ -988,8 +988,6 @@
--R Type: Polynomial HexadecimalExpansion
--E 7
)spool
-
-)spool
)lisp (bye)
@
<<HexadecimalExpansion.help>>=
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <gi...@ax...> - 2007-09-14 14:26:37
|
changelog | 5 +++++
src/hyper/hyper.pamphlet | 9 ++++++---
src/include/sman.h1 | 1 +
src/sman/bookvol6.pamphlet | 30 ++++++++++++++++++++++++------
4 files changed, 36 insertions(+), 9 deletions(-)
New commits:
commit 3d0f0cb9e97a1d2a244f0cc72d322e642070b037
Author: Tim Daly <root@localhost.localdomain>
Date: Wed Sep 12 15:47:31 2007 -0400
allow )hd restart
|
|
From: <da...@us...> - 2007-09-14 14:22:10
|
Revision: 744
http://axiom.svn.sourceforge.net/axiom/?rev=744&view=rev
Author: daly
Date: 2007-09-14 07:22:14 -0700 (Fri, 14 Sep 2007)
Log Message:
-----------
20070914 jap adapt changes for )hd restart to Axiom sources
20070914 jap Jose Alfredo Portes <doy...@gm...>
20070914 wxh src/sman/bookvol6 enable restart of hyperdoc with )hd
20070914 wxh src/include/sman.h1 enable restart of hyperdoc with )hd
20070914 wxh src/hyper/hyper enable restart of hyperdoc with )hd
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/hyper/hyper.pamphlet
trunk/axiom/src/include/sman.h1
trunk/axiom/src/sman/bookvol6.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-14 03:25:54 UTC (rev 743)
+++ trunk/axiom/changelog 2007-09-14 14:22:14 UTC (rev 744)
@@ -1,3 +1,8 @@
+20070914 jap adapt changes for )hd restart to Axiom sources
+20070914 jap Jose Alfredo Portes <doy...@gm...>
+20070914 wxh src/sman/bookvol6 enable restart of hyperdoc with )hd
+20070914 wxh src/include/sman.h1 enable restart of hyperdoc with )hd
+20070914 wxh src/hyper/hyper enable restart of hyperdoc with )hd
20070913 tpd src/input/Makefile schaum1.input added
20070913 tpd src/input/schaum1.input added
20070909 tpd src/algebra/newton.spad included in fffg.spad
Modified: trunk/axiom/src/hyper/hyper.pamphlet
===================================================================
--- trunk/axiom/src/hyper/hyper.pamphlet 2007-09-14 03:25:54 UTC (rev 743)
+++ trunk/axiom/src/hyper/hyper.pamphlet 2007-09-14 14:22:14 UTC (rev 744)
@@ -945,11 +945,14 @@
*/
if (open_server(MenuServerName) == -2) {
- fprintf(stderr, "(HyperDoc) Warning: Not connected to AXIOM Server!\n");
- MenuServerOpened = 0;
+ fprintf(stderr, "(HyperDoc) Warning: Not connected to AXIOM Server!\n");
+ MenuServerOpened = 0;
}
else
- MenuServerOpened = 1;
+ /* In order to allow hyperdoc restarts from the console we clean up
+ * the socket on exit */
+ atexit(&clean_socket);
+ MenuServerOpened = 1;
/*
Modified: trunk/axiom/src/include/sman.h1
===================================================================
--- trunk/axiom/src/include/sman.h1 2007-09-14 03:25:54 UTC (rev 743)
+++ trunk/axiom/src/include/sman.h1 2007-09-14 14:22:14 UTC (rev 744)
@@ -25,6 +25,7 @@
static void fork_Axiom(void);
static void start_the_Axiom(char * * envp);
static void clean_up_sockets(void);
+static void clean_hypertex_socket(void);
static void read_from_spad_io(int ptcNum);
static void read_from_manager(int ptcNum);
static void manage_spad_io(int ptcNum);
Modified: trunk/axiom/src/sman/bookvol6.pamphlet
===================================================================
--- trunk/axiom/src/sman/bookvol6.pamphlet 2007-09-14 03:25:54 UTC (rev 743)
+++ trunk/axiom/src/sman/bookvol6.pamphlet 2007-09-14 14:22:14 UTC (rev 744)
@@ -812,6 +812,8 @@
#define NadaDelShitsky 2
/* When a process dies start it up again */
#define DoItAgain 3
+/* When hypertex dies, clean its socket */
+#define CleanHypertexSocket 4
typedef struct spad_proc {
int proc_id; /* process id of child */
@@ -1405,7 +1407,8 @@
sprintf(prog, "%s -k -rv %s", HypertexProgram, VerifyRecordFile);
spawn_of_hell(prog, NadaDelShitsky);
}
- else spawn_of_hell(HypertexProgram, NadaDelShitsky);
+ /* If we restart hyperdoc from the axiom command prompt */
+ else spawn_of_hell(HypertexProgram, CleanHypertexSocket);
}
@
@@ -1515,8 +1518,18 @@
@
\subsection{clean\_up\_sockets}
+In order to be able to restart hyperdoc from the axiom command prompt
+we need to remove the socket for this server.
<<sman.cleanupsockets>>=
static void
+clean_hypertex_socket(void)
+{
+ char name[256];
+ sprintf(name, "%s%d", MenuServerName, server_num);
+ unlink(name);
+}
+
+static void
clean_up_sockets(void)
{
char name[256];
@@ -1526,8 +1539,7 @@
unlink(name);
sprintf(name, "%s%d", SessionIOName, server_num);
unlink(name);
- sprintf(name, "%s%d", MenuServerName, server_num);
- unlink(name);
+ clean_hypertex_socket();
}
@
@@ -1700,11 +1712,12 @@
stat = 0;
dead_baby = wait(&stat);
/* Check the value of dead_baby, since wait may have returned
- a pid but subsequently we have received a signal. Yeuch! */
+ a pid but subsequently we have received a signal. Yeuch!
+ In order to restart hyperdoc from the axiom command prompt
+ we no longer call clean_up_terminal */
if (dead_baby == -1 && death_signal) {
kill_all_children();
clean_up_sockets();
- clean_up_terminal();
sleep(2);
exit(0);
}
@@ -1728,10 +1741,12 @@
continue;
}
switch(proc->death_action) {
+ /* In order to restart hyperdoc from the axiom command prompt
+ we no longer call clean_up_terminal. Instead we've added a
+ case to just clean up the socket. */
case Die:
kill_all_children();
clean_up_sockets();
- clean_up_terminal();
sleep(2);
exit(0);
case NadaDelShitsky:
@@ -1739,6 +1754,9 @@
case DoItAgain:
spawn_of_hell(proc->command, DoItAgain);
break;
+ case CleanHypertexSocket:
+ clean_hypertex_socket();
+ break;
}
}
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <da...@us...> - 2007-09-14 03:25:55
|
Revision: 743
http://axiom.svn.sourceforge.net/axiom/?rev=743&view=rev
Author: daly
Date: 2007-09-13 20:25:54 -0700 (Thu, 13 Sep 2007)
Log Message:
-----------
20070913 tpd src/input/Makefile schaum1.input added
20070913 tpd src/input/schaum1.input added
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/input/Makefile.pamphlet
Added Paths:
-----------
trunk/axiom/src/input/schaum1.input.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-11 05:08:07 UTC (rev 742)
+++ trunk/axiom/changelog 2007-09-14 03:25:54 UTC (rev 743)
@@ -1,3 +1,5 @@
+20070913 tpd src/input/Makefile schaum1.input added
+20070913 tpd src/input/schaum1.input added
20070909 tpd src/algebra/newton.spad included in fffg.spad
20070909 tpd src/algebra/Makefile remove newton.spad (duplicate)
20070907 tpd src/algebra/acplot.spad fix PlaneAlgebraicCurvePlot.help NOISE
Modified: trunk/axiom/src/input/Makefile.pamphlet
===================================================================
--- trunk/axiom/src/input/Makefile.pamphlet 2007-09-11 05:08:07 UTC (rev 742)
+++ trunk/axiom/src/input/Makefile.pamphlet 2007-09-14 03:25:54 UTC (rev 743)
@@ -345,6 +345,7 @@
r21bugsbig.regress r21bugs.regress radff.regress radix.regress \
realclos.regress reclos.regress repa6.regress robidoux.regress \
roman.regress roots.regress ruleset.regress rules.regress \
+ schaum1.regress \
scherk.regress scope.regress segbind.regress seg.regress \
series2.regress series.regress sersolve.regress set.regress \
sincosex.regress sint.regress skew.regress slowint.regress \
@@ -601,7 +602,8 @@
${OUT}/radff.input ${OUT}/radix.input ${OUT}/realclos.input \
${OUT}/reclos.input ${OUT}/regset.input \
${OUT}/robidoux.input ${OUT}/roman.input ${OUT}/roots.input \
- ${OUT}/ruleset.input ${OUT}/rules.input ${OUT}/saddle.input \
+ ${OUT}/ruleset.input ${OUT}/rules.input ${OUT}/schaum1.input \
+ ${OUT}/saddle.input \
${OUT}/scherk.input ${OUT}/scope.input \
${OUT}/segbind.input ${OUT}/seg.input ${OUT}/series2.input \
${OUT}/series.input ${OUT}/sersolve.input ${OUT}/set.input \
@@ -879,6 +881,7 @@
${DOC}/robidoux.input.dvi ${DOC}/roman.input.dvi \
${DOC}/romnum.as.dvi ${DOC}/roots.input.dvi \
${DOC}/ruleset.input.dvi ${DOC}/rules.input.dvi \
+ ${DOC}/schaum1.input.dvi \
${DOC}/s01eaf.input.dvi ${DOC}/s13aaf.input.dvi \
${DOC}/s13acf.input.dvi ${DOC}/s13adf.input.dvi \
${DOC}/s14aaf.input.dvi ${DOC}/s14abf.input.dvi \
Added: trunk/axiom/src/input/schaum1.input.pamphlet
===================================================================
--- trunk/axiom/src/input/schaum1.input.pamphlet (rev 0)
+++ trunk/axiom/src/input/schaum1.input.pamphlet 2007-09-14 03:25:54 UTC (rev 743)
@@ -0,0 +1,1265 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/input schaum1.input}
+\author{Timothy Daly}
+\maketitle
+\eject
+\tableofcontents
+\eject
+\section{\cite{1}:14.59~~~~~$\displaystyle\int{\frac{dx}{ax+b}~dx}$}
+$$\int{\frac{dx}{ax+b}~dx}==\frac{1}{a}~\ln(ax+b)$$
+<<*>>=
+)spool schaum1.output
+)set message test on
+)set message auto off
+)clear all
+
+--S 1
+integrate(1/(a*x+b),x)
+--R
+--R log(a x + b)
+--R (1) ------------
+--R a
+--R Type: Union(Expression Integer,...)
+--E 1
+@
+\section{\cite{1}:14.60~~~~~$\displaystyle\int{\frac{x~dx}{ax+b}}$}
+$$\int{\frac{x~dx}{ax+b}}=\frac{x}{a}-\frac{b}{a^2}~\ln(ax+b)$$
+<<*>>=
+)clear all
+
+--S 2
+integrate(x/(a*x+b),x)
+--R
+--R
+--R - b log(a x + b) + a x
+--R (1) ----------------------
+--R 2
+--R a
+--R Type: Union(Expression Integer,...)
+--E 2
+@
+\section{\cite{1}:14.61~~~~~$\displaystyle\int{\frac{x^2~dx}{ax+b}}$}
+$$\int{\frac{x^2~dx}{ax+b}}=
+\frac{(ax+b)^2}{2a^3}-\frac{2b(ax+b)}{a^3}+\frac{b^2}{a^3}~\ln(ax+b)$$
+<<*>>=
+)clear all
+
+--S 3
+nn:=integrate(x^2/(a*x+b),x)
+--R
+--R 2 2 2
+--R 2b log(a x + b) + a x - 2a b x
+--R (1) -------------------------------
+--R 3
+--R 2a
+--R Type: Union(Expression Integer,...)
+--E 3
+@
+To see that these are the same answers we put the prior result over
+a common fraction:
+<<*>>=
+--S 4
+mm:=((a*x+b)^2-2*2*b*(a*x+b)+2*b^2*log(a*x+b))/(2*a^3)
+--R
+--R 2 2 2 2
+--R 2b log(a x + b) + a x - 2a b x - 3b
+--R (2) -------------------------------------
+--R 3
+--R 2a
+--R Type: Expression Integer
+--E 4
+@
+and we take their difference:
+<<*>>=
+--S 5
+pp:=mm-nn
+--R
+--R 2
+--R 3b
+--R (3) - ---
+--R 3
+--R 2a
+--R Type: Expression Integer
+--E 5
+@
+which is a constant with respect to x, and thus the constant C.
+<<*>>=
+--S 6
+D(pp,x)
+--R
+--R (4) 0
+--R Type: Expression Integer
+--E 6
+@
+Alternatively we can differentiate the answers with respect to x:
+<<*>>=
+--S 7
+D(nn,x)
+--R
+--R 2
+--R x
+--R (5) -------
+--R a x + b
+--R Type: Expression Integer
+--E 7
+@
+<<*>>=
+--S 8
+D(mm,x)
+--R
+--R 2
+--R x
+--R (6) -------
+--R a x + b
+--R Type: Expression Integer
+--E 8
+@
+and see that they are indeed the same.
+
+\section{\cite{1}:14.62~~~~~$\displaystyle\int{\frac{x^3~dx}{ax+b}}$}
+$$\int{\frac{x^3~dx}{ax+b}}=
+\frac{(ax+b)^3}{3a^4}-\frac{3b(ax+b)^2}{2a^4}+
+\frac{3b^2(ax+b)}{a^4}-\frac{b^3}{a^4}~\ln(ax+b)$$
+<<*>>=
+)clear all
+
+--S 9
+aa:=integrate(x^3/(a*x+b),x)
+--R
+--R 3 3 3 2 2 2
+--R - 6b log(a x + b) + 2a x - 3a b x + 6a b x
+--R (1) --------------------------------------------
+--R 4
+--R 6a
+--R Type: Union(Expression Integer,...)
+--E 9
+@
+and the book expression is:
+<<*>>=
+--S 10
+bb:=(a*x+b)^3/(3*a^4)-(3*b*(a*x+b)^2)/(2*a^4)+(3*b^2*(a*x+b))/a^4-(b^3/a^4)*log(a*x+b)
+--R
+--R 3 3 3 2 2 2 3
+--R - 6b log(a x + b) + 2a x - 3a b x + 6a b x + 11b
+--R (2) ---------------------------------------------------
+--R 4
+--R 6a
+--R Type: Expression Integer
+--E 10
+@
+
+The difference is a constant with respect to x:
+<<*>>=
+--S 11
+aa-bb
+--R
+--R 3
+--R 11b
+--R (3) - ----
+--R 4
+--R 6a
+--R Type: Expression Integer
+--E 11
+@
+
+If we differentiate each expression we see
+<<*>>=
+--S 12
+cc:=D(aa,x)
+--R
+--R 3
+--R x
+--R (4) -------
+--R a x + b
+--R Type: Expression Integer
+--E 12
+@
+<<*>>=
+--S 13
+dd:=D(bb,x)
+--R
+--R 3
+--R x
+--R (5) -------
+--R a x + b
+--R Type: Expression Integer
+--E 13
+@
+<<*>>=
+--S 14
+cc-dd
+--R
+--R (6) 0
+--R Type: Expression Integer
+--E 14
+@
+
+\section{\cite{1}:14.63~~~~~$\displaystyle\int{\frac{dx}{x~(ax+b)}}$}
+$$\int{\frac{dx}{x~(ax+b)}}=\frac{1}{b}~\ln\left(\frac{x}{ax+b}\right)$$
+<<*>>=
+)clear all
+
+--S 15
+ff:=integrate(1/(x*(a*x+b)),x)
+--R
+--R - log(a x + b) + log(x)
+--R (1) -----------------------
+--R b
+--R Type: Union(Expression Integer,...)
+--E 15
+@
+but we know that $$\log(a)-\log(b)=\log(\frac{a}{b})$$
+
+We can express this fact as a rule:
+<<*>>=
+--S 16
+logdiv:=rule(log(a)-log(b) == log(a/b))
+--R
+--R a
+--I (2) - log(b) + log(a) + %I == log(-) + %I
+--R b
+--R Type: RewriteRule(Integer,Integer,Expression Integer)
+--E 16
+@
+and use this rule to rewrite the logs into divisions:
+<<*>>=
+--S 17
+logdiv ff
+--R
+--R x
+--R log(-------)
+--R a x + b
+--R (3) ------------
+--R b
+--R Type: Expression Integer
+--E 17
+@
+so we can see the equivalence directly.
+
+\section{\cite{1}:14.64~~~~~$\displaystyle\int{\frac{dx}{x^2~(ax+b)}}$}
+$$\int{\frac{dx}{x^2~(ax+b)}}=
+-\frac{1}{bx}+\frac{a}{b^2}~\ln\left(\frac{ax+b}{x}\right)$$
+<<*>>=
+)clear all
+
+--S 18
+aa:=integrate(1/(x^2*(a*x+b)),x)
+--R
+--R a x log(a x + b) - a x log(x) - b
+--R (1) ---------------------------------
+--R 2
+--R b x
+--R Type: Union(Expression Integer,...)
+--E 18
+@
+
+The original form given in the book expands to:
+<<*>>=
+--S 19
+bb:=-1/(b*x)+a/b^2*log((a*x+b)/x)
+--R
+--R a x + b
+--R a x log(-------) - b
+--R x
+--R (2) --------------------
+--R 2
+--R b x
+--R Type: Expression Integer
+--E 19
+@
+
+We can define the following rule to expand log forms:
+<<*>>=
+--S 20
+divlog:=rule(log(a/b) == log(a) - log(b))
+--R
+--R a
+--R (3) log(-) == - log(b) + log(a)
+--R b
+--R Type: RewriteRule(Integer,Integer,Expression Integer)
+--E 20
+@
+and apply it to the book form:
+<<*>>=
+--S 21
+cc:= divlog bb
+--R
+--R a x log(a x + b) - a x log(x) - b
+--R (4) ---------------------------------
+--R 2
+--R b x
+--R Type: Expression Integer
+--E 21
+@
+and we can now see that the results are identical.
+<<*>>=
+--S 22
+aa-cc
+--R
+--R (5) 0
+--R Type: Expression Integer
+--E 22
+@
+
+\section{\cite{1}:14.65~~~~~$\displaystyle\int{\frac{dx}{x^3~(ax+b)}}$}
+$$\int{\frac{dx}{x^3~(ax+b)}}=
+\frac{2ax-b}{2b^2x^2}+\frac{a^2}{b^3}~\ln\left(\frac{x}{ax+b}\right)$$
+<<*>>=
+)clear all
+--S 23
+aa:=integrate(1/(x^3*(a*x+b)),x)
+--R
+--R 2 2 2 2 2
+--R - 2a x log(a x + b) + 2a x log(x) + 2a b x - b
+--R (1) -----------------------------------------------
+--R 3 2
+--R 2b x
+--R Type: Union(Expression Integer,...)
+--E 23
+@
+
+<<*>>=
+--S 24
+bb:=(2*a*x-b)/(2*b^2*x^2)+a^2/b^3*log(x/(a*x+b))
+--R
+--R 2 2 x 2
+--R 2a x log(-------) + 2a b x - b
+--R a x + b
+--R (2) -------------------------------
+--R 3 2
+--R 2b x
+--R Type: Expression Integer
+--E 24
+@
+
+<<*>>=
+--S 25
+divlog:=rule(log(a/b) == log(a) - log(b))
+--R
+--R a
+--R (3) log(-) == - log(b) + log(a)
+--R b
+--R Type: RewriteRule(Integer,Integer,Expression Integer)
+--E 25
+@
+
+<<*>>=
+--S 26
+cc:=divlog bb
+--R
+--R 2 2 2 2 2
+--R - 2a x log(a x + b) + 2a x log(x) + 2a b x - b
+--R (4) -----------------------------------------------
+--R 3 2
+--R 2b x
+--R Type: Expression Integer
+--E 26
+@
+
+<<*>>=
+--S 27
+cc-aa
+--R
+--R (5) 0
+--R Type: Expression Integer
+--E 27
+@
+
+\section{\cite{1}:14.66~~~~~$\displaystyle\int{\frac{dx}{(ax+b)^2}}$}
+$$\int{\frac{dx}{(ax+b)^2}}=\frac{-1}{a~(ax+b)}$$
+<<*>>=
+)clear all
+
+--S 28
+integrate(1/(a*x+b)^2,x)
+--R
+--R 1
+--R (1) - ---------
+--R 2
+--R a x + a b
+--R Type: Union(Expression Integer,...)
+--E 28
+@
+
+\section{\cite{1}:14.67~~~~~$\displaystyle\int{\frac{x~dx}{(ax+b)^2}}$}
+$$\int{\frac{x~dx}{(ax+b)^2}}=
+\frac{b}{a^2~(ax+b)}+\frac{1}{a^2}~\ln(ax+b)$$
+<<*>>=
+)clear all
+
+--S 29
+integrate(x/(a*x+b)^2,x)
+--R
+--R (a x + b)log(a x + b) + b
+--R (1) -------------------------
+--R 3 2
+--R a x + a b
+--R Type: Union(Expression Integer,...)
+--E 29
+@
+and the book form expands to:
+<<*>>=
+--S 30
+b/(a^2*(a*x+b))+(1/a^2)*log(a*x+b)
+--R
+--R (a x + b)log(a x + b) + b
+--R (2) -------------------------
+--R 3 2
+--R a x + a b
+--R Type: Expression Integer
+--E 30
+@
+
+\section{\cite{1}:14.68~~~~~$\displaystyle\int{\frac{x^2~dx}{(ax+b)^2}}$}
+$$\int{\frac{x^2~dx}{(ax+b)^2}}=
+\frac{ax+b}{a^3}-\frac{b^2}{a^3~(ax+b)}
+-\frac{2b}{a^3}~\ln(ax+b)$$
+<<*>>=
+)clear all
+
+--S 31
+aa:=integrate(x^2/(a*x+b)^2,x)
+--R
+--R 2 2 2 2
+--R (- 2a b x - 2b )log(a x + b) + a x + a b x - b
+--R (1) ------------------------------------------------
+--R 4 3
+--R a x + a b
+--R Type: Union(Expression Integer,...)
+--E 31
+@
+and the book expression expands into
+<<*>>=
+--S 32
+bb:=(a*x+b)/a^3-b^2/(a^3*(a*x+b))-((2*b)/a^3)*log(a*x+b)
+--R
+--R 2 2 2
+--R (- 2a b x - 2b )log(a x + b) + a x + 2a b x
+--R (2) --------------------------------------------
+--R 4 3
+--R a x + a b
+--R Type: Expression Integer
+--E 32
+@
+
+These two expressions differ by the constant
+<<*>>=
+--S 33
+aa-bb
+--R
+--R b
+--R (3) - --
+--R 3
+--R a
+--R Type: Expression Integer
+--E 33
+@
+
+These are the same integrands as can be shown by differentiation:
+<<*>>=
+--S 34
+D(aa,x)
+--R
+--R 2
+--R x
+--R (4) ------------------
+--R 2 2 2
+--R a x + 2a b x + b
+--R Type: Expression Integer
+--E 34
+@
+
+<<*>>=
+--S 35
+D(bb,x)
+--R
+--R 2
+--R x
+--R (5) ------------------
+--R 2 2 2
+--R a x + 2a b x + b
+--R Type: Expression Integer
+--E 35
+@
+
+\section{\cite{1}:14.69~~~~~$\displaystyle\int{\frac{x^3~dx}{(ax+b)^2}}$}
+$$\int{\frac{x^3~dx}{(ax+b)^2}}=
+\frac{(ax+b)^2}{2a^4}-\frac{3b(ax+b)}{a^4}+\frac{b^3}{a^4(ax+b)}
++\frac{3b^2}{a^4}~\ln(ax+b)$$
+<<*>>=
+)clear all
+
+--S 36
+aa:=integrate(x^3/(a*x+b)^2,x)
+--R
+--R 2 3 3 3 2 2 2 3
+--R (6a b x + 6b )log(a x + b) + a x - 3a b x - 4a b x + 2b
+--R (1) ----------------------------------------------------------
+--R 5 4
+--R 2a x + 2a b
+--R Type: Union(Expression Integer,...)
+--E 36
+@
+
+<<*>>=
+--S 37
+bb:=(a*x+b)^2/(2*a^4)-(3*b*(a*x+b))/a^4+b^3/(a^4*(a*x+b))+(3*b^2/a^4)*log(a*x+b)
+--R
+--R 2 3 3 3 2 2 2 3
+--R (6a b x + 6b )log(a x + b) + a x - 3a b x - 9a b x - 3b
+--R (2) ----------------------------------------------------------
+--R 5 4
+--R 2a x + 2a b
+--R Type: Expression Integer
+--E 37
+@
+
+<<*>>=
+--S 38
+aa-bb
+--R
+--R 2
+--R 5b
+--R (3) ---
+--R 4
+--R 2a
+--R Type: Expression Integer
+--E 38
+@
+
+<<*>>=
+--S 39
+cc:=D(aa,x)
+--R
+--R 3
+--R x
+--R (4) ------------------
+--R 2 2 2
+--R a x + 2a b x + b
+--R Type: Expression Integer
+--E 39
+@
+
+<<*>>=
+--S 40
+dd:=D(bb,x)
+--R
+--R 3
+--R x
+--R (5) ------------------
+--R 2 2 2
+--R a x + 2a b x + b
+--R Type: Expression Integer
+--E 40
+@
+
+<<*>>=
+--S 41
+cc-dd
+--R
+--R (6) 0
+--R Type: Expression Integer
+--E 41
+@
+
+\section{\cite{1}:14.70~~~~~$\displaystyle\int{\frac{dx}{x~(ax+b)^2}}$}
+$$\int{\frac{dx}{x~(ax+b)^2}}=
+\frac{1}{b~(ax+b)}+\frac{1}{b^2}~\ln\left(\frac{x}{ax+b}\right)$$
+<<*>>=
+)clear all
+
+--S 42
+aa:=integrate(1/(x*(a*x+b)^2),x)
+--R
+--R (- a x - b)log(a x + b) + (a x + b)log(x) + b
+--R (1) ---------------------------------------------
+--R 2 3
+--R a b x + b
+--R Type: Union(Expression Integer,...)
+--E 42
+@
+and the book says:
+<<*>>=
+--S 43
+bb:=(1/(b*(a*x+b))+(1/b^2)*log(x/(a*x+b)))
+--R
+--R x
+--R (a x + b)log(-------) + b
+--R a x + b
+--R (2) -------------------------
+--R 2 3
+--R a b x + b
+--R Type: Expression Integer
+--E 43
+@
+
+So we look at the divlog rule again:
+<<*>>=
+--S 44
+divlog:=rule(log(a/b) == log(a) - log(b))
+--R
+--R a
+--R (3) log(-) == - log(b) + log(a)
+--R b
+--R Type: RewriteRule(Integer,Integer,Expression Integer)
+--E 44
+@
+
+we apply it:
+<<*>>=
+--S 45
+cc:=divlog bb
+--R
+--R (- a x - b)log(a x + b) + (a x + b)log(x) + b
+--R (4) ---------------------------------------------
+--R 2 3
+--R a b x + b
+--R Type: Expression Integer
+--E 45
+@
+and we difference the two to find they are identical:
+<<*>>=
+--S 46
+cc-aa
+--R
+--R (5) 0
+--R Type: Expression Integer
+--E 46
+@
+
+\section{\cite{1}:14.71~~~~~$\displaystyle\int{\frac{dx}{x^2~(ax+b)^2}}$}
+$$\int{\frac{dx}{x^2~(ax+b)^2}}=
+\frac{-a}{b^2~(ax+b)}-\frac{1}{b^2~x}+
+\frac{2a}{b^3}~\ln\left(\frac{ax+b}{x}\right)$$
+<<*>>=
+)clear all
+
+--S 47
+aa:=integrate(1/(x^2*(a*x+b)^2),x)
+--R
+--R 2 2 2 2 2
+--R (2a x + 2a b x)log(a x + b) + (- 2a x - 2a b x)log(x) - 2a b x - b
+--R (1) ---------------------------------------------------------------------
+--R 3 2 4
+--R a b x + b x
+--R Type: Union(Expression Integer,...)
+--E 47
+@
+and the book says:
+<<*>>=
+--S 48
+bb:=(-a/(b^2*(a*x+b)))-(1/(b^2*x))+((2*a)/b^3)*log((a*x+b)/x)
+--R
+--R 2 2 a x + b 2
+--R (2a x + 2a b x)log(-------) - 2a b x - b
+--R x
+--R (2) ------------------------------------------
+--R 3 2 4
+--R a b x + b x
+--R Type: Expression Integer
+--E 48
+@
+which calls for our divlog rule:
+<<*>>=
+--S 49
+divlog:=rule(log(a/b) == log(a) - log(b))
+--R
+--R a
+--R (3) log(-) == - log(b) + log(a)
+--R b
+--R Type: RewriteRule(Integer,Integer,Expression Integer)
+--E 49
+@
+which we use to transform the result:
+<<*>>=
+--S 50
+cc:=divlog bb
+--R
+--R 2 2 2 2 2
+--R (2a x + 2a b x)log(a x + b) + (- 2a x - 2a b x)log(x) - 2a b x - b
+--R (4) ---------------------------------------------------------------------
+--R 3 2 4
+--R a b x + b x
+--R Type: Expression Integer
+--E 50
+@
+and we show they are identical:
+<<*>>=
+--S 51
+dd:=aa-cc
+--R
+--R (5) 0
+--R Type: Expression Integer
+--E 51
+@
+
+\section{\cite{1}:14.72~~~~~$\displaystyle\int{\frac{dx}{x^3~(ax+b)^2}}$}
+$$\int{\frac{dx}{x^3~(ax+b)^2}}=
+-\frac{(ax+b)^2}{2b^4x^2}+\frac{3a(ax+b)}{b^4x}-
+\frac{a^3x}{b^4(ax+b)}-\frac{3a^2}{b^4}~\ln\left(\frac{ax+b}{x}\right)$$
+<<*>>=
+)clear all
+
+--S 52
+aa:=integrate(1/(x^3*(a*x+b)^2),x)
+--R
+--R (1)
+--R 3 3 2 2 3 3 2 2 2 2
+--R (- 6a x - 6a b x )log(a x + b) + (6a x + 6a b x )log(x) + 6a b x
+--R +
+--R 2 3
+--R 3a b x - b
+--R /
+--R 4 3 5 2
+--R 2a b x + 2b x
+--R Type: Union(Expression Integer,...)
+--E 52
+@
+
+<<*>>=
+--S 53
+bb:=-(a*x+b)^2/(2*b^4*x^2)+(3*a*(a*x+b))/(b^4*x)-(a^3*x)/(b^4*(a*x+b))-((3*a^2)/b^4)*log((a*x+b)/x)
+--R
+--R 3 3 2 2 a x + b 3 3 2 2 2 3
+--R (- 6a x - 6a b x )log(-------) + 3a x + 9a b x + 3a b x - b
+--R x
+--R (2) ---------------------------------------------------------------
+--R 4 3 5 2
+--R 2a b x + 2b x
+--R Type: Expression Integer
+--E 53
+@
+
+<<*>>=
+--S 54
+divlog:=rule(log(a/b) == log(a) - log(b))
+--R
+--R a
+--R (3) log(-) == - log(b) + log(a)
+--R b
+--R Type: RewriteRule(Integer,Integer,Expression Integer)
+--E 54
+@
+
+<<*>>=
+--S 55
+cc:=divlog bb
+--R
+--R (4)
+--R 3 3 2 2 3 3 2 2 3 3
+--R (- 6a x - 6a b x )log(a x + b) + (6a x + 6a b x )log(x) + 3a x
+--R +
+--R 2 2 2 3
+--R 9a b x + 3a b x - b
+--R /
+--R 4 3 5 2
+--R 2a b x + 2b x
+--R Type: Expression Integer
+--E 55
+@
+
+<<*>>=
+--S 56
+cc-aa
+--R
+--R 2
+--R 3a
+--R (5) ---
+--R 4
+--R 2b
+--R Type: Expression Integer
+--E 56
+@
+
+<<*>>=
+--S 57
+dd:=D(aa,x)
+--R
+--R 1
+--R (6) ---------------------
+--R 2 5 4 2 3
+--R a x + 2a b x + b x
+--R Type: Expression Integer
+--E 57
+@
+
+<<*>>=
+--S 58
+ee:=D(bb,x)
+--R
+--R 1
+--R (7) ---------------------
+--R 2 5 4 2 3
+--R a x + 2a b x + b x
+--R Type: Expression Integer
+--E 58
+@
+
+<<*>>=
+--S 59
+dd-ee
+--R
+--R (8) 0
+--R Type: Expression Integer
+--E 59
+@
+
+\section{\cite{1}:14.73~~~~~$\displaystyle\int{\frac{dx}{(ax+b)^3}}$}
+$$\int{\frac{dx}{(ax+b)^3}}=\frac{-1}{2a(ax+b)^2}$$
+<<*>>=
+)clear all
+
+--S 60
+aa:=integrate(1/(a*x+b)^3,x)
+--R
+--R 1
+--R (1) - ----------------------
+--R 3 2 2 2
+--R 2a x + 4a b x + 2a b
+--R Type: Union(Expression Integer,...)
+--E 60
+@
+
+{\bf NOTE: }There is a missing factor of $1/a$ in the published book.
+This factor has been inserted here.
+<<*>>=
+--S 61
+bb:=-1/(2*a*(a*x+b)^2)
+--R
+--R 1
+--R (2) - ----------------------
+--R 3 2 2 2
+--R 2a x + 4a b x + 2a b
+--R Type: Fraction Polynomial Integer
+--E 61
+@
+
+<<*>>=
+--S 62
+aa-bb
+--R
+--R (3) 0
+--R Type: Expression Integer
+--E 62
+@
+
+\section{\cite{1}:14.74~~~~~$\displaystyle\int{\frac{x~dx}{(ax+b)^3}}$}
+$$\int{\frac{x~dx}{(ax+b)^3}}=
+\frac{-1}{a^2(ax+b)}+\frac{b}{2a^2(ax+b)^2}$$
+<<*>>=
+)clear all
+
+--S 63
+aa:=integrate(x/(a*x+b)^3,x)
+--R
+--R - 2a x - b
+--R (1) ----------------------
+--R 4 2 3 2 2
+--R 2a x + 4a b x + 2a b
+--R Type: Union(Expression Integer,...)
+--E 63
+@
+
+<<*>>=
+--S 64
+bb:=-1/(a^2*(a*x+b))+b/(2*a^2*(a*x+b)^2)
+--R
+--R - 2a x - b
+--R (2) ----------------------
+--R 4 2 3 2 2
+--R 2a x + 4a b x + 2a b
+--R Type: Fraction Polynomial Integer
+--E 64
+@
+
+<<*>>=
+--S 65
+aa-bb
+--R
+--R (3) 0
+--R Type: Expression Integer
+--E 65
+@
+
+\section{\cite{1}:14.75~~~~~$\displaystyle\int{\frac{x^2~dx}{(ax+b)^3}}$}
+$$\int{\frac{x^2~dx}{(ax+b)^3}}=
+\frac{2b}{a^3(ax+b)}-\frac{b^2}{2a^3(ax+b)^2}+
+\frac{1}{a^3}~\ln(ax+b)$$
+<<*>>=
+)clear all
+
+--S 66
+aa:=integrate(x^2/(a*x+b)^3,x)
+--R
+--R 2 2 2 2
+--R (2a x + 4a b x + 2b )log(a x + b) + 4a b x + 3b
+--R (1) -------------------------------------------------
+--R 5 2 4 3 2
+--R 2a x + 4a b x + 2a b
+--R Type: Union(Expression Integer,...)
+--E 66
+@
+
+<<*>>=
+--S 67
+bb:=(2*b)/(a^3*(a*x+b))-(b^2)/(2*a^3*(a*x+b)^2)+1/a^3*log(a*x+b)
+--R
+--R 2 2 2 2
+--R (2a x + 4a b x + 2b )log(a x + b) + 4a b x + 3b
+--R (2) -------------------------------------------------
+--R 5 2 4 3 2
+--R 2a x + 4a b x + 2a b
+--R Type: Expression Integer
+--E 67
+@
+
+<<*>>=
+--S 68
+aa-bb
+--R
+--R (3) 0
+--R Type: Expression Integer
+--E 68
+@
+
+\section{\cite{1}:14.76~~~~~$\displaystyle\int{\frac{x^3~dx}{(ax+b)^3}}$}
+$$\int{\frac{x^3~dx}{(ax+b)^3}}=
+\frac{x}{a^3}-\frac{3b^2}{a^4(ax+b)}+\frac{b^3}{2a^4(ax+b)^2}-
+\frac{3b}{a^4}~\ln(ax+b)$$
+<<*>>=
+)clear all
+--S 69
+aa:=integrate(x^3/(a*x+b)^3,x)
+--R
+--R (1)
+--R 2 2 2 3 3 3 2 2 2 3
+--R (- 6a b x - 12a b x - 6b )log(a x + b) + 2a x + 4a b x - 4a b x - 5b
+--R ------------------------------------------------------------------------
+--R 6 2 5 4 2
+--R 2a x + 4a b x + 2a b
+--R Type: Union(Expression Integer,...)
+--E 69
+@
+
+<<*>>=
+--S 70
+bb:=(x/a^3)-(3*b^2)/(a^4*(a*x+b))+b^3/(2*a^4*(a*x+b)^2)-(3*b)/a^4*log(a*x+b)
+--R
+--R (2)
+--R 2 2 2 3 3 3 2 2 2 3
+--R (- 6a b x - 12a b x - 6b )log(a x + b) + 2a x + 4a b x - 4a b x - 5b
+--R ------------------------------------------------------------------------
+--R 6 2 5 4 2
+--R 2a x + 4a b x + 2a b
+--R Type: Expression Integer
+--E 70
+@
+
+<<*>>=
+--S 71
+aa-bb
+--R
+--R (3) 0
+--R Type: Expression Integer
+--E 71
+@
+
+\section{\cite{1}:14.77~~~~~$\displaystyle\int{\frac{dx}{x(ax+b)^3}}$}
+$$\int{\frac{dx}{x(ax+b)^3}}=
+\frac{3}{2b(ax+b)^2}+\frac{2ax}{2b^2(ax+b)^2}-
+\frac{1}{b^3}*\ln\left(\frac{ax+b}{x}\right)$$
+
+{\bf NOTE: }The equation given in the book is wrong. This is correct.
+
+<<*>>=
+)clear all
+
+--S 72
+aa:=integrate(1/(x*(a*x+b)^3),x)
+--R
+--R (1)
+--R 2 2 2 2 2 2
+--R (- 2a x - 4a b x - 2b )log(a x + b) + (2a x + 4a b x + 2b )log(x)
+--R +
+--R 2
+--R 2a b x + 3b
+--R /
+--R 2 3 2 4 5
+--R 2a b x + 4a b x + 2b
+--R Type: Union(Expression Integer,...)
+--E 72
+@
+
+<<*>>=
+--S 73
+bb:=3/(2*b*(a*x+b)^2)+(2*a*x)/(2*b^2*(a*x+b)^2)-1/b^3*log((a*x+b)/x)
+--R
+--R 2 2 2 a x + b 2
+--R (- 2a x - 4a b x - 2b )log(-------) + 2a b x + 3b
+--R x
+--R (2) ---------------------------------------------------
+--R 2 3 2 4 5
+--R 2a b x + 4a b x + 2b
+--R Type: Expression Integer
+--E 73
+@
+
+<<*>>=
+--S 74
+divlog:=rule(log(a/b) == log(a) - log(b))
+--R
+--R a
+--R (3) log(-) == - log(b) + log(a)
+--R b
+--R Type: RewriteRule(Integer,Integer,Expression Integer)
+--E 74
+@
+
+<<*>>=
+--S 75
+cc:=divlog bb
+--R
+--R (4)
+--R 2 2 2 2 2 2
+--R (- 2a x - 4a b x - 2b )log(a x + b) + (2a x + 4a b x + 2b )log(x)
+--R +
+--R 2
+--R 2a b x + 3b
+--R /
+--R 2 3 2 4 5
+--R 2a b x + 4a b x + 2b
+--R Type: Expression Integer
+--E 75
+@
+
+<<*>>=
+--S 76
+aa-cc
+--R
+--R (5) 0
+--R Type: Expression Integer
+--E 76
+@
+
+\section{\cite{1}:14.78~~~~~$\displaystyle\int{\frac{dx}{x^2(ax+b)^3}}$}
+$$\int{\frac{dx}{x^2(ax+b)^3}}=
+\frac{-a}{2b^2(ax+b)^2}-\frac{2a}{b^3(ax+b)}-
+\frac{1}{b^3x}+\frac{3a}{b^4}~\ln\left(\frac{ax+b}{x}\right)$$
+<<*>>=
+)clear all
+
+--S 77
+aa:=integrate(1/(x^2*(a*x+b)^3),x)
+--R
+--R (1)
+--R 3 3 2 2 2
+--R (6a x + 12a b x + 6a b x)log(a x + b)
+--R +
+--R 3 3 2 2 2 2 2 2 3
+--R (- 6a x - 12a b x - 6a b x)log(x) - 6a b x - 9a b x - 2b
+--R /
+--R 2 4 3 5 2 6
+--R 2a b x + 4a b x + 2b x
+--R Type: Union(Expression Integer,...)
+--E 77
+@
+
+<<*>>=
+--S 78
+bb:=-a/(2*b^2*(a*x+b)^2)-(2*a)/(b^3*(a*x+b))-1/(b^3*x)+((3*a)/b^4)*log((a*x+b)/x)
+--R
+--R 3 3 2 2 2 a x + b 2 2 2 3
+--R (6a x + 12a b x + 6a b x)log(-------) - 6a b x - 9a b x - 2b
+--R x
+--R (2) ----------------------------------------------------------------
+--R 2 4 3 5 2 6
+--R 2a b x + 4a b x + 2b x
+--R Type: Expression Integer
+--E 78
+@
+
+<<*>>=
+--S 79
+divlog:=rule(log(a/b) == log(a) - log(b))
+--R
+--R a
+--R (3) log(-) == - log(b) + log(a)
+--R b
+--R Type: RewriteRule(Integer,Integer,Expression Integer)
+--E 79
+@
+
+<<*>>=
+--S 80
+cc:=divlog bb
+--R
+--R (4)
+--R 3 3 2 2 2
+--R (6a x + 12a b x + 6a b x)log(a x + b)
+--R +
+--R 3 3 2 2 2 2 2 2 3
+--R (- 6a x - 12a b x - 6a b x)log(x) - 6a b x - 9a b x - 2b
+--R /
+--R 2 4 3 5 2 6
+--R 2a b x + 4a b x + 2b x
+--R Type: Expression Integer
+--E 80
+@
+
+<<*>>=
+--S 81
+cc-aa
+--R
+--R (5) 0
+--R Type: Expression Integer
+--E 81
+@
+
+\section{\cite{1}:14.79~~~~~$\displaystyle\int{\frac{dx}{x^3(ax+b)^3}}$}
+$$\int{\frac{dx}{x^3(ax+b)^3}}=$$
+$$-\frac{1}{2bx^2(ax+b)^2}+
+\frac{2a}{b^2x(ax+b)^2}+
+\frac{9a^2}{b^3(ax+b)^2}+
+\frac{6a^3x}{b^4(ax+b)^2}-
+\frac{6a^2}{b^5}~\ln\left(\frac{ax+b}{x}\right)$$
+
+{\bf NOTE: }The equation given in the book is wrong. This is correct.
+
+<<*>>=
+)clear all
+
+--S 82
+aa:=integrate(1/(x^3*(a*x+b)^3),x)
+--R
+--R (1)
+--R 4 4 3 3 2 2 2
+--R (- 12a x - 24a b x - 12a b x )log(a x + b)
+--R +
+--R 4 4 3 3 2 2 2 3 3 2 2 2 3 4
+--R (12a x + 24a b x + 12a b x )log(x) + 12a b x + 18a b x + 4a b x - b
+--R /
+--R 2 5 4 6 3 7 2
+--R 2a b x + 4a b x + 2b x
+--R Type: Union(Expression Integer,...)
+--E 82
+@
+
+<<*>>=
+--S 83
+bb:=-1/(2*b*x^2*(a*x+b)^2)_
+ +(2*a)/(b^2*x*(a*x+b)^2)_
+ +(9*a^2)/(b^3*(a*x+b)^2)_
+ +(6*a^3*x)/(b^4*(a*x+b)^2)_
+ +(-6*a^2)/b^5*log((a*x+b)/x)
+--R
+--R (2)
+--R 4 4 3 3 2 2 2 a x + b 3 3 2 2 2
+--R (- 12a x - 24a b x - 12a b x )log(-------) + 12a b x + 18a b x
+--R x
+--R +
+--R 3 4
+--R 4a b x - b
+--R /
+--R 2 5 4 6 3 7 2
+--R 2a b x + 4a b x + 2b x
+--R Type: Expression Integer
+--E 83
+@
+<<*>>=
+--S 84
+cc:=aa-bb
+--R
+--R 2 2 2 a x + b
+--R - 6a log(a x + b) + 6a log(x) + 6a log(-------)
+--R x
+--R (3) -----------------------------------------------
+--R 5
+--R b
+--R Type: Expression Integer
+--E 84
+@
+
+<<*>>=
+--S 85
+divlog:=rule(log(a/b) == log(a) - log(b))
+--R
+--R a
+--R (4) log(-) == - log(b) + log(a)
+--R b
+--R Type: RewriteRule(Integer,Integer,Expression Integer)
+--E 85
+@
+
+<<*>>=
+--S 86
+divlog cc
+--R
+--R (5) 0
+--R Type: Expression Integer
+--E 86
+@
+
+\section{\cite{1}:14.80~~~~~$\displaystyle\int{(ax+b)^n~dx}$}
+$$\int{(ax+b)^n~dx}=
+\frac{(ax+b)^{n+1}}{(n+1)a}{\rm\ provided\ }n \ne -1$$
+<<*>>=
+)clear all
+--S 87
+aa:=integrate((a*x+b)^n,x)
+--R
+--R n log(a x + b)
+--R (a x + b)%e
+--R (1) -------------------------
+--R a n + a
+--R Type: Union(Expression Integer,...)
+--E 87
+@
+
+<<*>>=
+--S 88
+explog:=rule(%e^(n*log(x)) == x^n)
+--R
+--R n log(x) n
+--R (2) %e == x
+--R Type: RewriteRule(Integer,Integer,Expression Integer)
+--E 88
+@
+
+<<*>>=
+--S 89
+explog aa
+--R
+--R n
+--R (a x + b)(a x + b)
+--R (3) -------------------
+--R a n + a
+--R Type: Expression Integer
+--E 89
+@
+
+\section{\cite{1}:14.81~~~~~$\displaystyle\int{x(ax+b)^n~dx}$}
+$$\int{x(ax+b)^n~dx}=
+\frac{(ax+b)^{n+2}}{(n+2)a^2}-\frac{b(ax+b)^{n+1}}{(n+1)a^2}
+{\rm\ provided\ }n \ne -1,-2$$
+
+\section{\cite{1}:14.82~~~~~$\displaystyle\int{x^2(ax+b)^n~dx}$}
+$$\int{x^2(ax+b)^n~dx}=
+\frac{(ax+b)^{n+2}}{(n+3)a^3}-
+\frac{2b(ax+b)^{n+2}}{(n+2)a^3}+
+\frac{b^2(ax+b)^{n+1}}{(n+1)a^3}
+{\rm\ provided\ }n \ne -1,-2,-3$$
+
+<<*>>=
+)spool
+)lisp (bye)
+@
+
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} Spiegel, Murray R.
+{\sl Mathematical Handbook of Formulas and Tables}\\
+Schaum's Outline Series McGraw-Hill 1968 pp60-61
+\end{thebibliography}
+\end{document}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <gi...@ax...> - 2007-09-14 03:24:26
|
changelog | 2 +
src/input/Makefile.pamphlet | 5 +-
src/input/schaum1.input.pamphlet | 1265 ++++++++++++++++++++++++++++++++++++++
3 files changed, 1271 insertions(+), 1 deletions(-)
New commits:
commit ae3f97ff1fdff1a9e61b45c68cfb0121b67a0363
Author: Tim Daly <root@localhost.localdomain>
Date: Wed Sep 12 05:17:12 2007 -0400
20070913 tpd src/input/Makefile schaum1.input added
20070913 tpd src/input/schaum1.input added
|
|
From: <da...@us...> - 2007-09-11 05:08:05
|
Revision: 742
http://axiom.svn.sourceforge.net/axiom/?rev=742&view=rev
Author: daly
Date: 2007-09-10 22:08:07 -0700 (Mon, 10 Sep 2007)
Log Message:
-----------
add missing algebra files
Added Paths:
-----------
trunk/axiom/src/algebra/fffg.spad.pamphlet
trunk/axiom/src/algebra/mantepse.spad.pamphlet
trunk/axiom/src/algebra/rec.spad.pamphlet
trunk/axiom/src/algebra/ssolve.spad.pamphlet
Added: trunk/axiom/src/algebra/fffg.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/fffg.spad.pamphlet (rev 0)
+++ trunk/axiom/src/algebra/fffg.spad.pamphlet 2007-09-11 05:08:07 UTC (rev 742)
@@ -0,0 +1,666 @@
+\documentclass{article}
+\usepackage{axiom,amsthm}
+\newtheorem{ToDo}{ToDo}[section]
+\begin{document}
+\title{fffg.spad}
+\author{Martin Rubey}
+\maketitle
+\begin{abstract}
+ The packages defined in this file provide fast fraction free rational
+ interpolation algorithms.
+\end{abstract}
+\tableofcontents
+\section{package FAMR2 FiniteAbelianMonoidRingFunctions2}
+<<package FAMR2 FiniteAbelianMonoidRingFunctions2>>=
+)abbrev package FAMR2 FiniteAbelianMonoidRingFunctions2
+++ Author: Martin Rubey
+++ Description:
+++ This package provides a mapping function for \spadtype{FiniteAbelianMonoidRing}
+FiniteAbelianMonoidRingFunctions2(E: OrderedAbelianMonoid,
+ R1: Ring,
+ A1: FiniteAbelianMonoidRing(R1, E),
+ R2: Ring,
+ A2: FiniteAbelianMonoidRing(R2, E)) _
+ : Exports == Implementation where
+ Exports == with
+
+ map: (R1 -> R2, A1) -> A2
+ ++ \spad{map}(f, a) applies the map f to each coefficient in a. It is
+ ++ assumed that f maps 0 to 0
+
+ Implementation == add
+
+ map(f: R1 -> R2, a: A1): A2 ==
+ if zero? a then 0$A2
+ else
+ monomial(f leadingCoefficient a, degree a)$A2 + map(f, reductum a)
+
+@
+
+\section{package FFFG FractionFreeFastGaussian}
+<<package FFFG FractionFreeFastGaussian>>=
+)abbrev package FFFG FractionFreeFastGaussian
+++ Author: Martin Rubey
+++ Description:
+++ This package implements the interpolation algorithm proposed in Beckermann,
+++ Bernhard and Labahn, George, Fraction-free computation of matrix rational
+++ interpolants and matrix GCDs, SIAM Journal on Matrix Analysis and
+++ Applications 22.
+FractionFreeFastGaussian(D, V): Exports == Implementation where
+ D: Join(IntegralDomain, GcdDomain)
+ V: AbelianMonoidRing(D, NonNegativeInteger)
+
+ SUP ==> SparseUnivariatePolynomial
+
+ cFunction ==> (NonNegativeInteger, Vector SUP D) -> D
+
+ CoeffAction ==> (NonNegativeInteger, NonNegativeInteger, V) -> D
+
+ Exports == with
+
+ fffg: (List D, cFunction, List NonNegativeInteger) -> Matrix SUP D
+ ++ \spad{fffg} is the general algorithm as proposed by Beckermann and
+ ++ Labahn.
+ ++
+ ++ The second argument, c(k, M) computes c_k(M) as in Equation (2). Note
+ ++ that the information about f is therefore encoded in c.
+
+ interpolate: (List D, List D, NonNegativeInteger) -> Fraction SUP D
+ ++ \spad{interpolate(xlist, ylist, deg} returns the rational function with
+ ++ numerator degree at most \spad{deg} and denominator degree at most
+ ++ \spad{#xlist-deg-1} that interpolates the given points using
+ ++ fraction free arithmetic. Note that rational interpolation does not
+ ++ guarantee that all given points are interpolated correctly:
+ ++ unattainable points may make this impossible.
+
+@
+
+\begin{ToDo}
+ The following function could be moved to [[FFFGF]], parallel to
+ [[generalInterpolation]]. However, the reason for moving
+ [[generalInterpolation]] for fractions to a separate package was the need of
+ a generic signature, hence the extra argument [[VF]] to [[FFFGF]]. In the
+ special case of rational interpolation, this extra argument is not necessary,
+ since we are always returning a fraction of [[SUP]]s, and ignore [[V]]. In
+ fact, [[V]] is not needed for [[fffg]] itself, only if we want to specify a
+ [[CoeffAction]].
+
+ Thus, maybe it would be better to move [[fffg]] to a separate package?
+\end{ToDo}
+
+<<package FFFG FractionFreeFastGaussian>>=
+ interpolate: (List Fraction D, List Fraction D, NonNegativeInteger)
+ -> Fraction SUP D
+ ++ \spad{interpolate(xlist, ylist, deg} returns the rational function with
+ ++ numerator degree \spad{deg} that interpolates the given points using
+ ++ fraction free arithmetic.
+
+ generalInterpolation: (List D, CoeffAction,
+ Vector V, List NonNegativeInteger) -> Matrix SUP D
+ ++ \spad{generalInterpolation(C, CA, f, eta)} performs Hermite-Pade
+ ++ approximation using the given action CA of polynomials on the elements
+ ++ of f. The result is guaranteed to be correct up to order
+ ++ |eta|-1. Given that eta is a "normal" point, the degrees on the
+ ++ diagonal are given by eta. The degrees of column i are in this case
+ ++ eta + e.i - [1,1,...,1], where the degree of zero is -1.
+ ++
+ ++ The first argument C is the list of coefficients c_{k,k} in the
+ ++ expansion <x^k> z g(x) = sum_{i=0}^k c_{k,i} <x^i> g(x).
+ ++
+ ++ The second argument, CA(k, l, f), should return the coefficient of x^k
+ ++ in z^l f(x).
+
+ generalInterpolation: (List D, CoeffAction,
+ Vector V, NonNegativeInteger, NonNegativeInteger)
+ -> Stream Matrix SUP D
+ ++ \spad{generalInterpolation(C, CA, f, sumEta, maxEta)} applies
+ ++ \spad{generalInterpolation(C, CA, f, eta)} for all possible \spad{eta}
+ ++ with maximal entry \spad{maxEta} and sum of entries at most
+ ++ \spad{sumEta}.
+ ++
+ ++ The first argument C is the list of coefficients c_{k,k} in the
+ ++ expansion <x^k> z g(x) = sum_{i=0}^k c_{k,i} <x^i> g(x).
+ ++
+ ++ The second argument, CA(k, l, f), should return the coefficient of x^k
+ ++ in z^l f(x).
+
+ generalCoefficient: (CoeffAction, Vector V,
+ NonNegativeInteger, Vector SUP D) -> D
+ ++ \spad{generalCoefficient(action, f, k, p)} gives the coefficient of
+ ++ x^k in p(z)\dot f(x), where the action of z^l on a polynomial in x is
+ ++ given by action, i.e., action(k, l, f) should return the coefficient
+ ++ of x^k in z^l f(x).
+
+ ShiftAction: (NonNegativeInteger, NonNegativeInteger, V) -> D
+ ++ \spad{ShiftAction(k, l, g)} gives the coefficient of x^k in z^l g(x),
+ ++ where \spad{z*(a+b*x+c*x^2+d*x^3+...) = (b*x+2*c*x^2+3*d*x^3+...)}. In
+ ++ terms of sequences, z*u(n)=n*u(n).
+
+ ShiftC: NonNegativeInteger -> List D
+ ++ \spad{ShiftC} gives the coefficients c_{k,k} in the expansion <x^k> z
+ ++ g(x) = sum_{i=0}^k c_{k,i} <x^i> g(x), where z acts on g(x) by
+ ++ shifting. In fact, the result is [0,1,2,...]
+
+ DiffAction: (NonNegativeInteger, NonNegativeInteger, V) -> D
+ ++ \spad{DiffAction(k, l, g)} gives the coefficient of x^k in z^l g(x),
+ ++ where z*(a+b*x+c*x^2+d*x^3+...) = (a*x+b*x^2+c*x^3+...), i.e.,
+ ++ multiplication with x.
+
+ DiffC: NonNegativeInteger -> List D
+ ++ \spad{DiffC} gives the coefficients c_{k,k} in the expansion <x^k> z
+ ++ g(x) = sum_{i=0}^k c_{k,i} <x^i> g(x), where z acts on g(x) by
+ ++ shifting. In fact, the result is [0,0,0,...]
+
+ qShiftAction: (D, NonNegativeInteger, NonNegativeInteger, V) -> D
+ ++ \spad{qShiftAction(q, k, l, g)} gives the coefficient of x^k in z^l
+ ++ g(x), where z*(a+b*x+c*x^2+d*x^3+...) =
+ ++ (a+q*b*x+q^2*c*x^2+q^3*d*x^3+...). In terms of sequences,
+ ++ z*u(n)=q^n*u(n).
+
+ qShiftC: (D, NonNegativeInteger) -> List D
+ ++ \spad{qShiftC} gives the coefficients c_{k,k} in the expansion <x^k> z
+ ++ g(x) = sum_{i=0}^k c_{k,i} <x^i> g(x), where z acts on g(x) by
+ ++ shifting. In fact, the result is [1,q,q^2,...]
+
+ Implementation ==> add
+
+-------------------------------------------------------------------------------
+-- Shift Operator
+-------------------------------------------------------------------------------
+
+-- ShiftAction(k, l, f) is the CoeffAction appropriate for the shift operator.
+
+ ShiftAction(k: NonNegativeInteger, l: NonNegativeInteger, f: V): D ==
+ k**l*coefficient(f, k)
+
+
+ ShiftC(total: NonNegativeInteger): List D ==
+ [i::D for i in 0..total-1]
+
+-------------------------------------------------------------------------------
+-- q-Shift Operator
+-------------------------------------------------------------------------------
+
+-- q-ShiftAction(k, l, f) is the CoeffAction appropriate for the q-shift operator.
+
+ qShiftAction(q: D, k: NonNegativeInteger, l: NonNegativeInteger, f: V): D ==
+ q**(k*l)*coefficient(f, k)
+
+
+ qShiftC(q: D, total: NonNegativeInteger): List D ==
+ [q**i for i in 0..total-1]
+
+-------------------------------------------------------------------------------
+-- Differentiation Operator
+-------------------------------------------------------------------------------
+
+-- DiffAction(k, l, f) is the CoeffAction appropriate for the differentiation
+-- operator.
+
+ DiffAction(k: NonNegativeInteger, l: NonNegativeInteger, f: V): D ==
+ coefficient(f, (k-l)::NonNegativeInteger)
+
+
+ DiffC(total: NonNegativeInteger): List D ==
+ [0 for i in 1..total]
+
+-------------------------------------------------------------------------------
+-- general - suitable for functions f
+-------------------------------------------------------------------------------
+
+-- get the coefficient of z^k in the scalar product of p and f, the action
+-- being defined by coeffAction
+
+ generalCoefficient(coeffAction: CoeffAction, f: Vector V,
+ k: NonNegativeInteger, p: Vector SUP D): D ==
+ res: D := 0
+ for i in 1..#f repeat
+
+-- Defining a and b and summing only over those coefficients that might be
+-- nonzero makes a huge difference in speed
+ a := f.i
+ b := p.i
+ for l in minimumDegree b..degree b repeat
+ if not zero? coefficient(b, l)
+ then res := res + coefficient(b, l) * coeffAction(k, l, a)
+ res
+
+
+ generalInterpolation(C: List D, coeffAction: CoeffAction,
+ f: Vector V,
+ eta: List NonNegativeInteger): Matrix SUP D ==
+
+ c: cFunction := generalCoefficient(coeffAction, f,
+ (#1-1)::NonNegativeInteger, #2)
+ fffg(C, c, eta)
+
+
+
+-------------------------------------------------------------------------------
+-- general - suitable for functions f - trying all possible degree combinations
+-------------------------------------------------------------------------------
+
+@
+
+The following function returns the lexicographically next vector with
+non-negative components smaller than [[p]] with the same sum as [[v]].
+
+<<package FFFG FractionFreeFastGaussian>>=
+ nextVector!(p: NonNegativeInteger, v: List NonNegativeInteger)
+ : Union("failed", List NonNegativeInteger) ==
+ n := #v
+ pos := position(#1 < p, v)
+ zero? pos => return "failed"
+ if pos = 1 then
+ sum: Integer := v.1
+ for i in 2..n repeat
+ if v.i < p and sum > 0 then
+ v.i := v.i + 1
+ sum := sum - 1
+ for j in 1..i-1 repeat
+ if sum > p then
+ v.j := p
+ sum := sum - p
+ else
+ v.j := sum::NonNegativeInteger
+ sum := 0
+ return v
+ else sum := sum + v.i
+ return "failed"
+ else
+ v.pos := v.pos + 1
+ v.(pos-1) := (v.(pos-1) - 1)::NonNegativeInteger
+
+ v
+@
+
+The following function returns the stream of all possible degree vectors,
+beginning with [[v]], where the degree vectors are sorted in reverse
+lexicographic order. Furthermore, the entries are all less or equal to [[p]]
+and their sum equals the sum of the entries of [[v]]. We assume that the
+entries of [[v]] are also all less or equal to [[p]].
+
+<<package FFFG FractionFreeFastGaussian>>=
+ vectorStream(p: NonNegativeInteger, v: List NonNegativeInteger)
+ : Stream List NonNegativeInteger == delay
+ next := nextVector!(p, copy v)
+ (next case "failed") => empty()$Stream(List NonNegativeInteger)
+ cons(next, vectorStream(p, next))
+@
+
+[[vectorStream2]] skips every second entry of [[vectorStream]].
+
+<<package FFFG FractionFreeFastGaussian>>=
+ vectorStream2(p: NonNegativeInteger, v: List NonNegativeInteger)
+ : Stream List NonNegativeInteger == delay
+ next := nextVector!(p, copy v)
+ (next case "failed") => empty()$Stream(List NonNegativeInteger)
+ next2 := nextVector!(p, copy next)
+ (next2 case "failed") => cons(next, empty())
+ cons(next2, vectorStream2(p, next2))
+@
+
+This version of [[generalInterpolation]] returns a stream of solutions, one for
+each possible degree vector. Thus, it only needs to apply the previously
+defined [[generalInterpolation]] to each degree vector. These are generated by
+[[vectorStream]] and [[vectorStream2]] respectively.
+
+If [[f]] consists of two elements only, we can skip every second degree vector:
+note that [[fffg]], and thus also [[generalInterpolation]], returns a matrix
+with [[#f]] columns, each corresponding to a solution of the interpolation
+problem. More precisely, the $i$\textsuperscript{th} column is a solution with
+degrees [[eta]]$-(1,1,\dots,1)+e_i$. Thus, in the case of $2\times 2$ matrices,
+[[vectorStream]] would produce solutions corresponding to $(d,0), (d-1,1);
+(d-1,1), (d-2, 2); (d-2, 2), (d-3,3)\dots$, i.e., every second matrix is
+redundant.
+
+Although some redundancy exists also for higher dimensional [[f]], the scheme
+becomes much more complicated, thus we did not implement it.
+
+<<package FFFG FractionFreeFastGaussian>>=
+ generalInterpolation(C: List D, coeffAction: CoeffAction,
+ f: Vector V,
+ sumEta: NonNegativeInteger,
+ maxEta: NonNegativeInteger)
+ : Stream Matrix SUP D ==
+
+ <<generate an initial degree vector>>
+
+ if #f = 2 then
+ map(generalInterpolation(C, coeffAction, f, #1),
+ cons(eta, vectorStream2(maxEta, eta)))
+ $StreamFunctions2(List NonNegativeInteger,
+ Matrix SUP D)
+ else
+ map(generalInterpolation(C, coeffAction, f, #1),
+ cons(eta, vectorStream(maxEta, eta)))
+ $StreamFunctions2(List NonNegativeInteger,
+ Matrix SUP D)
+@
+
+We need to generate an initial degree vector, being the minimal element in
+reverse lexicographic order, i.e., $m, m, \dots, m, k, 0, 0, \dots$, where $m$
+is [[maxEta]] and $k$ is the remainder of [[sumEta]] divided by
+[[maxEta]]. This is done by the following code:
+
+<<generate an initial degree vector>>=
+sum: Integer := sumEta
+entry: Integer
+eta: List NonNegativeInteger
+ := [(if sum < maxEta _
+ then (entry := sum; sum := 0) _
+ else (entry := maxEta; sum := sum - maxEta); _
+ entry::NonNegativeInteger) for i in 1..#f]
+@
+
+We want to generate all vectors with sum of entries being at most
+[[sumEta]]. Therefore the following is incorrect.
+<<BUG generate an initial degree vector>>=
+-- (sum > 0) => empty()$Stream(Matrix SUP D)
+@
+
+<<package FFFG FractionFreeFastGaussian>>=
+-------------------------------------------------------------------------------
+-- rational interpolation
+-------------------------------------------------------------------------------
+
+ interpolate(x: List Fraction D, y: List Fraction D, d: NonNegativeInteger)
+ : Fraction SUP D ==
+ gx := splitDenominator(x)$InnerCommonDenominator(D, Fraction D, _
+ List D, _
+ List Fraction D)
+ gy := splitDenominator(y)$InnerCommonDenominator(D, Fraction D, _
+ List D, _
+ List Fraction D)
+ r := interpolate(gx.num, gy.num, d)
+ elt(numer r, monomial(gx.den,1))/(gy.den*elt(denom r, monomial(gx.den,1)))
+
+
+ interpolate(x: List D, y: List D, d: NonNegativeInteger): Fraction SUP D ==
+-- berechne Interpolante mit Graden d und N-d-1
+ if (N := #x) ~= #y then
+ error "interpolate: number of points and values must match"
+ if N <= d then
+ error "interpolate: numerator degree must be smaller than number of data points"
+ c: cFunction := y.#1 * elt(#2.2, x.#1) - elt(#2.1, x.#1)
+ eta: List NonNegativeInteger := [d, (N-d)::NonNegativeInteger]
+ M := fffg(x, c, eta)
+
+ if zero?(M.(2,1)) then M.(1,2)/M.(2,2)
+ else M.(1,1)/M.(2,1)
+
+@
+Because of Lemma~5.3, [[M.1.(2,1)]] and [[M.1.(2,2)]] cannot both vanish,
+since [[eta_sigma]] is always $\sigma$-normal by Theorem~7.2 and therefore also
+para-normal, see Definition~4.2.
+
+Because of Lemma~5.1 we have that [[M.1.(*,2)]] is a solution of the
+interpolation problem, if [[M.1.(2,1)]] vanishes.
+
+<<package FFFG FractionFreeFastGaussian>>=
+-------------------------------------------------------------------------------
+-- fffg
+-------------------------------------------------------------------------------
+
+-- a major part of the time is spent here
+ recurrence(M: Matrix SUP D, lambda: NonNegativeInteger, m: NonNegativeInteger,
+ r: Vector D, d: D, z: SUP D, Ck: D, p: Vector D): Matrix SUP D ==
+
+ rLambda := qelt(r, lambda)
+ polyf := rLambda * (z - Ck::SUP D)
+
+ for i in 1..m repeat
+ Milambda := qelt(M, i, lambda)
+ newMilambda := polyf * Milambda
+
+
+-- update columns ~= lambda and calculate their sum
+ for l in 1..m | l ~= lambda repeat
+ rl := qelt(r, l)
+ Mil := ((qelt(M, i, l) * rLambda - Milambda * rl) exquo d)::SUP D
+ qsetelt!(M, i, l, Mil)
+
+ pl := qelt(p, l)
+ newMilambda := newMilambda - pl * Mil
+
+-- update column lambda
+
+ qsetelt!(M, i, lambda, (newMilambda exquo d)::SUP D)
+
+ M
+
+
+ fffg(C: List D, c: cFunction, eta: List NonNegativeInteger): Matrix SUP D ==
+
+-- eta is the vector of degrees. We compute M with degrees eta+e_i-1, i=1..m
+ z: SUP D := monomial(1, 1)
+ m: NonNegativeInteger := #eta
+ M: Matrix SUP D := scalarMatrix(m, 1)
+ d: D := 1
+ K: NonNegativeInteger := reduce(_+, eta)
+ etak: Vector NonNegativeInteger := zero(m)
+ r: Vector D := zero(m)
+ p: Vector D := zero(m)
+ Lambda: List Integer
+ lambdaMax: Integer
+ lambda: NonNegativeInteger
+
+ for k in 1..K repeat
+-- k = sigma+1
+
+ for l in 1..m repeat r.l := c(k, column(M, l))
+
+ Lambda := [eta.l-etak.l for l in 1..m | r.l ~= 0]
+
+-- if Lambda is empty, then M, d and etak remain unchanged. Otherwise, we look
+-- for the next closest para-normal point.
+
+ (empty? Lambda) => "iterate"
+
+ lambdaMax := reduce(max, Lambda)
+ lambda := 1
+ while eta.lambda-etak.lambda < lambdaMax or r.lambda = 0 repeat
+ lambda := lambda + 1
+
+-- Calculate leading coefficients
+
+ for l in 1..m | l ~= lambda repeat
+ if etak.l > 0 then
+ p.l := coefficient(M.(l, lambda), (etak.l-1)::NonNegativeInteger)
+ else
+ p.l := 0
+
+-- increase order and adjust degree constraints
+
+ M := recurrence(M, lambda, m, r, d, z, C.k, p)
+
+ d := r.lambda
+ etak.lambda := etak.lambda + 1
+
+ M
+
+@
+%$
+
+\section{package FFFGF FractionFreeFastGaussianFractions}
+<<package FFFGF FractionFreeFastGaussianFractions>>=
+)abbrev package FFFGF FractionFreeFastGaussianFractions
+++ Author: Martin Rubey
+++ Description:
+++ This package lifts the interpolation functions from
+++ \spadtype{FractionFreeFastGaussian} to fractions.
+FractionFreeFastGaussianFractions(D, V, VF): Exports == Implementation where
+ D: Join(IntegralDomain, GcdDomain)
+ V: FiniteAbelianMonoidRing(D, NonNegativeInteger)
+ VF: FiniteAbelianMonoidRing(Fraction D, NonNegativeInteger)
+
+ F ==> Fraction D
+
+ SUP ==> SparseUnivariatePolynomial
+
+ FFFG ==> FractionFreeFastGaussian
+
+ FAMR2 ==> FiniteAbelianMonoidRingFunctions2
+
+ cFunction ==> (NonNegativeInteger, Vector SUP D) -> D
+
+ CoeffAction ==> (NonNegativeInteger, NonNegativeInteger, V) -> D
+-- coeffAction(k, l, f) is the coefficient of x^k in z^l f(x)
+
+ Exports == with
+
+ generalInterpolation: (List D, CoeffAction, Vector VF, List NonNegativeInteger)
+ -> Matrix SUP D
+ ++ \spad{generalInterpolation(l, CA, f, eta)} performs Hermite-Pade
+ ++ approximation using the given action CA of polynomials on the elements
+ ++ of f. The result is guaranteed to be correct up to order
+ ++ |eta|-1. Given that eta is a "normal" point, the degrees on the
+ ++ diagonal are given by eta. The degrees of column i are in this case
+ ++ eta + e.i - [1,1,...,1], where the degree of zero is -1.
+
+ generalInterpolation: (List D, CoeffAction,
+ Vector VF, NonNegativeInteger, NonNegativeInteger)
+ -> Stream Matrix SUP D
+ ++ \spad{generalInterpolation(l, CA, f, sumEta, maxEta)} applies
+ ++ generalInterpolation(l, CA, f, eta) for all possible eta with maximal
+ ++ entry maxEta and sum of entries sumEta
+
+ Implementation == add
+
+ multiplyRows!(v: Vector D, M: Matrix SUP D): Matrix SUP D ==
+ n := #v
+ for i in 1..n repeat
+ for j in 1..n repeat
+ M.(i,j) := v.i*M.(i,j)
+
+ M
+
+ generalInterpolation(C: List D, coeffAction: CoeffAction,
+ f: Vector VF, eta: List NonNegativeInteger): Matrix SUP D ==
+ n := #f
+ g: Vector V := new(n, 0)
+ den: Vector D := new(n, 0)
+
+ for i in 1..n repeat
+ c := coefficients(f.i)
+ den.i := commonDenominator(c)$CommonDenominator(D, F, List F)
+ g.i := map(retract(#1*den.i)@D, f.i)
+ $FAMR2(NonNegativeInteger, Fraction D, VF, D, V)
+
+ M := generalInterpolation(C, coeffAction, g, eta)$FFFG(D, V)
+
+-- The following is necessary since I'm multiplying each row with a factor, not
+-- each column. Possibly I could factor out gcd den, but I'm not sure whether
+-- this is efficient.
+
+ multiplyRows!(den, M)
+
+ generalInterpolation(C: List D, coeffAction: CoeffAction,
+ f: Vector VF,
+ sumEta: NonNegativeInteger,
+ maxEta: NonNegativeInteger)
+ : Stream Matrix SUP D ==
+
+ n := #f
+ g: Vector V := new(n, 0)
+ den: Vector D := new(n, 0)
+
+ for i in 1..n repeat
+ c := coefficients(f.i)
+ den.i := commonDenominator(c)$CommonDenominator(D, F, List F)
+ g.i := map(retract(#1*den.i)@D, f.i)
+ $FAMR2(NonNegativeInteger, Fraction D, VF, D, V)
+
+ c: cFunction := generalCoefficient(coeffAction, g,
+ (#1-1)::NonNegativeInteger, #2)$FFFG(D, V)
+
+
+ MS: Stream Matrix SUP D
+ := generalInterpolation(C, coeffAction, g, sumEta, maxEta)$FFFG(D, V)
+
+-- The following is necessary since I'm multiplying each row with a factor, not
+-- each column. Possibly I could factor out gcd den, but I'm not sure whether
+-- this is efficient.
+
+ map(multiplyRows!(den, #1), MS)$Stream(Matrix SUP D)
+
+@
+
+
+\section{package NEWTON NewtonInterpolation}
+<<package NEWTON NewtonInterpolation>>=
+)abbrev package NEWTON NewtonInterpolation
+++ Description:
+++ This package exports Newton interpolation for the special case where the
+++ result is known to be in the original integral domain
+NewtonInterpolation F: Exports == Implementation where
+ F: IntegralDomain
+ Exports == with
+
+ newton: List F -> SparseUnivariatePolynomial F
+
+ ++ \spad{newton}(l) returns the interpolating polynomial for the values
+ ++ l, where the x-coordinates are assumed to be [1,2,3,...,n] and the
+ ++ coefficients of the interpolating polynomial are known to be in the
+ ++ domain F. I.e., it is a very streamlined version for a special case of
+ ++ interpolation.
+
+ Implementation == add
+
+ differences(yl: List F): List F ==
+ [y2-y1 for y1 in yl for y2 in rest yl]
+
+ z: SparseUnivariatePolynomial(F) := monomial(1,1)
+
+-- we assume x=[1,2,3,...,n]
+ newtonAux(k: F, fact: F, yl: List F): SparseUnivariatePolynomial(F) ==
+ if empty? rest yl
+ then ((yl.1) exquo fact)::F::SparseUnivariatePolynomial(F)
+ else ((yl.1) exquo fact)::F::SparseUnivariatePolynomial(F)
+ + (z-k::SparseUnivariatePolynomial(F)) _
+ * newtonAux(k+1$F, fact*k, differences yl)
+
+
+ newton yl == newtonAux(1$F, 1$F, yl)
+
+@
+%$
+
+\section{License}
+<<license>>=
+--Copyright (c) 2006-2007, Martin Rubey <Mar...@un...>
+--
+--Redistribution and use in source and binary forms, with or without
+--modification, are permitted provided that the following conditions are
+--met:
+--
+-- - Redistributions of source code must retain the above copyright
+-- notice, this list of conditions and the following disclaimer.
+--
+-- - Redistributions in binary form must reproduce the above copyright
+-- notice, this list of conditions and the following disclaimer in
+-- the documentation and/or other materials provided with the
+-- distribution.
+--
+--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+@
+
+<<*>>=
+<<license>>
+
+<<package FAMR2 FiniteAbelianMonoidRingFunctions2>>
+<<package FFFG FractionFreeFastGaussian>>
+<<package FFFGF FractionFreeFastGaussianFractions>>
+<<package NEWTON NewtonInterpolation>>
+@
+\end{document}
Added: trunk/axiom/src/algebra/mantepse.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/mantepse.spad.pamphlet (rev 0)
+++ trunk/axiom/src/algebra/mantepse.spad.pamphlet 2007-09-11 05:08:07 UTC (rev 742)
@@ -0,0 +1,2390 @@
+\documentclass{article}
+\usepackage{axiom,amsthm,amsmath}
+\newtheorem{ToDo}{ToDo}[section]
+
+\newcommand{\Axiom}{{\tt Axiom}}
+\newcommand{\Rate}{{\tt Rate}}
+\newcommand{\GFUN}{{\tt GFUN}}
+\begin{document}
+\title{mantepse.spad}
+\author{Martin Rubey}
+\maketitle
+\begin{abstract}
+ The packages defined in this file enable {\Axiom} to guess formulas for
+ sequences of, for example, rational numbers or rational functions, given the
+ first few terms. It extends and complements Christian Krattenthaler's
+ program \Rate\ and the relevant parts of Bruno Salvy and Paul Zimmermann's
+ \GFUN.
+\end{abstract}
+\tableofcontents
+\section{domain UFPS UnivariateFormalPowerSeries}
+<<domain UFPS UnivariateFormalPowerSeries>>=
+)abbrev domain UFPS UnivariateFormalPowerSeries
+UnivariateFormalPowerSeries(Coef: Ring) ==
+ UnivariateTaylorSeries(Coef, 'x, 0$Coef)
+
+@
+%$
+
+\section{package UFPS1 UnivariateFormalPowerSeriesFunctions}
+<<package UFPS1 UnivariateFormalPowerSeriesFunctions>>=
+)abbrev package UFPS1 UnivariateFormalPowerSeriesFunctions
+UnivariateFormalPowerSeriesFunctions(Coef: Ring): Exports == Implementation
+ where
+
+ UFPS ==> UnivariateFormalPowerSeries Coef
+
+ Exports == with
+
+ hadamard: (UFPS, UFPS) -> UFPS
+
+ Implementation == add
+
+ hadamard(f, g) ==
+ series map(#1*#2, coefficients f, coefficients g)
+ $StreamFunctions3(Coef, Coef, Coef)
+
+@
+%$
+
+\section{domain GOPT GuessOption}
+<<domain GOPT GuessOption>>=
+)abbrev domain GOPT GuessOption
+++ Author: Martin Rubey
+++ Description: GuessOption is a domain whose elements are various options used
+++ by \spadtype{Guess}.
+GuessOption(): Exports == Implementation where
+
+ Exports == SetCategory with
+
+ maxDerivative: Integer -> %
+ ++ maxDerivative(d) specifies the maximum derivative in an algebraic
+ ++ differential equation. maxDerivative(-1) specifies that the maximum
+ ++ derivative can be arbitrary. This option is expressed in the form
+ ++ \spad{maxDerivative == d}.
+
+ maxShift: Integer -> %
+ ++ maxShift(d) specifies the maximum shift in a recurrence
+ ++ equation. maxShift(-1) specifies that the maximum shift can be
+ ++ arbitrary. This option is expressed in the form \spad{maxShift == d}.
+
+ maxPower: Integer -> %
+ ++ maxPower(d) specifies the maximum degree in an algebraic differential
+ ++ equation. For example, the degree of (f'')^3 f' is 4. maxPower(-1)
+ ++ specifies that the maximum exponent can be arbitrary. This option is
+ ++ expressed in the form \spad{maxPower == d}.
+
+ homogeneous: Boolean -> %
+ ++ homogeneous(d) specifies whether we allow only homogeneous algebraic
+ ++ differential equations. This option is expressed in the form
+ ++ \spad{homogeneous == d}.
+
+ maxLevel: Integer -> %
+ ++ maxLevel(d) specifies the maximum number of recursion levels operators
+ ++ guessProduct and guessSum will be applied. maxLevel(-1) specifies
+ ++ that all levels are tried. This option is expressed in the form
+ ++ \spad{maxLevel == d}.
+
+ maxDegree: Integer -> %
+ ++ maxDegree(d) specifies the maximum degree of the coefficient
+ ++ polynomials in an algebraic differential equation or a recursion with
+ ++ polynomial coefficients. For rational functions with an exponential
+ ++ term, \spad{maxDegree} bounds the degree of the denominator
+ ++ polynomial.
+ ++ maxDegree(-1) specifies that the maximum
+ ++ degree can be arbitrary. This option is expressed in the form
+ ++ \spad{maxDegree == d}.
+
+ allDegrees: Boolean -> %
+ ++ allDegrees(d) specifies whether all possibilities of the degree vector
+ ++ - taking into account maxDegree - should be tried. This is mainly
+ ++ interesting for rational interpolation. This option is expressed in
+ ++ the form \spad{allDegrees == d}.
+
+ safety: NonNegativeInteger -> %
+ ++ safety(d) specifies the number of values reserved for testing any
+ ++ solutions found. This option is expressed in the form \spad{safety ==
+ ++ d}.
+
+ one: Boolean -> %
+ ++ one(d) specifies whether we are happy with one solution. This option
+ ++ is expressed in the form \spad{ d}.
+
+ debug: Boolean -> %
+ ++ debug(d) specifies whether we want additional output on the
+ ++ progress. This option is expressed in the form \spad{debug == d}.
+
+ functionName: Symbol -> %
+ ++ functionName(d) specifies the name of the function given by the
+ ++ algebraic differential equation or recurrence. This option is
+ ++ expressed in the form \spad{functionName == d}.
+
+ variableName: Symbol -> %
+ ++ variableName(d) specifies the variable used in by the algebraic
+ ++ differential equation. This option is expressed in the form
+ ++ \spad{variableName == d}.
+
+ indexName: Symbol -> %
+ ++ indexName(d) specifies the index variable used for the formulas. This
+ ++ option is expressed in the form \spad{indexName == d}.
+
+ displayAsGF: Boolean -> %
+ ++ displayAsGF(d) specifies whether the result is a generating function
+ ++ or a recurrence. This option should not be set by the user, but rather
+ ++ by the HP-specification.
+
+ option : (List %, Symbol) -> Union(Any, "failed")
+ ++ option() is not to be used at the top level;
+ ++ option determines internally which drawing options are indicated in
+ ++ a draw command.
+
+ option?: (List %, Symbol) -> Boolean
+ ++ option?() is not to be used at the top level;
+ ++ option? internally returns true for drawing options which are
+ ++ indicated in a draw command, or false for those which are not.
+
+ checkOptions: List % -> Void
+ ++ checkOptions checks whether an option is given twice
+
+ Implementation ==> add
+ import AnyFunctions1(Boolean)
+ import AnyFunctions1(Symbol)
+ import AnyFunctions1(Integer)
+ import AnyFunctions1(NonNegativeInteger)
+
+ Rep := Record(keyword: Symbol, value: Any)
+
+ maxLevel d == ["maxLevel"::Symbol, d::Any]
+ maxDerivative d == ["maxDerivative"::Symbol, d::Any]
+ maxShift d == maxDerivative d
+ maxDegree d == ["maxDegree"::Symbol, d::Any]
+ allDegrees d == ["allDegrees"::Symbol, d::Any]
+ maxPower d == ["maxPower"::Symbol, d::Any]
+ safety d == ["safety"::Symbol, d::Any]
+ homogeneous d == ["homogeneous"::Symbol, d::Any]
+ debug d == ["debug"::Symbol, d::Any]
+ one d == ["one"::Symbol, d::Any]
+ functionName d == ["functionName"::Symbol, d::Any]
+ variableName d == ["variableName"::Symbol, d::Any]
+ indexName d == ["indexName"::Symbol, d::Any]
+ displayAsGF d == ["displayAsGF"::Symbol, d::Any]
+
+ coerce(x:%):OutputForm == x.keyword::OutputForm = x.value::OutputForm
+ x:% = y:% == x.keyword = y.keyword and x.value = y.value
+
+ option?(l, s) ==
+ for x in l repeat
+ x.keyword = s => return true
+ false
+
+ option(l, s) ==
+ for x in l repeat
+ x.keyword = s => return(x.value)
+ "failed"
+
+ checkOptions l ==
+ if not empty? l then
+ if find((first l).keyword = #1.keyword, rest l) case "failed"
+ then checkOptions rest l
+ else error "GuessOption: Option specified twice"
+
+@
+
+\section{package GOPT0 GuessOptionFunctions0}
+<<package GOPT0 GuessOptionFunctions0>>=
+)abbrev package GOPT0 GuessOptionFunctions0
+++ Author: Martin Rubey
+++ Description: GuessOptionFunctions0 provides operations that extract the
+++ values of options for \spadtype{Guess}.
+GuessOptionFunctions0(): Exports == Implementation where
+
+ LGOPT ==> List GuessOption
+
+ Exports == SetCategory with
+
+ maxLevel: LGOPT -> Integer
+ ++ maxLevel returns the specified maxLevel or -1 as default.
+
+ maxPower: LGOPT -> Integer
+ ++ maxPower returns the specified maxPower or -1 as default.
+
+ maxDerivative: LGOPT -> Integer
+ ++ maxDerivative returns the specified maxDerivative or -1 as default.
+
+ maxShift: LGOPT -> Integer
+ ++ maxShift returns the specified maxShift or -1 as default.
+
+ maxDegree: LGOPT -> Integer
+ ++ maxDegree returns the specified maxDegree or -1 as default.
+
+ allDegrees: LGOPT -> Boolean
+ ++ allDegrees returns whether all possibilities of the degree vector
+ ++ should be tried, the default being false.
+
+ safety: LGOPT -> NonNegativeInteger
+ ++ safety returns the specified safety or 1 as default.
+
+ one: LGOPT -> Boolean
+ ++ one returns whether we need only one solution, default being true.
+
+ homogeneous: LGOPT -> Boolean
+ ++ homogeneous returns whether we allow only homogeneous algebraic
+ ++ differential equations, default being false
+
+ functionName: LGOPT -> Symbol
+ ++ functionName returns the name of the function given by the algebraic
+ ++ differential equation, default being f
+
+ variableName: LGOPT -> Symbol
+ ++ variableName returns the name of the variable used in by the
+ ++ algebraic differential equation, default being x
+
+ indexName: LGOPT -> Symbol
+ ++ indexName returns the name of the index variable used for the
+ ++ formulas, default being n
+
+ displayAsGF: LGOPT -> Boolean
+ ++ displayAsGF specifies whether the result is a generating function
+ ++ or a recurrence. This option should not be set by the user, but rather
+ ++ by the HP-specification, therefore, there is no default.
+
+ debug: LGOPT -> Boolean
+ ++ debug returns whether we want additional output on the progress,
+ ++ default being false
+
+ Implementation == add
+
+ maxLevel l ==
+ if (opt := option(l, "maxLevel" :: Symbol)) case "failed" then
+ -1
+ else
+ retract(opt :: Any)$AnyFunctions1(Integer)
+
+ maxDerivative l ==
+ if (opt := option(l, "maxDerivative" :: Symbol)) case "failed" then
+ -1
+ else
+ retract(opt :: Any)$AnyFunctions1(Integer)
+
+ maxShift l == maxDerivative l
+
+ maxDegree l ==
+ if (opt := option(l, "maxDegree" :: Symbol)) case "failed" then
+ -1
+ else
+ retract(opt :: Any)$AnyFunctions1(Integer)
+
+ allDegrees l ==
+ if (opt := option(l, "allDegrees" :: Symbol)) case "failed" then
+ false
+ else
+ retract(opt :: Any)$AnyFunctions1(Boolean)
+
+ maxPower l ==
+ if (opt := option(l, "maxPower" :: Symbol)) case "failed" then
+ -1
+ else
+ retract(opt :: Any)$AnyFunctions1(Integer)
+
+ safety l ==
+ if (opt := option(l, "safety" :: Symbol)) case "failed" then
+ 1$NonNegativeInteger
+ else
+ retract(opt :: Any)$AnyFunctions1(Integer)::NonNegativeInteger
+
+ one l ==
+ if (opt := option(l, "one" :: Symbol)) case "failed" then
+ true
+ else
+ retract(opt :: Any)$AnyFunctions1(Boolean)
+
+ debug l ==
+ if (opt := option(l, "debug" :: Symbol)) case "failed" then
+ false
+ else
+ retract(opt :: Any)$AnyFunctions1(Boolean)
+
+ homogeneous l ==
+ if (opt := option(l, "homogeneous" :: Symbol)) case "failed" then
+ false
+ else
+ retract(opt :: Any)$AnyFunctions1(Boolean)
+
+ variableName l ==
+ if (opt := option(l, "variableName" :: Symbol)) case "failed" then
+ "x" :: Symbol
+ else
+ retract(opt :: Any)$AnyFunctions1(Symbol)
+
+ functionName l ==
+ if (opt := option(l, "functionName" :: Symbol)) case "failed" then
+ "f" :: Symbol
+ else
+ retract(opt :: Any)$AnyFunctions1(Symbol)
+
+ indexName l ==
+ if (opt := option(l, "indexName" :: Symbol)) case "failed" then
+ "n" :: Symbol
+ else
+ retract(opt :: Any)$AnyFunctions1(Symbol)
+
+ displayAsGF l ==
+ if (opt := option(l, "displayAsGF" :: Symbol)) case "failed" then
+ error "GuessOption: displayAsGF not set"
+ else
+ retract(opt :: Any)$AnyFunctions1(Boolean)
+
+@
+
+\section{package GUESS Guess}
+<<package GUESS Guess>>=
+)abbrev package GUESS Guess
+++ Author: Martin Rubey
+++ Description: This package implements guessing of sequences. Packages for the
+++ most common cases are provided as \spadtype{GuessInteger},
+++ \spadtype{GuessPolynomial}, etc.
+Guess(F, S, EXPRR, R, retract, coerce): Exports == Implementation where
+ F: Field -- zB.: FRAC POLY PF 5
+-- in F wird interpoliert und gerechnet
+
+ S: GcdDomain
+
+-- in guessExpRat möchte ich die Wurzeln von Polynomen über F bestimmen. Wenn F
+-- ein Quotientenkörper ist, dann kann ich den Nenner loswerden.
+-- F wäre dann also QFCAT S
+
+ R: Join(OrderedSet, IntegralDomain) -- zB.: FRAC POLY INT
+
+-- die Ergebnisse werden aber in EXPRR dargestellt
+-- EXPRR: Join(ExpressionSpace, IntegralDomain,
+ EXPRR: Join(FunctionSpace Integer, IntegralDomain,
+ RetractableTo R, RetractableTo Symbol,
+ RetractableTo Integer, CombinatorialOpsCategory,
+ PartialDifferentialRing Symbol) with
+ _* : (%,%) -> %
+ _/ : (%,%) -> %
+ _*_* : (%,%) -> %
+ numerator : % -> %
+ denominator : % -> %
+ ground? : % -> Boolean
+
+ -- zB.: EXPR INT
+-- EXPR FRAC POLY INT ist ja verboten. Deshalb kann ich nicht einfach EXPR R
+-- verwenden
+-- EXPRR existiert, falls irgendwann einmal EXPR PF 5 unterstützt wird.
+
+
+-- das folgende möchte ich eigentlich loswerden.
+
+ retract: R -> F -- zB.: i+->i
+ coerce: F -> EXPRR -- zB.: i+->i
+-- Achtung EXPRR ~= EXPR R
+
+ LGOPT ==> List GuessOption
+ GOPT0 ==> GuessOptionFunctions0
+
+ NNI ==> NonNegativeInteger
+ PI ==> PositiveInteger
+ EXPRI ==> Expression Integer
+ GUESSRESULT ==> List Record(function: EXPRR, order: NNI)
+
+ UFPSF ==> UnivariateFormalPowerSeries F
+ UFPS1 ==> UnivariateFormalPowerSeriesFunctions
+
+ UFPSS ==> UnivariateFormalPowerSeries S
+
+ SUP ==> SparseUnivariatePolynomial
+
+ UFPSSUPF ==> UnivariateFormalPowerSeries SUP F
+
+ FFFG ==> FractionFreeFastGaussian
+ FFFGF ==> FractionFreeFastGaussianFractions
+
+-- CoeffAction
+ DIFFSPECA ==> (NNI, NNI, SUP S) -> S
+
+ DIFFSPECAF ==> (NNI, NNI, UFPSSUPF) -> SUP F
+
+ DIFFSPECAX ==> (NNI, Symbol, EXPRR) -> EXPRR
+
+-- the diagonal of the C-matrix
+ DIFFSPECC ==> NNI -> List S
+
+
+ HPSPEC ==> Record(guessStream: UFPSF -> Stream UFPSF,
+ degreeStream: Stream NNI,
+ testStream: UFPSSUPF -> Stream UFPSSUPF,
+ exprStream: (EXPRR, Symbol) -> Stream EXPRR,
+ A: DIFFSPECA,
+ AF: DIFFSPECAF,
+ AX: DIFFSPECAX,
+ C: DIFFSPECC)
+
+-- note that empty?(guessStream.o) has to return always. In other words, if the
+-- stream is finite, empty? should recognize it.
+
+ DIFFSPECN ==> EXPRR -> EXPRR -- eg.: i+->q^i
+
+ GUESSER ==> (List F, LGOPT) -> GUESSRESULT
+
+ FSUPS ==> Fraction SUP S
+ FSUPF ==> Fraction SUP F
+
+ V ==> OrderedVariableList(['a1, 'A])
+ POLYF ==> SparseMultivariatePolynomial(F, V)
+ FPOLYF ==> Fraction POLYF
+ FSUPFPOLYF ==> Fraction SUP FPOLYF
+ POLYS ==> SparseMultivariatePolynomial(S, V)
+ FPOLYS ==> Fraction POLYS
+ FSUPFPOLYS ==> Fraction SUP FPOLYS
+
+@
+
+The following should be commented out when not debugging, see also
+Section~\ref{sec:Hermite-Pade}.
+
+<<package GUESS Guess>>=
+ --@<<implementation: Guess - Hermite-Pade - Types for Operators>>
+@
+
+<<package GUESS Guess>>=
+ Exports == with
+
+ guess: List F -> GUESSRESULT
+ ++ \spad{guess l} applies recursively \spadfun{guessRec} and
+ ++ \spadfun{guessADE} to the successive differences and quotients of
+ ++ the list. Default options as described in
+ ++ \spadtype{GuessOptionFunctions0} are used.
+
+ guess: (List F, LGOPT) -> GUESSRESULT
+ ++ \spad{guess(l, options)} applies recursively \spadfun{guessRec}
+ ++ and \spadfun{guessADE} to the successive differences and quotients
+ ++ of the list. The given options are used.
+
+ guess: (List F, List GUESSER, List Symbol) -> GUESSRESULT
+ ++ \spad{guess(l, guessers, ops)} applies recursively the given
+ ++ guessers to the successive differences if ops contains the symbol
+ ++ guessSum and quotients if ops contains the symbol guessProduct to
+ ++ the list. Default options as described in
+ ++ \spadtype{GuessOptionFunctions0} are used.
+
+ guess: (List F, List GUESSER, List Symbol, LGOPT) -> GUESSRESULT
+ ++ \spad{guess(l, guessers, ops)} applies recursively the given
+ ++ guessers to the successive differences if ops contains the symbol
+ ++ \spad{guessSum} and quotients if ops contains the symbol
+ ++ \spad{guessProduct} to the list. The given options are used.
+
+ guessExpRat: List F -> GUESSRESULT
+ ++ \spad{guessExpRat l} tries to find a function of the form
+ ++ n+->(a+b n)^n r(n), where r(n) is a rational function, that fits
+ ++ l.
+
+ guessExpRat: (List F, LGOPT) -> GUESSRESULT
+ ++ \spad{guessExpRat(l, options)} tries to find a function of the
+ ++ form n+->(a+b n)^n r(n), where r(n) is a rational function, that
+ ++ fits l.
+
+ if F has RetractableTo Symbol and S has RetractableTo Symbol then
+
+ guessExpRat: Symbol -> GUESSER
+ ++ \spad{guessExpRat q} returns a guesser that tries to find a
+ ++ function of the form n+->(a+b q^n)^n r(q^n), where r(n) is a
+ ++ rational function, that fits l.
+
+ guessHP: (LGOPT -> HPSPEC) -> GUESSER
+ ++ \spad{guessHP f} constructs an operation that applies Hermite-Pade
+ ++ approximation to the series generated by the given function f.
+
+ guessADE: List F -> GUESSRESULT
+ ++ \spad{guessADE l} tries to find an algebraic differential equation
+ ++ for a generating function whose first Taylor coefficients are
+ ++ given by l, using the default options described in
+ ++ \spadtype{GuessOptionFunctions0}.
+
+ guessADE: (List F, LGOPT) -> GUESSRESULT
+ ++ \spad{guessADE(l, options)} tries to find an algebraic
+ ++ differential equation for a generating function whose first Taylor
+ ++ coefficients are given by l, using the given options.
+
+ guessAlg: List F -> GUESSRESULT
+ ++ \spad{guessAlg l} tries to find an algebraic equation for a
+ ++ generating function whose first Taylor coefficients are given by
+ ++ l, using the default options described in
+ ++ \spadtype{GuessOptionFunctions0}. It is equivalent to
+ ++ \spadfun{guessADE}(l, maxDerivative == 0).
+
+ guessAlg: (List F, LGOPT) -> GUESSRESULT
+ ++ \spad{guessAlg(l, options)} tries to find an algebraic equation
+ ++ for a generating function whose first Taylor coefficients are
+ ++ given by l, using the given options. It is equivalent to
+ ++ \spadfun{guessADE}(l, options) with \spad{maxDerivative == 0}.
+
+ guessHolo: List F -> GUESSRESULT
+ ++ \spad{guessHolo l} tries to find an ordinary linear differential
+ ++ equation for a generating function whose first Taylor coefficients
+ ++ are given by l, using the default options described in
+ ++ \spadtype{GuessOptionFunctions0}. It is equivalent to
+ ++ \spadfun{guessADE}\spad{(l, maxPower == 1)}.
+
+ guessHolo: (List F, LGOPT) -> GUESSRESULT
+ ++ \spad{guessHolo(l, options)} tries to find an ordinary linear
+ ++ differential equation for a generating function whose first Taylor
+ ++ coefficients are given by l, using the given options. It is
+ ++ equivalent to \spadfun{guessADE}\spad{(l, options)} with
+ ++ \spad{maxPower == 1}.
+
+ guessPade: (List F, LGOPT) -> GUESSRESULT
+ ++ \spad{guessPade(l, options)} tries to find a rational function
+ ++ whose first Taylor coefficients are given by l, using the given
+ ++ options. It is equivalent to \spadfun{guessADE}\spad{(l,
+ ++ maxDerivative == 0, maxPower == 1, allDegrees == true)}.
+
+ guessPade: List F -> GUESSRESULT
+ ++ \spad{guessPade(l, options)} tries to find a rational function
+ ++ whose first Taylor coefficients are given by l, using the default
+ ++ options described in \spadtype{GuessOptionFunctions0}. It is
+ ++ equivalent to \spadfun{guessADE}\spad{(l, options)} with
+ ++ \spad{maxDerivative == 0, maxPower == 1, allDegrees == true}.
+
+ guessRec: List F -> GUESSRESULT
+ ++ \spad{guessRec l} tries to find an ordinary difference equation
+ ++ whose first values are given by l, using the default options
+ ++ described in \spadtype{GuessOptionFunctions0}.
+
+ guessRec: (List F, LGOPT) -> GUESSRESULT
+ ++ \spad{guessRec(l, options)} tries to find an ordinary difference
+ ++ equation whose first values are given by l, using the given
+ ++ options.
+
+ guessPRec: (List F, LGOPT) -> GUESSRESULT
+ ++ \spad{guessPRec(l, options)} tries to find a linear recurrence
+ ++ with polynomial coefficients whose first values are given by l,
+ ++ using the given options. It is equivalent to
+ ++ \spadfun{guessRec}\spad{(l, options)} with \spad{maxPower == 1}.
+
+ guessPRec: List F -> GUESSRESULT
+ ++ \spad{guessPRec l} tries to find a linear recurrence with
+ ++ polynomial coefficients whose first values are given by l, using
+ ++ the default options described in
+ ++ \spadtype{GuessOptionFunctions0}. It is equivalent to
+ ++ \spadfun{guessRec}\spad{(l, maxPower == 1)}.
+
+ guessRat: (List F, LGOPT) -> GUESSRESULT
+ ++ \spad{guessRat(l, options)} tries to find a rational function
+ ++ whose first values are given by l, using the given options. It is
+ ++ equivalent to \spadfun{guessRec}\spad{(l, maxShift == 0, maxPower
+ ++ == 1, allDegrees == true)}.
+
+ guessRat: List F -> GUESSRESULT
+ ++ \spad{guessRat l} tries to find a rational function whose first
+ ++ values are given by l, using the default options described in
+ ++ \spadtype{GuessOptionFunctions0}. It is equivalent to
+ ++ \spadfun{guessRec}\spad{(l, maxShift == 0, maxPower == 1,
+ ++ allDegrees == true)}.
+
+ diffHP: LGOPT -> HPSPEC
+ ++ \spad{diffHP options} returns a specification for Hermite-Pade
+ ++ approximation with the differential operator
+
+ shiftHP: LGOPT -> HPSPEC
+ ++ \spad{shiftHP options} returns a specification for Hermite-Pade
+ ++ approximation with the shift operator
+
+ if F has RetractableTo Symbol and S has RetractableTo Symbol then
+ shiftHP: Symbol -> (LGOPT -> HPSPEC)
+ ++ \spad{shiftHP options} returns a specification for
+ ++ Hermite-Pade approximation with the $q$-shift operator
+
+ diffHP: Symbol -> (LGOPT -> HPSPEC)
+ ++ \spad{diffHP options} returns a specification for Hermite-Pade
+ ++ approximation with the $q$-dilation operator
+
+ guessRec: Symbol -> GUESSER
+ ++ \spad{guessRec q} returns a guesser that finds an ordinary
+ ++ q-difference equation whose first values are given by l, using
+ ++ the given options.
+
+ guessPRec: Symbol -> GUESSER
+ ++ \spad{guessPRec q} returns a guesser that tries to find
+ ++ a linear q-recurrence with polynomial coefficients whose first
+ ++ values are given by l, using the given options. It is
+ ++ equivalent to \spadfun{guessRec}\spad{(q)} with
+ ++ \spad{maxPower == 1}.
+
+ guessRat: Symbol -> GUESSER
+ ++ \spad{guessRat q} returns a guesser that tries to find a
+ ++ q-rational function whose first values are given by l, using
+ ++ the given options. It is equivalent to \spadfun{guessRec} with
+ ++ \spad{(l, maxShift == 0, maxPower == 1, allDegrees == true)}.
+
+ guessADE: Symbol -> GUESSER
+ ++ \spad{guessADE q} returns a guesser that tries to find an
+ ++ algebraic differential equation for a generating function whose
+ ++ first Taylor coefficients are given by l, using the given
+ ++ options.
+
+ --@<<debug exports: Guess>>
+@
+
+<<debug exports: Guess>>=
+termAsUFPSF: (UFPSF, List Integer, DIFFSPECS, DIFFSPEC1) -> UFPSF
+termAsUFPSF2: (UFPSF, List Integer, DIFFSPECS, DIFFSPEC1) -> UFPSF
+termAsEXPRR: (EXPRR, Symbol, List Integer, DIFFSPECX, DIFFSPEC1X) -> EXPRR
+termAsUFPSSUPF: (UFPSSUPF, List Integer, DIFFSPECSF, DIFFSPEC1F) -> UFPSSUPF
+termAsUFPSSUPF2: (UFPSSUPF, List Integer, DIFFSPECSF, DIFFSPEC1F) -> UFPSSUPF
+
+ShiftSXGF: (EXPRR, Symbol, NNI) -> EXPRR
+ShiftAXGF: (NNI, Symbol, EXPRR) -> EXPRR
+
+FilteredPartitionStream: LGOPT -> Stream List Integer
+
+guessInterpolate: (List SUP F, List NNI, HPSPEC) -> Matrix SUP S
+testInterpolant: (List SUP S, List F, List UFPSSUPF, List EXPRR, List EXPRR, _
+ NNI, HPSPEC, Symbol, BasicOperator, LGOPT) _
+ -> Union("failed", Record(function: EXPRR, order: NNI))
+
+@
+
+<<package GUESS Guess>>=
+ Implementation == add
+ <<implementation: Guess>>
+@
+
+<<implementation: Guess>>=
+
+-- We have to put this chunk at the beginning, because otherwise it will take
+-- very long to compile.
+
+<<implementation: Guess - guessExpRat - Order and Degree>>
+
+<<implementation: Guess - Utilities>>
+<<implementation: Guess - guessExpRat>>
+<<implementation: Guess - Hermite-Pade>>
+<<implementation: Guess - guess>>
+@
+
+\subsection{general utilities}
+
+<<implementation: Guess - Utilities>>=
+checkResult(res: EXPRR, n: Symbol, l: Integer, list: List F,
+ options: LGOPT): NNI ==
+ for i in l..1 by -1 repeat
+ den := eval(denominator res, n::EXPRR, (i-1)::EXPRR)
+ if den = 0 then return i::NNI
+ num := eval(numerator res, n::EXPRR, (i-1)::EXPRR)
+ if list.i ~= retract(retract(num/den)@R)
+ then return i::NNI
+ 0$NNI
+
+SUPS2SUPF(p: SUP S): SUP F ==
+ if F is S then
+ p pretend SUP(F)
+ else if F is Fraction S then
+ map(coerce(#1)$Fraction(S), p)
+ $SparseUnivariatePolynomialFunctions2(S, F)
+ else error "Type parameter F should be either equal to S or equal _
+ to Fraction S"
+
+@
+
+\subsection{guessing rational functions with an exponential term}
+
+<<implementation: Guess - guessExpRat>>=
+<<implementation: Guess - guessExpRat - Utilities>>
+<<implementation: Guess - guessExpRat - Main>>
+@
+
+\subsubsection{Utilities}
+
+\paragraph{conversion routines}
+
+<<implementation: Guess - guessExpRat - Utilities>>=
+F2FPOLYS(p: F): FPOLYS ==
+ if F is S then
+ p::POLYF::FPOLYF pretend FPOLYS
+ else if F is Fraction S then
+ numer(p)$Fraction(S)::POLYS/denom(p)$Fraction(S)::POLYS
+ else error "Type parameter F should be either equal to S or equal _
+ to Fraction S"
+
+MPCSF ==> MPolyCatFunctions2(V, IndexedExponents V,
+ IndexedExponents V, S, F,
+ POLYS, POLYF)
+
+SUPF2EXPRR(xx: Symbol, p: SUP F): EXPRR ==
+ zero? p => 0
+ (coerce(leadingCoefficient p))::EXPRR * (xx::EXPRR)**degree p
+ + SUPF2EXPRR(xx, reductum p)
+
+FSUPF2EXPRR(xx: Symbol, p: FSUPF): EXPRR ==
+ (SUPF2EXPRR(xx, numer p)) / (SUPF2EXPRR(xx, denom p))
+
+
+POLYS2POLYF(p: POLYS): POLYF ==
+ if F is S then
+ p pretend POLYF
+ else if F is Fraction S then
+ map(coerce(#1)$Fraction(S), p)$MPCSF
+ else error "Type parameter F should be either equal to S or equal _
+ to Fraction S"
+
+SUPPOLYS2SUPF(p: SUP POLYS, a1v: F, Av: F): SUP F ==
+ zero? p => 0
+ lc: POLYF := POLYS2POLYF leadingCoefficient(p)
+ monomial(retract(eval(lc, [index(1)$V, index(2)$V]::List V,
+ [a1v, Av])),
+ degree p)
+ + SUPPOLYS2SUPF(reductum p, a1v, Av)
+
+
+SUPFPOLYS2FSUPPOLYS(p: SUP FPOLYS): Fraction SUP POLYS ==
+ cden := splitDenominator(p)
+ $UnivariatePolynomialCommonDenominator(POLYS, FPOLYS,SUP FPOLYS)
+
+ pnum: SUP POLYS
+ := map(retract(#1 * cden.den)$FPOLYS, p)
+ $SparseUnivariatePolynomialFunctions2(FPOLYS, POLYS)
+ pden: SUP POLYS := (cden.den)::SUP POLYS
+
+ pnum/pden
+
+
+POLYF2EXPRR(p: POLYF): EXPRR ==
+ map(convert(#1)@Symbol::EXPRR, coerce(#1)@EXPRR, p)
+ $PolynomialCategoryLifting(IndexedExponents V, V,
+ F, POLYF, EXPRR)
+
+@
+
+\paragraph{mimicking $q$-anaologa}
+
+<<implementation: Guess - guessExpRat - Utilities>>=
+defaultD: DIFFSPECN
+defaultD(expr: EXPRR): EXPRR == expr
+
+-- applies n+->q^n or whatever DN is to i
+DN2DL: (DIFFSPECN, Integer) -> F
+DN2DL(DN, i) == retract(retract(DN(i::EXPRR))@R)
+
+<<implementation: Guess - guessExpRat - evalResultant>>
+<<implementation: Guess - guessExpRat - substitute>>
+@
+
+\subsubsection{reducing the degree of the polynomials}
+
+The degree of [[poly3]] is governed by $(a0+x_m a_1)^{x_m}$. Therefore, we
+substitute $A-x_m a1$ for $a$, which reduces the degree in $a_1$ by
+$x_m-x_{i+1}$.
+
+<<implementation: Guess - guessExpRat - substitute>>=
+p(xm: Integer, i: Integer, va1: V, vA: V, basis: DIFFSPECN): FPOLYS ==
+ vA::POLYS::FPOLYS + va1::POLYS::FPOLYS _
+ * F2FPOLYS(DN2DL(basis, i) - DN2DL(basis, xm))
+
+p2(xm: Integer, i: Symbol, a1v: F, Av: F, basis: DIFFSPECN): EXPRR ==
+ coerce(Av) + coerce(a1v)*(basis(i::EXPRR) - basis(xm::EXPRR))
+
+@
+
+<<not interesting implementation: Guess - guessExpRat - substitute>>=
+p(xm: Integer, i: Integer, va1: V, vA: V, basis: DIFFSPECN): FPOLYS ==
+ vA::POLYS::FPOLYS + (i-xm)*va1::POLYS::FPOLYS
+
+p2(xm: Integer, i: Symbol, a1v: F, Av: F, basis: DIFFSPECN): EXPRR ==
+ coerce(Av) + coerce(a1v)*(i::EXPRR - xm::EXPRR)
+
+@
+
+\subsubsection{Order and Degree}
+
+The following expressions for order and degree of the resultants [[res1]] and
+[[res2]] in [[guessExpRatAux]] were first guessed -- partially with the aid of
+[[guessRat]], and then proven to be correct.
+
+<<implementation: Guess - guessExpRat - Order and Degree>>=
+ord1(x: List Integer, i: Integer): Integer ==
+ n := #x - 3 - i
+ x.(n+1)*reduce(_+, [x.j for j in 1..n], 0) + _
+ 2*reduce(_+, [reduce(_+, [x.k*x.j for k in 1..j-1], 0) _
+ for j in 1..n], 0)
+
+ord2(x: List Integer, i: Integer): Integer ==
+ if zero? i then
+ n := #x - 3 - i
+ ord1(x, i) + reduce(_+, [x.j for j in 1..n], 0)*(x.(n+2)-x.(n+1))
+ else
+ ord1(x, i)
+
+deg1(x: List Integer, i: Integer): Integer ==
+ m := #x - 3
+ (x.(m+3)+x.(m+1)+x.(1+i))*reduce(_+, [x.j for j in 2+i..m], 0) + _
+ x.(m+3)*x.(m+1) + _
+ 2*reduce(_+, [reduce(_+, [x.k*x.j for k in 2+i..j-1], 0) _
+ for j in 2+i..m], 0)
+
+deg2(x: List Integer, i: Integer): Integer ==
+ m := #x - 3
+ deg1(x, i) + _
+ (x.(m+3) + reduce(_+, [x.j for j in 2+i..m], 0)) * _
+ (x.(m+2)-x.(m+1))
+
+@
+
+[[evalResultant]] evaluates the resultant of [[p1]] and [[p2]] at [[d-o+1]]
+points, so that we can recover it by interpolation.
+
+<<implementation: Guess - guessExpRat - evalResultant>>=
+evalResultant(p1: POLYS, p2: POLYS, o: Integer, d: Integer, va1: V, vA: V)_
+: List S ==
+ res: List S := []
+ d1 := degree(p1, va1)
+ d2 := degree(p2, va1)
+ lead: S
+ for k in 1..d-o+1 repeat
+ p1atk := univariate eval(p1, vA, k::S)
+ p2atk := univariate eval(p2, vA, k::S)
+
+@
+
+It may happen, that the leading coefficients of one or both of the polynomials
+changes, when we evaluate it at $k$. In this case, we need to correct this by
+multiplying with the corresponding power of the leading coefficient of the
+other polynomial.
+
+Consider the Sylvester matrix of the original polynomials. We want to evaluate
+it at $A=k$. If the first few leading coefficients of $p2$ vanish, the first
+few columns of the Sylvester matrix have triangular shape, with the leading
+coeffi...
[truncated message content] |
|
From: <gi...@ax...> - 2007-09-09 23:17:35
|
changelog | 2 +
src/algebra/Makefile.pamphlet | 3 +-
src/algebra/newton.spad.pamphlet | 46 --------------------------------------
3 files changed, 3 insertions(+), 48 deletions(-)
New commits:
commit b1815feeea1eab7906c7888561f120e686a72097
Author: Tim Daly <root@localhost.localdomain>
Date: Sat Sep 8 12:15:45 2007 -0400
20070909 tpd src/algebra/newton.spad included in fffg.spad
20070909 tpd src/algebra/Makefile remove newton.spad (duplicate)
|
|
From: <da...@us...> - 2007-09-09 23:17:26
|
Revision: 741
http://axiom.svn.sourceforge.net/axiom/?rev=741&view=rev
Author: daly
Date: 2007-09-09 16:17:28 -0700 (Sun, 09 Sep 2007)
Log Message:
-----------
20070909 tpd src/algebra/newton.spad included in fffg.spad
20070909 tpd src/algebra/Makefile remove newton.spad (duplicate)
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/algebra/Makefile.pamphlet
Removed Paths:
-------------
trunk/axiom/src/algebra/newton.spad.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-09 15:26:06 UTC (rev 740)
+++ trunk/axiom/changelog 2007-09-09 23:17:28 UTC (rev 741)
@@ -1,3 +1,5 @@
+20070909 tpd src/algebra/newton.spad included in fffg.spad
+20070909 tpd src/algebra/Makefile remove newton.spad (duplicate)
20070907 tpd src/algebra/acplot.spad fix PlaneAlgebraicCurvePlot.help NOISE
20070907 tpd src/algebra/Makefile make regression respect NOISE variable
20070907 tpd src/input/Makefile make regression respect NOISE variable
Modified: trunk/axiom/src/algebra/Makefile.pamphlet
===================================================================
--- trunk/axiom/src/algebra/Makefile.pamphlet 2007-09-09 15:26:06 UTC (rev 740)
+++ trunk/axiom/src/algebra/Makefile.pamphlet 2007-09-09 23:17:28 UTC (rev 741)
@@ -1265,7 +1265,7 @@
${OUTSRC}/mts.spad ${OUTSRC}/multfact.spad ${OUTSRC}/multpoly.spad \
${OUTSRC}/multsqfr.spad \
${OUTSRC}/naalgc.spad ${OUTSRC}/naalg.spad \
- ${OUTSRC}/newdata.spad ${OUTSRC}/newpoint.spad ${OUTSRC}/newton.spad \
+ ${OUTSRC}/newdata.spad ${OUTSRC}/newpoint.spad \
${OUTSRC}/newpoly.spad ${OUTSRC}/nlinsol.spad ${OUTSRC}/nlode.spad \
${OUTSRC}/npcoef.spad \
${OUTSRC}/nregset.spad \
@@ -1427,7 +1427,6 @@
${DOC}/multsqfr.spad.dvi \
${DOC}/naalgc.spad.dvi ${DOC}/naalg.spad.dvi ${DOC}/ndftip.as.dvi \
${DOC}/nepip.as.dvi ${DOC}/newdata.spad.dvi ${DOC}/newpoint.spad.dvi \
- ${DOC}/newton.spad.dvi \
${DOC}/newpoly.spad.dvi ${DOC}/nlinsol.spad.dvi ${DOC}/nlode.spad.dvi \
${DOC}/noptip.as.dvi ${DOC}/npcoef.spad.dvi ${DOC}/nqip.as.dvi \
${DOC}/nrc.as.dvi ${DOC}/nregset.spad.dvi ${DOC}/nsfip.as.dvi \
Deleted: trunk/axiom/src/algebra/newton.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/newton.spad.pamphlet 2007-09-09 15:26:06 UTC (rev 740)
+++ trunk/axiom/src/algebra/newton.spad.pamphlet 2007-09-09 23:17:28 UTC (rev 741)
@@ -1,46 +0,0 @@
-\documentclass{article}
-\usepackage{amsthm,amsmath,axiom}
-\begin{document}
-\title{newton.spad}
-\author{Martin Rubey}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\tableofcontents
-\section{package NEWTON NewtonInterpolation}
-<<package NEWTON NewtonInterpolation>>=
-)abbrev package NEWTON NewtonInterpolation
-++ Description:
-++ This package exports Newton interpolation for the special case where the
-++ result is known to be in the original integral domain
-NewtonInterpolation F: Exports == Implementation where
- F: IntegralDomain
- Exports == with
- newton: List F -> SparseUnivariatePolynomial F
-
- Implementation == add
-
- differences(yl: List F): List F ==
- [y2-y1 for y1 in yl for y2 in rest yl]
-
- z:SparseUnivariatePolynomial(F) := monomial(1,1)
-
--- we assume x=[1,2,3,...,n]
- newtonAux(k: F, fact: F, yl: List F): SparseUnivariatePolynomial(F) ==
- if empty? rest yl
- then ((yl.1) exquo fact)::F::SparseUnivariatePolynomial(F)
- else ((yl.1) exquo fact)::F::SparseUnivariatePolynomial(F)
- + (z-k::SparseUnivariatePolynomial(F)) _
- * newtonAux(k+1$F, fact*k, differences yl)
-
-
- newton yl == newtonAux(1$F, 1$F, yl)
-@
-<<*>>=
-<<package NEWTON NewtonInterpolation>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <gi...@ax...> - 2007-09-09 15:26:25
|
changelog | 3 ++
src/algebra/Makefile.pamphlet | 40 ++++++++++++++++++++++++++++++-------
src/algebra/acplot.spad.pamphlet | 6 +++++
src/input/Makefile.pamphlet | 33 +++++++++++++++++++++---------
4 files changed, 64 insertions(+), 18 deletions(-)
New commits:
commit 91a2e01df20b57ee32f05aa4fe054bd844e6e9a9
Author: Tim Daly <root@localhost.localdomain>
Date: Sat Sep 8 07:45:59 2007 -0400
20070907 tpd src/algebra/acplot.spad fix PlaneAlgebraicCurvePlot.help NOISE
20070907 tpd src/algebra/Makefile make regression respect NOISE variable
20070907 tpd src/input/Makefile make regression respect NOISE variable
|
|
From: <da...@us...> - 2007-09-09 15:26:06
|
Revision: 740
http://axiom.svn.sourceforge.net/axiom/?rev=740&view=rev
Author: daly
Date: 2007-09-09 08:26:06 -0700 (Sun, 09 Sep 2007)
Log Message:
-----------
20070907 tpd src/algebra/acplot.spad fix PlaneAlgebraicCurvePlot.help NOISE
20070907 tpd src/algebra/Makefile make regression respect NOISE variable
20070907 tpd src/input/Makefile make regression respect NOISE variable
Modified Paths:
--------------
trunk/axiom/changelog
trunk/axiom/src/algebra/Makefile.pamphlet
trunk/axiom/src/algebra/acplot.spad.pamphlet
trunk/axiom/src/input/Makefile.pamphlet
Modified: trunk/axiom/changelog
===================================================================
--- trunk/axiom/changelog 2007-09-08 04:44:16 UTC (rev 739)
+++ trunk/axiom/changelog 2007-09-09 15:26:06 UTC (rev 740)
@@ -1,3 +1,6 @@
+20070907 tpd src/algebra/acplot.spad fix PlaneAlgebraicCurvePlot.help NOISE
+20070907 tpd src/algebra/Makefile make regression respect NOISE variable
+20070907 tpd src/input/Makefile make regression respect NOISE variable
20070906 tpd Makefile int/sman/axiom command to target command for install
20070906 tpd src/sman/Makefile copy axiom command to int
20070906 tpd merge help files branch
Modified: trunk/axiom/src/algebra/Makefile.pamphlet
===================================================================
--- trunk/axiom/src/algebra/Makefile.pamphlet 2007-09-08 04:44:16 UTC (rev 739)
+++ trunk/axiom/src/algebra/Makefile.pamphlet 2007-09-09 15:26:06 UTC (rev 740)
@@ -2057,10 +2057,23 @@
\begin{verbatim}
%.regress: %.input
@ echo algebra regression testing $*
- @ rm -f $*.output
- @ echo ')read $*.input' | ${TESTSYS}
- @ echo ')lisp (regress "$*.output")' | ${TESTSYS} \
- | egrep -v '(Timestamp|Version)' | tee $*.regress
+ @ (cd ${MID} ; \
+ rm -f $*.output ; \
+ if [ -z "${NOISE}" ] ; then \
+ echo ')read $*.input' | ${TESTSYS} ; \
+ else \
+ echo ')read $*.input' | ${TESTSYS} >${TMP}/trace ; \
+ fi ; \
+ rm $*.input ; \
+ if [ -z "${NOISE}" ] ; then \
+ echo ')lisp (regress "$*.output")' | ${TESTSYS} \
+ | egrep -v '(Timestamp|Version)' | tee $*.regress ; \
+ else \
+ echo ')lisp (regress "$*.output")' | ${TESTSYS} \
+ | egrep -v '(Timestamp|Version)' > $*.regress ; \
+ fi ; \
+ fgrep "regression result" $*.regress )
+
\end{verbatim}
The input Makefile extracts {\bf algebra.regress} and then calls
make to process this file.
@@ -2191,10 +2204,21 @@
%.regress: %.input
@ echo algebra regression testing $*
- @ rm -f $*.output
- @ echo ')read $*.input' | ${TESTSYS}
- @ echo ')lisp (regress "$*.output")' | ${TESTSYS} \
- | egrep -v '(Timestamp|Version)' | tee $*.regress
+ @ rm -f $*.output
+ @ if [ -z "${NOISE}" ] ; then \
+ echo ')read $*.input' | ${TESTSYS} ; \
+ else \
+ echo ')read $*.input' | ${TESTSYS} >${TMP}/trace ; \
+ fi
+ @ rm $*.input
+ @ if [ -z "${NOISE}" ] ; then \
+ echo ')lisp (regress "$*.output")' | ${TESTSYS} \
+ | egrep -v '(Timestamp|Version)' | tee $*.regress ; \
+ else \
+ echo ')lisp (regress "$*.output")' | ${TESTSYS} \
+ | egrep -v '(Timestamp|Version)' > $*.regress ; \
+ fi
+ @ fgrep "regression result" $*.regress
all: ${REGRESS}
@echo algebra test cases complete.
Modified: trunk/axiom/src/algebra/acplot.spad.pamphlet
===================================================================
--- trunk/axiom/src/algebra/acplot.spad.pamphlet 2007-09-08 04:44:16 UTC (rev 739)
+++ trunk/axiom/src/algebra/acplot.spad.pamphlet 2007-09-09 15:26:06 UTC (rev 740)
@@ -218,6 +218,10 @@
\section{domain ACPLOT PlaneAlgebraicCurvePlot}
<<PlaneAlgebraicCurvePlot.input>>=
-- acplot.spad.pamphlet PlaneAlgebraicCurvePlot.input
+)spool PlaneAlgebraicCurvePlot.output
+)set message test on
+)set message auto off
+)clear all
--S 1 of 1
makeSketch(x+y,x,y,-1/2..1/2,-1/2..1/2)$ACPLOT
--R (1) ACPLOT
@@ -228,6 +232,8 @@
--R [-0.5,0.5]
--R Type: PlaneAlgebraicCurvePlot
--E 1
+)spool
+)lisp (bye)
@
<<PlaneAlgebraicCurvePlot.help>>=
====================================================================
Modified: trunk/axiom/src/input/Makefile.pamphlet
===================================================================
--- trunk/axiom/src/input/Makefile.pamphlet 2007-09-08 04:44:16 UTC (rev 739)
+++ trunk/axiom/src/input/Makefile.pamphlet 2007-09-09 15:26:06 UTC (rev 740)
@@ -383,8 +383,13 @@
@ echo ')set message auto off' >> tmp.input
@ echo ')read $*' >> tmp.input
@ echo ')lisp (bye)' >> tmp.input
- @ echo 'systemCommand "read tmp.input"' | ${TESTSYS} \
- | egrep -v '(Timestamp|Version)' | tee $*.output
+ @ if [ -z "${NOISE}" ] ; then \
+ echo 'systemCommand "read tmp.input"' | ${TESTSYS} \
+ | egrep -v '(Timestamp|Version)' | tee $*.output ; \
+ else \
+ echo 'systemCommand "read tmp.input"' | ${TESTSYS} \
+ | egrep -v '(Timestamp|Version)' > $*.output ; \
+ fi
@ rm tmp.input
@
@@ -402,15 +407,23 @@
<<regression tests>>=
%.regress: %.input
@ echo regression testing $*
- @ rm -f $*.output
- @ cp ${IN}/$*.input.pamphlet ${MID}
- @ echo ')read $*.input' | ${TESTSYS}
- @ rm ${MID}/$*.input.pamphlet
- @ rm $*.input
- @ echo ')lisp (regress "$*.output")' | ${TESTSYS} \
- | egrep -v '(Timestamp|Version)' | tee $*.regress
+ @ (cd ${MID} ; \
+ rm -f $*.output ; \
+ if [ -z "${NOISE}" ] ; then \
+ echo ')read $*.input' | ${TESTSYS} ; \
+ else \
+ echo ')read $*.input' | ${TESTSYS} >${TMP}/trace ; \
+ fi ; \
+ rm $*.input ; \
+ if [ -z "${NOISE}" ] ; then \
+ echo ')lisp (regress "$*.output")' | ${TESTSYS} \
+ | egrep -v '(Timestamp|Version)' | tee $*.regress ; \
+ else \
+ echo ')lisp (regress "$*.output")' | ${TESTSYS} \
+ | egrep -v '(Timestamp|Version)' > $*.regress ; \
+ fi ; \
+ fgrep "regression result" $*.regress )
-
@
\section{The Makefile}
This Makefile will be notangled into the [[input]] subdirectory.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|