diff --git a/books/bookvol10.2.pamphlet b/books/bookvol10.2.pamphlet
index 300a1b4..606738e 100644
--- a/books/bookvol10.2.pamphlet
+++ b/books/bookvol10.2.pamphlet
@@ -13392,7 +13392,7 @@ inverse matrix [[j**i for i in 0..4] for j in 1..5]
 --R MatrixCategory(R: Ring,Row: FiniteLinearAggregate t#1,Col: FiniteLinearAggregate t#1)  is a category constructor
 --R Abbreviation for MatrixCategory is MATCAT 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.2.spad.pamphlet to see algebra source code for MATCAT 
+--R Issue )edit bookvol10.2.pamphlet to see algebra source code for MATCAT 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?*? : (Row,%) -> Row                  ?*? : (%,Col) -> Col
diff --git a/books/bookvol10.3.pamphlet b/books/bookvol10.3.pamphlet
index d509cbe..c1e57d5 100644
--- a/books/bookvol10.3.pamphlet
+++ b/books/bookvol10.3.pamphlet
@@ -1638,7 +1638,7 @@ b := generator()$Sae
 --R Any  is a domain constructor
 --R Abbreviation for Any is ANY 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for ANY 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for ANY 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?=? : (%,%) -> Boolean                any : (SExpression,None) -> %
@@ -2198,7 +2198,7 @@ latex a
 --R ArrayStack S: SetCategory  is a domain constructor
 --R Abbreviation for ArrayStack is ASTACK 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for ASTACK 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for ASTACK 
 --R
 --R------------------------------- Operations --------------------------------
 --R arrayStack : List S -> %              bag : List S -> %
@@ -19344,7 +19344,7 @@ latex a
 --R Dequeue S: SetCategory  is a domain constructor
 --R Abbreviation for Dequeue is DEQUEUE 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for DEQUEUE 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for DEQUEUE 
 --R
 --R------------------------------- Operations --------------------------------
 --R back : % -> S                         bag : List S -> %
@@ -42296,7 +42296,7 @@ coerce a
 --R Heap S: OrderedSet  is a domain constructor
 --R Abbreviation for Heap is HEAP 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for HEAP 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for HEAP 
 --R
 --R------------------------------- Operations --------------------------------
 --R bag : List S -> %                     copy : % -> %
@@ -60279,7 +60279,7 @@ have to be switched by swapping names.
 --R MathMLFormat  is a domain constructor
 --R Abbreviation for MathMLFormat is MMLFORM 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for MMLFORM 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for MMLFORM 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?=? : (%,%) -> Boolean                coerce : OutputForm -> String
@@ -65447,7 +65447,7 @@ sample()$NOTTING(PF(1783))
 --R NottinghamGroup F: FiniteFieldCategory  is a domain constructor
 --R Abbreviation for NottinghamGroup is NOTTING 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for NOTTING 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for NOTTING 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?*? : (%,%) -> %                      ?**? : (%,Integer) -> %
@@ -75385,7 +75385,7 @@ listBranches(refined)
 --R PlaneAlgebraicCurvePlot  is a domain constructor
 --R Abbreviation for PlaneAlgebraicCurvePlot is ACPLOT 
 --R This constructor is not exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for ACPLOT 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for ACPLOT 
 --R
 --R------------------------------- Operations --------------------------------
 --R coerce : % -> OutputForm              refine : (%,DoubleFloat) -> %
@@ -81566,7 +81566,7 @@ latex a
 --R Queue S: SetCategory  is a domain constructor
 --R Abbreviation for Queue is QUEUE 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for QUEUE 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for QUEUE 
 --R
 --R------------------------------- Operations --------------------------------
 --R back : % -> S                         bag : List S -> %
@@ -97527,7 +97527,7 @@ latex a
 --R Stack S: SetCategory  is a domain constructor
 --R Abbreviation for Stack is STACK 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for STACK 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for STACK 
 --R
 --R------------------------------- Operations --------------------------------
 --R bag : List S -> %                     copy : % -> %
@@ -102372,7 +102372,7 @@ in addition we need to add a line defining [[PI2]] in [[formatPlex]]:
 --R TexFormat  is a domain constructor
 --R Abbreviation for TexFormat is TEX 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for TEX 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for TEX 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?=? : (%,%) -> Boolean                coerce : OutputForm -> %
@@ -110460,7 +110460,7 @@ x^2*a
 --R UnivariateSkewPolynomial(x: Symbol,R: Ring,sigma: Automorphism R,delta: (R -> R))  is a domain constructor
 --R Abbreviation for UnivariateSkewPolynomial is OREUP 
 --R This constructor is not exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for OREUP 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for OREUP 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?*? : (R,%) -> %                      ?*? : (%,R) -> %
diff --git a/books/bookvol10.4.pamphlet b/books/bookvol10.4.pamphlet
index c748eda..6703c8d 100644
--- a/books/bookvol10.4.pamphlet
+++ b/books/bookvol10.4.pamphlet
@@ -3699,22 +3699,22 @@ credits()
 --RRichard Anderson       George Andrews         S.J. Atkins
 --RHenry Baker            Stephen Balzac         Yurij Baransky
 --RDavid R. Barton        Gerald Baumgartner     Gilbert Baumslag
---RMichael Becker         Jay Belanger           David Bindel
---RFred Blair             Vladimir Bondarenko    Mark Botch
---RAlexandre Bouyer       Peter A. Broadbery     Martin Brock
---RManuel Bronstein       Stephen Buchwald       Florian Bundschuh
---RLuanne Burns           William Burge
+--RMichael Becker         Nelson H. F. Beebe     Jay Belanger
+--RDavid Bindel           Fred Blair             Vladimir Bondarenko
+--RMark Botch             Alexandre Bouyer       Peter A. Broadbery
+--RMartin Brock           Manuel Bronstein       Stephen Buchwald
+--RFlorian Bundschuh      Luanne Burns           William Burge
 --RQuentin Carpent        Robert Caviness        Bruce Char
 --ROndrej Certik          Cheekai Chin           David V. Chudnovsky
 --RGregory V. Chudnovsky  Josh Cohen             Christophe Conil
 --RDon Coppersmith        George Corliss         Robert Corless
 --RGary Cornell           Meino Cramer           Claire Di Crescenzo
 --RDavid Cyganski
---RTimothy Daly Sr.       Timothy Daly Jr.       James H. Davenport
---RDidier Deshommes       Michael Dewar          Jean Della Dora
---RGabriel Dos Reis       Claire DiCrescendo     Sam Dooley
---RLionel Ducos           Lee Duhem              Martin Dunstan
---RBrian Dupee            Dominique Duval
+--RNathaniel Daly         Timothy Daly Sr.       Timothy Daly Jr.
+--RJames H. Davenport     Didier Deshommes       Michael Dewar
+--RJean Della Dora        Gabriel Dos Reis       Claire DiCrescendo
+--RSam Dooley             Lionel Ducos           Lee Duhem
+--RMartin Dunstan         Brian Dupee            Dominique Duval
 --RRobert Edwards         Heow Eide-Goodman      Lars Erickson
 --RRichard Fateman        Bertfried Fauser       Stuart Feldman
 --RJohn Fletcher          Brian Ford             Albrecht Fortenbacher
@@ -3765,8 +3765,8 @@ credits()
 --RJonathan Steinbach     Fabio Stumbo           Christine Sundaresan
 --RRobert Sutor           Moss E. Sweedler       Eugene Surowitz
 --RMax Tegmark            James Thatcher         Balbir Thomas
---RMike Thomas            Dylan Thurston         Barry Trager
---RThemos T. Tsikas
+--RMike Thomas            Dylan Thurston         Steve Toleque
+--RBarry Trager           Themos T. Tsikas
 --RGregory Vanuxem
 --RBernhard Wall          Stephen Watt           Jaap Weel
 --RJuergen Weiss          M. Weller              Mark Wegman
@@ -3795,7 +3795,7 @@ summary()
 --R ApplicationProgramInterface  is a package constructor
 --R Abbreviation for ApplicationProgramInterface is API 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.4.spad.pamphlet to see algebra source code for API 
+--R Issue )edit bookvol10.4.pamphlet to see algebra source code for API 
 --R
 --R------------------------------- Operations --------------------------------
 --R credits : () -> Void                  getDomains : Symbol -> Set Symbol
@@ -129365,7 +129365,7 @@ map(f,q)
 --R QuaternionCategoryFunctions2(QR: QuaternionCategory R,R: CommutativeRing,QS: QuaternionCategory S,S: CommutativeRing)  is a package constructor
 --R Abbreviation for QuaternionCategoryFunctions2 is QUATCT2 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.4.spad.pamphlet to see algebra source code for QUATCT2 
+--R Issue )edit bookvol10.4.pamphlet to see algebra source code for QUATCT2 
 --R
 --R------------------------------- Operations --------------------------------
 --R map : ((R -> S),QR) -> QS            
@@ -145829,7 +145829,7 @@ definingPolynomial %phi1
 --R TransSolvePackage R: Join(OrderedSet,EuclideanDomain,RetractableTo Integer,LinearlyExplicitRingOver Integer,CharacteristicZero)  is a package constructor
 --R Abbreviation for TransSolvePackage is SOLVETRA 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.4.spad.pamphlet to see algebra source code for SOLVETRA 
+--R Issue )edit bookvol10.4.pamphlet to see algebra source code for SOLVETRA 
 --R
 --R------------------------------- Operations --------------------------------
 --R solve : Expression R -> List Equation Expression R
diff --git a/books/bookvol5.pamphlet b/books/bookvol5.pamphlet
index 2c97918..1a58ea4 100644
--- a/books/bookvol5.pamphlet
+++ b/books/bookvol5.pamphlet
@@ -36681,7 +36681,7 @@ constructor abbreviation to pamphlet file name.
   (setq |$constructorList| nil) ;; affects buildLibdb
   (setq *sourcefiles* (build-name-to-pamphlet-hash 
     (concatenate 'string (|getEnv| "AXIOM") 
-      "/../../src/algebra/*.spad.pamphlet")))
+      "/../../src/algebra/*.pamphlet")))
   (|buildLibdb|)
   (|dbSplitLibdb|)
 ; (|dbAugmentConstructorDataTable|)
diff --git a/changelog b/changelog
index d72155e..1f33580 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,31 @@
+20100404 tpd src/axiom-website/patches.html 20100404.01.tpd.patch
+20100404 tpd src/algebra/nsfip.as removed
+20100404 tpd src/algebra/nrc.as removed
+20100404 tpd src/algebra/nqip.as removed
+20100404 tpd src/algebra/noptip.as removed
+20100404 tpd src/algebra/nepip.as removed
+20100404 tpd src/algebra/ndftip.as removed
+20100404 tpd src/algebra/mlift.spad.jhd removed
+20100404 tpd src/algebra/iviews.as removed
+20100404 tpd src/algebra/invutils.as removed
+20100404 tpd src/algebra/invtypes.as removed
+20100404 tpd src/algebra/invrender.as removed
+20100404 tpd src/algebra/invnode.as removed
+20100404 tpd src/algebra/interval.as removed
+20100404 tpd src/algebra/herm.as removed
+20100404 tpd src/algebra/ffrac.as removed
+20100404 tpd src/algebra/axtimer.as removed
+20100404 tpd src/algebra/Makefile remove unused .as.pamphlet files
+20100404 tpd src/input/unittest1.input fix regression tests
+20100404 tpd src/input/gstbl.input fix regression tests
+20100404 tpd src/input/grpthry.input fix regression tests
+20100404 tpd src/input/chtheorem.input fix regression tests
+20100404 tpd src/input/biquat.input fix regression tests
+20100404 tpd books/bookvol10.4 fix regression tests
+20100404 tpd books/bookvol10.3 fix regression tests
+20100404 tpd books/bookvol10.2 fix regression tests
+20100404 tpd books/bookvol5 change .spad.pamphlet to just .pamphlet
+20100404 tpd src/Makefile change .spad.pamphlet to .pamphlet
 20100403 tpd src/axiom-website/patches.html 20100403.02.tpd.patch
 20100403 tpd Makefile clean books from src/algebra
 20100403 tpd src/axiom-website/patches.html 20100403.01.tpd.patch
diff --git a/src/Makefile.pamphlet b/src/Makefile.pamphlet
index eb05c33..8e3b75a 100644
--- a/src/Makefile.pamphlet
+++ b/src/Makefile.pamphlet
@@ -499,10 +499,10 @@ ${SRC}/algebra/Makefile: ${SRC}/algebra/Makefile.pamphlet
                 ${SRC}/algebra/Makefile.pamphlet
 	@( cd algebra ; ${DOCUMENT} ${NOISE} Makefile ; \
         cp Makefile.dvi ${MNT}/${SYS}/doc/src/algebra.Makefile.dvi ; \
-	cp ${SPD}/books/bookvol10.2.pamphlet bookvol10.2.spad.pamphlet ; \
-	cp ${SPD}/books/bookvol10.3.pamphlet bookvol10.3.spad.pamphlet ; \
-	cp ${SPD}/books/bookvol10.4.pamphlet bookvol10.4.spad.pamphlet ; \
-	cp ${SPD}/books/bookvol10.5.pamphlet bookvol10.5.spad.pamphlet ; \
+	cp ${SPD}/books/bookvol10.2.pamphlet . ; \
+	cp ${SPD}/books/bookvol10.3.pamphlet . ; \
+	cp ${SPD}/books/bookvol10.4.pamphlet . ; \
+	cp ${SPD}/books/bookvol10.5.pamphlet . ; \
 	echo 30a extracting findAlgebraFiles from \
                  ${SRC}/algebra/Makefile.pamphlet ; \
 	${TANGLE} -t8 -RfindAlgebraFiles Makefile.pamphlet \
diff --git a/src/algebra/Makefile.pamphlet b/src/algebra/Makefile.pamphlet
index 3f2928e..385c70d 100644
--- a/src/algebra/Makefile.pamphlet
+++ b/src/algebra/Makefile.pamphlet
@@ -15798,38 +15798,6 @@ it to interpret special characters, in our case newlines.
 SHELL=bash
 
 @
-\subsection{The ALDORFILES list}
-<<environment>>=
-
-ALDORFILES= \
- ${MID}/axtimer.as \
- ${MID}/ffrac.as \
- ${MID}/herm.as \
- ${MID}/interval.as \
- ${MID}/invnode.as ${MID}/invrender.as \
- ${MID}/invtypes.as ${MID}/invutils.as \
- ${MID}/iviews.as \
- ${MID}/ndftip.as \
- ${MID}/nepip.as \
- ${MID}/noptip.as ${MID}/nqip.as \
- ${MID}/nrc.as ${MID}/nsfip.as 
-
-@
-\subsection{The DOCFILES list}
-<<environment>>=
-
-DOCFILES= \
- ${DOC}/herm.as.dvi \
- ${DOC}/interval.as.dvi \
- ${DOC}/invnode.as.dvi ${DOC}/invrender.as.dvi \
- ${DOC}/invtypes.as.dvi ${DOC}/invutils.as.dvi  \
- ${DOC}/iviews.as.dvi \
- ${DOC}/ndftip.as.dvi \
- ${DOC}/nepip.as.dvi  \
- ${DOC}/noptip.as.dvi ${DOC}/nqip.as.dvi \
- ${DOC}/nrc.as.dvi  ${DOC}/nsfip.as.dvi 
-
-@
 \section{The Makefile Stanzas}
 A [[spad]] pamphlet can contain many Axiom [[categories]], [[domains]], and
 [[packages]]. 
@@ -15986,10 +15954,10 @@ ${MID}/%.o: ${MID}/%.lsp
 @
 <<genericSPADfiles>>=
 
-${OUTSRC}/%.spad: ${IN}/%.spad.pamphlet
-	@ echo tangling $*.spad.pamphlet to $*.spad
+${OUTSRC}/%.spad: ${IN}/%.pamphlet
+	@ echo tangling $*.pamphlet to $*.spad
 	@(cd ${OUTSRC} ; \
-	 ${TANGLE} ${IN}/$*.spad.pamphlet >$*.spad )
+	 ${TANGLE} ${IN}/$*.pamphlet >$*.spad )
 
 @
 <<genericDOCfiles>>=
@@ -16061,19 +16029,19 @@ means that the pattern is an extended regular expression) because
 extended regular expressions allows the use of alternatives
 written as (domain|package|category). Thus the command
 \begin{verbatim}
- grep -E '@<<(domain|package|category) .*>>=' *.spad.pamphlet 
+ grep -E '@<<(domain|package|category) .*>>=' *.pamphlet 
 \end{verbatim}
 will scan the algebra files looking for special chunknames.
 Axiom's chunk names are written in a stylized form so that each
 algebra chunk name begins with one of those three symbols. Thus 
-in zerodim.spad.pamphlet the LexTriangularPackage chunkname is:
+in bookvol10.3.pamphlet the LexTriangularPackage chunkname is:
 \begin{verbatim}
 @<<package LEXTRIPK LexTriangularPackage>>
 \end{verbatim}
 so this grep will generate an output line, prefixed by the filename
 that looks like:
 \begin{verbatim}
-zerodim.spad.pamphlet:@<<package LEXTRIPK LexTriangularPackage>>=
+bookvol10.3.pamphlet:@<<package LEXTRIPK LexTriangularPackage>>=
 \end{verbatim}
 There can be many lines of output per pamphlet file, one for
 each domain, package and category cod chunk contained in the file.
@@ -16084,7 +16052,7 @@ Step 2 is an [[awk]] command line.
 
 <<findSpadFiles>>=
 
-grep -E '@<<(domain|package|category) .*>>=' *.spad.pamphlet | sort | uniq | \
+grep -E '@<<(domain|package|category) .*>>=' *.pamphlet | sort | uniq | \
 awk -F: '{
   chunk=substr($2,3,length($2)-5);
   split(chunk,part," ");
@@ -16129,15 +16097,15 @@ the case above this would resolve to [[\${MID}/LEXTRIPK.spad]].
 
 For the line given above it outputs the following:
 \begin{verbatim}
-${MID}/LEXTRIPK.spad: ${IN}/zerodim.spad.pamphlet
-	${TANGLE} -R"package LEXTRIPK LexTriangularPackage" ${IN}/zerodim.spad.pamphlet > ${MID}/LEXTRIPK.spad
+${MID}/LEXTRIPK.spad: ${IN}/bookvol10.3.pamphlet
+	${TANGLE} -R"package LEXTRIPK LexTriangularPackage" ${IN}/bookvol10.3.pamphlet > ${MID}/LEXTRIPK.spad
 \end{verbatim}
 
 \subsection{Find the algebra bootstrap code}
 Step 3 works like step 1 above except that we are looking for
 chunk names that have the "BOOTSTRAP" string. The output will look like:
 \begin{verbatim}
-vector.spad.pamphlet:@<<VECTOR.lsp BOOTSTRAP>>=
+bookvol10.3.pamphlet:@<<VECTOR.lsp BOOTSTRAP>>=
 \end{verbatim}
 This output, which can consist of many lines per input file is piped
 into [[awk]].
@@ -16147,7 +16115,7 @@ For each of the above output lines we run an [[grep]] command:
 
 <<findBootstrapFiles>>=
 
-grep -E '@<<.*BOOTSTRAP>>=' *.spad.pamphlet | sort | uniq | \
+grep -E '@<<.*BOOTSTRAP>>=' *.pamphlet | sort | uniq | \
 awk -F: '{
   chunk=substr($2,3,length($2)-5);
   split(chunk,part," ");
@@ -16171,8 +16139,8 @@ ${MID}/VECTOR.lsp
 \end{verbatim}
 Finally we output two lines:
 \begin{verbatim}
-${MID}/vector.spad.pamphlet: ${IN}/vector.spad.pamphlet
-	${TANGLE} -R"VECTOR.lsp BOOTSTRAP" ${IN}/vector.spad.pamphlet >${MID}/VECTOR.lsp
+${MID}/vector.spad.pamphlet: ${IN}/bookvol10.3.pamphlet
+	${TANGLE} -R"VECTOR.lsp BOOTSTRAP" ${IN}/bookvol10.3.pamphlet >${MID}/VECTOR.lsp
 \end{verbatim}
 
 The first line is the stanza head and creates a dependence between
@@ -18216,7 +18184,7 @@ ${HELP}/ZeroDimensionalSolvePackage.help: ${BOOKS}/bookvol10.4.pamphlet
 <<layer22>>
 <<order>>
 
-all: src ${OUT}/libdb.text ${DOCFILES} ${SPADBIN}/index.html gloss
+all: src ${OUT}/libdb.text ${SPADBIN}/index.html gloss
 	@ echo 4302 finished ${IN}
 
 parallelhelp: syntaxhelp ${SPADHELP}
@@ -18257,8 +18225,6 @@ everything: lib db cmd gloss
 src:	${ORDER}
 	@ echo 4304 Building nrlibS from spad sources
 
-document: ${DOCFILES}
-
 <<genericRules>>
 <<spadhelp>>
 <<ps (DOC from SRC)>>
diff --git a/src/algebra/axtimer.as.pamphlet b/src/algebra/axtimer.as.pamphlet
deleted file mode 100644
index 1f046bb..0000000
--- a/src/algebra/axtimer.as.pamphlet
+++ /dev/null
@@ -1,191 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra axtimer.as}
-\author{The Axiom Team}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
---------------------------------------------------------------------------------
---
--- BasicMath: timer.as --- Objects for tracking time
---
--------------------------------------------------------------------------------- 
---
-
--- ToDo: Write so that start!(x) ; start!(x) ; stop!(x); stop!(x) works.
-
-#include "axiom.as"
-
-Z ==> Integer;
-
-Duration ==> Record(hours: Z, mins: Z, seconds: Z, milliseconds: Z);
-
-+++ Timer
-+++ History: 22/5/94 Original version by MB.
-+++	      9/4/97 [Peter Broadbery] incorporated into Basicmath.
-+++	      7/8/97 [PAB] Hacked into axiom.
-+++ Timer is a type whose elements are stopwatch timers, which can be used
-+++ to time precisely various sections of code.
-+++ The precision can be up to 1 millisecond but depends on the operating system.
-+++ The times returned are CPU times used by the process that created the timer.
-
-Timer: BasicType with {
-	HMS:    Z -> Duration;
-		++ Returns (h, m, s, u) where n milliseconds is equal
-		++ to h hours, m minutes, s seconds and u milliseconds.
-	read:   %  -> Z;
-		++ Reads the timer t without stopping it.
-		++ Returns the total accumulated time in milliseconds by all
-		++ the start/stop cycles since t was created or last reset.
-		++ If t is running, the time since the last start is added in,
-		++ and t is not stopped or affected.
-	reset!: %  -> %;
-		++ Resets the timer t to 0 and stops it if it is running.
-		++ Returns the timer t after it is reset.
-	start!: %  -> Z;
-		++ Starts or restarts t, without resetting it to 0,
-		++ It has no effect on t if it is already running.
-		++ Returns 0 if t was already running, the absolute time at which
-		++ the start/restart was done otherwise.
-	stop!:  %  -> Z;
-		++ Stops t without resetting it to 0.
-		++ It has no effect on t if it is not running.
-		++ Returns the elapsed time in milliseconds since the last time t
-		++ was restarted, 0 if t was not running.
-	timer:  () -> %;
-		++ Creates a timer, set to 0 and stopped.
-		++ Returns the timer that has been created.
-
-	coerce: % -> OutputForm;	
-#if 0
-	gcTimer: () -> %;
-		++ Returns the system garbage collection timer.
-		++ Do not use for other purposes!
-#endif
-} == add {
-        -- time     = total accumulated time since created or reset
-        -- start    = absolute time of last start
-        -- running? = true if currently running, false if currently stopped
-	Rep ==> Record(time:Z, start:Z, running?:Boolean);
-	import {
-		BOOT_:_:GET_-INTERNAL_-RUN_-TIME: () -> Integer;
-	} from Foreign Lisp;
-	cpuTime(): Integer == BOOT_:_:GET_-INTERNAL_-RUN_-TIME();
-
-	import from Rep, Z, Boolean;
-
-	timer():% == per [0, 0, false];
-
-	read(t:%):Z == {
-		rec := rep t;
-		ans := rec.time;
-		if rec.running? then ans := ans + cpuTime() - rec.start;
-		ans
-	}
-
-	stop!(t:%):Z == {
-		local ans:Z;
-		rec := rep t;
-		if rec.running? then {
-			ans := cpuTime() - rec.start;
-			rec.time := rec.time + ans;
-			rec.running? := false
-		}
-		else ans := 0;
-		ans
-	}
-
-	start!(t:%):Z == {
-		local ans:Z;
-		rec := rep t;
-		if not(rec.running?) then {
-			rec.start := ans := cpuTime();
-			rec.running? := true;
-		}
-		else ans := 0;
-		ans
-	}
-
-	reset!(t:%):% == {
-		rec := rep t;
-		rec.time := rec.start := 0;
-		rec.running? := false;
-		t
-	}
-
-	HMS(m:Z): Duration == {
-		import from Record(quotient: Integer, remainder: Integer);
-		(h, m) := explode divide(m, 3600000);
-		(m, s) := explode divide(m, 60000);
-		(s, l) := explode divide(s, 1000);
-		[h, m, s, l]
-	}
-
-#if 0
-	import {
-		gcTimer: () -> Pointer;
-	} from Foreign C;
-	
-	gcTimer(): % == (gcTimer())@Pointer pretend %;
-#endif
-	coerce(x: %): OutputForm == {
-		import from List OutputForm;
-		assign(name: String, val: OutputForm): OutputForm == {
-			import from Symbol;
-			blankSeparate [outputForm(name::Symbol), outputForm("="::Symbol), val];
-		}
-		state: Symbol := coerce(if rep(x).running? then "on" else "off");
-		bracket [outputForm("Timer:"::Symbol), 
-			 commaSeparate [assign("state", outputForm state),
-				        assign("value", (read x)::OutputForm)]];
-	}
-
-	(a: %) = (b: %): Boolean == error "No equality for Timers";
-}
-
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/ffrac.as.pamphlet b/src/algebra/ffrac.as.pamphlet
deleted file mode 100644
index 8eccc8a..0000000
--- a/src/algebra/ffrac.as.pamphlet
+++ /dev/null
@@ -1,204 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra ffrac.as}
-\author{Michael Richardson}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\begin{verbatim}
-
--- FormalFraction
-
--- N.B. ndftip.as inlines this, must be recompiled if this is.
-
--- To test:
--- sed '1,/^#if NeverAssertThis/d;/#endif/d' < ffrac.as > ffrac.input    
--- axiom
--- )set nag host <some machine running nagd>
--- )r ffrac.input
-
-\end{verbatim}
-\section{FormalFraction}
-<<FormalFraction>>=
-
-#include "axiom.as"
-
-FFRAC ==> FormalFraction ;
- 
-OF    ==> OutputForm ;
-SC    ==> SetCategory ;
-FRAC  ==> Fraction ;
-ID    ==> IntegralDomain ;
-
-+++ Author: M.G. Richardson
-+++ Date Created: 1996 Jan. 23
-+++ Date Last Updated:
-+++ Basic Functions:
-+++ Related Constructors: Fraction
-+++ Also See:
-+++ AMS Classifications:
-+++ Keywords:
-+++ References:
-+++ Description:
-+++ This type represents formal fractions - that is, pairs displayed as
-+++ fractions with no simplification.
-+++
-+++ If the elements of the pair have a type X which is an integral
-+++ domain, a FFRAC X can be coerced to a FRAC X, provided that this
-+++ is a valid type.  A FRAC X can always be coerced to a FFRAC X.
-+++ If the type of the elements is a Field, a FFRAC X can be coerced
-+++ to X.
-+++
-+++ Formal fractions are used to return results from numerical methods
-+++ which determine numerator and denominator separately, to enable
-+++ users to inspect these components and recognise, for example,
-+++ ratios of very small numbers as potentially indeterminate.
-
-FormalFraction(X : SC) : SC with {
-
--- Could generalise further to allow numerator and denominator to be of
--- different types, X and Y, both SCs.  "Left as an exercise."
-
-  / : (X,X) -> % ;
-++ / forms the formal quotient of two items.
-
-  numer : % -> X ;
-++ numer returns the numerator of a FormalFraction.
-
-  denom : % -> X ;
-++ denom returns the denominator of a FormalFraction.
-
-  if X has ID then {
-  
-    coerce : % -> FRAC(X pretend ID) ;
-++ coerce x converts a FormalFraction over an IntegralDomain to a
-++ Fraction over that IntegralDomain.
-
-    coerce : FRAC(X pretend ID) -> % ;
-++ coerce converts a Fraction to a FormalFraction.
-
-  }
-
-  if X has Field then coerce : % -> (X pretend Field) ;
-
-} == add {
-
-  import from Record(num : X, den : X) ;
-  
-  Rep == Record(num : X, den : X) ; -- representation
-
-  ((x : %) = (y : %)) : Boolean ==
-    ((rep(x).num = rep(y).num) and (rep(x).den = rep(y).den)) ;
-
-  ((n : X)/(d : X)) : % == per(record(n,d)) ;
-
-  coerce(r : %) : OF == (rep(r).num :: OF) / (rep(r).den :: OF) ;
-
-  numer(r : %) : X == rep(r).num ;
-
-  denom(r : %) : X == rep(r).den ;
-
-  if X has ID then {
- 
-    coerce(r : %) : FRAC(X pretend ID)
-      == ((rep(r).num)/(rep(r).den)) @ (FRAC(X pretend ID)) ;
-
-    coerce(x : FRAC(X pretend ID)) : % == x pretend % ; 
-
-  }
-
-  if X has Field then coerce(r : %) : (X pretend Field)
-    == ((rep(r).num)/(rep(r).den)) $ (X pretend Field) ;
-
-}
-
-#if NeverAssertThis
-
-)lib ffrac
-
-f1 : FormalFraction Integer
-f1 := 6/3
-
---       6
---       -
---       3
-
-f2 := (3.6/2.4)$FormalFraction Float
-
---       3.6
---       ---
---       2.4
-
-numer f1
-
---       6
-
-denom f2
-
---       2.4
-
-f1 :: FRAC INT
-
---       2
-
-% :: FormalFraction Integer      
-
---       2
---       -
---       1
-
-f2 :: Float
-
---       1.5
-
-output "End of tests"
-
-#endif
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
-<<FormalFraction>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/herm.as.pamphlet b/src/algebra/herm.as.pamphlet
deleted file mode 100644
index 81915ba..0000000
--- a/src/algebra/herm.as.pamphlet
+++ /dev/null
@@ -1,369 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra herm.as}
-\author{Michael Richardson}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\begin{verbatim}
--- N.B. ndftip.as inlines this, must be recompiled if this is.
-
--- To test:
--- sed -ne '1,/^#if NeverAssertThis/d;/#endif/d;p' < herm.as > herm.input    
--- axiom
--- )set nag host <some machine running nagd>
--- )r herm.input
-\end{verbatim}
-\section{PackedHermitianSequence}
-<<PackedHermitianSequence>>=
-#include "axiom.as"
-
-INT ==> Integer ;
-NNI ==> NonNegativeInteger ;
-PHS ==> PackedHermitianSequence ;
- 
-+++ Author: M.G. Richardson
-+++ Date Created: 1995 Nov. 24
-+++ Date Last Updated:
-+++ Basic Functions:
-+++ Related Constructors: Vector
-+++ Also See:
-+++ AMS Classifications:
-+++ Keywords:
-+++ References:
-+++ Description:
-+++ This type represents packed Hermitian sequences - that is, complex
-+++ sequences s, whose tails ("rest s", in Axiom terms) are conjugate to
-+++ themselves reversed - by real sequences, in the "standard" manner;
-+++ in this, the real parts of the elements of the first "half" of the
-+++ tail are stored there and the imaginary parts are stored in reverse
-+++ order in the second "half" of the tail.
-+++ (If the tail has an odd number of elements, its middle element is
-+++ real and is stored unchanged.  The "halves" mentioned above then
-+++ refer to the elements before and after the middle, respectively.)
-
-PackedHermitianSequence(R : CommutativeRing) : LinearAggregate(R) with{
-
-  pHS : List R -> % ;
-++ pHS(l) converts the list l to a packedHermitianSequence.
-
-  expand : % -> Vector Complex R ;
-++ expand(h) converts the packedHermitianSequence h to a Hermitian
-++ sequence (a complex vector).
-
-  packHS : Vector Complex R -> % ;
-++ packHS(v) checks that the complex vector v represents a Hermitian
-++ sequences and, if so, converts it to a packedHermitianSequence;
-++ otherwise an error message is printed.
-
-  conjHerm : % -> % ;
-++ conjHerm(h) returns the packedHermitianSequence which represents the
-++ Hermitian sequence conjugate to that represented by h.
-
- coerce : % -> OutputForm -- shouldn't need this, should be inherited
-			  -- from Vector.
-
-} == Vector(R) add {
-  Rep ==> Vector R;
-  import from Rep;
-  import from INT ;
-  import from R ;
-  import from Vector R ;
-  import from Complex R ;
-  import from Vector Complex R ;
-  import from ErrorFunctions ;
-  import from String ;
-  import from List String ;
-
-  local (..)(a:INT,b:INT):Generator INT == {
-    generate {
-      t := a ;
-      while (t <= b) repeat {
-        yield t ;
-        t := t + 1 ;
-        }
-      }
-    }
-
-  pHS(l : List R) : % == (vector l) pretend PHS R ;
-
-  expand(h : %) : Vector Complex R == {
- 
-    local len : NNI ;
-    local nvals, npairs, n1, n2 : INT ;
-    local fullh : Vector Complex R ;
-
-    {
-    len := # h ;
-    nvals := len pretend INT ; -- pretend since :: fails
-    npairs := (nvals - 1) quo 2 ;
-    fullh := new(len, 0) ;
-    (nvals = 0) => () ;
-    fullh.1 := complex(h.1,0) ;
-    (nvals = 1) => () ;
-    fullh.(npairs+2) := complex(h.(npairs+2),0) ; -- need if even length
-                                                  -- (not worth testing)
-    for j in 1 .. npairs repeat {
-      n1 := j + 1 ;
-      n2 := nvals - j + 1 ;
-      fullh.n1 := complex(h.n1, h.n2) ;
-      fullh.n2 := complex(h.n1, -(h.n2)) ;
-      }
-
-    }
-      
-    fullh
-
-  }
-
-  packHS(v : Vector Complex R) : % == {
-
-    local len : NNI ;
-    local nonhs : String ==  "The argument of packHS is not Hermitian" ;
-    local nvals, testprs, n1, n2 : INT ;
-    local hpacked : Vector R ;
-    local v1, v2 : Complex R ;
-    local r1, i1, r2, i2 : R ;
-
-    {
-    len := # v ;
-    nvals := len pretend INT ; -- pretend since :: fails
-    testprs := nvals quo 2 ;
-    hpacked := new(len, 0) ;
-    (nvals = 0) => () ;
-    if imag(v.1) ~= 0 
-    then error [nonhs, " - the first element must be real."]
-    else {
-      hpacked.1 := real(v.1) ;
-      (nvals = 1) => () ;
-      for j in 1 .. testprs repeat {
-        n1 := j + 1 ;
-        n2 := nvals - j + 1 ;
-        v1 := v.n1 ;
-        v2 := v.n2 ;
-        r1 := real v1 ;
-        i1 := imag v1 ;
-        r2 := real v2 ;
-        i2 := imag v2 ;
-        if r1 ~= r2 or i1 ~= -i2
-        then if n1 = n2
-	     then error [nonhs,
-			 " - element ",
-			 string(n1),
-			 " must be real to be self-conjugate."]
-             else error [nonhs,
-                         " - elements ",
-			 string(n1),
-			 " and ",
-			 string(n2),
-			 " are not conjugate."]
-        else {
-	  hpacked.n2 := i1 ; -- This order means that when the tail of v
-	  hpacked.n1 := r1 ; -- has odd length, the (real part) of its
-			     -- middle element ends up in that position.
-	  }
-	}
-      }
-    }
-
-    hpacked pretend %
-
-  }
-
-  local set!(x: %, i: INT, v: R): () == {
-	(rep x).i := v;
-  }
-  conjHerm(h : %) : %  == {
-
-    local len : NNI ;
-    local nvals, npairs : INT ;
-    local ch : % ;
-
-    ch := copy h ;
-    len := # h ;
-    (len < 3) => ch ; -- these Hermitian sequences are self-conjugate.
-    nvals := len pretend INT ; -- pretend since :: fails
-    npairs := (nvals - 1) quo 2 ;
-    for j in (nvals - npairs + 1) .. nvals repeat ch.j := - h.j ;
-    ch
-
-  }
-   
-  import from List OutputForm ;
-  
-  coerce(h : %) : OutputForm ==
-    bracket commaSeparate [
-       qelt(h, k) :: OutputForm for k in minIndex h .. maxIndex h]
-
-}
-
-#if NeverAssertThis
-
-)lib herm
-
-h0 := pHS([] :: List INT)
-
---       []
-
-h1 := pHS [1]
-
---       [1]
-
-h2 := pHS [1,2]
-
---       [1,2]
-
-h3 := pHS [1,2,3]
-
---       [1,2,3]
-
-h4 := pHS [1,2,3,4]
-
---       [1,2,3,4]
-
-h5 := pHS [1,2,3,4,5]
-
---       [1,2,3,4,5]
-
-
-f0 := expand h0
-
---       []
-
-f1 := expand h1
-
---       [1]
-
-f2 := expand h2
-
---       [1,2]
-
-f3 := expand h3
-
---       [1,2 + 3%i,2 - 3%i]
-
-f4 := expand h4
-
---       [1,2 + 4%i,3,2 - 4%i]
-
-f5 := expand h5
-
---       [1,2 + 5%i,3 + 4%i,3 - 4%i,2 - 5%i]
-      
-packHS f0
-
---       []
-
-packHS f1
-
---       [1]
-
-packHS f2
-
---       [1,2]
-
-packHS f3
-
---       [1,2,3]
-
-packHS f4
-
---       [1,2,3,4]
-
-packHS f5
-
---       [1,2,3,4,5]
-
-packHS vector[%i,3,3,3]
-
--- Error signalled from user code:
---    The argument of packHS is not Hermitian - the first element must
---    be real.
-
-packHS vector [1, 3, 5, 7]
-
--- Error signalled from user code:
---    The argument of packHS is not Hermitian - elements 2 and 4 are 
---    not conjugate.
-
-packHS [1, 3, %i, 3]
-
--- Error signalled from user code:
---    The argument of packHS is not Hermitian - element 3 must be real
---    to be self-conjugate.
-
-conjHerm h0
-
---       []
-
-conjHerm h1
-
---       [1]
-
-conjHerm h2
-
---       [1,2]
-
-conjHerm h3
-
---       [1,2,- 3]
-
-conjHerm h4
-
---       [1,2,3,- 4]
-
-conjHerm h5
-
---       [1,2,3,- 4,- 5]
-
-output "End of tests"
-
-#endif
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
-<<PackedHermitianSequence>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/interval.as.pamphlet b/src/algebra/interval.as.pamphlet
deleted file mode 100644
index 765d732..0000000
--- a/src/algebra/interval.as.pamphlet
+++ /dev/null
@@ -1,564 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra interval.as}
-\author{Mike Dewar}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{IntervalCategory}
-<<IntervalCategory>>=
-#include "axiom.as"
-
-+++ Author: Mike Dewar
-+++ Date Created: November 1996
-+++ Date Last Updated:
-+++ Basic Functions:
-+++ Related Constructors: 
-+++ Also See:
-+++ AMS Classifications:
-+++ Keywords:
-+++ References:
-+++ Description:
-+++ This category is an implementation of interval arithmetic and transcendental
-+++ functions over intervals.
-
-FUNCAT ==> Join(FloatingPointSystem,TranscendentalFunctionCategory);
-
-define IntervalCategory(R:FUNCAT): Category ==
- Join(GcdDomain, OrderedSet, TranscendentalFunctionCategory, RadicalCategory,
-      RetractableTo(Integer))
- with {
-  approximate;
-  interval : (R,R) -> %;
-    ++ interval(inf,sup) creates a new interval, either \axiom{[inf,sup]} if
-    ++ \axiom{inf <= sup} or \axiom{[sup,in]} otherwise.
-  qinterval : (R,R) -> %;
-    ++ qinterval(inf,sup) creates a new interval \axiom{[inf,sup]}, without
-    ++ checking the ordering on the elements.
-  interval : R -> %;
-    ++ interval(f) creates a new interval around f.
-  interval : Fraction Integer -> %;
-    ++ interval(f) creates a new interval around f.
-  inf : % -> R;
-    ++ inf(u) returns the infinum of \axiom{u}.
-  sup : % -> R;
-    ++ sup(u) returns the supremum of \axiom{u}.
-  width : % -> R;
-    ++ width(u) returns \axiom{sup(u) - inf(u)}.
-  positive? : % -> Boolean;
-    ++ positive?(u) returns \axiom{true} if every element of u is positive,
-    ++ \axiom{false} otherwise.
-  negative? : % -> Boolean;
-    ++ negative?(u) returns \axiom{true} if every element of u is negative,
-    ++ \axiom{false} otherwise.
-  contains? : (%,R) -> Boolean;
-    ++ contains?(i,f) returns true if \axiom{f} is contained within the interval
-    ++ \axiom{i}, false otherwise.
-}
-
-@
-\section{Interval}
-<<Interval>>=
-+++ Author: Mike Dewar
-+++ Date Created: November 1996
-+++ Date Last Updated:
-+++ Basic Functions:
-+++ Related Constructors: 
-+++ Also See:
-+++ AMS Classifications:
-+++ Keywords:
-+++ References:
-+++ Description:
-+++ This domain is an implementation of interval arithmetic and transcendental
-+++ functions over intervals.
-
-Interval(R:FUNCAT): IntervalCategory(R) == add {
-
-  import from Integer;
-  import from R;
-
-  Rep ==> Record(Inf:R, Sup:R);
-
-  import from Rep;
-
-  local roundDown(u:R):R == 
-    if zero?(u) then float(-1,-(bits() pretend Integer));
-                else float(mantissa(u) - 1,exponent(u));
-
-  local roundUp(u:R):R   == 
-    if zero?(u) then float(1, -(bits()) pretend Integer);
-                else float(mantissa(u) + 1,exponent(u));
-
-  -- Sometimes the float representation does not use all the bits (e.g. when
-  -- representing an integer in software using arbitrary-length Integers as
-  -- your mantissa it is convenient to keep them exact).  This function 
-  -- normalises things so that rounding etc. works as expected.  It is only
-  -- called when creating new intervals.
-  local normaliseFloat(u:R):R == 
-    if zero? u then u else {
-    m : Integer := mantissa u;
-    b : Integer := bits() pretend Integer;
-    l : Integer := length(m);
-    if (l < b) then {
-      BASE : Integer := base()$R pretend Integer;
-      float(m*BASE**((b-l) pretend PositiveInteger),exponent(u)-b+l);
-    }
-    else
-      u;
-  }
-
-  interval(i:R,s:R):% == {
-    i > s =>  per [roundDown normaliseFloat s,roundUp normaliseFloat i];
-    per [roundDown normaliseFloat i,roundUp normaliseFloat s];
-  }
-
-  interval(f:R):% == { 
-    zero?(f) => 0;
-    one?(f)  => 1;
-    -- This next part is necessary to allow e.g. mapping between Expressions:
-    -- AXIOM assumes that Integers stay as Integers!
-    import from Union(value1:Integer,failed:'failed');
-    fnew : R := normaliseFloat f;
-    retractIfCan(f)@Union(value1:Integer,failed:'failed') case value1 =>
-      per [fnew,fnew];
-    per [roundDown fnew, roundUp fnew];
-  }
-
-  qinterval(i:R,s:R):% ==
-    per [roundDown normaliseFloat i,roundUp normaliseFloat s];
-
-  local exactInterval(i:R,s:R):% == per [i,s];
-  local exactSupInterval(i:R,s:R):% == per [roundDown i,s];
-  local exactInfInterval(i:R,s:R):% == per [i,roundUp s];
-
-  inf(u:%):R == (rep u).Inf;
-  sup(u:%):R == (rep u).Sup;
-  width(u:%):R == (rep u).Sup - (rep u).Inf;
-
-  contains?(u:%,f:R):Boolean == (f > inf(u)) and (f < sup(u));
-
-  positive?(u:%):Boolean == inf(u) > 0;
-  negative?(u:%):Boolean == sup(u) < 0;
-
-  (<)(a:%,b:%):Boolean ==
-    if inf(a) < inf(b) then
-      true
-    else if inf(a) > inf(b) then
-      false
-    else
-      sup(a) < sup(b);
-
-  (+)(a:%,b:%):% == {
-    -- A couple of blatent hacks to preserve the Ring Axioms!
-    if zero?(a) then return(b) else if zero?(b) then return(a);
-    if a=b then return qinterval(2*inf(a),2*sup(a));
-    qinterval(inf(a) + inf(b), sup(a) + sup(b));
-  }
-
-  (-)(a:%,b:%):% ==  {
-    if zero?(a) then return(-b) else if zero?(b) then return(a);
-    if a=b then 0 else qinterval(inf(a) - sup(b), sup(a) - inf(b));
-  }
-
-  (*)(a:%,b:%):% == {
-    -- A couple of blatent hacks to preserve the Ring Axioms!
-    if one?(a) then return(b) else if one?(b) then return(a);
-    if zero?(a) then return(0) else if zero?(b) then return(0);
-    prods : List R :=  sort [inf(a)*inf(b),sup(a)*sup(b),
-                             inf(a)*sup(b),sup(a)*inf(b)];
-    qinterval(first prods, last prods);
-  }
-
-  (*)(a:Integer,b:%):% == {
-    if (a > 0) then 
-      qinterval(a*inf(b),a*sup(b));
-    else if (a < 0) then
-      qinterval(a*sup(b),a*inf(b));
-    else
-      0;
-  }
-
-  (*)(a:PositiveInteger,b:%):% == qinterval(a*inf(b),a*sup(b));
-
-  (**)(a:%,n:PositiveInteger):% == {
-    contains?(a,0) and zero?((n pretend Integer) rem 2) =>
-      interval(0,max(inf(a)**n,sup(a)**n)); 
-    interval(inf(a)**n,sup(a)**n);
-  }
-
-  (^) (a:%,n:PositiveInteger):% ==  {
-    contains?(a,0) and zero?((n pretend Integer) rem 2) => 
-      interval(0,max(inf(a)**n,sup(a)**n)); 
-    interval(inf(a)**n,sup(a)**n);
-  }
-
-  (-)(a:%):% == exactInterval(-sup(a),-inf(a));
-
-  (=)(a:%,b:%):Boolean == (inf(a)=inf(b)) and (sup(a)=sup(b));
-  (~=)(a:%,b:%):Boolean == (inf(a)~=inf(b)) or (sup(a)~=sup(b));
-
-  1:% == {one : R := normaliseFloat 1; per([one,one])};
-  0:% == per([0,0]);
-
-  recip(u:%):Union(value1:%,failed:'failed') == {
-   contains?(u,0) => [failed];
-   vals:List R := sort[1/inf(u),1/sup(u)];
-   [qinterval(first vals, last vals)];
-  }
-
-  unit?(u:%):Boolean == contains?(u,0);
-
-  exquo(u:%,v:%):Union(value1:%,failed:'failed') == {
-   contains?(v,0) => [failed];
-   one?(v) => [u];
-   u=v => [1];
-   u=-v => [-1];
-   vals:List R := sort[inf(u)/inf(v),inf(u)/sup(v),sup(u)/inf(v),sup(u)/sup(v)];
-   [qinterval(first vals, last vals)];
-  }
-
-  gcd(u:%,v:%):% == 1;
-
-  coerce(u:Integer):% == {
-    ur := normaliseFloat(u::R);
-    exactInterval(ur,ur);
-  }
-
-  interval(u:Fraction Integer):% == {
-    import { log2 : % -> %;
-             coerce : Integer -> %;
-             retractIfCan : % -> Union(value1:Integer,failed:'failed');}
-    from Float;
-    flt := u::R;
-
-    -- Test if the representation in R is exact
-    --den := denom(u)::Float;
-    local bin : Union(value1:Integer,failed:'failed');
-    bin := retractIfCan(log2(denom(u)::Float));
-    bin case value1 and length(numer u)$Integer < (bits() pretend Integer) => {
-      flt := normaliseFloat flt;
-      exactInterval(flt,flt);
-    }
-
-    qinterval(flt,flt);
-  }
-
-  retractIfCan(u:%):Union(value1:Integer,failed:'failed') == {
-    not zero? width(u) => [failed];
-    retractIfCan inf u;
-  }
-
-  retract(u:%):Integer == {
-    not zero? width(u) =>
-      error "attempt to retract a non-Integer interval to an Integer";
-    retract inf u;
-  }
-
-  coerce(u:%):OutputForm ==
-    bracket([coerce inf(u), coerce sup(u)]$List(OutputForm));
-
-  characteristic():NonNegativeInteger == 0;
-
-
-  -- Explicit export from TranscendentalFunctionCategory
-  pi():% == qinterval(pi(),pi());
-
-  -- From ElementaryFunctionCategory
-  log(u:%):% == {
-    positive?(u) => qinterval(log inf u, log sup u);
-    error "negative logs in interval";
-  }
-
-  exp(u:%):% == qinterval(exp inf u, exp sup u);
-
-  (**)(u:%,v:%):% == {
-    zero?(v) => if zero?(u) then error "0**0 is undefined" else 1;
-    one?(u)  => 1;
-    expts : List R :=  sort [inf(u)**inf(v),sup(u)**sup(v),
-                             inf(u)**sup(v),sup(u)**inf(v)];
-    qinterval(first expts, last expts);
-  }
-
-  -- From TrigonometricFunctionCategory
-
-  -- This function checks whether an interval contains a value of the form
-  -- `offset + 2 n pi'.
-  local hasTwoPiMultiple(offset:R,Pi:R,i:%):Boolean == {
-    import from Integer;
-    next : Integer := retract ceiling( (inf(i) - offset)/(2*Pi) );
-    contains?(i,offset+2*next*Pi);
-  }
-
-  -- This function checks whether an interval contains a value of the form
-  -- `offset + n pi'.
-  local hasPiMultiple(offset:R,Pi:R,i:%):Boolean == {
-    import from Integer;
-    next : Integer := retract ceiling( (inf(i) - offset)/Pi );
-    contains?(i,offset+next*Pi);
-  }
-
-  sin(u:%):% == {
-    import from Integer;
-    Pi : R := pi();
-    hasOne? : Boolean := hasTwoPiMultiple(Pi/(2::R),Pi,u);
-    hasMinusOne? : Boolean := hasTwoPiMultiple(3*Pi/(2::R),Pi,u);
-
-    if hasOne? and hasMinusOne? then 
-      exactInterval(-1,1);
-    else {
-      vals : List R := sort [sin inf u, sin sup u];
-      if hasOne? then
-        exactSupInterval(first vals, 1);
-      else if hasMinusOne? then
-        exactInfInterval(-1,last vals);
-      else
-        qinterval(first vals, last vals);
-    }
-  }
-
-  cos(u:%):% == {
-    Pi : R := pi();
-    hasOne? : Boolean := hasTwoPiMultiple(0,Pi,u);
-    hasMinusOne? : Boolean := hasTwoPiMultiple(Pi,Pi,u);
-
-    if hasOne? and hasMinusOne? then 
-      exactInterval(-1,1);
-    else {
-      vals : List R := sort [cos inf u, cos sup u];
-      if hasOne? then
-        exactSupInterval(first vals, 1);
-      else if hasMinusOne? then
-        exactInfInterval(-1,last vals);
-      else
-        qinterval(first vals, last vals);
-    }
-  }
-
-  tan(u:%):% == {
-    Pi : R := pi();
-    if width(u) > Pi then
-      error "Interval contains a singularity"
-    else {
-      -- Since we know the interval is less than pi wide, monotonicity implies
-      -- that there is no singularity.  If there is a singularity on a endpoint
-      -- of the interval the user will see the error generated by R.
-      lo : R := tan inf u; 
-      hi : R := tan sup u;
-
-      lo > hi => error "Interval contains a singularity";
-      qinterval(lo,hi);
-    }
-  }
-
-  csc(u:%):% == {
-    Pi : R := pi();
-    if width(u) > Pi then
-      error "Interval contains a singularity"
-    else {
-      import from Integer;
-      -- singularities are at multiples of Pi
-      if hasPiMultiple(0,Pi,u) then error "Interval contains a singularity";
-      vals : List R := sort [csc inf u, csc sup u];
-      if hasTwoPiMultiple(Pi/(2::R),Pi,u) then 
-        exactInfInterval(1,last vals);
-      else if hasTwoPiMultiple(3*Pi/(2::R),Pi,u) then
-        exactSupInterval(first vals,-1);
-      else
-        qinterval(first vals, last vals);
-    }
-  }
-
-  sec(u:%):% == {
-    Pi : R := pi();
-    if width(u) > Pi then
-      error "Interval contains a singularity"
-    else {
-      import from Integer;
-      -- singularities are at Pi/2 + n Pi
-      if hasPiMultiple(Pi/(2::R),Pi,u) then
-        error "Interval contains a singularity";
-      vals : List R := sort [sec inf u, sec sup u];
-      if hasTwoPiMultiple(0,Pi,u) then 
-        exactInfInterval(1,last vals);
-      else if hasTwoPiMultiple(Pi,Pi,u) then
-        exactSupInterval(first vals,-1);
-      else
-        qinterval(first vals, last vals);
-    }
-  }
-
-
-  cot(u:%):% == {
-    Pi : R := pi();
-    if width(u) > Pi then
-      error "Interval contains a singularity"
-    else {
-      -- Since we know the interval is less than pi wide, monotonicity implies
-      -- that there is no singularity.  If there is a singularity on a endpoint
-      -- of the interval the user will see the error generated by R.
-      hi : R := cot inf u; 
-      lo : R := cot sup u;
-
-      lo > hi => error "Interval contains a singularity";
-      qinterval(lo,hi);
-    }
-  }
-
-  -- From ArcTrigonometricFunctionCategory
-
-  asin(u:%):% == {
-    lo : R := inf(u);
-    hi : R := sup(u);
-    if (lo < -1) or (hi > 1) then error "asin only defined on the region -1..1";
-    qinterval(asin lo,asin hi);
-  }
-
-  acos(u:%):% == {
-    lo : R := inf(u);
-    hi : R := sup(u);
-    if (lo < -1) or (hi > 1) then error "acos only defined on the region -1..1";
-    qinterval(acos hi,acos lo);
-  }
-
-  atan(u:%):% == qinterval(atan inf u, atan sup u);
-
-  acot(u:%):% == qinterval(acot sup u, acot inf u);
-
-  acsc(u:%):% == {
-    lo : R := inf(u);
-    hi : R := sup(u);
-    if ((lo <= -1) and (hi >= -1)) or ((lo <= 1) and (hi >= 1)) then
-      error "acsc not defined on the region -1..1";
-    qinterval(acsc hi, acsc lo);
-  }
-
-  asec(u:%):% == {
-    lo : R := inf(u);
-    hi : R := sup(u);
-    if ((lo < -1) and (hi > -1)) or ((lo < 1) and (hi > 1)) then
-      error "asec not defined on the region -1..1";
-    qinterval(asec lo, asec hi);
-  }
-
-  -- From HyperbolicFunctionCategory
-
-  tanh(u:%):% == qinterval(tanh inf u, tanh sup u);
-
-  sinh(u:%):% == qinterval(sinh inf u, sinh sup u);
-
-  sech(u:%):% == {
-    negative? u => qinterval(sech inf u, sech sup u);
-    positive? u => qinterval(sech sup u, sech inf u);
-    vals : List R := sort [sech inf u, sech sup u];
-    exactSupInterval(first vals,1);
-  }
-
-  cosh(u:%):% == {
-    negative? u => qinterval(cosh sup u, cosh inf u);
-    positive? u => qinterval(cosh inf u, cosh sup u);
-    vals : List R := sort [cosh inf u, cosh sup u];
-    exactInfInterval(1,last vals);
-  }
-
-  csch(u:%):% == {
-    contains?(u,0) => error "csch: singularity at zero";
-    qinterval(csch sup u, csch inf u);
-  }
-
-  coth(u:%):% == {
-    contains?(u,0) => error "coth: singularity at zero";
-    qinterval(coth sup u, coth inf u);
-  }
-
-  -- From ArcHyperbolicFunctionCategory
-
-  acosh(u:%):% == {
-    inf(u)<1 => error "invalid argument: acosh only defined on the region 1..";
-    qinterval(acosh inf u, acosh sup u);
-  }
-
-  acoth(u:%):% == {
-    lo : R := inf(u);
-    hi : R := sup(u);
-    if ((lo <= -1) and (hi >= -1)) or ((lo <= 1) and (hi >= 1)) then
-      error "acoth not defined on the region -1..1";
-    qinterval(acoth hi, acoth lo);
-  }
-
-  acsch(u:%):% == {
-    contains?(u,0) => error "acsch: singularity at zero";
-    qinterval(acsch sup u, acsch inf u);
-  }
-
-  asech(u:%):% == {
-    lo : R := inf(u);
-    hi : R := sup(u);
-    if  (lo <= 0) or (hi > 1) then 
-      error "asech only defined on the region 0 < x <= 1";
-    qinterval(asech hi, asech lo);
-  }
-
-  asinh(u:%):% == qinterval(asinh inf u, asinh sup u);
-
-  atanh(u:%):% == {
-    lo : R := inf(u);
-    hi : R := sup(u);
-    if  (lo <= -1) or (hi >= 1) then 
-      error "atanh only defined on the region -1 < x < 1";
-    qinterval(atanh lo, atanh hi);
-  }
-
-  -- From RadicalCategory
-  (**)(u:%,n:Fraction Integer):% == interval(inf(u)**n,sup(u)**n);
-  
-}
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
-<<IntervalCategory>>
-<<Interval>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/invnode.as.pamphlet b/src/algebra/invnode.as.pamphlet
deleted file mode 100644
index d051753..0000000
--- a/src/algebra/invnode.as.pamphlet
+++ /dev/null
@@ -1,340 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra invnode.as}
-\author{The Axiom Team}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{IVNodeCategory}
-<<IVNodeCategory>>=
-#include "axiom"
-
-POINT ==> Point DoubleFloat;
-
-local DF2S;
-DF2S(u:DoubleFloat):String == {
-    STRINGIMAGE ==> VMLISP_:_:STRINGIMAGE;
-    import { STRINGIMAGE : DoubleFloat -> String} from Foreign Lisp;
-    STRINGIMAGE(u);
-}
-
-+++ Category  of all open inventor node types
-+++ Uses IVObject as a 'generic' value.
-define IVNodeCategory: Category == SetCategory with {
-	quickWrite: (TextFile, %) 	-> ();
-		++ Quick version. Not guaranteed to terminate
-	children:      % 		-> List IVNodeObject;
-	addChild!:    (%, IVNodeObject) -> ();
-	fields:        %		-> List IVField;
-	className:     %		-> String;	
-	coerce:	       %		-> IVNodeObject;
-	default {
-		import from Symbol;
-		quickWrite(out: TextFile, node: %): () == {
-			write!(out, className(node));
-			write!(out, " {");
-			writeLine!(out);
-			import from List IVField;
-			import from IVValue;
-			for field in fields node repeat {
-				write!(out, string name field);
-				write!(out, " ");
-				invWrite(out, value field);
-			}
-			writeLine!(out, "}");
-		}
-		coerce(x: %): IVNodeObject == 
-			make(% pretend IVNodeCategory, x);
-
-		coerce(x: %): OutputForm == {
-			import from String;
-			coerce className x;
-		}
-	}
-}
-
-@
-\section{IVLeafNodeCategory}
-<<IVLeafNodeCategory>>=
-+++ Category for leaves --- just adds a few defaults to make life 
-+++ easy.
-define IVLeafNodeCategory: Category == IVNodeCategory with {
-	default {
-		children(v: %): List IVNodeObject == [];
-		addChild!(v: %, new: IVNodeObject): () == 
-			error "can't add child to a leaf";
-	}
-}
-
-@
-\section{IVNodeObject}
-<<IVNodeObject>>=
--- virtual functions for fun and profit...
-IVNodeObject: IVNodeCategory with {
-	make:     (T: IVNodeCategory, T) -> %;
-	coerce:   (T: IVNodeCategory, %) -> T;
-	uniqueID: % -> Integer;
-} == add {
-	Rep ==> Record(NT: IVNodeCategory, val: NT, idx: Integer);
-	import from Rep;
-	default z: Integer;
-
-	local iCount: Integer := 0;
-	local explode: (o: %) -> (NodeType: IVNodeCategory, NodeType);
-
-	uniqueID(o: %): Integer == rep(o).idx;
-
-	explode(o: %): (NodeType: IVNodeCategory, v: NodeType) == {
-		(NT, val, id) == explode rep o;
-		(NT, val);
-	}
-
-	make(T: IVNodeCategory, val: T):   % == {
-		free iCount := iCount + 1;
-		per [T, val, iCount];
-	}
-	coerce(T: IVNodeCategory, val: %): T == {
-		(type, v, id) == explode rep val;
-		v pretend T;
-	}
-
-	-- The '0' functions are needed to turn non-constants 
-	-- (eg. fn return values) -- into constants.
-	children(v: %): List IVNodeObject == {
-		children0(NodeType: IVNodeCategory, val: NodeType): 
-							List IVNodeObject == 
-			children val;
-		children0 explode v;
-	}
-
-	fields(v: %): List IVField == {
-		fields0(NodeType: IVNodeCategory, val: NodeType): List IVField == 
-			fields val;
-		fields0 explode v;
-	}
-
-	className(v: %): String == {
-		name0(NodeType: IVNodeCategory, val: NodeType): String == 
-			className(val)$NodeType;
-		name0 explode v;
-	}
-
-	addChild!(v: %, child: %): () == { 
-		addChild0!(NodeType: IVNodeCategory, val: NodeType): () == 
-			addChild!(val, child);
-		addChild0! explode v;
-	}
-
-	-- BasicType stuff
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality on ivobject";
-}
-
-@
-\section{IVNodeConnection}
-<<IVNodeConnection>>=
-IVNodeConnection: with {
-	bracket: (IVNodeObject, Symbol) -> %;
-	field: % -> Symbol;
-	node:  % -> IVNodeObject;
-} == add {
-	Rep ==> Record(o: IVNodeObject, f: Symbol);
-	import from Rep;
-
-	[o: IVNodeObject, f: Symbol]: % == per [o,f];
-	field(c: %): Symbol == rep(c).f;
-	node(c: %): IVNodeObject == rep(c).o;
-}
-
-@
-\section{IVValue}
-<<IVValue>>=
-IVValue: BasicType with {
-	DECL(T, fld, flg) ==> {
-		coerce: % -> T;
-		flg: % -> Boolean;
-		fld: T -> %;
-	}
-	DECL(DoubleFloat,        float,     float?);
-	DECL(IVNodeObject,       node, 	    node?);
-	DECL(Boolean,            bool, 	    bool?);
-	DECL(SingleInteger,      int, 	    int?);
-	DECL(String, 	     	 string,    string?);
-	DECL(Symbol,	         symbol,    symbol?);
-	DECL(POINT,		 point,     point?);
-	DECL(List DoubleFloat,   floatlist, floatlist?);
-	DECL(List SingleInteger, intlist,   intlist?);
-	DECL(List POINT,	 pointlist, pointlist?);
-	DECL(IVNodeConnection,   connect,   connect?);
-
-	invWrite: (TextFile, %) -> ();
-} == add {
-	Rep ==> Union(  float:     DoubleFloat,
-			node:      IVNodeObject,
-			bool:      Boolean,
-			int:       SingleInteger,
-			string:    String,
-			symbol:    Symbol,
-			point:	   POINT,
-			intlist:   List SingleInteger,
-			floatlist: List DoubleFloat,
-			pointlist: List POINT,
-			connect:   IVNodeConnection
-		     );
-	import from Rep;
-	
-	Accessor(T, fld, flg) ==> { 
-		coerce(x: %):  	 T == rep(x).fld;
-		flg(x: %): Boolean == rep(x) case fld;
-		fld(x: T): %       == per [x, fld];
-	}
-	Accessor(DoubleFloat,        float, 	float?);
-	Accessor(IVNodeObject,       node, 	node?);
-	Accessor(Boolean,            bool, 	bool?);
-	Accessor(SingleInteger,      int, 	int?);
-	Accessor(String, 	     string,	string?);
-	Accessor(Symbol,	     symbol,    symbol?);
-	Accessor(POINT,		     point,     point?);
-	Accessor(List DoubleFloat,   floatlist, floatlist?);
-	Accessor(List SingleInteger, intlist,   intlist?);
-	Accessor(List POINT,	     pointlist, pointlist?);
-	Accessor(IVNodeConnection,   connect,   connect?);
-
-	local ppoint(out: TextFile, val: POINT, dim: Integer): () == {
-		for i in 1..dim repeat {
-			write!(out, DF2S(val.(i::Integer)));
-			write!(out, " ");
-		}
-	}
-	invWrite(out: TextFile, val: %): () == {
-		import from Float, Integer;
-		float? val => {
-			writeLine!(out, 
-				   convert(convert(val::DoubleFloat)$Float));
-		}
-		node? val or connect? val => {
-			error "Sorry, can't write a node here";
-			--writeLine!(out, val::IVNodeObject);
-		}
-		bool? val => {
-			writeLine!(out, 
-				   if val::Boolean then "true" else "false");
-		}
-		int? val => {
-			writeLine!(out, 
-				 convert(convert(val::SingleInteger)@Integer));
-		}
-		string? val => {
-			writeLine!(out, val::String);
-		}
-		symbol? val => {
-			writeLine!(out, string(val::Symbol));
-		}
-		point? val => {
-			ppoint(out, rep(val).point, 3);
-			writeLine!(out, "");
-		}
-		floatlist? val => {
-			write!(out, "[ ");
-			for fl in val::List DoubleFloat repeat {
-				write!(out,convert(convert(fl)$Float));
-				write!(out, ", ");
-			}
-			writeLine!(out, "]");
-		}
-		intlist? val => {
-			write!(out, "[ ");
-			for i in val::List SingleInteger repeat {
-				write!(out,convert(convert(i)@Integer));
-				write!(out, ", ");
-			}
-			writeLine!(out, "]");
-		}
-		pointlist? val => {
-			write!(out, "[ ");
-			for p in val::List POINT repeat {
-				ppoint(out, p, 3);
-				writeLine!(out, ",");
-			}
-			writeLine!(out, "]");
-		}
-		never
-	}
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality for values";
-}
-
-@
-\section{IVField}
-<<IVField>>=
-IVField: BasicType with {
-	new: (Symbol,IVValue) -> %;
-	name:  % -> Symbol;
-	value: % -> IVValue;
-} == add {
-	Rep ==> Record(name: Symbol, v: IVValue);
-	import from Rep;
-
-	new(name: Symbol, val: IVValue): % == per [name, val];
-	name(f: %):  Symbol == rep(f).name;
-	value(f: %): IVValue == rep(f).v;
-
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality for values";
-}
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
-<<IVNodeCategory>>
-<<IVLeafNodeCategory>>
-<<IVNodeObject>>
-<<IVNodeConnection>>
-<<IVValue>>
-<<IVField>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/invrender.as.pamphlet b/src/algebra/invrender.as.pamphlet
deleted file mode 100644
index a107b9b..0000000
--- a/src/algebra/invrender.as.pamphlet
+++ /dev/null
@@ -1,172 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra invrender.as}
-\author{The Axiom Team}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{RenderTools}
-<<RenderTools>>=
-#include "axiom"
-
-POINT ==> Point DoubleFloat;
-NNI   ==> NonNegativeInteger;
-SI    ==> SingleInteger;
-
-RenderTools: with {
-	renderToFile!:  (FileName, ThreeSpace DoubleFloat) -> ();
-	makeSceneGraph: (ThreeSpace DoubleFloat) 	   -> IVNodeObject;
-} == add {
-	import from IVUtilities;
-
-	renderToFile!(f: FileName, space: ThreeSpace DoubleFloat): () == {
-		root := makeSceneGraph(space);
-		write!(f, root);
-	}
-
-	local makePts: (lpts: List POINT, 
-			indicies: List List List NonNegativeInteger) -> 
-				(List POINT, List DoubleFloat);
-
-	local makePts(lp: List POINT, indicies: List List List NonNegativeInteger): 
-				(List POINT, List DoubleFloat) == {
-		local colorIdx: Integer;
-		indexList := concat concat indicies;
-		coordpts := lp;
-		if (# first lp = 4) then colorIdx := 4 else colorIdx := 3;
-		colors := [pt.colorIdx for pt in lp];
-		(coordpts, colors)
-	}
-
-	local makeBaseColor(l: List DoubleFloat): IVBaseColor == {
-		-- This works by interpolating between blue and green (via cyan).
-		-- There may well be better ways...
-		import from POINT;
-		import from List POINT;
-		import from DoubleFloat;
-		import from List DoubleFloat;
-		low  := 10000.0;
-		high := -10000.0;
-		for df in l repeat {
-			if low  > df then low  := df;
-			if high < df then high := df;
-		}
-		if (high = low) then high := high + 1.0;
-		new [ point([0, (df - low)/(high - low), (high - df)/(high - low)])
-			for df in l]
-	}
-	makeSceneGraph(space: ThreeSpace DoubleFloat): IVNodeObject == {
-		import from List ThreeSpace DoubleFloat;
-		import from List List List NNI;
-		import from Integer;
-		import from Symbol;
-		import from IVValue;
-		check(space);
-		lpts := lp(space);
-	      	indicies := lllip(space);
-		root: IVSeparator  := new();
-		(coordpts, colorvalues) := makePts(lpts, indicies);
-		coords: IVCoordinate3 := new coordpts;
-		colors: IVBaseColor   := makeBaseColor(colorvalues);
-		addChild!(root, coerce coords);
-		addChild!(root, coerce colors);
-		binding: IVBasicNode := make "MaterialBinding";
-		addField!(binding, "value", symbol "PER__VERTEX");
-		addChild!(root, coerce binding);
-		offset: NNI := 0;
-	      	for ss in components space
-	      	for index in indicies repeat {
-			local coordIndex: List NNI;
-			default i: Integer;
-	        	closedCurve? ss => {
-        			n: Integer := (#(index.1))::Integer;
-          			coordIndex := 
-					[offset+coerce(i) for i in 0..n::Integer];
-				-- Close the curve
-          			setlast!(coordIndex,offset);
-          			curve :  IVIndexedLineSet := new coordIndex;
-          			addChild!(root, coerce curve);
-          			offset := offset+n::NNI;
-			}
-	        	curve? ss => {
-	          		n := (#(index.1))::Integer;
-	          		coordIndex := 
-					[offset+coerce(i) for i in 0..(n-1)];
-	          		curve :  IVIndexedLineSet := new coordIndex;
-	          		addChild!(root, coerce curve);
-	          		offset := offset+n::NNI;
-			}
-	        	polygon? ss => {
-	          		vertices := #(index.1) + #(index.2);
-	          		face : IVFaceSet := new(vertices::SI,offset::SI);
-	          		addChild!(root, coerce face);
-	          		offset := offset+vertices;
-			}
-	        	mesh? ss => {
-	          		xStep: SingleInteger := (#index)::SingleInteger;
-	          		yStep: SingleInteger := (#(first index))::SingleInteger;
-	          		quadMesh : IVQuadMesh := 
-						new(xStep,yStep,offset::SingleInteger);
-	          		addChild!(root, coerce quadMesh);
-	          		offset := offset+coerce(xStep*yStep);
-			}
-	        	point? ss => {
-	          		pt : IVPointSet := new(offset::SingleInteger, 
-						       1$SingleInteger);
-	          		addChild!(root, coerce pt);
-	          		offset := offset+1;
-			}
-	        	error "Unrecognised SubSpace component";
-		}
-		coerce root;
-	}
-}
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
-<<RenderTools>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/invtypes.as.pamphlet b/src/algebra/invtypes.as.pamphlet
deleted file mode 100644
index 5ada870..0000000
--- a/src/algebra/invtypes.as.pamphlet
+++ /dev/null
@@ -1,302 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra invtypes.as}
-\author{The Axiom Team}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{IVSimpleInnerNode}
-<<IVSimpleInnerNode>>=
-#include "axiom"
-
-import from IVValue, Symbol;
-
-POINT ==> Point DoubleFloat;
-NNI ==> NonNegativeInteger;
-
-IVSimpleInnerNode: with {
-	new: 	   ()                -> %;
-	addChild!: (%, IVNodeObject) -> ();
-	children:   %                -> List IVNodeObject;
-	fields:     %                -> List IVField;
-
-	=: (%, %) -> Boolean;
-
-} == add {
-	Rep ==> Record(lst: List IVNodeObject);
-	import from Rep;
-
-	new(): % == per [[]];
-	addChild!(v: %, new: IVNodeObject): () == {
-		rep(v).lst := concat!(rep(v).lst, new);
-	}
-
-	children(v: %): List IVNodeObject == rep(v).lst;
-
-	fields(node: %): List IVField == [];
-
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality on IVInnerNodes";
-}
-
-@
-\section{IVSeparator}
-<<IVSeparator>>=
-IVSeparator: IVNodeCategory with {
-	new: () -> %;
-} == IVSimpleInnerNode add {
-	className(v: %): String == "Separator";
-}
-
-@
-\section{IVGroup}
-<<IVGroup>>=
-IVGroup: IVNodeCategory with {
-	new: () -> %;
-} == IVSimpleInnerNode add {
-	className(v: %): String == "Group";
-}
-
-@
-\section{IVCoordinate3}
-<<IVCoordinate3>>=
-IVCoordinate3: IVLeafNodeCategory with {
-	new: List POINT -> %;
-} == add {
-	Rep ==> List POINT;
-	className(x: %): String == "Coordinate3";
-
-	new(l: List POINT): % == per l;
-	fields(v: %): List IVField == [ new("point", pointlist rep v)];
-	
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality on IVCoord3";
-}
-
-@
-\section{IVCoordinate4}
-<<IVCoordinate4>>=
-IVCoordinate4: IVLeafNodeCategory with {
-	new: List POINT -> %;
-} == add {
-	Rep ==> List POINT;
-	import from Rep;
-
-	className(x: %): String == "Coordinate4";
-
-	new(l: List POINT): % == per l;
-	fields(v: %): List IVField == [ new("point", pointlist rep v)];
-	
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality on IVCoord4";
-}
-
-@
-\section{IVQuadMesh}
-<<IVQuadMesh>>=
-IVQuadMesh: IVLeafNodeCategory with {
-	new: (SingleInteger, SingleInteger, SingleInteger) -> %;
-} == add {
-	Rep ==> Record( rowc: SingleInteger, 
-			colc: SingleInteger, 
-			start: SingleInteger);
-	import from Rep;
-
-	className(x: %): String == "QuadMesh";
-	
-	new(rc: SingleInteger, cc: SingleInteger, start: SingleInteger): % == 
-		per [rc, cc, start];
-
-	fields(v: %): List IVField == [
-		new("verticesPerColumn", int rep(v).colc),
-		new("verticesPerRow",    int rep(v).rowc),
-		new("startIndex", 	 int rep(v).start)
-	];
-
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality on IVQuadMesh";
-}
-
-@
-\section{IVBaseColor}
-<<IVBaseColor>>=
-IVBaseColor: IVLeafNodeCategory with {
-	new: List POINT -> %;
-} == add {
-	Rep ==> List POINT;
-	import from Rep;
-
-	className(x: %): String == "BaseColor";
-
-	new(l: List POINT): % == per l;
-	fields(v: %): List IVField == [ new("rgb", pointlist rep v) ];
-	
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality on IVBaseColor";
-}
-
-@
-\section{IVIndexedLineSet}
-<<IVIndexedLineSet>>=
-IVIndexedLineSet: IVLeafNodeCategory with {
-	new: List NNI -> %;
-	new: List SingleInteger -> %;
-} == add {
-	Rep ==> List SingleInteger;
-	import from Rep;
-
-	className(x: %): String == "IndexedLineSet";
-
-	new(l: List SingleInteger): % == per l;
-	new(l: List NNI): 	    % == new [ coerce n for n in l];
-
-	fields(v: %): List IVField == [ new("points", intlist rep v) ];
-	
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality on IVBaseColor";
-}
-
-@
-\section{IVFaceSet}
-<<IVFaceSet>>=
-IVFaceSet: IVLeafNodeCategory with {
-	new: (SingleInteger, SingleInteger) -> %;
-} == add {
-	Rep ==> Record(startIndex: SingleInteger, numVertices: SingleInteger);
-	import from Rep;
-
-	className(x: %): String == "FaceSet";
-
-	new(x: SingleInteger, y: SingleInteger): % == per [x,y];
-	fields(v: %): List IVField == [ 
-		new("numVertices", int rep(v).numVertices),
-		new("startIndex",  int rep(v).startIndex)
-	];
-	
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality on IVFaceSet";
-}
-
-@
-\section{IVPointSet}
-<<IVPointSet>>=
-IVPointSet: IVLeafNodeCategory with {
-	new: (SingleInteger, SingleInteger) -> %;
-} == add {
-	Rep ==> Record(startIndex: SingleInteger, numPoints: SingleInteger);
-	import from Rep;
-
-	className(x: %): String == "PointSet";
-
-	new(x: SingleInteger, y: SingleInteger): % == per [x,y];
-
-	fields(v: %): List IVField == [ 
-		new("numPoints", int rep(v).numPoints),
-		new("startIndex",  int rep(v).startIndex)
-	];
-	
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality on IVFaceSet";
-}
-
-@
-\section{IVBasicNode}
-<<IVBasicNode>>=
-IVBasicNode: IVNodeCategory with {
-	make: String -> %;
-	addField!: (%, IVField) -> ();
-	addField!: (%, Symbol, IVValue) -> ();
-} == add {
-	Rep ==> Record(class: String,
-		       kids: List IVNodeObject, 
-		       fields: List IVField);
-	import from Rep, IVField;
-
-	make(name: String): % == per [name, [], []];
-
-	className(node: %): String == rep(node).class;
-	children(node: %): List IVNodeObject == rep(node).kids;
-	fields(node: %): List IVField == rep(node).fields;
-
-	addField!(node: %, fld: IVField): () == {
-		rep(node).fields := cons(fld, rep(node).fields);
-	}
-
-	addChild!(node: %, kid: IVNodeObject): () == {
-		rep(node).kids := cons(kid, rep(node).kids);
-	}
-	
-	addField!(node: %, sym: Symbol, val: IVValue): () ==
-		addField!(node, new(sym, val));
-
-	--
-	sample: % == % pretend %;
-	(=)(a: %, b: %): Boolean == error "no equality on IVBasicNode";
-}
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
-<<IVSimpleInnerNode>>
-<<IVSeparator>>
-<<IVGroup>>
-<<IVCoordinate3>>
-<<IVCoordinate4>>
-<<IVQuadMesh>>
-<<IVBaseColor>>
-<<IVIndexedLineSet>>
-<<IVFaceSet>>
-<<IVPointSet>>
-<<IVBasicNode>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/invutils.as.pamphlet b/src/algebra/invutils.as.pamphlet
deleted file mode 100644
index e375133..0000000
--- a/src/algebra/invutils.as.pamphlet
+++ /dev/null
@@ -1,172 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra invutils.as}
-\author{The Axiom Team}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{IVUtilities}
-<<IVUtilities>>=
-#include "axiom"
-
-IVUtilities: with {
-	walkTree: ((IVNodeObject, Boolean) -> Boolean, 
-		   (IVNodeObject, Boolean) -> Boolean, 
-		   IVNodeObject, Boolean) -> ();
-	write!:   (TextFile, IVNodeObject) -> ();
-	write!:   (FileName, IVNodeObject) -> ();
-} == add {
-	-- walk a tree from 'root', and call f on each node.
-	-- nodesOnly will stop the recursion finding subnodes within
-	-- fields.
-	-- preFn is a function that takes a node, a flag indicating if the 
-	-- node has already been traversed.  It returns a flag
- 	-- indicating if the traversal should descend the node.
-	walkTree(preFn:  (IVNodeObject, Boolean) -> Boolean, 
-		 postFn: (IVNodeObject, Boolean) -> Boolean, 
-		 root: 	    IVNodeObject, 
-		 nodesOnly: Boolean): () == {
-		tab: Table(Integer, IVNodeObject) := table();
-		innerWalk(node: IVNodeObject): () == {
-			import from List IVNodeObject;
-			import from List IVField;
-			import from IVValue;
-			present := key?(uniqueID node, tab);
-			not preFn(node, present) => return;
-			tab.(uniqueID node) := node;
-			for child in children node repeat
-				innerWalk(child);
-			if not nodesOnly then {
-				for fld in fields node repeat {
-					import from IVNodeConnection;
-					connect? value fld =>
-						innerWalk(node coerce value fld);
-					node? value fld =>
-						innerWalk(coerce value fld);
-				}
-			}
-			postFn(node, false);
-		}
-		innerWalk(root);
-	}
-
-	write!(out: TextFile, root: IVNodeObject): () == {
-		import from Boolean;
-		names: Table(Integer, IVNodeObject) := table();
-
-		getName(node: IVNodeObject): String == {
-			import from Integer;
-		convert uniqueID node;
-		}
-
-		doNamingVisit(node: IVNodeObject, flag: Boolean): Boolean == {
-			if flag then names.(uniqueID node) := node;
-			flag
-		}
-		writeNodeHeader(node: IVNodeObject): () == {
-			present := key?(uniqueID node, names);
-			if present then {
-				write!(out, "DEF ");
-				write!(out, getName node);
-			}
-		}
-		doPrintingVisit(node: IVNodeObject, 
-				flag: Boolean): Boolean == {
-			if flag then {
-				write!(out, "USE ");
-				write!(out, getName node);
-				return false;
-			}
-			write!(out, className(node));
-			writeLine!(out, " {");
-			import from List IVField, Symbol;
-			for field in fields node repeat {
-				import from IVNodeConnection;
-				val: IVValue := value field;
-				write!(out, string name field);
-				write!(out, " ");
-				node? val => {
-					walkTree(doPrintingVisit,
-						 doFinalPrint,
-						 coerce val, false);
-				}
-				connect? val => {
-					walkTree(doPrintingVisit,
-						 doFinalPrint, 
-						 node coerce val, false);
-					write!(out, ".");
-					writeLine!(out,
-						   string field coerce val);
-				}
-				-- simple case:
-				invWrite(out, value field);
-			}
-			return true;
-		}
-
-		doFinalPrint(node: IVNodeObject, x: Boolean): Boolean == {
-			writeLine!(out, "}");
-			true;
-		}
-		doNothing(node: IVNodeObject, x: Boolean): Boolean == x;
-
-		writeLine!(out, "#Inventor V2.0 ascii");
-		walkTree(doNamingVisit, doNothing, root, true);
-		walkTree(doPrintingVisit, doFinalPrint, root, false);
-	}
-
-	write!(file: FileName, root:IVNodeObject): () == {
-		out: TextFile := open(file, "output");
-		write!(out, root);
-		close!(out);
-	}	
-}
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
-<<IVUtilities>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/iviews.as.pamphlet b/src/algebra/iviews.as.pamphlet
deleted file mode 100644
index 37a09af..0000000
--- a/src/algebra/iviews.as.pamphlet
+++ /dev/null
@@ -1,330 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra iviews.as}
-\author{The Axiom Team}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{InventorDataSink}
-<<InventorDataSink>>=
-#include "axiom"
-#assert Real
-
-NNI 	==> NonNegativeInteger;
-SI 	==> SingleInteger;
-DF 	==> DoubleFloat;
-POINT 	==> Point DF;
-SPACE3  ==> ThreeSpace DoubleFloat;
-
-DefaultSize ==> 65530;
-
-Value ==> Symbol;
-
-InventorDataSink: with {
-	CoercibleTo OutputForm;
-	new:      () -> %;
-	dispose!: % -> ();
-
-	put!:     (%, SI) -> ();
-	put!:     (%, DF) -> ();
-	put!: (%, String) -> ();
-
-	vstart!: (%, 'int,float', SI) -> ();
-	vput!:  (%, SI) -> ();
-	vput!:  (%, DF) -> ();
-
-	lstart!:   % -> ();
-	lend!:     % -> ();
-	export from 'int,float'
-} == add {
-#if Real
-	-- No rep (we cheat!)
-	import from SI;
-	valOf(x) ==> x pretend Value;
-	default sink: %;
-	import {
-		LISP_:_:GR_-GET_-MEM_-AREA: SI -> %;
-		LISP_:_:GR_-KILL_-MEM_-AREA: % -> ();
-		LISP_:_:GR_-PUT_-ITEM:	    (%, Value) -> ();
-		LISP_:_:GR_-PUT_-LSTART:	     % -> ();
-		LISP_:_:GR_-PUT_-LEND:	             % -> ();
-		LISP_:_:GR_-INIT_-VECTOR: (%, Value, Value) -> %;
-		LISP_:_:GR_-ADD_-TO_-VECTOR:     (%, Value) -> %;
-	} from Foreign Lisp;
-
-	new(): % == LISP_:_:GR_-GET_-MEM_-AREA(DefaultSize);
-	dispose!(sink): () == LISP_:_:GR_-KILL_-MEM_-AREA(sink);
-
-	put!(sink, si: SI): ()     == LISP_:_:GR_-PUT_-ITEM(sink, valOf(si));
-	put!(sink, st: String): () == LISP_:_:GR_-PUT_-ITEM(sink, valOf(st));
-	put!(sink, fl: DF): ()     == LISP_:_:GR_-PUT_-ITEM(sink, valOf(fl));
-
-	vstart!(sink, type: 'int,float', sz: SI): () == {
-		local sym: Symbol;
-		if type = int then 
-			sym := coerce("integer");
-		else 
-			sym := coerce("float");
-		LISP_:_:GR_-INIT_-VECTOR(sink, valOf(sym), valOf(sz));
-	}
-
-	vput!(sink, si: SI): () == 
-		LISP_:_:GR_-ADD_-TO_-VECTOR(sink, valOf(si));
-	vput!(sink, df: DF): () == 
-		LISP_:_:GR_-ADD_-TO_-VECTOR(sink, valOf(df));
-
-	lstart!(sink): () == LISP_:_:GR_-PUT_-LSTART sink;
-	lend!(sink): () == LISP_:_:GR_-PUT_-LEND sink;
-
-	coerce(sink): OutputForm == {
-		[outputForm "aSink"]
-	}
-#else
-	Rep ==> Record(count: NonNegativeInteger);
-	import from Rep, NNI;
-	default sink: %;
-	coerce(sink): OutputForm == {
-		import from List OutputForm;
-		bracket [outputForm "Sink: ", 
-			 outputForm coerce rep(sink).count];
-	}
-
-	local addn!(sink, n: NNI): () == 
-		rep(sink).count := rep(sink).count + n;
-	new(): % == per [0];
-	dispose!(sink): () == dispose! rep sink;
-	
-	put!(sink, n: SI): () == addn!(sink, 1 + 4);
-	put!(sink, f: DF): () == addn!(sink, 1 + 4);
-	put!(sink, s: String): () == {
-		addn!(sink, #s + 1 + 1);
-	}
-	
-	vstart!(sink, type: 'int, float', n: SI): () == {
-		addn!(sink, 1 + n::NNI*4);
-	}
-
-	vput!(sink, n: SI): () == {};
-	vput!(sink, n: DF): () == {};
-
-	lstart!(sink): () == addn!(sink, 1);
-	lend!(sink): () == addn!(sink, 1);
-
-#endif
-}
-
-@
-\section{InventorViewPort}
-<<InventorViewPort>>=
-InventorViewPort: with {
-	new:      () -> %;
-	new: ThreeSpace DoubleFloat -> %;
-	addData!: (%, InventorDataSink) -> %;
-	addData!: (%, ThreeSpace DoubleFloat) -> %;
-} == add {
-#if Real
-	import {
-		LISP_:_:GR_-MAKE_-VIEW: (SI) -> %;
-		LISP_:_:GR_-SET_-DATA:  (%, InventorDataSink) -> ();
-	} from Foreign Lisp;
-	import from SingleInteger;
-
-	new(): % == LISP_:_:GR_-MAKE_-VIEW(0);
-
-	new(space: ThreeSpace DoubleFloat): % == {
-		import from InventorDataSink;
-		import from InventorRenderPackage;
-		view: % := new();
-		addData!(view, space);
-		view
-	}
-
-	addData!(view: %, data: InventorDataSink): % == {
-		LISP_:_:GR_-SET_-DATA(view, data);
-		view;
-	}
-
-	addData!(view: %, space: ThreeSpace DoubleFloat): % == {
-		import from InventorRenderPackage;
-		sink: InventorDataSink := new();
-		render(sink, space, cartesian$CoordinateSystems(DoubleFloat));
-		addData!(view, sink);
-		view
-	}
-
-#else
-	Rep ==> SingleInteger;
-	import from Rep;
-	new(): % == per 1;
-	new(x: ThreeSpace DoubleFloat): % == per 2;
-	addData!(view: %, data: InventorDataSink): % == view;
-#endif
-	
-}
-
-@
-\section{InventorRenderPackage}
-<<InventorRenderPackage>>=
-InventorRenderPackage: with {
-	render: (InventorDataSink, ThreeSpace DoubleFloat, POINT->POINT) -> ();
-} == add {
-	default sink: InventorDataSink;
-	default space: ThreeSpace DoubleFloat;
-	default transform: POINT->POINT;
-	import from SI;
-
-	local put!(sink, dims: UniversalSegment SI, 
-		   lp: List Point DoubleFloat,
-		   f: Point DoubleFloat -> Point DoubleFloat): () == {
-		import from NNI, Integer;
-		i : SI := 0;
-		for x in dims repeat i:= i+1;
-		vstart!(sink, float, i*(coerce #lp));
-		for p in lp repeat {
-			p1 := f(p);
-			for idx in dims repeat 
-				vput!(sink, p1.(idx::Integer));
-		}
-	}
-
-	local put!(sink, lp: List SI): () == {
-		import from NNI;
-		vstart!(sink, int, coerce #lp);
-		for p in lp repeat {
-			vput!(sink, p);
-		}
-	}
-
-	local putPoints!(sink, transform, 
-			 lpts: List POINT, indexList: List NNI): () == {
-		import from Integer;
-		if not sorted? indexList 
-		then {
-			-- not nice!
-			lst: List POINT := [];
-			for idx in indexList repeat 
-				lst := cons(lpts.(coerce idx), lst);
-			lpts := reverse! lst;
-		}
-		put!(sink, 1..3, lpts, transform);
-      		if (# first lpts) = 4 
-		then {
-			put!(sink, "Colours");
-			put!(sink, 4..4, lpts, transform);
-		}
-	}
-	render(sink, space, transform): () == {
-	 	default ss: SPACE3;
-	      	default i: NNI;
-		import from List POINT;
-		import from List List List NNI;
-		import from List List NNI;
-		import from List SPACE3;
-		import from SingleInteger;
-		put!(sink, "ThreeDScene");
-      		-- Get the point data
-	      	check(space);
-      		indices := lllip(space);
-      		lpts := lp(space);
-      		indexList := concat concat indices;
-		put!(sink, "Points");
-		putPoints!(sink, transform, lpts, indexList);
-      		offset : SI := 0;
-		lstart!(sink);
-		for ss in components(space) for index in indices repeat {
-			closedCurve? ss => {
-				put!(sink, "closedCurve");
-				n: SI := coerce #(first index);
-				put!(sink, offset);
-				put!(sink, n);
-				offset := offset + n;
-			}
-			curve? ss=> {
-				put!(sink, "curve");
-				n: SI := coerce #(first index);
-				put!(sink, offset);
-				put!(sink, n);
-				offset := offset + n;
-			}
-			polygon? ss => {
-				local vertices: SI;
-				put!(sink, "polygon");
-				vertices := coerce(#(first index)
-						   + #(first rest index));
-				put!(sink, offset);
-				put!(sink, vertices);
-				offset := offset+vertices;
-			}
-			mesh? ss=> {
-				local xStep, yStep: SI;
-				put!(sink, "mesh");
-				xStep := coerce #index;
-          			yStep := coerce #(first index);
-				put!(sink, offset);
-				put!(sink, xStep);
-				put!(sink, yStep);
-				offset := offset+xStep*yStep;
-			}
-			point? ss => {
-				put!(sink, "points");
-				put!(sink, offset);
-				put!(sink, 1);
-				offset := offset+1;
-			}
-			error "Unrecognised SubSpace component";
-		}
-		lend!(sink);
-	}
-	
-}
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
-<<InventorDataSink>>
-<<InventorViewPort>>
-<<InventorRenderPackage>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/mlift.spad.jhd.pamphlet b/src/algebra/mlift.spad.jhd.pamphlet
deleted file mode 100644
index 637d1b3..0000000
--- a/src/algebra/mlift.spad.jhd.pamphlet
+++ /dev/null
@@ -1,272 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra mlift.spad.jhd}
-\author{Patrizia Gianni}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{package MLIFT MultivariateLifting}
-<<package MLIFT MultivariateLifting>>=
-)abbrev package MLIFT MultivariateLifting
-++ Author : P.Gianni.
-++ Description: 
-++ This package provides the functions for the multivariate "lifting", using
-++ an algorithm of Paul Wang.
-++ This package will work for every euclidean domain R which has property
-++ F, i.e. there exists a factor operation in \spad{R[x]}.
- 
-MultivariateLifting(E,OV,R,P) : C == T
- where
-  OV        :   OrderedSet
-  E         :   OrderedAbelianMonoidSup
-  R         :   EuclideanDomain  -- with property "F"
-  Z         ==> Integer
-  BP        ==> SparseUnivariatePolynomial R
-  P         :   PolynomialCategory(R,E,OV)
-  SUP       ==> SparseUnivariatePolynomial P
-  NNI       ==> NonNegativeInteger
-  Term      ==> Record(expt:NNI,pcoef:P)
-  VTerm     ==> List Term
-  Table     ==> Vector List BP
-  L         ==> List
- 
-  C == with
-    corrPoly:      (SUP,L OV,L R,L NNI,L SUP,Table,R) -> Union(L SUP,"failed")
-    lifting:       (SUP,L OV,L BP,L R,L P,L NNI,R) -> Union(L SUP,"failed")
-    lifting1:      (SUP,L OV,L SUP,L R,L P,L VTerm,L NNI,Table,R) ->
-                                                Union(L SUP,"failed")
- 
-  T == add
-    GenExEuclid(R,BP)
-    NPCoef(BP,E,OV,R,P)
-    IntegerCombinatoricFunctions(Z)
-
-    SUPF2   ==> SparseUnivariatePolynomialFunctions2
-
-    DetCoef ==> Record(deter:L SUP,dterm:L VTerm,nfacts:L BP,
-                       nlead:L P)
- 
-              ---   local functions   ---
-    normalDerivM    :    (P,Z,OV)     ->  P
-    normalDeriv     :     (SUP,Z)     ->  SUP
-    subslead        :     (SUP,P)     ->  SUP
-    subscoef        :   (SUP,L Term)  ->  SUP
-    maxDegree       :   (SUP,OV)      ->  NonNegativeInteger
- 
- 
-    corrPoly(m:SUP,lvar:L OV,fval:L R,ld:L NNI,flist:L SUP,
-             table:Table,pmod:R):Union(L SUP,"failed") ==
-      --  The correction coefficients are evaluated recursively.
-      --   Extended Euclidean algorithm for the multivariate case.
- 
-      -- the polynomial is univariate  --
-      #lvar=0 =>
-        lp:=solveid(map(ground,m)$SUPF2(P,R),pmod,table)
-        if lp case "failed" then return "failed"
-        lcoef:= [map(coerce,mp)$SUPF2(R,P) for mp in lp::L BP]
- 
- 
-      diff,ddiff,pol,polc:SUP
-      listpolv,listcong:L SUP
-      deg1:NNI:= ld.first
-      np:NNI:= #flist
-      a:P:= fval.first ::P
-      y:OV:=lvar.first
-      lvar:=lvar.rest
-      listpolv:L SUP := [map(eval(#1,y,a),f1) for f1 in flist]
-      um:=map(eval(#1,y,a),m)
-      flcoef:=corrPoly(um,lvar,fval.rest,ld.rest,listpolv,table,pmod)
-      if flcoef case "failed" then return "failed"
-      else lcoef:=flcoef :: L SUP
-      listcong:=[*/[flist.i for i in 1..np | i^=l] for l in 1..np]
-      polc:SUP:= (monomial(1,y,1) - a)::SUP
-      pol := 1$SUP
-      diff:=m- +/[lcoef.i*listcong.i for i in 1..np]
-      for l in 1..deg1 repeat
-        if diff=0 then return lcoef
-        pol := pol*polc
-        (ddiff:= map(eval(normalDerivM(#1,l,y),y,a),diff)) = 0 => "next l"
-        fbeta := corrPoly(ddiff,lvar,fval.rest,ld.rest,listpolv,table,pmod)
-        if fbeta case "failed" then return "failed"
-        else beta:=fbeta :: L SUP
-        lcoef := [lcoef.i+beta.i*pol  for i in 1..np]
-        diff:=diff- +/[listcong.i*beta.i for i in 1..np]*pol
-      lcoef
- 
- 
- 
-    lifting1(m:SUP,lvar:L OV,plist:L SUP,vlist:L R,tlist:L P,_
-      coeflist:L VTerm,listdeg:L NNI,table:Table,pmod:R) :Union(L SUP,"failed") ==
-    -- The factors of m (multivariate) are determined ,
-    -- We suppose to know the true univariate factors
-    -- some coefficients are determined
-      conglist:L SUP:=empty()
-      nvar : NNI:= #lvar
-      pol,polc:P
-      mc,mj:SUP
-      testp:Boolean:= (not empty?(tlist))
-      lalpha : L SUP := empty()
-      tlv:L P:=empty()
-      subsvar:L OV:=empty()
-      subsval:L R:=empty()
-      li:L OV := lvar
-      ldeg:L NNI:=empty()
-      clv:L VTerm:=empty()
-      --j =#variables, i=#factors
-      for j in 1..nvar repeat
-        x  := li.first
-        li := rest li
-        conglist:= plist
-        v := vlist.first
-        vlist := rest vlist
-        degj := listdeg.j
-        ldeg := cons(degj,ldeg)
-        subsvar:=cons(x,subsvar)
-        subsval:=cons(v,subsval)
- 
-      --substitute the determined coefficients
-        if testp then
-          if j<nvar then
-            tlv:=[eval(p,li,vlist) for p in tlist]
-            clv:=[[[term.expt,eval(term.pcoef,li,vlist)]$Term
-                   for term in clist] for clist  in coeflist]
-          else (tlv,clv):=(tlist,coeflist)
-          plist :=[subslead(p,lcp) for p in plist for lcp in tlv]
-          if not(empty? coeflist) then
-            plist:=[subscoef(tpol,clist)
-                   for tpol in plist for clist in clv]
-        mj := map(eval(#1,li,vlist),m)  --m(x1,..,xj,aj+1,..,an
-        polc := x::P - v::P  --(xj-aj)
-        pol:= 1$P
-      --Construction of Rik, k in 1..right degree for xj+1
-        for k in 1..degj repeat  --I can exit before
-         pol := pol*polc
-         mc := */[term for term in plist]-mj
-         if mc=0 then leave "next var"
-         --Modulus Dk
-         mc:=map(normalDerivM(#1,k,x),mc)
-         (mc := map(eval(#1,[x],[v]),mc))=0 => "next k"
-         flalpha:=corrPoly(mc,subsvar.rest,subsval.rest,
-                          ldeg.rest,conglist,table,pmod)
-         if flalpha case "failed" then return "failed"
-         else lalpha:=flalpha :: L SUP
-         plist:=[term-alpha*pol for term in plist for alpha in lalpha]
-        for term in plist repeat degj:=degj-maxDegree(term,x)
-        degj ^= 0 => return "failed"
-      plist
-        --There are not extraneous factors
- 
-    maxDegree(um:SUP,x:OV):NonNegativeInteger ==
-       ans:NonNegativeInteger:=0
-       while um ^= 0 repeat
-          ans:=max(ans,degree(leadingCoefficient um,x))
-          um:=reductum um
-       ans
-
-    lifting(um:SUP,lvar:L OV,plist:L BP,vlist:L R,
-            tlist:L P,listdeg:L NNI,pmod:R):Union(L SUP,"failed") ==
-    -- The factors of m (multivariate) are determined, when the
-    --  univariate true factors are known and some coefficient determined
-      nplist:List SUP:=[map(coerce,pp)$SUPF2(R,P) for pp in plist]
-      empty? tlist =>
-        table:=tablePow(degree um,pmod,plist)
-        table case "failed" => error "Table construction failed in MLIFT"
-        lifting1(um,lvar,nplist,vlist,tlist,empty(),listdeg,table,pmod)
-      ldcoef:DetCoef:=npcoef(um,plist,tlist)
-      if not empty?(listdet:=ldcoef.deter) then
-        if #listdet = #plist  then return listdet
-        plist:=ldcoef.nfacts
-        nplist:=[map(coerce,pp)$SUPF2(R,P) for pp in plist]
-        um:=(um exquo */[pol for pol in listdet])::SUP
-        tlist:=ldcoef.nlead
-        tab:=tablePow(degree um,pmod,plist.rest)
-      else tab:=tablePow(degree um,pmod,plist)
-      tab case "failed" => error "Table construction failed in MLIFT"
-      table:Table:=tab
-      ffl:=lifting1(um,lvar,nplist,vlist,tlist,ldcoef.dterm,listdeg,table,pmod)
-      if ffl case "failed" then return "failed"
-      append(listdet,ffl:: L SUP)
-
-    -- normalDerivM(f,m,x) = the normalized (divided by m!) m-th
-    -- derivative with respect to x of the multivariate polynomial f
-    normalDerivM(g:P,m:Z,x:OV) : P ==
-     multivariate(normalDeriv(univariate(g,x),m),x)
-
-    normalDeriv(f:SUP,m:Z) : SUP ==
-     (n1:Z:=degree f) < m => 0$SUP
-     n1=m => leadingCoefficient f :: SUP
-     k:=binomial(n1,m)
-     ris:SUP:=0$SUP
-     n:Z:=n1
-     while n>= m repeat
-       while n1>n repeat
-         k:=(k*(n1-m)) quo n1
-         n1:=n1-1
-       ris:=ris+monomial(k*leadingCoefficient f,(n-m)::NNI)
-       f:=reductum f
-       n:=degree f
-     ris
-
-    subslead(m:SUP,pol:P):SUP ==
-      dm:NNI:=degree m
-      monomial(pol,dm)+reductum m
- 
-    subscoef(um:SUP,lterm:L Term):SUP ==
-      dm:NNI:=degree um
-      new:=monomial(leadingCoefficient um,dm)
-      for k in dm-1..0 by -1 repeat
-        i:NNI:=k::NNI
-        empty?(lterm) or lterm.first.expt^=i =>
-                                new:=new+monomial(coefficient(um,i),i)
-        new:=new+monomial(lterm.first.pcoef,i)
-        lterm:=lterm.rest
-      new
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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 MLIFT MultivariateLifting>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/ndftip.as.pamphlet b/src/algebra/ndftip.as.pamphlet
deleted file mode 100644
index 3c9c0ea..0000000
--- a/src/algebra/ndftip.as.pamphlet
+++ /dev/null
@@ -1,1174 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra ndftip.as}
-\author{Michael Richardson}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{NagDiscreteFourierTransformInterfacePackage}
-<<NagDiscreteFourierTransformInterfacePackage>>=
-+++ Author: M.G. Richardson
-+++ Date Created: 1995 Dec. 08
-+++ Date Last Updated:
-+++ Basic Functions:
-+++ Related Constructors:
-+++ Also See:
-+++ AMS Classifications:
-+++ Keywords:
-+++ References:
-+++ Description:
-+++ This package provides Axiom-like interfaces to the NAG
-+++ Finite Fourier Transform routines in the NAGlink.
-
-NagDiscreteFourierTransformInterfacePackage: with {
-
-  nagDFT : VDF -> VCDF ;                                --  test  1
-
-++ nagDFT(seq) calculates the discrete Fourier transform of a sequence
-++ of real data values 
-#if saturn
-++ $x_{1} \ldots x_{n}$
-#else
-++ \spad{x[1] .. x[n]}
-#endif
-++ supplied in the vector seq.
-++ Note that the definition used for the discrete Fourier transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{-i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(x[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06EAF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06eaf.
-
-  nagDFT : VCDF -> VCDF ;                               --  test  3
-
-++ nagDFT(seq) calculates the discrete Fourier transform of a sequence
-++ of complex data values 
-#if saturn
-++ $z_{1} \ldots z_{n}$
-#else
-++ \spad{z[1] .. z[n]}
-#endif
-++ supplied in the vector seq.
-++ Note that the definition used for the discrete Fourier transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{-i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(z[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06ECF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06ecf.
-
-  nagDFT : PHSDF -> VDF ;                               --  test  7
-
-++ nagDFT(seq) calculates the discrete Fourier transform of a Hermitian
-++ sequence of complex data values,
-#if saturn
-++ $z_{1} \ldots z_{n}$
-#else
-++ \spad{z[1] .. z[n]}
-#endif
-++ supplied in the PackedHermitianSequence seq.
-++ Note that the definition used for the discrete Fourier transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{-i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(z[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06EBF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06ebf.
-
-  nagDFT : LVDF -> LVCDF ;                               --  test 10, 19
-
-++ nagDFT(seqs) calculates the discrete Fourier transform of each of a
-++ list of sequences of real data values 
-#if saturn
-++ $x_{1} \ldots x_{n}$
-#else
-++ \spad{x[1] .. x[n]}
-#endif
-++ supplied in the list of vectors, seqs.
-++ Note that the definition used for the discrete Fourier transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{-i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(x[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06FPF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06fpf.
-
-  nagDFT : LVCDF -> LVCDF ;                             --  test 16
-
-++ nagDFT(seqs) calculates the discrete Fourier transform of each of a
-++ list of sequences of complex data values 
-#if saturn
-++ $z_{1} \ldots z_{n}$
-#else
-++ \spad{z[1] .. z[n]}
-#endif
-++ supplied in the list of vectors, seqs.
-++ Note that the definition used for the discrete Fourier transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{-i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(z[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06FRF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06frf.
-
-  nagDFT : LPHSDF -> LVDF ;                              --  test 12, 21
-
-++ nagDFT(seq) calculates the discrete Fourier transform of a each of a
-++ list of Hermitian sequences of complex data values,
-#if saturn
-++ $z_{1} \ldots z_{n}$
-#else
-++ \spad{z[1] .. z[n]}
-#endif
-++ supplied in the List PackedHermitianSequence, seq.
-++ Note that the definition used for the discrete Fourier transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{-i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(z[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06FQF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06fqf.
-
-  nagInverseDFT : VDF -> VCDF ;                         --  test  8
-
-++ nagInverseDFT(seq) calculates the inverse discrete Fourier
-++ transform of a sequence of real data values 
-#if saturn
-++ $x_{1} \ldots x_{n}$
-#else
-++ \spad{x[1] .. x[n]}
-#endif
-++ supplied in the vector seq.
-++ Note that the definition used for the inverse discrete Fourier
-++ transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(x[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06EAF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06eaf.
-
-  nagInverseDFT : VCDF -> VCDF ;                         --  test  2,  4
-
-++ nagInverseDFT(seq) calculates the inverse discrete Fourier
-++ transform of a sequence of complex data values 
-#if saturn
-++ $z_{1} \ldots z_{n}$
-#else
-++ \spad{z[1] .. z[n]}
-#endif
-++ supplied in the vector seq.
-++ Note that the definition used for the inverse discrete Fourier
-++ transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(z[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06ECF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06ecf.
-
-  nagInverseDFT : PHSDF -> VDF ;                        --  test  6
-
-++ nagInverseDFT(seq) calculates the inverse discrete Fourier transform
-++ of a Hermitian sequence of complex data values 
-#if saturn
-++ $z_{1} \ldots z_{n}$
-#else
-++ \spad{z[1] .. z[n]}
-#endif
-++ supplied in the PackedHermitianSequence seq.
-++ Note that the definition used for the inverse discrete Fourier
-++ transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(z[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06EBF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06ebf.
-
-  nagInverseDFT : LVDF -> LVCDF ;                       --  test 13
-
-++ nagInverseDFT(seqs) calculates the inverse discrete Fourier
-++ transform of each of a list of sequences of real data values 
-#if saturn
-++ $x_{1} \ldots x_{n}$
-#else
-++ \spad{x[1] .. x[n]}
-#endif
-++ supplied in the list of vectors, seqs.
-++ Note that the definition used for the inverse discrete Fourier
-++ transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(x[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06FPF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06fpf.
-
-  nagInverseDFT : LVCDF -> LVCDF ;                       --  test 11, 17
-
-++ nagInverseDFT(seqs) calculates the inverse discrete Fourier
-++ transform of each of a list of sequences of complex data values 
-#if saturn
-++ $z_{1} \ldots z_{n}$
-#else
-++ \spad{z[1] .. z[n]}
-#endif
-++ supplied in the list of vectors, seqs.
-++ Note that the definition used for the inverse discrete Fourier
-++ transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(z[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06FRF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06frf.
-
-  nagInverseDFT : LPHSDF -> LVDF ;                      --  test 15
-
-++ nagInverseDFT(seqs) calculates the inverse discrete Fourier transform
-++ of each of a list of Hermitian sequences of complex data values 
-#if saturn
-++ $z_{1} \ldots z_{n}$
-#else
-++ \spad{z[1] .. z[n]}
-#endif
-++ supplied in the List PackedHermitianSequence, seqs.
-++ Note that the definition used for the inverse discrete Fourier
-++ transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} z_{j} e^{i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(z[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06FQF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06fqf.
-
-  nagHermitianDFT : VDF -> PHSDF ;                      --  test  5
-
-++ nagHermitianDFT(seq) calculates the discrete Fourier transform, in
-++ packed Hermitian form, of a sequence of real data values 
-#if saturn
-++ $x_{1} \ldots x_{n}$
-#else
-++ \spad{x[1] .. x[n]}
-#endif
-++ supplied in the vector seq.
-++ Note that the definition used for the discrete Fourier transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{-i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(x[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06EAF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06eaf.
-
-  nagHermitianDFT : LVDF -> LPHSDF ;                     --  test 14, 20
-
-++ nagHermitianDFT(seqs) calculates the discrete Fourier transform, in
-++ packed Hermitian form, of each of a list of sequences of real data
-++ values 
-#if saturn
-++ $x_{1} \ldots x_{n}$
-#else
-++ \spad{x[1] .. x[n]}
-#endif
-++ supplied in the list of vectors, seqs.
-++ Note that the definition used for the discrete Fourier transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{-i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(x[j]*%e^(-i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06FPF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06fpf.
-
-  nagHermitianInverseDFT : VDF -> PHSDF ;               --  test  9
-
-++ nagHermitianInverseDFT(seq) calculates the inverse discrete Fourier
-++ transform, in packed Hermitian form, of a sequence of real data
-++ values 
-#if saturn
-++ $x_{1} \ldots x_{n}$
-#else
-++ \spad{x[1] .. x[n]}
-#endif
-++ supplied in the vector seq.
-++ Note that the definition used for the inverse discrete Fourier
-++ transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(x[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06EAF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06eaf.
-
-  nagHermitianInverseDFT : LVDF -> LPHSDF ;             --  test 18
-
-++ nagHermitianInverseDFT(seqs) calculates the inverse discrete Fourier
-++ transform, in packed Hermitian form, of each of a list of sequences
-++ of real data values 
-#if saturn
-++ $x_{1} \ldots x_{n}$
-#else
-++ \spad{x[1] .. x[n]}
-#endif
-++ supplied in the list of vectors, seqs.
-++ Note that the definition used for the inverse discrete Fourier
-++ transform is
-#if saturn
-++ \[ \frac{1}{\sqrt{n} \sum_{j=0}^{n-1} x_{j} e^{i \frac{2 \pi j k}{n}
-++ \qquad k = 0 \ldots  n - 1 \]
-#else
-++ \spad{1/sqrt(n)*sum(x[j]*%e^(i*2*%pi*j*k/n), j=0..(n-1))} for
-++ \spad{k=0..(n-1)}.
-#endif
-++ The numerical calculation is performed by the NAG routine C06FPF.
-++
-++ For more detailed information, please consult the NAG
-++ manual via the Browser page for the operation c06fpf.
-
-} == add {
-
-  import from AnyFunctions1 MDF ;
-  import from CDF;
-  import from ErrorFunctions ;
-  import from LLDF ;
-  import from MCDF ;
-  import from MDF ;
-  import from NagResultChecks ;
-  import from NagSeriesSummationPackage ;
-  import from PHSDF;
-  import from STRG ;
-  import from List STRG ;
-  import from Symbol ;
-  import from VDF ;
-
-  local (..)(a:INT,b:INT):Generator INT == {
-    generate {
-      t := a ;
-      while (t <= b) repeat {
-        yield t ;
-        t := t + 1 ;
-        }
-      }
-    }
-
-  local ipIfail : INT := -1 ;
-
--- First, the functions corresponding to single NAGlink calls of C06E
--- routines (single vector transforms):
-
--- c06eaf:
-
-  nagHermitianDFT(seq : VDF) : PHSDF ; == {
-    local lseq : INT ;
-
-    lseq := ((# seq)@NNI) pretend INT ; -- @ to eliminate SI possibility
-    row(checkMxDF(c06eaf(lseq,matrix [members seq],ipIfail),
-		   "x",
-		    "C06EAF"),
-			      1)
-				 pretend PHSDF
-  }
-
--- c06ebf:
-
-  nagDFT(seq : PHSDF) : VDF == {
-    local lseq : INT ;
-
-    lseq := ((# seq)@NNI) pretend INT ; -- @ to eliminate SI possibility
-    row(checkMxDF(c06ebf(lseq,matrix [members seq],ipIfail),
-		   "x",
-		    "C06EBF"),
-			      1)
-  }
-
--- c06ecf:
-
-  nagDFT(seq : VCDF) : VCDF == {
-    local nseq : NNI ;
-    local lseq : INT ;
-    local rvec, ivec : VDF ;
-    local cvec : VCDF ;
-    local c06ecfResult : RSLT ;
-
-    nseq := # seq ;
-    lseq := nseq pretend INT ;
-    rvec := new(nseq,0) ;
-    ivec := new(nseq,0) ;
-    for i in 1..lseq repeat {
-      rvec(i) := real seq(i) ;
-      ivec(i) := imag seq(i) ;
-    }
-    c06ecfResult := c06ecf(lseq,
-			    matrix [members rvec],
-			     matrix [members ivec],
-			      ipIfail) ;
-    rvec := row(checkMxDF(c06ecfResult,"x","C06ECF"),1) ;
-    ivec := row((retract(c06ecfResult."y") @ MDF),1) ;
-    cvec := new(nseq,0) ;
-    for i in 1..lseq repeat cvec(i) := complex(rvec(i),ivec(i)) ;
-    cvec
-  }
-
--- inverse transforms, in terms of these and functions from PHS:
-
-  nagInverseDFT(seq : PHSDF) : VDF == nagDFT conjHerm seq ;
-
-  nagHermitianInverseDFT(seq : VDF) : PHSDF
-			== conjHerm nagHermitianDFT seq ;
-
-  nagInverseDFT(seq : VCDF) : VCDF == {
-    local nseq : NNI ;
-    local lseq : INT ;
-    local rvec, ivec : VDF ;
-    local cvec : VCDF ;
-    local c06ecfResult : RSLT ;
-
-    nseq := # seq ;
-    lseq := nseq pretend INT ;
-    rvec := new(nseq,0) ;
-    ivec := new(nseq,0) ;
-    for i in 1..lseq repeat {
-      rvec(i) := real seq(i) ;
-      ivec(i) := - imag seq(i) ;
-    }
-    c06ecfResult := c06ecf(lseq,
-			    matrix [members rvec],
-			     matrix [members ivec],
-			      ipIfail) ;
-    rvec := row(checkMxDF(c06ecfResult,"x","C06ECF"),1) ;
-    ivec := row((retract(c06ecfResult."y") @ MDF),1) ;
-    cvec := new(nseq,0) ;
-    for i in 1..lseq repeat cvec(i) := complex(rvec(i), - ivec(i)) ;
-    cvec
-  }
-
--- "Full form" equivalents of c06eaf and inverse:
-
-  nagDFT(seq : VDF) : VCDF == expand nagHermitianDFT seq ;
-
-  nagInverseDFT(seq : VDF) : VCDF == expand nagHermitianInverseDFT seq ;
-
-
--- Next, the functions corresponding to single NAGlink calls of C06F
--- routines (multiple vector transforms):
-
--- basic routines:
-
--- c06fpf
-
-  nagHermitianDFT(seqs : LVDF) : LPHSDF ; == {
-
-    local nr, nc : NNI ;
-    local inr, inc : INT ;
-    local seqMat, trig, result : MDF ;
-    local nextSeq : PHSDF ;
-    local hermDFTs : LPHSDF ;
-
-    nr := # seqs ;
-    inr := nr pretend INT ;
-    nc := # (seqs.1) ;
-    inc := nc pretend INT ;
-    seqMat := new(nr,nc,0) ;
-    for j in 1 .. inc repeat seqMat(1,j) := (seqs.1).j ;
-    for i in 2 .. inr repeat
-      if (# seqs.i) ~= nc 
-      then error ["The data sequences in nagHermitianDFT must all",
-		  " have the same length.  ",
-		  "The length of sequence 1 is ",
-		  string(inc),
-		  "that of sequence ",
-		  string(i pretend INT),
-		  " is ",
-		  string((# seqs.i)@NNI pretend INT), -- @ avoids SI
-		  "."]
-      else for j in 1 .. inc repeat seqMat(i,j) := (seqs.i).j ;
-    trig := new(1@NNI,2*nc,0) ;
-    result :=
-      checkMxDF(c06fpf(inr,inc,"i",seqMat,trig,ipIfail),"x","C06FPF") ;
-    hermDFTs := [] ;
-    for i in inr .. 1  by -1 repeat {
-      nextSeq := new(nc,0) ;
-      for j in 1 .. inc repeat nextSeq(j) := result(1,(j-1)*inr + i) ;
-      hermDFTs := cons(nextSeq,hermDFTs) ;
-    }
-    hermDFTs
-  }
-
--- c06fqf
-
-  nagDFT(seqs : LPHSDF) : LVDF == {
-    
-    local nr, nc : NNI ;
-    local inr, inc : INT ;
-    local seqMat, trig, result : MDF ;
-    local nextSeq : VDF ;
-    local dfts : LVDF ;
-  
-    nr := # seqs ;
-    inr := nr pretend INT ;
-    nc := # (seqs.1) ;
-    inc := nc pretend INT ;
-    seqMat := new(nr,nc,0) ;
-    for j in 1 .. inc repeat seqMat(1,j) := (seqs.1).j ;
-    for i in 2 .. inr repeat
-      if (# seqs.i) ~= nc 
-      then error ["The data sequences in nagDFT must all",
-                  " have the same length.  ",
-                  "The length of sequence 1 is ",
-                  string(inc),
-                  "that of sequence ",
-                  string(i pretend INT),
-                  " is ",
-                  string((# seqs.i)@NNI pretend INT), -- @ avoids SI
-                  "."]
-      else for j in 1 .. inc repeat seqMat(i,j) := (seqs.i).j ;
-    trig := new(1@NNI,2*nc,0) ;
-    result :=
-      checkMxDF(c06fqf(inr,inc,"i",seqMat,trig,ipIfail),"x","C06FQF") ;
-    dfts := [] ;
-    for i in inr .. 1 by -1 repeat {
-      nextSeq := new(nc,0) ;
-      for j in 1 .. inc repeat nextSeq(j) := result(1,(j-1)*inr + i) ;
-      dfts := cons(nextSeq,dfts) ;
-    }
-    dfts
-  }
-
--- c06frf
-
-  nagDFT(seqs : LVCDF) : LVCDF == {
-
-    local nr, nc : NNI ;
-    local inr, inc : INT ;
-    local trig, rMat, iMat : MDF ;
-    local result : RSLT ;
-    local nextSeq : VCDF ;
-    local dfts : LVCDF ;
-
-    nr := # seqs ;
-    inr := nr pretend INT ;
-    nc := # (seqs.1) ;
-    inc := nc pretend INT ;
-    rMat := new(nr,nc,0) ;
-    iMat := new(nr,nc,0) ;
-    for j in 1 .. inc repeat {
-      rMat(1,j) := real((seqs.1).j) ;
-      iMat(1,j) := imag((seqs.1).j) ;
-    }
-    for i in 2 .. inr repeat {
-      if (# seqs.i) ~= nc
-      then error ["The data sequences in nagDFT must all",
-		  " have the same length.  ",
-		  "The length of sequence 1 is ",
-		  string(inc),
-		  "that of sequence ",
-		  string(i pretend INT),
-		  " is ",
-		  string((# seqs.i)@NNI pretend INT), -- @ avoids SI
-		  "."]
-      else for j in 1 .. inc repeat {
-	rMat(i,j) := real((seqs.i).j) ;
-	iMat(i,j) := imag((seqs.i).j) ;
-      }
-    }
-    trig := new(1@NNI,2*nc,0) ;
-    result := c06frf(inr,inc,"i",rMat,iMat,trig,ipIfail) ;
-    rMat := checkMxDF(result, "x", "C06FRF") ;
-    iMat := retract(result."y") @ MDF ;
-    dfts := [] ;
-    for i in inr .. 1 by -1 repeat {
-      nextSeq := new(nc,0) ;
-      for j in 1 .. inc repeat
-	nextSeq(j) := complex(rMat(1,(j-1)*inr+i),iMat(1,(j-1)*inr+i)) ;
-      dfts := cons(nextSeq,dfts) ;
-    }
-    dfts
-  }
-
--- inverse transforms, in terms of these and functions from PHS:
-
-  nagInverseDFT(seqs : LVCDF) : LVCDF == {
-
-    local nr, nc : NNI ;
-    local inr, inc : INT ;
-    local conjSeq : VCDF ;
-    local temp, invdfts : LVCDF ;
-
-    nr := # seqs ;
-    inr := nr pretend INT ;
-    temp := [] ;
-    for i in inr .. 1 by -1 repeat {
-      nc := #(seqs.i) ;
-      inc := nc pretend INT ;
-      conjSeq := new(nc,0) ;
-      for j in 1 .. inc repeat
-        conjSeq(j) := conjugate((seqs.i).j) ;
-      temp := cons(conjSeq,temp) ;
-    }
-    temp := nagDFT temp ;
-    invdfts := [] ;
-    for i in inr .. 1 by -1 repeat {
-      conjSeq := new(nc,0) ;
-      for j in 1 .. inc repeat -- know inc is constant after nagDFT call
-	conjSeq(j) := conjugate((temp.i).j) ;
-      invdfts := cons(conjSeq,invdfts) ;
-    }
-    invdfts
-  }
-
-  nagInverseDFT(seqs : LPHSDF) : LVDF == {
-    local nr : NNI ;
-    local inr : INT ;
-    local conjSeqs : LPHSDF ;
-  
-    nr := # seqs ;
-    inr := nr pretend INT ;
-    conjSeqs := [] ;
-    for i in inr .. 1 by -1 repeat
-      conjSeqs := cons(conjHerm(seqs.i),conjSeqs) ;
-    nagDFT conjSeqs ;
-  }
-
-  nagHermitianInverseDFT(seqs : LVDF) : LPHSDF == {
-    local nr : NNI ;
-    local inr : INT ;
-    local conjSeqs, invSeqs : LPHSDF ;
-  
-    nr := # seqs ;
-    inr := nr pretend INT ;
-    conjSeqs := nagHermitianDFT seqs ;
-    invSeqs := [] ;
-    for i in inr .. 1 by -1 repeat
-      invSeqs := cons(conjHerm(conjSeqs.i),invSeqs) ;
-    invSeqs
-  }
-
--- "Full form" equivalents of c06fpf and inverse:
-
-  nagDFT(seqs : LVDF) : LVCDF == {
-  
-    local nr : NNI ;
-    local inr : INT ;
-    local hermdfts : LPHSDF ;
-    local dfts : LVCDF ;
-
-    nr := # seqs ;
-    inr := nr pretend INT ;
-    hermdfts := nagHermitianDFT seqs ;
-    dfts := [] ;
-    for i in inr .. 1 by -1 repeat
-      dfts := cons(expand(hermdfts.i),dfts) ;
-    dfts
-  }
-
-  nagInverseDFT(seqs : LVDF) : LVCDF == {
-    local nr : NNI ;
-    local inr : INT ;
-    local hermdfts : LPHSDF ;
-    local invdfts : LVCDF ;
-
-    nr := # seqs ;
-    inr := nr pretend INT ;
-    hermdfts := nagHermitianDFT seqs ;
-    invdfts := [] ;
-    for i in inr .. 1 by -1 repeat
-      invdfts := cons(expand conjHerm(hermdfts.i),invdfts) ;
-    invdfts
-  }
-    
-}
-
-#if NeverAssertThis
-
--- Note that the conversions of results from DoubleFloat to Float
--- will become unnecessary if outputGeneral is extended to apply to
--- DoubleFloat quantities.  Those results not converted will, of
--- course, then be displayed to 6 s.f.
-
-)lib nrc
-)lib herm
-)lib ndftip
-
-outputGeneral 6
-
-seqA := [0.34907,0.54890,0.74776,0.94459,1.1385,1.3285,1.5137];
-
-seqB := [0.34907 - 0.37168*%i,  _
-         0.54890 - 0.35669*%i,  _
-         0.74776 - 0.31175*%i,  _
-         0.94459 - 0.23702*%i,  _
-         1.13850 - 0.13274*%i,  _
-         1.32850 + 0.00074*%i,  _
-         1.51370 + 0.16298*%i];
-
-hseqC : PackedHermitianSequence DoubleFloat
-hseqC := packHS [0.34907,        _
-                 0.54890 + %i*1.51370,  _
-                 0.74776 + %i*1.32850,  _
-                 0.94459 + %i*1.13850,  _
-                 0.94459 - %i*1.13850,  _
-                 0.74776 - %i*1.32850,  _
-                 0.54890 - %i*1.51370];
-
-seqsD : List Vector DoubleFloat;
-seqsD := [vector [0.3854, 0.6772, 0.1138, 0.6751, 0.6362, 0.1424], _
-          vector [0.5417, 0.2983, 0.1181, 0.7255, 0.8638, 0.8723], _
-          vector [0.9172, 0.0644, 0.6037, 0.6430, 0.0428, 0.4815]];
-
-seqsE : List PackedHermitianSequence DoubleFloat;
-seqsE := [pHS [0.3854, 0.6772, 0.1138, 0.6751, 0.6362, 0.1424], _
-          pHS [0.5417, 0.2983, 0.1181, 0.7255, 0.8638, 0.8723], _
-          pHS [0.9172, 0.0644, 0.6037, 0.6430, 0.0428, 0.4815]];
-
-seqsF : List Vector Complex DoubleFloat
-seqsF := [vector [0.3854 + 0.5417*%i, 0.6772 + 0.2983*%i,   _
-                  0.1138 + 0.1181*%i, 0.6751 + 0.7255*%i,   _
-                  0.6362 + 0.8638*%i, 0.1424 + 0.8723*%i],  _
-          vector [0.9172 + 0.9089*%i, 0.0644 + 0.3118*%i,   _
-                  0.6037 + 0.3465*%i, 0.6430 + 0.6198*%i,   _
-                  0.0428 + 0.2668*%i, 0.4815 + 0.1614*%i],  _
-          vector [0.1156 + 0.6214*%i, 0.0685 + 0.8681*%i,   _
-                  0.2060 + 0.7060*%i, 0.8630 + 0.8652*%i,   _
-                  0.6967 + 0.9190*%i, 0.2792 + 0.3355*%i]];
-
--- test  1
-
-dftA := nagDFT seqA;
-dftA :: Vector Complex Float :: Matrix Complex Float
-                             -- Matrix to force display as a column,
-                             -- Float to allow outputGeneral to work.
-
---       +         2.48361         +
---       |                         |
---       |- 0.265985 + 0.530898 %i |
---       |                         |
---       |- 0.257682 + 0.202979 %i |
---       |                         |
---       |- 0.256363 + 0.0580623 %i|
---       |                         |
---       |- 0.256363 - 0.0580623 %i|
---       |                         |
---       |- 0.257682 - 0.202979 %i |
---       |                         |
---       +- 0.265985 - 0.530898 %i +
-
--- test  2
-
-nagInverseDFT dftA :: Vector Float
- 
---       [0.34907,0.5489,0.74776,0.94459,1.1385,1.3285,1.5137]
-
--- test  3
-
-dftB := nagDFT seqB;
-dftB :: Vector Complex Float :: Matrix Complex Float
-
---       +  2.48361 - 0.471004 %i  +
---       |                         |
---       | - 0.5518 + 0.496841 %i  |
---       |                         |
---       |- 0.367113 + 0.0975621 %i|
---       |                         |
---       |- 0.287669 - 0.0586476 %i|
---       |                         |
---       |- 0.225057 - 0.174772 %i |
---       |                         |
---       |- 0.148251 - 0.308396 %i |
---       |                         |
---       + 0.0198297 - 0.564956 %i +
- 
--- test  4
- 
-(nagInverseDFT dftB) :: Vector Complex Float :: Matrix Complex Float
- 
---       +0.34907 - 0.37168 %i+
---       |                    |
---       |0.5489 - 0.35669 %i |
---       |                    |
---       |0.74776 - 0.31175 %i|
---       |                    |
---       |0.94459 - 0.23702 %i|
---       |                    |
---       |1.1385 - 0.13274 %i |
---       |                    |
---       |1.3285 + 0.00074 %i |
---       |                    |
---       +1.5137 + 0.16298 %i +
-
--- test  5
-
-hdftA := nagHermitianDFT seqA;
-(expand hdftA) :: Vector Complex Float :: Matrix Complex Float
-
---       +         2.48361         +
---       |                         |
---       |- 0.265985 + 0.530898 %i |
---       |                         |
---       |- 0.257682 + 0.202979 %i |
---       |                         |
---       |- 0.256363 + 0.0580623 %i|
---       |                         |
---       |- 0.256363 - 0.0580623 %i|
---       |                         |
---       |- 0.257682 - 0.202979 %i |
---       |                         |
---       +- 0.265985 - 0.530898 %i +
- 
--- test  6
- 
-(nagInverseDFT hdftA) :: Vector Float
-
---       [0.34907,0.5489,0.74776,0.94459,1.1385,1.3285,1.5137]
-
--- test  7
-
-dftC := nagDFT hseqC;
-dftC :: Vector Float
- 
--- [1.82616,1.86862,- 0.017503,0.502001,- 0.598725,- 0.0314404,- 2.62557]
-
--- test  8
-
-(nagInverseDFT dftC) :: Vector Complex Float
- 
--- [0.34907, 0.5489 + 1.5137 %i, 0.74776 + 1.3285 %i, 0.94459 + 1.1385 %i,
---  0.94459 - 1.1385 %i, 0.74776 - 1.3285 %i, 0.5489 - 1.5137 %i]
-
--- test  9
-
-nagHermitianInverseDFT dftC
- 
--- [0.34907000000000005, 0.54889999999999983, 0.74775999999999987,
---  0.94459000000000004, 1.1385000000000003, 1.3284999999999998,
---  1.5136999999999998]
-
--- test 10:
-
-dftsD := nagDFT seqsD;
-
-dftsD :: List Vector Complex Float
- 
--- [
---   [1.07373, - 0.104062 - 0.00438406 %i, 0.112554 - 0.373777 %i, - 0.146684,
---    0.112554 + 0.373777 %i, - 0.104062 + 0.00438406 %i]
---   ,
-
---   [1.39609, - 0.0365178 + 0.466584 %i, 0.077955 - 0.0607051 %i, - 0.152072,
---    0.077955 + 0.0607051 %i, - 0.0365178 - 0.466584 %i]
---   ,
-
---   [1.12374, 0.0914068 - 0.050841 %i, 0.393551 + 0.345775 %i, 0.153011,
---    0.393551 - 0.345775 %i, 0.0914068 + 0.050841 %i]
---   ]
-
--- test 11:
-
-invdftsD := nagInverseDFT dftsD ;
-invdftsD :: List Vector Complex Float
- 
--- [[0.3854,0.6772,0.1138,0.6751,0.6362,0.1424],
---  [0.5417,0.2983,0.1181,0.7255,0.8638,0.8723],
---  [0.9172,0.0644,0.6037,0.643,0.0428,0.4815]]
-
--- test 12:
-
-dftsE := nagDFT seqsE;
-dftsE :: List Vector Float
- 
--- [[1.0788,0.662291,- 0.239146,- 0.578284,0.459192,- 0.438816],
---  [0.857321,1.22614,0.353348,- 0.222169,0.341327,- 1.22908],
---  [1.18245,0.262509,0.674406,0.552278,0.0539906,- 0.478963]]
-
--- test 13:
-
-invdftsE := nagInverseDFT dftsE;
-invdftsE :: List Vector Complex Float
- 
--- [
---   [0.3854, 0.6772 + 0.1424 %i, 0.1138 + 0.6362 %i, 0.6751,
---    0.1138 - 0.6362 %i, 0.6772 - 0.1424 %i]
---   ,
-
---   [0.5417, 0.2983 + 0.8723 %i, 0.1181 + 0.8638 %i, 0.7255,
---    0.1181 - 0.8638 %i, 0.2983 - 0.8723 %i]
---   ,
-
---   [0.9172, 0.0644 + 0.4815 %i, 0.6037 + 0.0428 %i, 0.643,
---    0.6037 - 0.0428 %i, 0.0644 - 0.4815 %i]
---   ]
-
--- test 14:
-
-hdftsD := nagHermitianDFT seqsD;
-map(expand,hdftsD) :: List Vector Complex Float
- 
--- [
---   [1.07373, - 0.104062 - 0.00438406 %i, 0.112554 - 0.373777 %i, - 0.146684,
---    0.112554 + 0.373777 %i, - 0.104062 + 0.00438406 %i]
---   ,
-
---   [1.39609, - 0.0365178 + 0.466584 %i, 0.077955 - 0.0607051 %i, - 0.152072,
---    0.077955 + 0.0607051 %i, - 0.0365178 - 0.466584 %i]
---   ,
-
---   [1.12374, 0.0914068 - 0.050841 %i, 0.393551 + 0.345775 %i, 0.153011,
---    0.393551 - 0.345775 %i, 0.0914068 + 0.050841 %i]
---   ]
-
--- test 15:
-
-(nagInverseDFT hdftsD) :: List Vector Float
- 
--- [[0.3854,0.6772,0.1138,0.6751,0.6362,0.1424],
---  [0.5417,0.2983,0.1181,0.7255,0.8638,0.8723],
---  [0.9172,0.0644,0.6037,0.643,0.0428,0.4815]]
-
--- test 16:
-
-dftsF := nagDFT seqsF;
-dftsF :: List Vector Complex Float
- 
--- [
---   [1.07373 + 1.39609 %i, - 0.570647 - 0.0409019 %i, 0.173259 - 0.295822 %i,
---    - 0.146684 - 0.152072 %i, 0.0518489 + 0.451732 %i,
---    0.362522 - 0.0321337 %i]
---   ,
-
---   [1.12374 + 1.06765 %i, 0.172759 + 0.0385858 %i, 0.418548 + 0.748083 %i,
---    0.153011 + 0.17522 %i, 0.368555 + 0.0565331 %i, 0.0100542 + 0.140268 %i]
---   ,
-
---   [0.909985 + 1.76167 %i, - 0.305418 + 0.0624335 %i,
---    0.407884 - 0.0694786 %i, - 0.078547 + 0.0725049 %i,
---    - 0.119334 + 0.128511 %i, - 0.531409 - 0.433531 %i]
---   ]
-
--- test 17:
-
-invdftsF := nagInverseDFT dftsF ;
-invdftsF :: List Vector Complex Float
- 
--- [
---   [0.3854 + 0.5417 %i, 0.6772 + 0.2983 %i, 0.1138 + 0.1181 %i,
---    0.6751 + 0.7255 %i, 0.6362 + 0.8638 %i, 0.1424 + 0.8723 %i]
---   ,
-
---   [0.9172 + 0.9089 %i, 0.0644 + 0.3118 %i, 0.6037 + 0.3465 %i,
---    0.643 + 0.6198 %i, 0.0428 + 0.2668 %i, 0.4815 + 0.1614 %i]
---   ,
-
---   [0.1156 + 0.6214 %i, 0.0685 + 0.8681 %i, 0.206 + 0.706 %i,
---    0.863 + 0.8652 %i, 0.6967 + 0.919 %i, 0.2792 + 0.3355 %i]
---   ]
-
--- test 18:
-
-nagHermitianInverseDFT dftsE
- 
--- [
---   [0.38540000000000013, 0.67720000000000025, 0.11380000000000001,
---    0.67510000000000014, 0.63620000000000021, 0.14240000000000003]
---   ,
-
---   [0.54170000000000018, 0.29830000000000012, 0.1181, 0.72550000000000014,
---    0.86380000000000023, 0.87230000000000019]
---   ,
-
---   [0.91720000000000035, 0.064399999999999999, 0.60370000000000024,
---    0.64300000000000013, 0.042799999999999991, 0.48150000000000015]
---   ]
-
--- error tests:
-
--- test 19:
-
-nagDFT [vector [0.3854 + 0.5417*%i, 0.6772 + 0.2983*%i,   _
-                0.1138 + 0.1181*%i, 0.6751 + 0.7255*%i,   _
-                0.6362 + 0.8638*%i, 0.1424 + 0.8723*%i],  _
-        vector [0.1156 + 0.6214*%i, 0.0685 + 0.8681*%i,   _
-                0.6967 + 0.9190*%i, 0.2792 + 0.3355*%i]]
- 
--- Error signalled from user code:
---    The data sequences in nagDFT must all have the same length. The 
---    length of sequence 1 is 6 that of sequence 2 is 4.
-
--- test 20:
-
-nagHermitianDFT [vector [0.3854, 0.6751, 0.6362, 0.1424], _
-                 vector [0.5417, 0.7255, 0.8638, 0.8723], _
-                 vector [0.9172, 0.0428, 0.4815]]
- 
--- Error signalled from user code:
---    The data sequences in nagHermitianDFT must all have the same 
---    length. The length of sequence 1 is 4 that of sequence 3 is 3.
-
--- test 21:
-
-badSeqs : List PackedHermitianSequence DoubleFloat
-badSeqs := [pHS [0.3854, 0.1138, 0.6751, 0.6362, 0.1424],         _
-            pHS [0.5417, 0.2983, 0.1181, 0.7255, 0.8638, 0.8723], _
-            pHS [0.9172, 0.0644, 0.6037, 0.6430, 0.0428, 0.4815]];
- 
-nagDFT badSeqs
- 
--- Error signalled from user code:
---    The data sequences in nagDFT must all have the same length. The 
---    length of sequence 1 is 5 that of sequence 2 is 6.
-
-outputGeneral()
-
-output "End of tests"
-
-#endif
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
--- To test:
--- sed -ne '1,/^#if NeverAssertThis/d;/#endif/d;p' < ndftip.as > ndftip.input
--- axiom
--- )set nag host <some machine running nagd>
--- )r ndftip.input
-
-#unassert saturn
-
-#include "axiom.as"
-
-DF     ==> DoubleFloat ;
-CDF    ==> Complex DoubleFloat ;
-LDF    ==> List DoubleFloat ;
-LLDF   ==> List LDF ;
-VDF    ==> Vector DoubleFloat ;
-LVDF   ==> List VDF ;
-VCDF   ==> Vector Complex DoubleFloat ;
-LVCDF  ==> List VCDF ;
-MDF    ==> Matrix DoubleFloat ;
-MCDF   ==> Matrix Complex DoubleFloat ;
-INT    ==> Integer ;
-NNI    ==> NonNegativeInteger ;
-RSLT   ==> Result ;
-STRG   ==> String ;
-PHSDF  ==> PackedHermitianSequence DF;
-LPHSDF ==> List PackedHermitianSequence DF;
-
-<<NagDiscreteFourierTransformInterfacePackage>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/nepip.as.pamphlet b/src/algebra/nepip.as.pamphlet
deleted file mode 100644
index 4d36177..0000000
--- a/src/algebra/nepip.as.pamphlet
+++ /dev/null
@@ -1,626 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra nepip.as}
-\author{Michael Richardson}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{NagEigenInterfacePackage}
-<<NagEigenInterfacePackage>>=
-+++ Author: M.G. Richardson
-+++ Date Created: 1996 January 12
-+++ Date Last Updated:
-+++ Basic Functions:
-+++ Related Constructors:
-+++ Also See:
-+++ AMS Classifications:
-+++ Keywords:
-+++ References:
-+++ Description:
-+++ This package provides Axiom-like interfaces to the NAG generalised
-+++ eigenvalue and eigenvector routines in the NAGlink.
-
-DF     ==> DoubleFloat ;
-CDF    ==> Complex DoubleFloat ;
-FFCDF  ==> FormalFraction Complex DoubleFloat ;
-LFFCDF ==> List FormalFraction Complex DoubleFloat ;
-LDF    ==> List DoubleFloat ;
-LCDF   ==> List Complex DoubleFloat ;
-LLDF   ==> List LDF ;
-VDF    ==> Vector DoubleFloat ;
-LVDF   ==> List VDF ;
-VCDF   ==> Vector Complex DoubleFloat ;
-LVCDF  ==> List VCDF ;
-MDF    ==> Matrix DoubleFloat ;
-MCDF   ==> Matrix Complex DoubleFloat ;
-INT    ==> Integer ;
-NNI    ==> NonNegativeInteger ;
-RCD    ==> Record ;
-RSLT   ==> Result ;
-STRG   ==> String ;
-UNNRES ==> Union(a:LDF,b:LFFCDF) ; -- a & b are dummy tags
-RURSLV ==> RCD(eigenvalues : UNNRES, eigenvectors : LVCDF) ;
-
-NagEigenInterfacePackage: with {
-
-  nagEigenvalues : (MDF,MDF,DF) -> UNNRES ;
-
-++ nagEigenvalues(A,B,eps) returns a list of the eigenvalues
-#if saturn
-++ $ \lambda $
-#else
-++ \spad{l}
-#endif
-++ of the system
-#if saturn
-++ $ A x = \lambda B x $
-#else
-++ \spad{A*x = l*B*x}
-#endif
-++ 
-++ The numerical calculation is performed by one of the NAG routines
-++ F02ADF and F02BJF, depending on the the form of \spad{A} and B.
-++ The result is of type Union(List DoubleFloat, List FormalFraction
-++ Complex DoubleFloat), the first branch resulting from F02ADF and
-++ the second from F02BJF.  Note that in the latter case values should
-++ be inspected for numerically small numerators and denominators,
-++ ratios of which may be in effect indeterminate, before the result is
-++ converted to List Complex DoubleFloat.
-++
-++ The parameter eps, if positive, defines a tolerance to be used in
-++ recognising negligable matrix elements when F02BJF is called; setting
-++ this may result in faster execution with less accuracy.
-++
-++ For more detailed information, please consult the NAG manual
-++ via the Browser pages for the operations f02adf and f02bjf.
-
-  nagEigenvalues : (MDF,MDF) -> UNNRES ;
-
-++ nagEigenvalues(A,B) returns a list of the eigenvalues
-#if saturn
-++ $ \lambda $
-#else
-++ \spad{l}
-#endif
-++ of the system
-#if saturn
-++ $ A x = \lambda B x $
-#else
-++ \spad{A*x = l*B*x}
-#endif
-++ 
-++ The numerical calculation is performed by one of the NAG routines
-++ F02ADF and F02BJF, depending on the the form of \spad{A} and B.
-++ The result is of type Union(List DoubleFloat, List FormalFraction
-++ Complex DoubleFloat), the first branch resulting from F02ADF and
-++ the second from F02BJF.  Note that in the latter case values should
-++ be inspected for numerically small numerators and denominators,
-++ ratios of which may be in effect indeterminate, before the result is
-++ converted to List Complex DoubleFloat.
-++
-++ For more detailed information, please consult the NAG manual
-++ via the Browser pages for the operations f02adf and f02bjf.
-
-  nagEigenvectors : (MDF,MDF,DF) -> RURSLV ;
-
-++ nagEigenvectors(A,B,eps) returns a record consisting of a list of the
-++ eigenvalues
-#if saturn
-++ $ \lambda $
-#else
-++ \spad{l}
-#endif
-++ and a list of the corresponding eigenvectors of the system
-#if saturn
-++ $ A x = \lambda B x $
-#else
-++ \spad{A*x = l*B*x}
-#endif
-++ where
-#if saturn
-++ $A$ and $B$
-#else 
-++ \spad{A} and B
-#endif
-
-#if saturn
-++ $B$
-#else 
-++ B
-#endif
-++ is positive-definite.
-++ 
-++ The numerical calculation is performed by one of the NAG routines
-++ F02AEF and F02BJF, depending on the the form of \spad{A} and B.
-++ The first component of the result, \spad{eigenvalues},
-++ is of type Union(List DoubleFloat, List FormalFraction
-++ Complex DoubleFloat), the first branch resulting from F02AEF and
-++ the second from F02BJF.  Note that in the latter case values should
-++ be inspected for numerically small numerators and denominators,
-++ ratios of which may be in effect indeterminate, before the result is
-++ converted to List Complex DoubleFloat.
-++
-++ The parameter eps, if positive, defines a tolerance to be used in
-++ recognising negligable matrix elements when F02BJF is called; setting
-++ this may result in faster execution with less accuracy.
-++
-++ For more detailed information, please consult the NAG manual
-++ via the Browser pages for the operations f02aef and f02bjf.
-
-  nagEigenvectors : (MDF,MDF) -> RURSLV ;
-
-++ nagEigenvectors(A,B) returns a record consisting of a list of the
-++ eigenvalues
-#if saturn
-++ $ \lambda $
-#else
-++ \spad{l}
-#endif
-++ and a list of the corresponding eigenvectors of the system
-#if saturn
-++ $ A x = \lambda B x $
-#else
-++ \spad{A*x = l*B*x}
-#endif
-++ where
-#if saturn
-++ $A$ and $B$
-#else 
-++ \spad{A} and B
-#endif
-
-#if saturn
-++ $B$
-#else 
-++ B
-#endif
-++ is positive-definite.
-++ 
-++ The numerical calculation is performed by one of the NAG routines
-++ F02AEF and F02BJF, depending on the the form of \spad{A} and B.
-++ The first component of the result, \spad{eigenvalues},
-++ is of type Union(List DoubleFloat, List FormalFraction
-++ Complex DoubleFloat), the first branch resulting from F02AEF and
-++ the second from F02BJF.  Note that in the latter case values should
-++ be inspected for numerically small numerators and denominators,
-++ ratios of which may be in effect indeterminate, before the result is
-++ converted to List Complex DoubleFloat.
-++
-++ For more detailed information, please consult the NAG manual
-++ via the Browser pages for the operations f02aef and f02bjf.
-
-} == add {
-
-  import from AnyFunctions1 INT ;
-  import from AnyFunctions1 MDF ;
-  import from CDF;
-  import from ErrorFunctions ;
-  import from MDF ;
-  import from NagResultChecks ;
-  import from NagEigenPackage ;
-  import from List STRG ;
-  import from Symbol ;
-  import from VDF ;
-  import from Boolean ;
-  import from Result ;
-
-  local (..)(a:INT,b:INT):Generator INT == {
-    generate {
-      t := a ;
-      while (t <= b) repeat {
-        yield t ;
-        t := t + 1 ;
-        }
-      }
-    }
-
-  local ipIfail : INT := -1 ;
-
-  -- First, some local functions:
-
-  f02bjfEigVals(A : MDF, B : MDF, orderAB : INT, eps : DF) : LFFCDF == {
-
-  -- orderAB is the common order of the square matrices A and B.
-
-    local f02bjfResult : RSLT ;
-    local numR, numI, den : LDF ;
-
-    f02bjfResult := f02bjf(orderAB,orderAB,orderAB,eps,
-			   false,orderAB,A,B,ipIfail) ;
-    den := entries(row(checkMxDF(f02bjfResult, "beta", "F02BJF"), 1)) ;
-    numR := entries(row(retract(f02bjfResult."alfr") @ MDF, 1)) ;
-    numI := entries(row(retract(f02bjfResult."alfi") @ MDF, 1)) ;
-
-    [ (complex(r,i)/complex(d,0@DF))$FFCDF for r in numR
-					   for i in numI
-					    for d in den ]
-
-  }
- 
- 
-  f02bjfEigVecs(A : MDF, B : MDF, orderAB : INT, eps : DF) : RURSLV == {
-
-  -- orderAB is the common order of the square matrices A and B.
-
-    local size : NNI ;
-    local f02bjfResult : RSLT ;
-    local numR, numI, den : LDF ;
-    local eVals : UNNRES ;
-    local vecMat : MDF ;
-    local eVecs : LVCDF ;
-    local j : INT ;
-    local thisVec, leftVec : VCDF ;
-
-    size := orderAB pretend NNI ;
-
-    f02bjfResult := f02bjf(orderAB,orderAB,orderAB,eps,
-			   true,orderAB,A,B,ipIfail) ;
-
-    den := entries(row(checkMxDF(f02bjfResult, "beta", "F02BJF"), 1)) ;
-    numR := entries(row(retract(f02bjfResult."alfr") @ MDF, 1)) ;
-    numI := entries(row(retract(f02bjfResult."alfi") @ MDF, 1)) ;
-    vecMat := retract(f02bjfResult."v") @ MDF ;
-
-                                             -- outer [] for union type:
-    eVals := [[(complex(r,i)/complex(d,0@DF))$FFCDF for r in numR
-					             for i in numI
-					              for d in den]] ;
-
-    eVecs := [] ;
-    j := orderAB ;
-    while j > 0 repeat {
-      if numI.j ~= 0$DF then {
-	if j = 1 or numI.(j-1) = 0$DF then
-	  error("nagEigenvectors",
-	      "Inconsistent results returned from NAG routine F02BJF") ;
-        thisVec := new(size,0) ;
-	leftVec := new(size,0) ;
-	for i in 1 .. orderAB repeat {
-          thisVec.i := complex(vecMat(i,j-1),-vecMat(i,j)) ;
-	  leftVec.i := complex(vecMat(i,j-1),vecMat(i,j)) ;
-	}
-	eVecs := cons(leftVec,cons(thisVec,eVecs)) ;
-	j := j - 2;
-      }
-      else {
-	thisVec := new(size,0) ;
-	for i in 1 .. orderAB repeat
-	  thisVec.i := complex(vecMat(i,j),0@DF) ;
-        eVecs := cons(thisVec,eVecs) ;
-	j := j - 1 ;
-      }
-    }
-
-    [eVals,eVecs]
-
-  }
-  
-  
-  nagError(routine : STRG, opIfail : INT) : Exit ==
-    error ["An error was detected when calling the NAG Library routine ",
-           routine,
-           ".  The error number (IFAIL value) is ",
-           string(opIfail),
-           ", please consult the NAG manual via the Browser for",
-	   " diagnostic information."] ;
-
-  -- Then the exported functions:
-  
-  nagEigenvalues(A : MDF, B : MDF, eps : DF) : UNNRES  == {
-
-  -- Strategy: if either matrix is asymmetric, use F02BJF, otherwise
-  -- try F02ADF in case B is positive-definite.
-  -- If F02ADF gives IFAIL=1 (should happen quickly if at all),
-  -- not positive-definite, use less efficient F02BJF.
-
-    local rA, rB, cA, cB : NNI ;
-    local orderAB, opIfail: INT ;
-    local vals : UNNRES ;
-
-    rA := nrows A ;
-    rB := nrows B ;
-    
-    if rA ~= rB
-    then error("nagEigenvalues",
-		"the two matrices supplied are of different sizes.") ;
-    orderAB := rA pretend INT ;
-    
-    if symmetric?(A) and symmetric?(B) then {
-      f02adfResult := f02adf(orderAB,orderAB,orderAB,A,B,ipIfail) ;
-      opIfail := retract(f02adfResult."ifail") @ INT ;
-      if zero? opIfail then              -- using [] to give union type:
-        vals := [entries(row(retract(f02adfResult."r") @ MDF,1))] ;
-      else if opIfail = 1 then
-	vals := [f02bjfEigVals(A,B,orderAB,eps)]
-      else
-	nagError("F02BJF",opIfail) ;
-    }
-    else {
-      cA := ncols A ;
-      if cA ~= rA then
-	error("nagEigenvalues",
-	       "the first matrix supplied is not square") ;
-      cB := ncols B ;
-      if cB ~= rB then
-	error("nagEigenvalues",
-	       "the second matrix supplied is not square") ;
-      vals := [f02bjfEigVals(A,B,orderAB,eps)] ;
-    }
-
-  vals
-
-  }
-
-  
-  nagEigenvalues(A : MDF, B : MDF) : UNNRES
-		== nagEigenvalues(A,B,0@DF) ;
-
-
-  nagEigenvectors(A : MDF, B : MDF, eps : DF) : RURSLV == {
-
-  -- Strategy: if either matrix is asymmetric, use F02BJF, otherwise
-  -- try F02AEF in case B is positive-definite.
-  -- If F02AEF gives IFAIL=1 (should happen quickly if at all),
-  -- not positive-definite, use less efficient F02BJF.
-
-    local rA, rB, cA, cB : NNI ;
-    local orderAB, opIfail, j : INT ;
-    local eVals : UNNRES ;
-    local eVecs : LVCDF ;
-    local vecMat : MDF ;
-    local thisVec : VCDF ;
-    local f02aefResult : RSLT ;
-    local result : RURSLV ;
-
-    rA := nrows A ;
-    rB := nrows B ;
-    
-    if rA ~= rB
-    then error("nagEigenvectors",
-		"the two matrices supplied are of different sizes.") ;
-    orderAB := rA pretend INT ;
-    
-    if symmetric?(A) and symmetric?(B) then {
-      f02aefResult := f02aef(orderAB,orderAB,orderAB,
-			     orderAB,A,B,ipIfail) ;
-      opIfail := retract(f02aefResult."ifail") @ INT ;
-      if zero? opIfail then {
-	-- using [] to give union type:
-        eVals := [entries(row(retract(f02aefResult."r") @ MDF,1))] ;
-	vecMat := retract(f02aefResult."v") @ MDF ;
-        eVecs := [] ;
-        j := orderAB ;
-        while j > 0 repeat {
-          thisVec := new(rA,0) ;
-          for i in 1 .. orderAB repeat
-          thisVec.i := complex(vecMat(i,j),0@DF) ;
-          eVecs := cons(thisVec,eVecs) ;
-          j := j - 1 ;
-        }
-	result := [eVals,eVecs]
-      }
-      else if opIfail = 1 then
-	result := f02bjfEigVecs(A,B,orderAB,eps)
-      else
-	nagError("F02BJF",opIfail) ;
-    }
-    else {
-      cA := ncols A ;
-      if cA ~= rA then
-	error("nagEigenvectors",
-	       "the first matrix supplied is not square") ;
-      cB := ncols B ;
-      if cB ~= rB then
-	error("nagEigenvectors",
-	       "the second matrix supplied is not square") ;
-      result := f02bjfEigVecs(A,B,orderAB,eps) ;
-    }
-
-  result
-
-  }
-
-  
-  nagEigenvectors(A : MDF, B : MDF) : RURSLV
-		 == nagEigenvectors(A,B,0@DF) ;
-
-}
-
-#if NeverAssertThis
-
--- Note that the conversions of results from DoubleFloat to Float
--- will become unnecessary if outputGeneral is extended to apply to
--- DoubleFloat quantities.
-
-)lib nrc
-)lib ffrac
-)lib nepip
-
-outputGeneral 5
-
-mA1 := matrix [[ 0.5 ,   1.5 ,   6.6 ,   4.8],  _
-               [ 1.5 ,   6.5 ,  16.2 ,   8.6],  _
-               [ 6.6 ,  16.2 ,  37.6 ,   9.8],  _
-               [ 4.8 ,   8.6 ,   9.8 , -17.1]];
-
-mB1 := matrix[[ 1 ,  3 ,   4 ,  1],  _
-              [ 3 , 13 ,  16 , 11],  _
-              [ 4 , 16 ,  24 , 18],  _
-              [ 1 , 11 ,  18 , 27]];
-
-mA2 := matrix [[ 3.9 ,  12.5 , -34.5 ,  -0.5],  _
-               [ 4.3 ,  21.5 , -47.5 ,   7.5],  _
-               [ 4.3 ,  21.5 , -43.5 ,   3.5],  _
-               [ 4.4 ,  26.0 , -46.0 ,   6.0]];
-
-mB2 := matrix[[ 1 , 2 , -3 , 1],  _
-              [ 1 , 3 , -5 , 4],  _
-              [ 1 , 3 , -4 , 3],  _
-              [ 1 , 3 , -4 , 4]];
-
-nagEigenvalues(mA1,mB1) :: List Float
-
---       [- 3.0,- 1.0,2.0,4.0]
-
-vv1 := nagEigenvectors(mA1,mB1);
-(vv1.eigenvalues) :: List Float
-
---       [- 3.0,- 1.0,2.0,4.0]
-
-(vv1.eigenvectors) :: List Vector Complex Float
-
--- [[- 4.35,0.05,1.0,- 0.5], [- 2.05,0.15,0.5,- 0.5], [- 3.95,0.85,0.5,- 0.5],
---  [2.65,0.05,- 1.0,0.5]]
-
-nagEigenvalues(mA2,mB2)
-
--- all components are O(1) or more so:
-
-% :: List Complex Float
-
---       [2.0,3.0 + 4.0 %i,3.0 - 4.0 %i,4.0]
-
-vv2 := nagEigenvectors(mA2,mB2);
-vv2.eigenvalues
-
--- all components are O(1) or more so:
-
-% :: List Complex Float
-
---       [2.0,3.0 + 4.0 %i,3.0 - 4.0 %i,4.0]
-
-vv2.eigenvectors :: List Vector Complex Float
-
--- [[0.99606,0.0056917,0.062609,0.062609],
---
---   [0.94491, 0.18898 + 0.26077 E -14 %i, 0.11339 - 0.15119 %i,
---    0.11339 - 0.15119 %i]
---   ,
---
---   [0.94491, 0.18898 - 0.26077 E -14 %i, 0.11339 + 0.15119 %i,
---    0.11339 + 0.15119 %i]
---   ,
---  [0.98752,0.010972,- 0.032917,0.15361]]
-
--- The same call with eps=0.0001:
-
-vv2a := nagEigenvectors(mA2,mB2,0.0001);
-vv2a.eigenvalues :: List Complex Float
-
---       [1.9989,3.0003 + 3.9994 %i,3.0003 - 3.9994 %i,4.0]
-
-vv2a.eigenvectors :: List Vector Complex Float
-
--- [[0.99605,0.0057355,0.062656,0.062656],
---
---   [0.94491, 0.18899 - 0.000048882 %i, 0.11336 - 0.15119 %i,
---    0.11336 - 0.15119 %i]
---   ,
---
---   [0.94491, 0.18899 + 0.000048882 %i, 0.11336 + 0.15119 %i,
---    0.11336 + 0.15119 %i]
---   ,
---  [0.98751,0.011031,- 0.032912,0.15367]]
-
-mB1(1,1) := -1;
-
--- The next test should fail on F02ADF then call F02BJF:
-
-nagEigenvalues(mA1,mB1)
-
--- all components are O(1) or more so:
-
-% :: List Complex Float
-
---       [3.5016,- 1.5471,0.041212 + 0.21738 %i,0.041212 - 0.21738 %i]
-
--- Similarly, this should fail on F02AEF then call F02BJF:
-
-vv3 := nagEigenvectors(mA1,mB1);
-vv3.eigenvalues
-
--- all components are O(1) or more so:
-
-% :: List Complex Float
-
---       [3.5016,- 1.5471,0.041212 + 0.21738 %i,0.041212 - 0.21738 %i]
-
-vv3.eigenvectors :: List Vector Complex Float
-
---  [[- 0.034577,0.63045,- 0.75202,0.1892],
---   [0.17876,- 0.73845,0.047413,0.64845],
---
---    [0.80838, - 0.00095133 + 0.47557 %i, - 0.20354 - 0.21737 %i,
---     0.15404 + 0.089179 %i]
---    ,
---
---    [0.80838, - 0.00095133 - 0.47557 %i, - 0.20354 + 0.21737 %i,
---     0.15404 - 0.089179 %i]
---   ]
-
-outputGeneral()
-
-output "End of tests"
-
-#endif
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
--- NagEigenProblemInterfacePackage
-
--- To test:
--- sed '1,/^#if NeverAssertThis/d;/#endif/d' < nepip.as > nepip.input
--- axiom
--- )set nag host <some machine running nagd>
--- )r nepip.input
-
-#unassert saturn
-
-#include "axiom.as"
-
-<<NagEigenInterfacePackage>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/noptip.as.pamphlet b/src/algebra/noptip.as.pamphlet
deleted file mode 100644
index 150087b..0000000
--- a/src/algebra/noptip.as.pamphlet
+++ /dev/null
@@ -1,241 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra noptip.as}
-\author{Michael Richardson}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{NagOptimisationInterfacePackage}
-<<NagOptimisationInterfacePackage>>=
-+++ Author: M.G. Richardson
-+++ Date Created: 1996 Feb. 01
-+++ Date Last Updated:
-+++ Basic Functions:
-+++ Related Constructors:
-+++ Also See:
-+++ AMS Classifications:
-+++ Keywords:
-+++ References:
-+++ Description:
-+++ This package provides Axiom-like interfaces to some of the NAG
-+++ optimisation routines in the NAGlink.
-
-NagOptimisationInterfacePackage: with {
-
-  nagMin : (EDF,LEQPDF) -> LEQPDF ;
-
-} == add {
-
-  import from MINT ;
-  import from BOOL ;
-  import from LLDF ;
-  import from VDF ;
-  import from PDF ;
-  import from LEQEDF ;
-  import from MDF ;
-  import from EDF ;
-  import from FL ;
-  import from SMBL ;
-  import from A49 ;
-  import from A55 ;
-  import from U49 ;
-  import from U55 ;
-  import from NagOptimisationPackage ;
-  import from OF ;
-  import from LOF ;
-  import from LLOF ;
-  import from ListFunctions2(INT,OF) ;
-  import from NagResultChecks ;
-  import from EF2DFFL ;
-
-  local (..)(a:INT,b:INT):Generator INT == {
-    generate {
-      t := a ;
-      while (t <= b) repeat {
-        yield t ;
-        t := t + 1 ;
-        }
-      }
-    }
-
-  -- to avoid unrecognised versions of U49 type for e04dgf:
-  e04dgflocal := e04dgf$NagOptimisationPackage pretend
-   ((INT,DF,DF,INT,DF,BOOL, DF,DF,INT,INT,INT,INT,MDF,INT, U49)->RSLT) ;
-  
-  nagMin(objective:EDF,startList:LEQPDF) : LEQPDF == {
-
-  -- Note that, as objective is an EDF, subst and eval
-  -- for this have as 2nd parameters LEQEDFs.
-
-    local nv : INT ;
-    local substList : LEQEDF ;
-    local indxOb : EF ;
-    local startVals : LDF ;
-    local startListXDF : LEQEDF ;
-    local startFVal : DF ;
-    local e04dgfResult : RSLT ;
-    local location : LDF  ;
-
-     
-     nv := ((# startList)@NNI pretend INT) ; -- @ avoids SI
-
-     substList := [lhs(startList.i)::EDF
-	             = (script("x"::SMBL,[[i::OF]]@LLOF))::EDF
-		                                for i in 1..nv] ;
-     -- [x=x[1], y=x[2], etc.]
-
-     indxOb := map(convert$Float,subst(objective,substList)) ;
-     -- objective function as an EF with x[i]'s, as required by A49
-
-     startVals := [retract(rhs(startList.i))@DF for i in 1..nv] ;
-
-     startListXDF := [lhs(startList.i)::EDF = rhs(startList.i)::EDF
-					            for i in 1..nv] ;
-     startFVal := ground(eval(objective,startListXDF))::DF ; 
-     startFVal := startFVal * 1.015625 ;
-
--- Note that there appears to be a problem running the standard NAG
--- example on Suns with an exact value for startFVal.  It looks as if
--- this causes too large a stepsize, perhaps due to exception code
--- being obeyed in the Fortran.  Until this is fixed, using the above
--- slightly perturbed value (adding 1/64) seems to avoid the problem.
-
-     e04dgfResult := e04dgflocal(
-		nv,                        -- No.vbls.
-                                           --
-                                           -- "optional" params next:
-                                           --
-                startFVal,                 -- es(timated obj've fn val)
-                -1.0,                      -- fun:
-                -1,                        -- it:
-                -1.0,                      -- lin:
-                false,                     -- list:
-                -1.0,                      -- ma:
-                -2.0,                      -- opt: made < fun for safety
-                0,                         -- pr:
-                -1,                        -- sta:
-                -1,                        -- sto:
-                -1,                        -- ve:
-                                           --
-                matrix [startVals],        -- initial position estimate
-                -1,                        -- IFAIL
-                [retract(indxOb)@A49]@U49  -- objective function
-                ) ;
-  
-     location := entries(row(checkMxDF(e04dgfResult,"x","E04DGF"),1)) ;
-
-     [ ((retract(lhs(startList.i))@SMBL)::PDF
-	  = (location.i)::PDF)@EQPDF for i in 1..nv ]
-  
-  }
-  
-}
-
-#if NeverAssertThis
-
--- Note that the conversions of results from DoubleFloat to Float
--- will become unnecessary if outputGeneral is extended to apply to
--- DoubleFloat quantities.
-
-)lib nrc
-
-outputGeneral 5
-
-f := %e^x*(4*x^2 + 2*y^2 + 4*x*y + 2*y + 1);
-start := [x=-1.0, y=1.0];
-nagMin(f,start) :: List Equation Polynomial Float
-
---       [x= 0.5,y= - 1.0]
-
-outputGeneral()
-
-output "End of tests"
-
-#endif
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
--- To test:
--- sed '1,/^#if NeverAssertThis/d;/#endif/d' < noptip.as > noptip.input
--- axiom
--- )set nag host <some machine running nagd>
--- )r noptip.input
-
-#unassert saturn
-
-#include "axiom.as"
-
-INT     ==> Integer ;
-NNI     ==> NonNegativeInteger ;
-MINT    ==> Matrix INT ;
-DF      ==> DoubleFloat ;
-EDF     ==> Expression DF ;
-EQEDF   ==> Equation EDF ;
-LEQEDF  ==> List EQEDF ;
-LDF     ==> List DF ;
-LLDF    ==> List LDF ;
-VDF     ==> Vector DF ;
-MDF     ==> Matrix DF ;
-PDF     ==> Polynomial DF ;
-EQPDF   ==> Equation PDF ;
-LEQPDF  ==> List EQPDF ;
-FL      ==> Float ;
-EF      ==> Expression FL ;
-BOOL    ==> Boolean ;
-A49     ==> Asp49("OBJFUN") ;
-A55     ==> Asp55("CONFUN") ;
-U49     ==> Union(fn: FileName, fp: A49) ;
-U55     ==> Union(fn: FileName, fp: A55) ;
-SMBL    ==> Symbol ;
-RSLT    ==> Result ;
-OF      ==> OutputForm ;
-LOF     ==> List OF ;
-LLOF    ==> List LOF ;
-EF2DFFL ==> ExpressionFunctions2(DF,FL) ;
-
-<<NagOptimisationInterfacePackage>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/nqip.as.pamphlet b/src/algebra/nqip.as.pamphlet
deleted file mode 100644
index e3e8d96..0000000
--- a/src/algebra/nqip.as.pamphlet
+++ /dev/null
@@ -1,231 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra nqip.as}
-\author{Michael Richardson}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{NagQuadratureInterfacePackage}
-<<NagQuadratureInterfacePackage>>=
-+++ Author: M.G. Richardson
-+++ Date Created: 1995 Dec. 07
-+++ Date Last Updated:
-+++ Basic Functions:
-+++ Related Constructors:
-+++ Also See:
-+++ AMS Classifications:
-+++ Keywords:
-+++ References:
-+++ Description:
-+++ This package provides Axiom-like interfaces to some of the NAG
-+++ quadrature (numerical integration) routines in the NAGlink.
-
-NagQuadratureInterfacePackage: with {
-
-  nagPolygonIntegrate : (LDF, LDF) ->
-			        RCD(integral : DF, errorEstimate : DF) ;
-
-  ++ nagPolygonIntegrate(xlist,ylist) evaluates the definite integral
-#if saturn
-  ++ $\int_{x_{1}}^{x_{n}}y(x) \, dx$
-#else
-  ++ integrate(y(x), x=x[1]..x[n])
-#endif
-  ++ where the numerical value of the function \spad{y} is specified at
-  ++ the \spad{n} distinct points
-#if saturn
-  ++ $x_{1}, x_{2}, \ldots , x_{n}$.
-#else
-  ++ x[1], x[2] ... x[n].
-#endif
-  ++ The \spad{x} and \spad{y} values are specified in the lists
-  ++ \spad{xlist} and \spad{ylist}, respectively; the \spad{xlist}
-  ++ values must form a strictly monotonic sequence of four or more
-  ++ points.
-  ++ The calculation is performed by the NAG routine D01GAF.
-  ++
-  ++ An estimate of the numerical error in the calculation is also
-  ++ returned; however, by choosing unrepresentative data points to
-  ++ approximate the function it is possible to achieve an arbitrarily
-  ++ large difference between the true integral and the value
-  ++ calculated.
-  ++ For more detailed information, please consult the NAG
-  ++ manual via the Browser page for the operation d01gaf.
-
-  nagPolygonIntegrate : MDF -> RCD(integral : DF, errorEstimate : DF) ;
-
-
-} == add {
-
-  import from NagIntegrationPackage ;
-  import from NagResultChecks ;
-  import from AnyFunctions1 DF ;
-  import from STRG ;
-  import from List STRG ;
-  import from Symbol ;
-  import from LLDF ;
-  import from VDF ;
-  import from MDF ;
-  import from ErrorFunctions ;
-
-  local ipIfail : INT := -1 ;
-  local d01gafError : DF := 0 ;
-
-  nagPolygonIntegrate(xlist : LDF, ylist : LDF)
-		     : RCD(integral : DF, errorEstimate : DF) == {
- 
-    local lx, ly : INT ;
-    local d01gafResult : RSLT ;
-
-    lx := (# xlist) pretend INT ;
-    ly := (# ylist) pretend INT ;
-    if lx ~= ly
-    then error ["The lists supplied to nagPolygonIntegrate are of ",
-		"different lengths: ",
-		string(lx),
-		" and ",
-		string(ly),
-		"."]
-    else {
-      d01gafResult := d01gaf(matrix [xlist],matrix [ylist],lx,ipIfail) ;
-      [checkResult(d01gafResult,"ans","D01GAF"),
-       retract(d01gafResult."er") @ DF]
-    }
-  }
-
-  nagPolygonIntegrate(coords : MDF)
-		     : RCD(integral : DF, errorEstimate : DF) ==
-    if (ncols(coords) pretend INT) ~= 2
-    then error ["Please supply the coordinate matrix in ",
-		"nagPolygonIntegrate with each row consisting of ",
-		"a single x-y pair."]
-    else nagPolygonIntegrate(members column(coords,1),
-		             members column(coords,2)) ;
-
-}
-
-#if NeverAssertThis
-
--- Note that the conversions of results from DoubleFloat to Float
--- will become unnecessary if outputGeneral is extended to apply to
--- DoubleFloat quantities.
-
-)lib nrc
-)lib nqip
-
-outputGeneral 5
-
-xvals := [0.00,0.04,0.08,0.12,0.22,0.26,0.30,0.38,0.39,0.42,0.45, 
-               0.46,0.60,0.68,0.72,0.73,0.83,0.85,0.88,0.90,1.00];
-
-yvals := [4.0000,3.9936,3.9746,3.9432,3.8135,3.7467,3.6697,3.4943,
-                 3.4719,3.4002,3.3264,3.3017,2.9412,2.7352,2.6344,
-                        2.6094,2.3684,2.3222,2.2543,2.2099,2.0000];
-
-result := nagPolygonIntegrate(xvals,yvals);
-result.integral :: Float             
-
---       3.1414
-
-result.errorEstimate :: Float        
-
---       - 0.000025627
-
-coords := transpose matrix [xvals, yvals];
-result := nagPolygonIntegrate coords;
-result.integral :: Float             
-
---       3.1414
-
-result.errorEstimate :: Float        
-
---       - 0.000025627
-
-nagPolygonIntegrate([1,2,3],[1,2,3,4])
- 
---   Error signalled from user code:
---      The lists supplied to nagPolygonIntegrate are of different 
---      lengths: 3 and 4.
-
-nagPolygonIntegrate([[1,2,3],[4,5,6]])
-
---   Error signalled from user code:
---      Please supply the coordinate matrix in nagPolygonIntegrate with
---      each row consisting of single a x-y pair.
-
-outputGeneral()
-
-output "End of tests"
-
-#endif
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
--- NagQuadratureInterfacePackage
-
--- To test:
--- sed -ne '1,/^#if NeverAssertThis/d;/#endif/d;p' < nqip.as > nqip.input
--- axiom
--- )set nag host <some machine running nagd>
--- )r nqip.input
- 
-#unassert saturn
-
-#include "axiom.as"
-
-DF   ==> DoubleFloat ;
-LDF  ==> List DoubleFloat ;
-LLDF ==> List LDF ;
-VDF  ==> Vector DoubleFloat ;
-MDF  ==> Matrix DoubleFloat ;
-INT  ==> Integer ;
-RCD  ==> Record ;
-RSLT ==> Result ;
-STRG ==> String ;
-
-<<NagQuadratureInterfacePackage>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/nrc.as.pamphlet b/src/algebra/nrc.as.pamphlet
deleted file mode 100644
index d317958..0000000
--- a/src/algebra/nrc.as.pamphlet
+++ /dev/null
@@ -1,132 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra nrc.as}
-\author{Michael Richardson}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{NagResultChecks}
-<<NagResultChecks>>=
-+++ Author: M.G. Richardson
-+++ Date Created: 1995 Dec. 06
-+++ Date Last Updated:
-+++ Basic Functions:
-+++ Related Constructors:
-+++ Also See:
-+++ AMS Classifications:
-+++ Keywords:
-+++ References:
-+++ Description:
-
-NagResultChecks: with {
-
-  checkResult   : (RSLT, STRG, STRG) -> DF ;
-  checkCxResult : (RSLT, STRG, STRG) -> CDF ;
-  checkMxCDF    : (RSLT, STRG, STRG) -> MCDF ;
-  checkMxDF     : (RSLT, STRG, STRG) -> MDF ;
-
-} == add {
-
-  import from DF ;
-  import from SMBL ;
-  import from INT ;
-  import from AnyFunctions1 INT ;
-  import from AnyFunctions1 DF ;
-  import from AnyFunctions1 CDF ;
-  import from AnyFunctions1 MDF ;
-  import from AnyFunctions1 MCDF ;
-  import from ErrorFunctions ;
-  import from STRG ;
-  import from List STRG ;
-
-  checkResult(returnValue : RSLT, returnKey : STRG, routine : STRG) : DF ==
-    if not zero?(retract(returnValue."ifail") @ INT)
-       then nagError(routine, retract(returnValue."ifail") @ INT)
-       else retract(returnValue.(returnKey::SMBL)) @ DF ;
-  
-  checkCxResult(returnValue : RSLT, returnKey : STRG, routine : STRG) : CDF ==
-    if not zero?(retract(returnValue."ifail") @ INT)
-       then nagError(routine, retract(returnValue."ifail") @ INT)
-       else retract(returnValue.(returnKey::SMBL)) @ CDF ;
-  
-  checkMxDF(returnValue : RSLT, returnKey : STRG, routine : STRG) : MDF ==
-    if not zero?(retract(returnValue."ifail") @ INT)
-       then nagError(routine, retract(returnValue."ifail") @ INT)
-       else retract(returnValue.(returnKey::SMBL)) @ MDF ;
-  
-  checkMxCDF(returnValue : RSLT, returnKey : STRG, routine : STRG) : MCDF ==
-    if not zero?(retract(returnValue."ifail") @ INT)
-       then nagError(routine, retract(returnValue."ifail") @ INT)
-       else retract(returnValue.(returnKey::SMBL)) @ MCDF ;
-  
-  nagError(routine : STRG, opIfail : INT) : Exit ==
-    error ["An error was detected when calling the NAG Library routine ",
-           routine,
-           ".  The error number (IFAIL value) is ",
-           string(opIfail),
-           ", please consult the NAG manual via the Browser for",
-	   " diagnostic information."] ;
-}
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
--- N.B. ndftip.as, nqip.as and nsfip.as inline from nrc
--- and must be recompiled if this is.
- 
-#include "axiom.as"
-
-DF   ==> DoubleFloat ;
-CDF  ==> Complex DoubleFloat ;
-MDF  ==> Matrix DoubleFloat ;
-MCDF ==> Matrix Complex DoubleFloat ;
-INT  ==> Integer ;
-RSLT ==> Result ;
-SMBL ==> Symbol ;
-STRG ==> String ;
-
-<<NagResultChecks>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/algebra/nsfip.as.pamphlet b/src/algebra/nsfip.as.pamphlet
deleted file mode 100644
index d6bbba9..0000000
--- a/src/algebra/nsfip.as.pamphlet
+++ /dev/null
@@ -1,1223 +0,0 @@
-\documentclass{article}
-\usepackage{axiom}
-\begin{document}
-\title{\$SPAD/src/algebra nsfip.as}
-\author{Michael Richardson}
-\maketitle
-\begin{abstract}
-\end{abstract}
-\eject
-\tableofcontents
-\eject
-\section{NagSpecialFunctionsInterfacePackage}
-<<NagSpecialFunctionsInterfacePackage>>=
-+++ Author: M.G. Richardson
-+++ Date Created: 1995 Nov. 27
-+++ Date Last Updated:
-+++ Basic Functions:
-+++ Related Constructors:
-+++ Also See:
-+++ AMS Classifications:
-+++ Keywords:
-+++ References:
-+++ Description:
-+++ This package provides Axiom-like interfaces to those of the NAG
-+++ special functions in the NAGlink for which no equivalent
-+++ functionality is transparently present in Axiom.
-
-NagSpecialFunctionsInterfacePackage: with {
-
-  nagExpInt : DF -> DF ;
-
-  ++ nagExpInt calculates an approximation to the exponential integral,
-  ++ \spad{E1}, defined by
-#if saturn
-  ++ \[E_{1}(x) = \int_{x}^{\infty }\frac{e^{-u}}{u}\,du\]
-#else
-  ++ \spad{E1(x) = integrate(1/u*%e^u, u=x..%infinity)}
-#endif
-  ++ using the NAG routine S13AAF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s13aaf.
-
-  nagSinInt : DF -> DF ;
-
-  ++ nagSinInt calculates an approximation to the sine integral,
-  ++ \spad{Si}, defined by
-#if saturn
-  ++ \[{\rm Si} (x) = \int_{0}^{x}\frac{\sin u}{u}\,du\]
-#else
-  ++ \spad{Si(x) = integrate(1/u*sin(u), u=0..x)}
-#endif
-  ++ using the NAG routine S13ADF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s13adf.
-
-  nagCosInt : DF -> DF ;
-
-  ++ nagCosInt calculates an approximation to the cosine integral,
-  ++ \spad{Ci}, defined by
-#if saturn
-  ++ \[{\rm Ci} (x) =
-  ++ \gamma + \ln x+ \int_{0}^{x}\frac{\cos u- 1}{u}\,du\]
-#else
-  ++ \spad{Ci(x) = gamma + log x + integrate(1/u*cos(u), u=0..x)}
-  ++ where \spad{gamma} is Euler's constant,
-#endif
-  ++ using the NAG routine S13ACF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s13acf.
-
-  nagIncompleteGammaP : (DF, DF) -> DF ; -- to machine precision
-
-  ++ nagIncompleteGammaP evaluates the incomplete gamma function
-  ++ \spad{P}, defined by 
-#if saturn
-  ++ \[P(a,x) & = & \frac{1}{\Gamma(a)}\int_{0}^{x}t^{a-1}e^{-t}\,dt\]
-#else
-  ++ \spad{P(a,x) = 1/Gamma(a)*integrate(t^(a-1)%e^(-t),t=0..x)}
-#endif
-  ++ to machine precision, using the NAG routine S14BAF.
-
-  nagIncompleteGammaP : (DF, DF, DF) -> DF ;
-
-  ++ nagIncompleteGammaP(a,x,tol) evaluates the incomplete gamma
-  ++ function \spad{P}, defined by
-#if saturn
-  ++ \[P(a,x) & = & \frac{1}{\Gamma(a)}\int_{0}^{x}t^{a-1}e^{-t}\,dt\]
-#else
-  ++ \spad{P(a,x) = 1/Gamma(a)*integrate(t^(a-1)%e^(-t),t=0..x)}
-#endif
-  ++ to a relative accuracy \spad{tol}, using the NAG routine S14BAF.
-
-  nagIncompleteGammaQ : (DF, DF) -> DF ;
-
-  ++ nagIncompleteGammaQ evaluates the incomplete gamma function
-  ++ \spad{Q}, defined by 
-#if saturn
-  ++ \[Q(a,x)&=&\frac{1}{\Gamma(a)}\int_{x}^{\infty}t^{a-1}e^{-t}\,dt\]
-#else
-  ++ \spad{Q(a,x) = 1/Gamma(a)*integrate(t^(a-1)%e^(-t),t=x..%infinity)}
-#endif
-  ++ to machine precision, using the NAG routine S14BAF.
-
-  nagIncompleteGammaQ : (DF, DF, DF) -> DF ;
-
-  ++ nagIncompleteGammaQ(a,x,tol) evaluates the incomplete gamma
-  ++ function \spad{Q}, defined by
-#if saturn
-  ++ \[Q(a,x)&=&\frac{1}{\Gamma(a)}\int_{x}^{\infty}t^{a-1}e^{-t}\,dt\]
-#else
-  ++ \spad{Q(a,x) = 1/Gamma(a)*integrate(t^(a-1)%e^(-t),t=x..%infinity)}
-#endif
-  ++ to a relative accuracy \spad{tol}, using the NAG routine S14BAF.
-
-  nagErf : DF -> DF ;
-
-  ++ nagErf calculates an approximation to the error function,
-  ++ \spad{erf}, defined by
-#if saturn
-  ++ \[{\rm erf}\, x = \frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}\,dt\]
-#else
-  ++ \spad{erf(x) = 2/sqrt(\%pi)*integrate(\%e^(-t^2),t=0..x)}
-#endif
-  ++ using the NAG routine S15AEF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s15aef.
-
-  nagErfC : DF -> DF ;
-
-  ++ nagErfC calculates an approximation to the complementary error
-  ++ function \spad{erfc}, defined by
-#if saturn
-  ++ \[{\rm erfc}\,x =
-  ++       \frac{2} {\sqrt{\pi}}\int_{x}^{\infty}e^{-t^{2}}\,dt\]
-#else
-  ++ \spad{erfc(x) = 2/sqrt(%pi)*integrate(%e^(-t^2),t=x..%infinity)}
-#endif
-  ++ using the NAG routine S15ADF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s15adf.
-
-  nagDAiryAi : DF -> DF ;
-
-  ++ nagDAiryAi calculates an approximation to \spad{Ai'}, the
-  ++ derivative of the Airy function \spad{Ai}, using the NAG routine
-  ++ S17AJF.
-
-  nagDAiryAi : CDF -> CDF ;
-
-  ++ nagDAiryAi calculates an approximation to \spad{Ai'}, the
-  ++ derivative of the Airy function \spad{Ai}, using the NAG routine
-  ++ S17DGF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s17dgf.
-
-  nagDAiryBi : DF -> DF ;
-
-  ++ nagDAiryBi calculates an approximation to \spad{Bi'}, the
-  ++ derivative of the Airy function \spad{Bi}, using the NAG routine
-  ++ S17AKF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s17akf.
-
-  nagDAiryBi : CDF -> CDF ;
-
-  ++ nagDAiryBi calculates an approximation to \spad{Bi'}, the
-  ++ derivative of the Airy function \spad{Bi}, using the NAG routine
-  ++ S17DHF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s17dhf.
-
-  nagScaledDAiryAi : CDF -> CDF ;
-
-  ++ nagDAiryAi(z) calculates an approximation to \spad{Ai'(z)}, the
-  ++ derivative of the Airy function \spad{Ai(z)}, with the result
-  ++ scaled by a factor 
-#if saturn
-  ++ $e^{2z\sqrt{z}/ 3}$
-#else
-  ++ \spad{%e^(2*z*sqrt(z)/3)}
-#endif
-  ++ using the NAG routine S17DGF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s17dgf.
-
-  nagScaledDAiryBi : CDF -> CDF ;
-
-  ++ nagDAiryBi(z) calculates an approximation to \spad{Bi'(z)}, the
-  ++ derivative of the Airy function \spad{Bi(z)}, with the result
-  ++ scaled by a factor 
-#if saturn
-  ++ $e^{2z\sqrt{z}/ 3}$
-#else
-  ++ \spad{%e^(2*z*sqrt(z)/3)}
-#endif
-  ++ using the NAG routine S17DHF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s17dhf.
-
-  nagHankelH1 : (DF, CDF, INT) -> MCDF ;
-
-  ++ nagHankelH1(nu,z,n) calculates an approximation to a sequence of n
-  ++ values of the Hankel function
-#if saturn
-  ++ $H_{\nu + k}^{(1)}(z)$
-#else
-  ++ \spad{H1(nu + k, z)}
-#endif
-  ++ for non-negative nu and \spad{k = 0,1 ... n-1}, using the NAG
-  ++ routine S17DLF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s17dlf.
-
-  nagHankelH2 : (DF, CDF, INT) -> MCDF ;
-
-  ++ nagHankelH2(nu,z,n) calculates an approximation to a sequence of n
-  ++ values of the Hankel function
-#if saturn
-  ++ $H_{\nu + k}^{(2)}(z)$
-#else
-  ++ \spad{H2(nu + k, z)}
-#endif
-  ++ for non-negative nu and \spad{k = 0,1 ... n-1}, using the NAG
-  ++ routine S17DLF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s17dlf.
-
-  nagScaledHankelH1 : (DF, CDF, INT) -> MCDF ;
-
-  ++ nagHankelH1(nu,z,n) calculates an approximation to a sequence of n
-  ++ values of the Hankel function
-#if saturn
-  ++ $H_{\nu + k}^{(1)}(z)$
-#else
-  ++ \spad{H1(nu + k, z)}
-#endif
-  ++ for non-negative nu and \spad{k = 0,1 ... n-1}, with the result
-  ++ scaled by a factor
-#if saturn
-  ++ $e^{-iz}
-#else
-  ++ \spad{%e^(-%i*z)}
-#endif
-  ++ using the NAG routine S17DLF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s17dlf.
-
-  nagScaledHankelH2 : (DF, CDF, INT) -> MCDF ;
-
-  ++ nagHankelH2(nu,z,n) calculates an approximation to a sequence of n
-  ++ values of the Hankel function
-#if saturn
-  ++ $H_{\nu + k}^{(2)}(z)$
-#else
-  ++ \spad{H2(nu + k, z)}
-#endif
-  ++ for non-negative nu and \spad{k = 0,1 ... n-1}, with the result
-  ++ scaled by a factor
-#if saturn
-  ++ $e^{iz}
-#else
-  ++ \spad{%e^(%i*z)}
-#endif
-  ++ using the NAG routine S17DLF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s17dlf.
-
-  nagKelvinBer : DF -> DF ;
-
-  ++ nagKelvinBer calculates an approximation to the Kelvin function
-  ++ \spad{ber}, using the NAG routine S19AAF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s19aaf.
-
-  nagKelvinBei : DF -> DF ;
-
-  ++ nagKelvinBei calculates an approximation to the Kelvin function
-  ++ \spad{bei}, using the NAG routine S19ABF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s19abf.
-
-  nagKelvinKer : DF -> DF ;
-
-  ++ nagKelvinKer calculates an approximation to the Kelvin function
-  ++ \spad{ker}, using the NAG routine S19ACF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s19acf.
-
-  nagKelvinKei : DF -> DF ;
-
-  ++ nagKelvinKei calculates an approximation to the Kelvin function
-  ++ \spad{kei}, using the NAG routine S19ADF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s19adf.
-
-  nagFresnelS : DF -> DF ;
-
-  ++ nagFresnelS calculates an approximation to the Fresnel integral
-  ++ \spad{S}, defined by
-#if saturn
-  ++ \[S(x) = \int_{0}^{x}\sin\left(\frac{\pi}{2}t^{2}\right)\,dt\]
-#else
-  ++ \spad{S(x) = integrate(sin(%pi/2*t^2),t=0..x)}
-#endif
-  ++ using the NAG routine S20ACF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s20acf.
-
-  nagFresnelC : DF -> DF ;
-
-  ++ nagFresnelC calculates an approximation to the Fresnel integral
-  ++ \spad{C}, defined by
-#if saturn
-  ++ \[C(x) = \int_{0}^{x}\cos\left(\frac{\pi}{2}t^{2}\right)\,dt\]
-#else
-  ++ \spad{C(x) = integrate(cos(%pi/2*t^2),t=0..x)}
-#endif
-  ++ using the NAG routine S20ADF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s20adf.
-
-  nagEllipticIntegralRC : (DF, DF) -> DF ;
-
-  ++ nagEllipticIntegralRC(x,y) calculates an approximation to the
-  ++ elementary (degenerate symmetrised elliptic) integral
-#if saturn
-  ++ \[R_{C}(x,y) =
-  ++       \frac{1}{2}\int_{0}^{\infty}\frac{dt}{\sqrt{t+x}(t+y)}\]
-#else
-  ++ \spad{RC(x,y) = 1/2*integrate(1/(sqrt(t+x)*(t+y)),t=0..\infinity)}
-#endif
-  ++ using the NAG routine S21BAF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s21baf.
-
-  nagEllipticIntegralRF : (DF, DF, DF) -> DF ;
-
-  ++ nagEllipticIntegralRF(x,y,z) calculates an approximation to the
-  ++ symmetrised elliptic integral of the first kind,
-#if saturn
-  ++ \[R_{F}(x, y, z) =
-  ++     \frac{1}{2}\int_{0}^{\infty}\frac{dt}{\sqrt{(t+x)(t+y)(t+z)}}\]
-#else
-  ++ \spad{RF(x,y,z) =
-  ++       1/2*integrate(1/sqrt((t+x)*(t+y)*(t+z)),t=0..\infinity)}
-#endif
-  ++ using the NAG routine S21BBF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s21bbf.
-
-  nagEllipticIntegralRD : (DF, DF, DF) -> DF ;
-
-  ++ nagEllipticIntegralRD(x,y,z) calculates an approximation to the
-  ++ symmetrised elliptic integral of the second kind, 
-#if saturn
-  ++ \[R_{D}(x, y, z) =
-  ++ \frac{3}{2}\int_{0}^{\infty}\frac{dt}{\sqrt{(t+x)(t+y)(t+z)^{3}}}\]
-#else
-  ++ \spad{RD(x,y,z) =
-  ++       1/2*integrate(1/sqrt((t+x)*(t+y)*(t+z)^3),t=0..\infinity)}
-#endif
-  ++ using the NAG routine S21BCF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s21bcf.
-
-  nagEllipticIntegralRJ : (DF, DF, DF, DF) -> DF ;
-
-  ++ nagEllipticIntegralRJ(x,y,z,rho) calculates an approximation to
-  ++ the symmetrised elliptic integral of the third kind,
-#if saturn
-  ++ \[R_{J}(x, y, z, \rho ) =
-  ++       \frac{3}{2}\int_{0}^{\infty}
-  ++             \frac{dt}{(t+\rho)\sqrt{(t+x)(t+y)(t+z)}}\]
-#else
-  ++ \spad{RJ(x,y,z,rho) =
-  ++ 3/2*integrate(1/((t+rho)*sqrt((t+x)*(t+y)*(t+z))),t=0..\infinity))u}
-#endif
-  ++ using the NAG routine S21BDF.
-  ++ For detailed information on the accuracy, please consult the NAG
-  ++ manual via the Browser page for the operation s21bdf.
-
-} == add {
-
-  import from NagSpecialFunctionsPackage ;
-  import from NagResultChecks ;
-
-  local ipIfail : Integer := -1 ;
-
-  nagExpInt(x : DF) : DF ==
-    checkResult(s13aaf(x,ipIfail), "s13aafResult", "S13AAF") ;
-
-  nagCosInt(x : DF) : DF ==
-    checkResult(s13acf(x,ipIfail), "s13acfResult", "S13ACF") ;
-
-  nagSinInt(x : DF) : DF ==
-    checkResult(s13adf(x,ipIfail), "s13adfResult", "S13ADF") ;
-
-  nagIncompleteGammaP(a : DF, x : DF) : DF ==
-    checkResult(s14baf(a,x,0.0,ipIfail), "p", "S14BAF") ;
-
-  nagIncompleteGammaP(a : DF, x : DF, tol : DF) : DF ==
-    checkResult(s14baf(a,x,tol,ipIfail), "p", "S14BAF") ;
-
-  nagIncompleteGammaQ(a : DF, x : DF) : DF ==
-    checkResult(s14baf(a,x,0.0,ipIfail), "q", "S14BAF") ;
-
-  nagIncompleteGammaQ(a : DF, x : DF, tol : DF) : DF ==
-    checkResult(s14baf(a,x,tol,ipIfail), "q", "S14BAF") ;
-
-  nagErfC(x : DF) : DF ==
-    checkResult(s15adf(x,ipIfail), "s15adfResult", "S15ADF") ;
-
-  nagErf(x : DF) : DF ==
-    checkResult(s15aef(x,ipIfail), "s15aefResult", "S15AEF") ;
-
-  nagDAiryAi(x : DF) : DF ==
-    checkResult(s17ajf(x,ipIfail), "s17ajfResult", "S17AJF") ;
-
-  nagDAiryAi(z : CDF) : CDF ==
-    checkCxResult(s17dgf("d",z,"u",ipIfail), "ai", "S17DGF") ;
-
-  nagDAiryBi(x : DF) : DF ==
-    checkResult(s17akf(x,ipIfail), "s17akfResult", "S17AKF") ;
-
-  nagDAiryBi(z : CDF) : CDF ==
-    checkCxResult(s17dhf("d",z,"u",ipIfail), "bi", "S17DHF") ;
-
-  nagScaledDAiryAi(z : CDF) : CDF ==
-    checkCxResult(s17dgf("d",z,"s",ipIfail), "ai", "S17DGF") ;
-
-  nagScaledDAiryBi(z : CDF) : CDF ==
-    checkCxResult(s17dhf("d",z,"s",ipIfail), "bi", "S17DHF") ;
-
-  nagHankelH1(order : DF, z : CDF, n : INT) : Matrix CDF == 
-    checkMxCDF(s17dlf(1,order,z,n,"u",ipIfail), "cy", "S17DLF") ;
-
-  nagHankelH2(order : DF, z : CDF, n : INT) : Matrix CDF ==
-    checkMxCDF(s17dlf(2,order,z,n,"u",ipIfail), "cy", "S17DLF") ;
-
-  nagScaledHankelH1(order : DF, z : CDF, n : INT) : Matrix CDF ==
-    checkMxCDF(s17dlf(1,order,z,n,"s",ipIfail), "cy", "S17DLF") ;
-
-  nagScaledHankelH2(order : DF, z : CDF, n : INT) : Matrix CDF ==
-    checkMxCDF(s17dlf(2,order,z,n,"s",ipIfail), "cy", "S17DLF") ;
-
-  nagKelvinBer(x : DF) : DF ==
-    checkResult(s19aaf(x,ipIfail), "s19aafResult", "S19AAF") ;
-
-  nagKelvinBei(x : DF) : DF ==
-    checkResult(s19abf(x,ipIfail), "s19abfResult", "S19ABF") ;
-
-  nagKelvinKer(x : DF) : DF ==
-    checkResult(s19acf(x,ipIfail), "s19acfResult", "S19ACF") ;
-
-  nagKelvinKei(x : DF) : DF ==
-    checkResult(s19adf(x,ipIfail), "s19adfResult", "S19ADF") ;
-
-  nagFresnelS(x : DF) : DF ==
-    checkResult(s20acf(x,ipIfail), "s20acfResult", "S20ACF") ;
-
-  nagFresnelC(x : DF) : DF ==
-    checkResult(s20adf(x,ipIfail), "s20adfResult", "S20ADF") ;
-
-  nagEllipticIntegralRC(x : DF, y : DF) : DF ==
-    checkResult(s21baf(x,y,ipIfail), "s21bafResult", "S21BAF") ;
-
-  nagEllipticIntegralRF(x : DF, y : DF, z : DF) : DF ==
-    checkResult(s21bbf(x,y,z,ipIfail), "s21bbfResult", "S21BBF") ;
-
-  nagEllipticIntegralRD(x : DF, y : DF, z : DF) : DF ==
-    checkResult(s21bcf(x,y,z,ipIfail), "s21bcfResult", "S21BCF") ;
-
-  nagEllipticIntegralRJ(x : DF, y : DF, z : DF, rho : DF) : DF ==
-    checkResult(s21bdf(x,y,z,rho,ipIfail), "s21bdfResult", "S21BDF") ;
-}
-
-#if NeverAssertThis
-
--- Note that the conversions of Results from DoubleFloat to Float
--- will become unnecessary if outputGeneral is extended to apply to
--- DoubleFloat quantities.
-
-)lib nrc
-)lib nsfip
-
-outputGeneral 4
-
--- DF here means DoubleFloat.
--- Results converted to Float as outputGeneral not working on DF.
-
--- nagExpInt : DF -> DF ;
-
-nagExpInt(2) :: Float
-
---       0.0489
-
-nagExpInt(-1) :: Float
-
--- ** ABNORMAL EXIT from NAG Library routine S13AAF: IFAIL =     1
--- ** NAG soft failure - control returned
--- 
---   Error signalled from user code:
---      An error was detected when calling the NAG Library routine 
---      S13AAF. The error number (IFAIL value) is 1, please consult the 
---      NAG manual via the Browser for diagnostic information.
-
--- nagSinInt : DF -> DF ;
-
-nagSinInt(0) :: Float
-
---       0.0
-
-nagSinInt(0.2) :: Float
-
---       0.1996
-
-nagSinInt(0.4) :: Float
-
---       0.3965
-
-nagSinInt(0.6) :: Float
-
---       0.5881
-
-nagSinInt(0.8) :: Float
-
---       0.7721
-
-nagSinInt(1) :: Float
-
---       0.9461
-
--- nagCosInt : DF -> DF ;
-
-nagCosInt(0.2) :: Float
-
---       - 1.042
-
-nagCosInt(0.4) :: Float
-
---       - 0.3788
-
-nagCosInt(0.6) :: Float
-
---       - 0.02227
-
-nagCosInt(0.8) :: Float
-
---       0.1983
-
-nagCosInt(1) :: Float
-
---       0.3374
-
--- nagIncompleteGammaP : (DF, DF) -> DF ; (to machine precision)
-
-nagIncompleteGammaP(2,3) :: Float
-
---       0.8009
-
-nagIncompleteGammaP(7,1) :: Float
-
---       0.00008324
-
-nagIncompleteGammaP(0.5,99) :: Float
-
---       1.0
-
-nagIncompleteGammaP(20,21) :: Float
-
---       0.6157
-
-nagIncompleteGammaP(21,20) :: Float
-
---       0.4409
-
--- nagIncompleteGammaP : (DF, DF, DF) -> DF ; (to specified precision)
-
-nagIncompleteGammaP(7,1,0.1) :: Float
-
---       0.00008313
-
--- nagIncompleteGammaQ : (DF, DF) -> DF ; (to machine precision)
-
-nagIncompleteGammaQ(2,3) :: Float
-
---       0.1991
-
-nagIncompleteGammaQ(7,1) :: Float
-
---       0.9999
-
-nagIncompleteGammaQ(0.5,99) :: Float
-
---       0.5705 E -44
-
-nagIncompleteGammaQ(20,21) :: Float
-
---       0.3843
-
-nagIncompleteGammaQ(21,20) :: Float
-
---       0.5591
-
-nagIncompleteGammaQ(25,14) :: Float
-
---       0.995
-
--- nagIncompleteGammaQ : (DF, DF, DF) -> DF ; (to specified precision)
-
-nagIncompleteGammaQ(25,14,0.1) :: Float
-
---       0.9953
-
--- nagErf : DF -> DF ;
-
-nagErf(-6) :: Float
-
---       - 1.0
-
-nagErf(-4.5) :: Float
-
---       - 1.0
-
-nagErf(-1) :: Float
-
---       - 0.8427
-
-nagErf(1) :: Float
-
---       0.8427
-
-nagErf(4.5) :: Float
-
---       1.0
-
-nagErf(6) :: Float
-
---       1.0
-
--- nagErfC : DF -> DF ;
-
-nagErfC(-10) :: Float
-
---       2.0
-
-nagErfC(-1) :: Float
-
---       1.843
-
-nagErfC(0) :: Float
-
---       1.0
-
-nagErfC(1) :: Float
-
---       0.1573
-
-nagErfC(15) :: Float
-
---       0.7213 E -99
-
--- nagDAiryAi : DF -> DF ;
-
-nagDAiryAi(-10) :: Float
-
---       0.9963
-
-nagDAiryAi(-1) :: Float
-
---       - 0.01016
-
-nagDAiryAi(0) :: Float
-
---       - 0.2588
-
-nagDAiryAi(1) :: Float
-
---       - 0.1591
-
-nagDAiryAi(5) :: Float
-
---       - 0.0002474
-
-nagDAiryAi(10) :: Float
-
---       - 0.3521 E -9
-
-nagDAiryAi(20) :: Float
-
---       - 0.7586 E -26
-
--- nagDAiryAi : CDF -> CDF ;
-
-nagDAiryAi(0.3+0.4*%i) :: Complex Float
-
---       - 0.2612 + 0.03848 %i
-
--- nagDAiryBi : DF -> DF ;
-
-nagDAiryBi(-10) :: Float
-
---       0.1194
-
-nagDAiryBi(-1) :: Float
-
---       0.5924
-
-nagDAiryBi(0) :: Float
-
---       0.4483
-
-nagDAiryBi(1) :: Float
-
---       0.9324
-
-nagDAiryBi(5) :: Float
-
---       1436.0
-
-nagDAiryBi(10) :: Float
-
---       0.1429 E 10
-
-nagDAiryBi(20) :: Float
-
---       0.9382 E 26
-
--- nagDAiryBi : CDF -> CDF ;
-
-nagDAiryBi(0.3+0.4*%i) :: Complex Float
-
---       0.4093 + 0.07966 %i
-
--- nagScaledDAiryAi : CDF -> CDF ;
-
-nagScaledDAiryAi(0.3+0.4*%i) :: Complex Float
-
---       - 0.2744 - 0.02356 %i
-
--- nagScaledDAiryBi : CDF -> CDF ;
-
-nagScaledDAiryBi(0.3+0.4*%i) :: Complex Float
-
---       0.3924 + 0.07638 %i
-
--- nagHankelH1 : (DF, CDF, Int) -> List CDF ;
-
-nagHankelH1(0,0.3+0.4*%i,2) :: Matrix Complex Float
-
---       [0.3466 - 0.5588 %i  - 0.7912 - 0.8178 %i]
-
-nagHankelH1(2.3,2,2) :: Matrix Complex Float
-
---       [0.2721 - 0.7398 %i  0.08902 - 1.412 %i]
-
-nagHankelH1(2.12,-1,2) :: Matrix Complex Float
-
---       [- 0.7722 - 1.693 %i  2.601 + 6.527 %i]
-
--- nagHankelH2 : (DF, CDF, Int) -> List CDF ;
-
-nagHankelH2(6,3.1-1.6*%i,2) :: Matrix Complex Float
-
---       [- 1.371 - 1.28 %i  - 1.491 - 5.993 %i]
-
--- nagScaledHankelH1 : (DF, CDF, Int) -> List CDF ;
-
-nagScaledHankelH1(0,0.3+0.4*%i,2) :: Matrix Complex Float
-
---       [0.2477 - 0.9492 %i  - 1.488 - 0.8166 %i]
-
--- nagScaledHankelH2 : (DF, CDF, Int) -> List CDF ;
-
-nagScaledHankelH2(6,3.1-1.6*%i,2) :: Matrix Complex Float
-
---       [7.05 + 6.052 %i  8.614 + 29.35 %i]
-
--- nagKelvinBer : DF -> DF ;
-
-nagKelvinBer(0.1) :: Float
-
---       1.0
-
-nagKelvinBer(1) :: Float
-
---       0.9844
-
-nagKelvinBer(2.5) :: Float
-
---       0.4
-
-nagKelvinBer(5) :: Float
-
---       - 6.23
-
-nagKelvinBer(10) :: Float
-
---       138.8
-
-nagKelvinBer(15) :: Float
-
---       - 2967.0
-
-nagKelvinBer(60) :: Float
-
--- ** ABNORMAL EXIT from NAG Library routine S19AAF: IFAIL =     1
--- ** NAG soft failure - control returned
--- 
---   Error signalled from user code:
---      An error was detected when calling the NAG Library routine 
---      S19AAF. The error number (IFAIL value) is 1, please consult the 
---      NAG manual via the Browser for diagnostic information.
-
-nagKelvinBer(-1) :: Float
-
---       0.9844
-
--- nagKelvinBei : DF -> DF ;
-
-nagKelvinBei(0.1) :: Float
-
---       0.0025
-
-nagKelvinBei(1) :: Float
-
---       0.2496
-
-nagKelvinBei(2.5) :: Float
-
---       1.457
-
-nagKelvinBei(5) :: Float
-
---       0.116
-
-nagKelvinBei(10) :: Float
-
---       56.37
-
-nagKelvinBei(15) :: Float
-
---       - 2953.0
-
-nagKelvinBei(60) :: Float
-
--- ** ABNORMAL EXIT from NAG Library routine S19ABF: IFAIL =     1
--- ** NAG soft failure - control returned
--- 
---   Error signalled from user code:
---      An error was detected when calling the NAG Library routine 
---      S19ABF. The error number (IFAIL value) is 1, please consult the 
---      NAG manual via the Browser for diagnostic information.
-
-nagKelvinBei(-1) :: Float
-
---       0.2496
-
--- nagKelvinKer : DF -> DF ;
-
-nagKelvinKer(0) :: Float
-
--- ** ABNORMAL EXIT from NAG Library routine S19ACF: IFAIL =     2
--- ** NAG soft failure - control returned
--- 
---   Error signalled from user code:
---      An error was detected when calling the NAG Library routine 
---      S19ACF. The error number (IFAIL value) is 2, please consult the 
---      NAG manual via the Browser for diagnostic information.
-
-nagKelvinKer(0.1) :: Float
-
---       2.42
-
-nagKelvinKer(1) :: Float
-
---       0.2867
-
-nagKelvinKer(2.5) :: Float
-
---       - 0.06969
-
-nagKelvinKer(5) :: Float
-
---       - 0.01151
-
-nagKelvinKer(10) :: Float
-
---       0.0001295
-
-nagKelvinKer(15) :: Float
-
---       - 0.1514 E -7
-
-nagKelvinKer(1100) :: Float
-
--- ** ABNORMAL EXIT from NAG Library routine S19ACF: IFAIL =     1
--- ** NAG soft failure - control returned
--- 
---   Error signalled from user code:
---      An error was detected when calling the NAG Library routine 
---      S19ACF. The error number (IFAIL value) is 1, please consult the 
---      NAG manual via the Browser for diagnostic information.
-
-nagKelvinKer(-1) :: Float
-
--- ** ABNORMAL EXIT from NAG Library routine S19ACF: IFAIL =     2
--- ** NAG soft failure - control returned
--- 
---   Error signalled from user code:
---      An error was detected when calling the NAG Library routine 
---      S19ACF. The error number (IFAIL value) is 2, please consult the 
---      NAG manual via the Browser for diagnostic information.
-
--- nagKelvinKei : DF -> DF ;
-
-nagKelvinKei(0) :: Float
-
---       - 0.7854
-
-nagKelvinKei(0.1) :: Float
-
---       - 0.7769
-
-nagKelvinKei(1) :: Float
-
---       - 0.495
-
-nagKelvinKei(2.5) :: Float
-
---       - 0.1107
-
-nagKelvinKei(5) :: Float
-
---       0.01119
-
-nagKelvinKei(10) :: Float
-
---       - 0.0003075
-
-nagKelvinKei(15) :: Float
-
---       0.000007963
-
-nagKelvinKei(1100) :: Float
-
--- ** ABNORMAL EXIT from NAG Library routine S19ADF: IFAIL =     1
--- ** NAG soft failure - control returned
--- 
---   Error signalled from user code:
---      An error was detected when calling the NAG Library routine 
---      S19ADF. The error number (IFAIL value) is 1, please consult the 
---      NAG manual via the Browser for diagnostic information.
-
-nagKelvinKei(-1) :: Float
-
--- ** ABNORMAL EXIT from NAG Library routine S19ADF: IFAIL =     2
--- ** NAG soft failure - control returned
--- 
---   Error signalled from user code:
---      An error was detected when calling the NAG Library routine 
---      S19ADF. The error number (IFAIL value) is 2, please consult the 
---      NAG manual via the Browser for diagnostic information.
-
-
--- nagFresnelS : DF -> DF ;
-
-nagFresnelS(0) :: Float
-
---       0.0
-
-nagFresnelS(0.5) :: Float
-
---       0.06473
-
-nagFresnelS(1) :: Float
-
---       0.4383
-
-nagFresnelS(2) :: Float
-
---       0.3434
-
-nagFresnelS(4) :: Float
-
---       0.4205
-
-nagFresnelS(5) :: Float
-
---       0.4992
-
-nagFresnelS(6) :: Float
-
---       0.447
-
-nagFresnelS(8) :: Float
-
---       0.4602
-
-nagFresnelS(10) :: Float
-
---       0.4682
-
-nagFresnelS(-1) :: Float
-
---       - 0.4383
-
-nagFresnelS(1000) :: Float
-
---       0.4997
-
--- nagFresnelC : DF -> DF ;
-
-nagFresnelC(0) :: Float
-
---       0.0
-
-nagFresnelC(0.5) :: Float
-
---       0.4923
-
-nagFresnelC(1) :: Float
-
---       0.7799
-
-nagFresnelC(2) :: Float
-
---       0.4883
-
-nagFresnelC(4) :: Float
-
---       0.4984
-
-nagFresnelC(5) :: Float
-
---       0.5636
-
-nagFresnelC(6) :: Float
-
---       0.4995
-
-nagFresnelC(8) :: Float
-
---       0.4998
-
-nagFresnelC(10) :: Float
-
---       0.4999
-
-nagFresnelC(-1) :: Float
-
---       - 0.7799
-
-nagFresnelC(1000) :: Float
-
---       0.5
-
--- nagEllipticIntegralRC : (DF, DF) -> DF ;
-
-nagEllipticIntegralRC(0.5,1) :: Float
-
---       1.111
-
-nagEllipticIntegralRC(1,1) :: Float
-
---       1.0
-
-nagEllipticIntegralRC(1.5,1) :: Float
-
---       0.9312
-
--- nagEllipticIntegralRD : (DF, DF, DF) -> DF ;
-
-nagEllipticIntegralRD(0.5,0.5,1) :: Float
-
---       1.479
-
-nagEllipticIntegralRD(0.5,1,1) :: Float
-
---       1.211
-
-nagEllipticIntegralRD(0.5,1.5,1) :: Float
-
---       1.061
-
-nagEllipticIntegralRD(1,1,1) :: Float
-
---       1.0
-
-nagEllipticIntegralRD(1,1.5,1) :: Float
-
---       0.8805
-
-nagEllipticIntegralRD(1.5,1.5,1) :: Float
-
---       0.7775
-
--- nagEllipticIntegralRF : (DF, DF, DF) -> DF ;
-
-nagEllipticIntegralRF(0.5,1,1.5) :: Float
-
---       1.028
-
-nagEllipticIntegralRF(1,1.5,2) :: Float
-
---       0.826
-
-nagEllipticIntegralRF(1.5,2,2.5) :: Float
-
---       0.7116
-
--- nagEllipticIntegralRJ : (DF, DF, DF, DF) -> DF ;
-
-nagEllipticIntegralRJ(0.5,0.5,0.5,2) :: Float
-
---       1.118
-
-nagEllipticIntegralRJ(0.5,0.5,1,2) :: Float
-
---       0.9221
-
-nagEllipticIntegralRJ(0.5,0.5,1.5,2) :: Float
-
---       0.8115
-
-nagEllipticIntegralRJ(0.5,1,1,2) :: Float
-
---       0.7671
-
-nagEllipticIntegralRJ(0.5,1,1.5,2) :: Float
-
---       0.6784
-
-nagEllipticIntegralRJ(0.5,1.5,1.5,2) :: Float
-
---       0.6017
-
-nagEllipticIntegralRJ(1,1,1,2) :: Float
-
---       0.6438
-
-nagEllipticIntegralRJ(1,1,1.5,2) :: Float
-
---       0.5722
-
-nagEllipticIntegralRJ(1,1.5,1.5,2) :: Float
-
---       0.5101
-
-nagEllipticIntegralRJ(1.5,1.5,1.5,2) :: Float
-
---       0.4561
-
-outputGeneral()
-
-output "End of tests"
-
-#endif
-
-@
-\section{License}
-<<license>>=
---Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
---All rights reserved.
---
---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.
---
---    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
---      names of its contributors may be used to endorse or promote products
---      derived from this software without specific prior written permission.
---
---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>>
-
--- NagSpecialFunctionsInterfacePackage
-
--- To test:
--- sed -ne '1,/^#if NeverAssertThis/d;/#endif/d;p' < nsfip.as > nsfip.input
--- axiom
--- )set nag host <some machine running nagd>
--- )r nsfip.input
-
-#unassert saturn
-
-#include "axiom.as"
-
-DF   ==> DoubleFloat ;
-CDF  ==> Complex DoubleFloat ;
-MCDF ==> Matrix Complex DoubleFloat ;
-INT  ==> Integer ;
-RSLT ==> Result ;
-SMBL ==> Symbol ;
-STRG ==> String ;
-
-<<NagSpecialFunctionsInterfacePackage>>
-@
-\eject
-\begin{thebibliography}{99}
-\bibitem{1} nothing
-\end{thebibliography}
-\end{document}
diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html
index f55d217..fdac097 100644
--- a/src/axiom-website/patches.html
+++ b/src/axiom-website/patches.html
@@ -2611,5 +2611,7 @@ books/bookvol10.5 add Volume 10.5 Axiom Numerics<br/>
 books/bookvol10.5 BLAS1 regress, help, and function documentation<br/>
 <a href="patches/20100403.02.tpd.patch">20100403.02.tpd.patch</a>
 Makefile clean books from src/algebra<br/>
+<a href="patches/20100404.01.tpd.patch">20100404.01.tpd.patch</a>
+books/bookvol5 change .spad.pamphlet to just .pamphlet<br/>
  </body>
 </html>
diff --git a/src/input/biquat.input.pamphlet b/src/input/biquat.input.pamphlet
index 659ef15..fdb4aff 100644
--- a/src/input/biquat.input.pamphlet
+++ b/src/input/biquat.input.pamphlet
@@ -660,7 +660,7 @@ So let's find out more about this domain:
 --R RewriteRule(Base: SetCategory,R: Join(Ring,PatternMatchable Base,OrderedSet,ConvertibleTo Pattern Base),F: Join(FunctionSpace R,PatternMatchable Base,ConvertibleTo Pattern Base))  is a domain constructor
 --R Abbreviation for RewriteRule is RULE 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for RULE 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for RULE 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?=? : (%,%) -> Boolean                coerce : Equation F -> %
diff --git a/src/input/chtheorem.input.pamphlet b/src/input/chtheorem.input.pamphlet
index ca27405..5137490 100644
--- a/src/input/chtheorem.input.pamphlet
+++ b/src/input/chtheorem.input.pamphlet
@@ -79,16 +79,16 @@ M:=matrix([[random()$D for i in 1..4] for j in 1..4])
 --R 
 --R
 --R        +      3     2                                           3     2+
---R        |    %A  + %A                0                1        %A  + %A |
+--I        |    %A  + %A                0                1        %A  + %A |
 --R        |                                                               |
 --R        |  3     2                  2                             2     |
---R        |%A  + %A  + %A + 1       %A  + %A            1         %A  + 1 |
+--I        |%A  + %A  + %A + 1       %A  + %A            1         %A  + 1 |
 --R   (2)  |                                                               |
 --R        |      3     2            3                                     |
---R        |    %A  + %A           %A  + %A + 1          0         %A + 1  |
+--I        |    %A  + %A           %A  + %A + 1          0         %A + 1  |
 --R        |                                                               |
 --R        |                      3     2             2             3     2|
---R        +        0           %A  + %A  + %A + 1  %A  + %A + 1  %A  + %A +
+--I        +        0           %A  + %A  + %A + 1  %A  + %A + 1  %A  + %A +
 --R         Type: Matrix FiniteFieldExtensionByPolynomial(PrimeField 2,?**4+?+1)
 --E 2
 
@@ -97,7 +97,7 @@ p:=characteristicPolynomial(M,y)
 --R 
 --R
 --R         4      2       3      2           2            3     2
---R   (3)  y  + (%A  + %A)y  + (%A  + %A + 1)y  + %A y + %A  + %A  + %A
+--I   (3)  y  + (%A  + %A)y  + (%A  + %A + 1)y  + %A y + %A  + %A  + %A
 --R     Type: Polynomial FiniteFieldExtensionByPolynomial(PrimeField 2,?**4+?+1)
 --E 3
 
@@ -116,48 +116,48 @@ sp:=map(z+->(squareMatrix$SM)diagonalMatrix([z,z,z,z]),p)
 --R
 --R   (5)
 --R          +  2                                   +
---R          |%A  + %A     0         0         0    |
+--I          |%A  + %A     0         0         0    |
 --R          |                                      |
 --R          |            2                         |
---R      4   |   0      %A  + %A     0         0    | 3
+--I      4   |   0      %A  + %A     0         0    | 3
 --R     y  + |                                      |y
 --R          |                      2               |
---R          |   0         0      %A  + %A     0    |
+--I          |   0         0      %A  + %A     0    |
 --R          |                                      |
 --R          |                                2     |
---R          +   0         0         0      %A  + %A+
+--I          +   0         0         0      %A  + %A+
 --R   + 
 --R     +  2                                                   +
---R     |%A  + %A + 1       0             0             0      |
+--I     |%A  + %A + 1       0             0             0      |
 --R     |                                                      |
 --R     |                2                                     |
---R     |     0        %A  + %A + 1       0             0      | 2
+--I     |     0        %A  + %A + 1       0             0      | 2
 --R     |                                                      |y
 --R     |                              2                       |
---R     |     0             0        %A  + %A + 1       0      |
+--I     |     0             0        %A  + %A + 1       0      |
 --R     |                                                      |
 --R     |                                            2         |
---R     +     0             0             0        %A  + %A + 1+
+--I     +     0             0             0        %A  + %A + 1+
 --R   + 
---R     +%A  0   0   0 +
+--I     +%A  0   0   0 +
 --R     |              |
---R     |0   %A  0   0 |
+--I     |0   %A  0   0 |
 --R     |              |y
---R     |0   0   %A  0 |
+--I     |0   0   %A  0 |
 --R     |              |
---R     +0   0   0   %A+
+--I     +0   0   0   %A+
 --R   + 
 --R     +  3     2                                                     +
---R     |%A  + %A  + %A        0               0               0       |
+--I     |%A  + %A  + %A        0               0               0       |
 --R     |                                                              |
 --R     |                  3     2                                     |
---R     |      0         %A  + %A  + %A        0               0       |
+--I     |      0         %A  + %A  + %A        0               0       |
 --R     |                                                              |
 --R     |                                  3     2                     |
---R     |      0               0         %A  + %A  + %A        0       |
+--I     |      0               0         %A  + %A  + %A        0       |
 --R     |                                                              |
 --R     |                                                  3     2     |
---R     +      0               0               0         %A  + %A  + %A+
+--I     +      0               0               0         %A  + %A  + %A+
 --RType: Polynomial SquareMatrix(4,FiniteFieldExtensionByPolynomial(PrimeField 2,?**4+?+1))
 --E 5
 
@@ -166,16 +166,16 @@ sm:=squareMatrix(M)$SM
 --R 
 --R
 --R        +      3     2                                           3     2+
---R        |    %A  + %A                0                1        %A  + %A |
+--I        |    %A  + %A                0                1        %A  + %A |
 --R        |                                                               |
 --R        |  3     2                  2                             2     |
---R        |%A  + %A  + %A + 1       %A  + %A            1         %A  + 1 |
+--I        |%A  + %A  + %A + 1       %A  + %A            1         %A  + 1 |
 --R   (6)  |                                                               |
 --R        |      3     2            3                                     |
---R        |    %A  + %A           %A  + %A + 1          0         %A + 1  |
+--I        |    %A  + %A           %A  + %A + 1          0         %A + 1  |
 --R        |                                                               |
 --R        |                      3     2             2             3     2|
---R        +        0           %A  + %A  + %A + 1  %A  + %A + 1  %A  + %A +
+--I        +        0           %A  + %A  + %A + 1  %A  + %A + 1  %A  + %A +
 --RType: SquareMatrix(4,FiniteFieldExtensionByPolynomial(PrimeField 2,?**4+?+1))
 --E 6
 
@@ -394,16 +394,16 @@ M:=matrix([[random()$D for i in 1..4] for j in 1..4])
 --R 
 --R
 --R        +                      2            2       3     2         +
---R        |        %A          %A  + %A     %A      %A  + %A  + %A + 1|
+--I        |        %A          %A  + %A     %A      %A  + %A  + %A + 1|
 --R        |                                                           |
 --R        |  3     2             3     2    3              3          |
---R        |%A  + %A  + %A + 1  %A  + %A   %A  + %A       %A  + %A     |
+--I        |%A  + %A  + %A + 1  %A  + %A   %A  + %A       %A  + %A     |
 --R   (2)  |                                                           |
 --R        |  3     2             3          3                         |
---R        |%A  + %A  + %A + 1  %A  + %A   %A  + %A          %A        |
+--I        |%A  + %A  + %A + 1  %A  + %A   %A  + %A          %A        |
 --R        |                                                           |
 --R        |     3                 2           2         3     2       |
---R        +   %A  + %A + 1      %A  + 1     %A        %A  + %A  + %A  +
+--I        +   %A  + %A + 1      %A  + 1     %A        %A  + %A  + %A  +
 --R                                                Type: Matrix FiniteField(2,4)
 --E 23
 
@@ -412,7 +412,7 @@ p:=characteristicPolynomial(M,y)
 --R 
 --R
 --R         4      3       3      3       2         3
---R   (3)  y  + (%A  + %A)y  + (%A  + %A)y  + y + %A  + 1
+--I   (3)  y  + (%A  + %A)y  + (%A  + %A)y  + y + %A  + 1
 --R                                            Type: Polynomial FiniteField(2,4)
 --E 24
 
@@ -430,40 +430,40 @@ sp:=map(z+->(squareMatrix$SM)diagonalMatrix([z,z,z,z]),p)
 --R
 --R   (5)
 --R          +  3                                   +
---R          |%A  + %A     0         0         0    |
+--I          |%A  + %A     0         0         0    |
 --R          |                                      |
 --R          |            3                         |
---R      4   |   0      %A  + %A     0         0    | 3
+--I      4   |   0      %A  + %A     0         0    | 3
 --R     y  + |                                      |y
 --R          |                      3               |
---R          |   0         0      %A  + %A     0    |
+--I          |   0         0      %A  + %A     0    |
 --R          |                                      |
 --R          |                                3     |
---R          +   0         0         0      %A  + %A+
+--I          +   0         0         0      %A  + %A+
 --R   + 
 --R     +  3                                   +
---R     |%A  + %A     0         0         0    |
+--I     |%A  + %A     0         0         0    |
 --R     |                                      |
 --R     |            3                         |
---R     |   0      %A  + %A     0         0    | 2
+--I     |   0      %A  + %A     0         0    | 2
 --R     |                                      |y  + y
 --R     |                      3               |
---R     |   0         0      %A  + %A     0    |
+--I     |   0         0      %A  + %A     0    |
 --R     |                                      |
 --R     |                                3     |
---R     +   0         0         0      %A  + %A+
+--I     +   0         0         0      %A  + %A+
 --R   + 
 --R     +  3                               +
---R     |%A  + 1     0        0        0   |
+--I     |%A  + 1     0        0        0   |
 --R     |                                  |
 --R     |           3                      |
---R     |   0     %A  + 1     0        0   |
+--I     |   0     %A  + 1     0        0   |
 --R     |                                  |
 --R     |                    3             |
---R     |   0        0     %A  + 1     0   |
+--I     |   0        0     %A  + 1     0   |
 --R     |                                  |
 --R     |                             3    |
---R     +   0        0        0     %A  + 1+
+--I     +   0        0        0     %A  + 1+
 --R                            Type: Polynomial SquareMatrix(4,FiniteField(2,4))
 --E 26
 
@@ -472,16 +472,16 @@ sm:=squareMatrix(M)$SM
 --R 
 --R
 --R        +                      2            2       3     2         +
---R        |        %A          %A  + %A     %A      %A  + %A  + %A + 1|
+--I        |        %A          %A  + %A     %A      %A  + %A  + %A + 1|
 --R        |                                                           |
 --R        |  3     2             3     2    3              3          |
---R        |%A  + %A  + %A + 1  %A  + %A   %A  + %A       %A  + %A     |
+--I        |%A  + %A  + %A + 1  %A  + %A   %A  + %A       %A  + %A     |
 --R   (6)  |                                                           |
 --R        |  3     2             3          3                         |
---R        |%A  + %A  + %A + 1  %A  + %A   %A  + %A          %A        |
+--I        |%A  + %A  + %A + 1  %A  + %A   %A  + %A          %A        |
 --R        |                                                           |
 --R        |     3                 2           2         3     2       |
---R        +   %A  + %A + 1      %A  + 1     %A        %A  + %A  + %A  +
+--I        +   %A  + %A + 1      %A  + 1     %A        %A  + %A  + %A  +
 --R                                       Type: SquareMatrix(4,FiniteField(2,4))
 --E 27
 
diff --git a/src/input/grpthry.input.pamphlet b/src/input/grpthry.input.pamphlet
index 273b15c..4228770 100644
--- a/src/input/grpthry.input.pamphlet
+++ b/src/input/grpthry.input.pamphlet
@@ -122,7 +122,7 @@ member? ( y , g2 )
 --R PermutationGroup S: SetCategory  is a domain constructor
 --R Abbreviation for PermutationGroup is PERMGRP 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for PERMGRP 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for PERMGRP 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?<? : (%,%) -> Boolean                ?<=? : (%,%) -> Boolean
diff --git a/src/input/gstbl.input.pamphlet b/src/input/gstbl.input.pamphlet
index 1693367..9265425 100644
--- a/src/input/gstbl.input.pamphlet
+++ b/src/input/gstbl.input.pamphlet
@@ -23,14 +23,8 @@
 --S 1 of 1
 patrons: GeneralSparseTable(String, Integer, KeyedAccessFile(Integer), 0) := table() ;
 --R 
---R 
---RDaly Bug
---R   >> Error detected within library code:
---R   File is not readable
---I   "kaf1414.sdata"
---R
---R   Continuing to read the file...
 --R
+--R           Type: GeneralSparseTable(String,Integer,KeyedAccessFile Integer,0)
 --E 1
 )spool 
 )lisp (bye)
diff --git a/src/input/unittest1.input.pamphlet b/src/input/unittest1.input.pamphlet
index c9e6d5f..cb07077 100644
--- a/src/input/unittest1.input.pamphlet
+++ b/src/input/unittest1.input.pamphlet
@@ -473,7 +473,7 @@ Unit test the user level commands
 --R SquareMatrix(ndim: NonNegativeInteger,R: Ring)  is a domain constructor
 --R Abbreviation for SquareMatrix is SQMATRIX 
 --R This constructor is not exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for SQMATRIX 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for SQMATRIX 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?*? : (R,%) -> %                      ?*? : (%,R) -> %
@@ -577,7 +577,7 @@ Unit test the user level commands
 --R SquareMatrix(ndim: NonNegativeInteger,R: Ring)  is a domain constructor
 --R Abbreviation for SquareMatrix is SQMATRIX 
 --R This constructor is exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for SQMATRIX 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for SQMATRIX 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?*? : (R,%) -> %                      ?*? : (%,R) -> %
@@ -705,7 +705,7 @@ Unit test the user level commands
 --R SquareMatrix(ndim: NonNegativeInteger,R: Ring)  is a domain constructor
 --R Abbreviation for SquareMatrix is SQMATRIX 
 --R This constructor is not exposed in this frame.
---R Issue )edit bookvol10.3.spad.pamphlet to see algebra source code for SQMATRIX 
+--R Issue )edit bookvol10.3.pamphlet to see algebra source code for SQMATRIX 
 --R
 --R------------------------------- Operations --------------------------------
 --R ?*? : (R,%) -> %                      ?*? : (%,R) -> %
