diff --git a/books/bookvol10.4.pamphlet b/books/bookvol10.4.pamphlet
index a72d13e..e26e653 100644
--- a/books/bookvol10.4.pamphlet
+++ b/books/bookvol10.4.pamphlet
@@ -3641,6 +3641,120 @@ AnyFunctions1(S:Type): with
@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{package API ApplicationProgramInterface}
+<>=
+)sys rm -f ApplicationProgramInterface.output
+)spool ApplicationProgramInterface.output
+)set message test on
+)set message auto off
+)clear all
+--S 1 of 3
+getDomains 'Collection
+--R
+--R (1)
+--R {AssociationList, Bits, CharacterClass, DataList, EqTable, FlexibleArray,
+--R GeneralPolynomialSet, GeneralSparseTable, GeneralTriangularSet, HashTable,
+--R IndexedBits, IndexedFlexibleArray, IndexedList, IndexedOneDimensionalArray,
+--R IndexedString, IndexedVector, InnerTable, KeyedAccessFile, Library, List,
+--R ListMultiDictionary, Multiset, OneDimensionalArray, Point, PrimitiveArray,
+--R RegularChain, RegularTriangularSet, Result, RoutinesTable, Set,
+--R SparseTable, SquareFreeRegularTriangularSet, Stream, String, StringTable,
+--R Table, Vector, WuWenTsunTriangularSet}
+--R Type: Set Symbol
+--E 1
+
+--S 2 of 3
+difference(getDomains 'IndexedAggregate,getDomains 'Collection)
+--R
+--R (2)
+--R {DirectProduct, DirectProductMatrixModule, DirectProductModule,
+--R HomogeneousDirectProduct, OrderedDirectProduct,
+--R SplitHomogeneousDirectProduct}
+--R Type: Set Symbol
+--E 2
+
+--S 3 of 3
+)show ApplicationProgramInterface
+--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.pamphlet to see algebra source code for API
+--R
+--R------------------------------- Operations --------------------------------
+--R getDomains : Symbol -> Set Symbol
+--R
+--E 3
+)spool
+)lisp (bye)
+@
+<>=
+====================================================================
+ApplicationProgramInterface examples
+====================================================================
+
+The ApplicationProgramInterface exposes Axiom internal functions
+which might be useful for understanding, debugging, or creating
+tools.
+
+The getDomains function takes the name of a category and returns
+a set of domains which inherit from that category:
+
+ getDomains 'Collection
+
+ {AssociationList, Bits, CharacterClass, DataList, EqTable, FlexibleArray,
+ GeneralPolynomialSet, GeneralSparseTable, GeneralTriangularSet, HashTable,
+ IndexedBits, IndexedFlexibleArray, IndexedList, IndexedOneDimensionalArray,
+ IndexedString, IndexedVector, InnerTable, KeyedAccessFile, Library, List,
+ ListMultiDictionary, Multiset, OneDimensionalArray, Point, PrimitiveArray,
+ RegularChain, RegularTriangularSet, Result, RoutinesTable, Set,
+ SparseTable, SquareFreeRegularTriangularSet, Stream, String, StringTable,
+ Table, Vector, WuWenTsunTriangularSet}
+ Type: Set Symbol
+
+This can be used to form the set-difference of two categories:
+
+ difference(getDomains 'IndexedAggregate, getDomains 'Collection)
+
+ {DirectProduct, DirectProductMatrixModule, DirectProductModule,
+ HomogeneousDirectProduct, OrderedDirectProduct,
+ SplitHomogeneousDirectProduct}
+ Type: Set Symbol
+
+@
+\pagehead{ApplicationProgramInterface}{API}
+\pagepic{ps/v104applicationprograminterface.ps}{API}{1.00}
+
+{\bf Exports:}\\
+\begin{tabular}{ll}
+\end{tabular}
+
+<>=
+)abbrev package API ApplicationProgramInterface
+++ Author: Timothy Daly, Martin Rubey
+++ Date Created: 3 March 2009
+++ Date Last Updated: 3 March 2009
+++ Description: This package contains useful functions that
+++ expose Axiom system internals
+ApplicationProgramInterface(): Exports == Implementation where
+ Exports ==> with
+ getDomains : Symbol -> Set Symbol
+ ++ getDomains takes a category and returns the list of domains
+ ++ that have that category
+ ++
+ ++X getDomains 'IndexedAggregate
+ Implementation ==> add
+ getDomains(cat:Symbol):Set(Symbol) ==
+ set [symbol car first destruct a _
+ for a in (destruct domainsOf(cat,NIL$Lisp)$Lisp)::List(SExpression)]
+
+@
+<>=
+"API" [color="#FF4488",href="bookvol10.4.pdf#nameddest=APPRULE"]
+"ORDSET" [color="#4488FF",href="bookvol10.2.pdf#nameddest=ORDSET"]
+"API" -> "ORDSET"
+
+@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{package APPRULE ApplyRules}
\pagehead{ApplyRules}{APPRULE}
\pagepic{ps/v104applyrules.ps}{APPRULE}{1.00}
@@ -63433,6 +63547,685 @@ NagPartialDifferentialEquationsPackage(): Exports == Implementation where
@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{package NAGC02 NagPolynomialRootsPackage}
+<>=
+
+ C02(3NAG) Foundation Library (12/10/92) C02(3NAG)
+
+ C02 -- Zeros of Polynomials Introduction -- C02
+ Chapter C02
+ Zeros of Polynomials
+
+ 1. Scope of the Chapter
+
+ This chapter is concerned with computing the zeros of a
+ polynomial with real or complex coefficients.
+
+ 2. Background to the Problems
+
+ Let f(z) be a polynomial of degree n with complex coefficients
+ a :
+ i
+
+ n n-1 n-2
+ f(z)==a z +a z +a z +...+a z+a , a /=0.
+ 0 1 2 n-1 n 0
+
+ A complex number z is called a zero of f(z) (or equivalently a
+ 1
+ root of the equation f(z)=0), if:
+
+ f(z )=0.
+ 1
+
+ If z is a zero, then f(z) can be divided by a factor (z-z ):
+ 1 1
+
+ f(z)=(z-z )f (z) (1)
+ 1 1
+
+ where f (z) is a polynomial of degree n-1. By the Fundamental
+ 1
+ Theorem of Algebra, a polynomial f(z) always has a zero, and so
+ the process of dividing out factors (z-z ) can be continued until
+ i
+ we have a complete factorization of f(z)
+
+ f(z)==a (z-z )(z-z )...(z-z ).
+ 0 1 2 n
+
+ Here the complex numbers z ,z ,...,z are the zeros of f(z); they
+ 1 2 n
+ may not all be distinct, so it is sometimes more convenient to
+ write:
+
+ m m m
+ 1 2 k
+ f(z)==a (z-z ) (z-z ) ...(z-z ) , k<=n,
+ 0 1 2 k
+
+ with distinct zeros z ,z ,...,z and multiplicities m >=1. If
+ 1 2 k i
+ m =1, z is called a single zero, if m >1, z is called a
+ i i i i
+ multiple or repeated zero; a multiple zero is also a zero of the
+ derivative of f(z).
+
+ If the coefficients of f(z) are all real, then the zeros of f(z)
+ are either real or else occur as pairs of conjugate complex
+ numbers x+iy and x-iy. A pair of complex conjugate zeros are the
+ 2
+ zeros of a quadratic factor of f(z), (z +rz+s), with real
+ coefficients r and s.
+
+ Mathematicians are accustomed to thinking of polynomials as
+ pleasantly simple functions to work with. However the problem of
+ numerically computing the zeros of an arbitrary polynomial is far
+ from simple. A great variety of algorithms have been proposed, of
+ which a number have been widely used in practice; for a fairly
+ comprehensive survey, see Householder [1]. All general algorithms
+ are iterative. Most converge to one zero at a time; the
+ corresponding factor can then be divided out as in equation (1)
+ above - this process is called deflation or, loosely, dividing
+ out the zero - and the algorithm can be applied again to the
+ polynomial f (z). A pair of complex conjugate zeros can be
+ 1
+ divided out together - this corresponds to dividing f(z) by a
+ quadratic factor.
+
+ Whatever the theoretical basis of the algorithm, a number of
+ practical problems arise: for a thorough discussion of some of
+ them see Peters and Wilkinson [2] and Wilkinson [3]. The most
+ elementary point is that, even if z is mathematically an exact
+ 1
+ zero of f(z), because of the fundamental limitations of computer
+ arithmetic the computed value of f(z ) will not necessarily be
+ 1
+ exactly 0.0. In practice there is usually a small region of
+ values of z about the exact zero at which the computed value of
+ f(z) becomes swamped by rounding errors. Moreover in many
+ algorithms this inaccuracy in the computed value of f(z) results
+ in a similar inaccuracy in the computed step from one iterate to
+ the next. This limits the precision with which any zero can be
+ computed. Deflation is another potential cause of trouble, since,
+ in the notation of equation (1), the computed coefficients of
+ f (z) will not be completely accurate, especially if z is not an
+ 1 1
+ exact zero of f(z); so the zeros of the computed f (z) will
+ 1
+ deviate from the zeros of f(z).
+
+ A zero is called ill-conditioned if it is sensitive to small
+ changes in the coefficients of the polynomial. An ill-conditioned
+ zero is likewise sensitive to the computational inaccuracies just
+ mentioned. Conversely a zero is called well-conditioned if it is
+ comparatively insensitive to such perturbations. Roughly speaking
+ a zero which is well separated from other zeros is well-
+ conditioned, while zeros which are close together are ill-
+ conditioned, but in talking about 'closeness' the decisive factor
+ is not the absolute distance between neighbouring zeros but their
+ ratio: if the ratio is close to 1 the zeros are ill-conditioned.
+ In particular, multiple zeros are ill-conditioned. A multiple
+ zero is usually split into a cluster of zeros by perturbations in
+ the polynomial or computational inaccuracies.
+
+ 2.1. References
+
+ [1] Householder A S (1970) The Numerical Treatment of a Single
+ Nonlinear Equation. McGraw-Hill.
+
+ [2] Peters G and Wilkinson J H (1971) Practical Problems Arising
+ in the Solution of Polynomial Equations. J. Inst. Maths
+ Applics. 8 16--35.
+
+ [3] Wilkinson J H (1963) Rounding Errors in Algebraic Processes,
+ Chapter 2. HMSO.
+
+ 3. Recommendations on Choice and Use of Routines
+
+ 3.1. Discussion
+
+ Two routines are available: C02AFF for polynomials with complex
+ coefficients and C02AGF for polynomials with real coefficients.
+
+ C02AFF and C02AGF both use a variant of Laguerre's Method due to
+ BT Smith to calculate each zero until the degree of the deflated
+ polynomial is less than 3, whereupon the remaining zeros are
+ obtained using the 'standard' closed formulae for a quadratic or
+ linear equation.
+
+ The accuracy of the roots will depend on how ill-conditioned they
+ are. Peters and Wilkinson [2] describe techniques for estimating
+ the errors in the zeros after they have been computed.
+
+ 3.2. Index
+
+ Zeros of a complex polynomial C02AFF
+ Zeros of a real polynomial C02AGF
+
+
+ C02 -- Zeros of Polynomials Contents -- C02
+ Chapter C02
+
+ Zeros of Polynomials
+
+ C02AFF All zeros of complex polynomial, modified Laguerre method
+
+ C02AGF All zeros of real polynomial, modified Laguerre method
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C02AFF(3NAG) Foundation Library (12/10/92) C02AFF(3NAG)
+
+ C02 -- Zeros of Polynomials C02AFF
+ C02AFF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C02AFF finds all the roots of a complex polynomial equation,
+ using a variant of Laguerre's Method.
+
+ 2. Specification
+
+ SUBROUTINE C02AFF (A, N, SCALE, Z, W, IFAIL)
+ INTEGER N, IFAIL
+ DOUBLE PRECISION A(2,N+1), Z(2,N), W(4*(N+1))
+ LOGICAL SCALE
+
+ 3. Description
+
+ The routine attempts to find all the roots of the nth degree
+ complex polynomial equation
+
+ n n-1 n-2
+ P(z)=a z +a z +a z +...+a z+a =0.
+ 0 1 2 n-1 n
+
+ The roots are located using a modified form of Laguerre's Method,
+ originally proposed by Smith [2].
+
+ The method of Laguerre [3] can be described by the iterative
+ scheme
+
+ -n*P(z )
+ k
+ L(z )=z -z = ----------------,
+ k k+1 k
+ P'(z )+- /H(z )
+ k \/ k
+
+ 2
+ where H(z )=(n-1)*[(n-1)*(P'(z )) -n*P(z )P''(z )], and z is
+ k k k k 0
+ specified.
+
+ The sign in the denominator is chosen so that the modulus of the
+ Laguerre step at z , viz. |L(z )|, is as small as possible. The
+ k k
+ method can be shown to be cubically convergent for isolated roots
+ (real or complex) and linearly convergent for multiple roots.
+ The routine generates a sequence of iterates z , z , z ,..., such
+ 1 2 3
+ that |P(z )|<|P(z )| and ensures that z +L(z ) 'roughly'
+ k+1 k k+1 k+1
+ lies inside a circular region of radius |F| about z known to
+ k
+ contain a zero of P(z); that is, |L(z )|<=|F|, where F denotes
+ k+1
+ the Fejer bound (see Marden [1]) at the point z . Following Smith
+ k
+ [2], F is taken to be min(B,1.445*n*R), where B is an upper bound
+ for the magnitude of the smallest zero given by
+
+ 1/n
+ B=1.0001*min(\/n*L(z ),|r |,|a /a | ),
+ k 1 n 0
+
+ r is the zero X of smaller magnitude of the quadratic equation
+ 1
+
+ 2
+ 2(P''(z )/(2*n*(n-1)))X +2(P'(z )/n)X+P(z )=0
+ k k k
+
+ and the Cauchy lower bound R for the smallest zero is computed
+ (using Newton's Method) as the positive root of the polynomial
+ equation
+
+ n n-1 n-2
+ |a |z +|a |z +|a |z +...+|a |z-|a |=0.
+ 0 1 2 n-1 n
+
+ Starting from the origin, successive iterates are generated
+ according to the rule z =z +L(z ) for k = 1,2,3,... and L(z )
+ k+1 k k k
+ is 'adjusted' so that |P(z )|<|P(z )| and |L(z )|<=|F|. The
+ k+1 k k+1
+ iterative procedure terminates if P(z ) is smaller in absolute
+ k+1
+ value than the bound on the rounding error in P(z ) and the
+ k+1
+ current iterate z =z is taken to be a zero of P(z). The
+ p k+1
+ ~
+ deflated polynomial P(z)=P(z)/(z-z ) of degree n-1 is then
+ p
+ formed, and the above procedure is repeated on the deflated
+ polynomial until n<3, whereupon the remaining roots are obtained
+ via the 'standard' closed formulae for a linear (n = 1) or
+ quadratic (n = 2) equation.
+
+ To obtain the roots of a quadratic polynomial, C02AHF(*) can be
+ used.
+
+ 4. References
+
+ [1] Marden M (1966) Geometry of Polynomials. Mathematical
+ Surveys. 3 Am. Math. Soc., Providence, RI.
+
+ [2] Smith B T (1967) ZERPOL: A Zero Finding Algorithm for
+ Polynomials Using Laguerre's Method. Technical Report.
+ Department of Computer Science, University of Toronto,
+ Canada.
+
+ [3] Wilkinson J H (1965) The Algebraic Eigenvalue Problem.
+ Clarendon Press.
+
+ 5. Parameters
+
+ 1: A(2,N+1) -- DOUBLE PRECISION array Input
+ On entry: if A is declared with bounds (2,0:N), then A(1,i)
+ and A(2,i) must contain the real and imaginary parts of a
+ i
+ n-i
+ (i.e., the coefficient of z ), for i=0,1,...,n.
+ Constraint: A(1,0) /= 0.0 or A(2,0) /= 0.0.
+
+ 2: N -- INTEGER Input
+ On entry: the degree of the polynomial, n. Constraint: N >=
+ 1.
+
+ 3: SCALE -- LOGICAL Input
+ On entry: indicates whether or not the polynomial is to be
+ scaled. See Section 8 for advice on when it may be
+ preferable to set SCALE = .FALSE. and for a description of
+ the scaling strategy. Suggested value: SCALE = .TRUE..
+
+ 4: Z(2,N) -- DOUBLE PRECISION array Output
+ On exit: the real and imaginary parts of the roots are
+ stored in Z(1,i) and Z(2,i) respectively, for i=1,2,...,n.
+
+ 5: W(4*(N+1)) -- DOUBLE PRECISION array Workspace
+
+ 6: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ If on entry IFAIL = 0 or -1, explanatory error messages are
+ output on the current error message unit (as defined by X04AAF).
+
+ IFAIL= 1
+ On entry A(1,0) = 0.0 and A(2,0) = 0.0,
+
+ or N < 1.
+
+ IFAIL= 2
+ The iterative procedure has failed to converge. This error
+ is very unlikely to occur. If it does, please contact NAG
+ immediately, as some basic assumption for the arithmetic has
+ been violated. See also Section 8.
+
+ IFAIL= 3
+ Either overflow or underflow prevents the evaluation of P(z)
+ near some of its zeros. This error is very unlikely to
+ occur. If it does, please contact NAG immediately. See also
+ Section 8.
+
+ 7. Accuracy
+
+ All roots are evaluated as accurately as possible, but because of
+ the inherent nature of the problem complete accuracy cannot be
+ guaranteed.
+
+ 8. Further Comments
+
+ If SCALE = .TRUE., then a scaling factor for the coefficients is
+ chosen as a power of the base B of the machine so that the
+ EMAX-P
+ largest coefficient in magnitude approaches THRESH = B .
+ Users should note that no scaling is performed if the largest
+ coefficient in magnitude exceeds THRESH, even if SCALE = .TRUE..
+ (For definition of B, EMAX and P see the Chapter Introduction
+ X02.)
+
+ However, with SCALE = .TRUE., overflow may be encountered when
+ the input coefficients a ,a ,a ,...,a vary widely in magnitude,
+ 0 1 2 n
+ (4*P)
+ particularly on those machines for which B overflows. In
+ such cases, SCALE should be set to .FALSE. and the coefficients
+ scaled so that the largest coefficient in magnitude does not
+ (EMAX-2*P)
+ exceed B .
+
+ Even so, the scaling strategy used in C02AFF is sometimes
+ insufficient to avoid overflow and/or underflow conditions. In
+ such cases, the user is recommended to scale the independent
+ variable (z) so that the disparity between the largest and
+ smallest coefficient in magnitude is reduced. That is, use the
+ routine to locate the zeros of the polynomial d*P(cz) for some
+ suitable values of c and d. For example, if the original
+ -100 100 20 -10
+ polynomial was P(z)=2 i+2 z , then choosing c=2 and
+ 100 20
+ d=2 , for instance, would yield the scaled polynomial i+z ,
+ which is well-behaved relative to overflow and underflow and has
+ 10
+ zeros which are 2 times those of P(z).
+
+ If the routine fails with IFAIL = 2 or 3, then the real and
+ imaginary parts of any roots obtained before the failure occurred
+ are stored in Z in the reverse order in which they were found.
+ Let n denote the number of roots found before the failure
+ R
+ occurred. Then Z(1,n) and Z(2,n) contain the real and imaginary
+ parts of the 1st root found, Z(1,n-1) and Z(2,n-1) contain the
+ real and imaginary parts of the 2nd root found, ..., Z(1,n ) and
+ R
+ Z(2,n ) contain the real and imaginary parts of the n th root
+ R R
+ found. After the failure has occurred, the remaining 2*(n-n )
+ R
+ elements of Z contain a large negative number (equal to
+
+
+ -1/(X02AMF().\/2)).
+
+ 9. Example
+
+ 5 4 3 2
+ To find the roots of the polynomial a z +a z +a z +a z +a z+a =0,
+ 0 1 2 3 4 5
+ where a =(5.0+6.0i), a =(30.0+20.0i), a =-(0.2+6.0i),
+ 0 1 2
+ a =(50.0+100000.0i), a =-(2.0-40.0i) and a =(10.0+1.0i).
+ 3 4 5
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C02AGF(3NAG) Foundation Library (12/10/92) C02AGF(3NAG)
+
+ C02 -- Zeros of Polynomials C02AGF
+ C02AGF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C02AGF finds all the roots of a real polynomial equation, using a
+ variant of Laguerre's Method.
+
+ 2. Specification
+
+ SUBROUTINE C02AGF (A, N, SCALE, Z, W, IFAIL)
+ INTEGER N, IFAIL
+ DOUBLE PRECISION A(N+1), Z(2,N), W(2*(N+1))
+ LOGICAL SCALE
+
+ 3. Description
+
+ The routine attempts to find all the roots of the nth degree real
+ polynomial equation
+
+ n n-1 n-2
+ P(z)=a z +a z +a z +...+a z+a =0.
+ 0 1 2 n-1 n
+
+ The roots are located using a modified form of Laguerre's Method,
+ originally proposed by Smith [2].
+
+ The method of Laguerre [3] can be described by the iterative
+ scheme
+
+ -n*P(z )
+ k
+ L(z )=z -z = ----------------,
+ k k+1 k
+ P'(z )+- /H(z )
+ k \/ k
+
+ 2
+ where H(z )=(n-1)*[(n-1)*(P'(z )) -n*P(z )P''(z )], and z is
+ k k k k 0
+ specified.
+
+ The sign in the denominator is chosen so that the modulus of the
+ Laguerre step at z , viz. |L(z )|, is as small as possible. The
+ k k
+ method can be shown to be cubically convergent for isolated roots
+ (real or complex) and linearly convergent for multiple roots.
+ The routine generates a sequence of iterates z , z , z ,..., such
+ 1 2 3
+ that |P(z +1)|<|P(z )| and ensures that z +L(z ) 'roughly'
+ k k k+1 k+1
+ lies inside a circular region of radius |F| about z known to
+ k
+ contain a zero of P(z); that is, |L(z )|<=|F|, where F denotes
+ k+1
+ the Fejer bound (see Marden [1]) at the point z . Following Smith
+ k
+ [2], F is taken to be min(B,1.445*n*R), where B is an upper bound
+ for the magnitude of the smallest zero given by
+
+ 1/n
+ B=1.0001*min(\/n*L(z ),|r |,|a /a | ),
+ k 1 n 0
+
+ r is the zero X of smaller magnitude of the quadratic equation
+ 1
+
+ 2
+ 2(P''(z )/(2*n*(n-1)))X +2(P'(z )/n)X+P(z )=0
+ k k k
+
+ and the Cauchy lower bound R for the smallest zero is computed
+ (using Newton's Method) as the positive root of the polynomial
+ equation
+
+ n n-1 n-2
+ |a |z +|a |z +|a |z +...+|a |z-|a |=0.
+ 0 1 2 n-1 n
+
+ Starting from the origin, successive iterates are generated
+ according to the rule z =z +L(z ) for k=1,2,3,... and L(z ) is
+ k+1 k k k
+ k+1 k k+1
+ iterative procedure terminates if P(z ) is smaller in absolute
+ k+1
+ value than the bound on the rounding error in P(z ) and the
+ k+1
+ current iterate z =z is taken to be a zero of P(z) (as is its
+ p k-1
+
+
+ conjugate z if z is complex). The deflated polynomial
+ p p
+ ~
+ P(z)=P(z)/(z-z ) of degree n-1 if z is real
+ p p
+ ~
+ (P(z)=P(z)/((z-z )(z-z )) of degree n-2 if z is complex) is then
+ p p p
+ formed, and the above procedure is repeated on the deflated
+ polynomial until n<3, whereupon the remaining roots are obtained
+ via the 'standard' closed formulae for a linear (n = 1) or
+ quadratic (n = 2) equation.
+
+ To obtain the roots of a quadratic polynomial, C02AJF(*) can be
+ used.
+
+ 4. References
+
+ [1] Marden M (1966) Geometry of Polynomials. Mathematical
+ Surveys. 3 Am. Math. Soc., Providence, RI.
+
+ [2] Smith B T (1967) ZERPOL: A Zero Finding Algorithm for
+ Polynomials Using Laguerre's Method. Technical Report.
+ Department of Computer Science, University of Toronto,
+ Canada.
+
+ [3] Wilkinson J H (1965) The Algebraic Eigenvalue Problem.
+ Clarendon Press.
+
+ 5. Parameters
+
+ 1: A(N+1) -- DOUBLE PRECISION array Input
+ On entry: if A is declared with bounds (0:N), then A(i)
+ n-i
+ must contain a (i.e., the coefficient of z ), for
+ i
+ i=0,1,...,n. Constraint: A(0) /= 0.0.
+
+ 2: N -- INTEGER Input
+ On entry: the degree of the polynomial, n. Constraint: N >=
+ 1.
+
+ 3: SCALE -- LOGICAL Input
+ On entry: indicates whether or not the polynomial is to be
+ scaled. See Section 8 for advice on when it may be
+ preferable to set SCALE = .FALSE. and for a description of
+ the scaling strategy. Suggested value: SCALE = .TRUE..
+
+ 4: Z(2,N) -- DOUBLE PRECISION array Output
+ On exit: the real and imaginary parts of the roots are
+ stored in Z(1,i) and Z(2,i) respectively, for i=1,2,...,n.
+ Complex conjugate pairs of roots are stored in consecutive
+ pairs of elements of Z; that is, Z(1,i+1) = Z(1,i) and
+ Z(2,i+1)=-Z(2,i).
+
+ 5: W(2*(N+1)) -- DOUBLE PRECISION array Workspace
+
+ 6: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ If on entry IFAIL = 0 or -1, explanatory error messages are
+ output on the current error message unit (as defined by X04AAF).
+
+ IFAIL= 1
+ On entry A(0) = 0.0,
+
+ or N < 1.
+
+ IFAIL= 2
+ The iterative procedure has failed to converge. This error
+ is very unlikely to occur. If it does, please contact NAG
+ immediately, as some basic assumption for the arithmetic has
+ been violated. See also Section 8.
+
+ IFAIL= 3
+ Either overflow or underflow prevents the evaluation of P(z)
+ near some of its zeros. This error is very unlikely to
+ occur. If it does, please contact NAG immediately. See also
+ Section 8.
+
+ 7. Accuracy
+
+ All roots are evaluated as accurately as possible, but because of
+ the inherent nature of the problem complete accuracy cannot be
+ guaranteed.
+
+ 8. Further Comments
+
+ If SCALE = .TRUE., then a scaling factor for the coefficients is
+ chosen as a power of the base B of the machine so that the
+ EMAX-P
+ largest coefficient in magnitude approaches THRESH = B .
+ Users should note that no scaling is performed if the largest
+ coefficient in magnitude exceeds THRESH, even if SCALE = .TRUE..
+ (For definition of B, EMAX and P see the Chapter Introduction
+ X02.)
+
+ However, with SCALE = .TRUE., overflow may be encountered when
+ the input coefficients a ,a ,a ,...,a vary widely in magnitude,
+ 0 1 2 n
+ (4*P)
+ particularly on those machines for which B overflows. In
+ such cases, SCALE should be set to .FALSE. and the coefficients
+ scaled so that the largest coefficient in magnitude does not
+ (EMAX-2*P)
+ exceed B .
+
+ Even so, the scaling strategy used in C02AGF is sometimes
+ insufficient to avoid overflow and/or underflow conditions. In
+ such cases, the user is recommended to scale the independent
+ variable (z) so that the disparity between the largest and
+ smallest coefficient in magnitude is reduced. That is, use the
+ routine to locate the zeros of the polynomial d*P(cz) for some
+ suitable values of c and d. For example, if the original
+ -100 100 20 -10
+ polynomial was P(z)=2 +2 z , then choosing c=2 and
+ 100 20
+ d=2 , for instance, would yield the scaled polynomial 1+z ,
+ which is well-behaved relative to overflow and underflow and has
+ 10
+ zeros which are 2 times those of P(z).
+
+ If the routine fails with IFAIL = 2 or 3, then the real and
+ imaginary parts of any roots obtained before the failure occurred
+ are stored in Z in the reverse order in which they were found.
+ Let n denote the number of roots found before the failure
+ R
+ occurred. Then Z(1,n) and Z(2,n) contain the real and imaginary
+ parts of the 1st root found, Z(1,n-1) and Z(2,n-1) contain the
+ real and imaginary parts of the 2nd root found, ..., Z(1,n ) and
+ R
+ Z(2,n ) contain the real and imaginary parts of the n th root
+ R R
+ found. After the failure has occurred, the remaining 2*(n-n )
+ R
+ elements of Z contain a large negative number (equal to
+
+
+ -1/(X02AMF().\/2)).
+
+ 9. Example
+
+ To find the roots of the 5th degree polynomial
+ 5 4 3 2
+ z +2z +3z +4z +5z+6=0.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+@
\pagehead{NagPolynomialRootsPackage}{NAGC02}
\pagepic{ps/v104nagpolynomialrootspackage.ps}{NAGC02}{1.00}
@@ -63525,6 +64318,831 @@ NagPolynomialRootsPackage(): Exports == Implementation where
@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{package NAGC05 NagRootFindingPackage}
+<>=
+
+ C05(3NAG) Foundation Library (12/10/92) C05(3NAG)
+
+ C05 -- Roots of One or More Transcendental Equations
+ Introduction -- C05
+ Chapter C05
+ Roots of One or More Transcendental Equations
+
+ 1. Scope of the Chapter
+
+ This chapter is concerned with the calculation of real zeros of
+ continuous real functions of one or more variables. (Complex
+ equations must be expressed in terms of the equivalent larger
+ system of real equations.)
+
+ 2. Background to the Problems
+
+ The chapter divides naturally into two parts.
+
+ 2.1. A Single Equation
+
+ The first deals with the real zeros of a real function of a
+ single variable f(x).
+
+ At present, there is only one routine with a simple calling
+ sequence. This routine assumes that the user can determine an
+ initial interval [a,b] within which the desired zero lies, that
+ is f(a)*f(b)<0, and outside which all other zeros lie. The
+ routine then systematically subdivides the interval to produce a
+ final interval containing the zero. This final interval has a
+ length bounded by the user's specified error requirements; the
+ end of the interval where the function has smallest magnitude is
+ returned as the zero. This routine is guaranteed to converge to a
+ simple zero of the function. (Here we define a simple zero as a
+ zero corresponding to a sign-change of the function.) The
+ algorithm used is due to Bus and Dekker.
+
+ 2.2. Systems of Equations
+
+ The routines in the second part of this chapter are designed to
+ solve a set of nonlinear equations in n unknowns
+
+
+ T
+ f (x)=0, i=1,2,...,n, x=(x ,x ,...,x ) (1)
+ i 1 2 n
+
+ where T stands for transpose.
+
+ It is assumed that the functions are continuous and
+ differentiable so that the matrix of first partial derivatives of
+ the functions, the Jacobian matrix J (x)=ddf /ddx evaluated at
+ ij i j
+ the point x, exists, though it may not be possible to calculate
+ it directly.
+
+ The functions f must be independent, otherwise there will be an
+ i
+ infinity of solutions and the methods will fail. However, even
+ when the functions are independent the solutions may not be
+ unique. Since the methods are iterative, an initial guess at the
+ solution has to be supplied, and the solution located will
+ usually be the one closest to this initial guess.
+
+ 2.3. References
+
+ [1] Gill P E and Murray W (1976) Algorithms for the Solution of
+ the Nonlinear Least-squares Problem. NAC 71 National
+ Physical Laboratory.
+
+ [2] More J J, Garbow B S and Hillstrom K E (1974) User Guide for
+ Minpack-1. ANL-80-74 Argonne National Laboratory.
+
+ [3] Ortega J M and Rheinboldt W C (1970) Iterative Solution of
+ Nonlinear Equations in Several Variables. Academic Press.
+
+ [4] Rabinowitz P (1970) Numerical Methods for Nonlinear
+ Algebraic Equations. Gordon and Breach.
+
+ 3. Recommendations on Choice and Use of Routines
+
+ 3.1. Zeros of Functions of One Variable
+
+ There is only one routine (C05ADF) for solving a single nonlinear
+ equation. This routine is designed for solving problems where the
+ function f(x) whose zero is to be calculated, can be coded as a
+ user-supplied routine.
+
+ C05ADF may only be used when the user can supply an interval
+ [a,b] containing the zero, that is f(a)*f(b)<0.
+
+ 3.2. Solution of Sets of Nonlinear Equations
+
+ The solution of a set of nonlinear equations
+
+ f (x ,x ,...,x )=0, i=1,2,...,n (2)
+ i 1 2 n
+
+ can be regarded as a special case of the problem of finding a
+ minimum of a sum of squares
+
+ m
+ / 2
+ s(x)= | [f (x ,x ,...,x )] (m>=n). (3)
+ / i 1 2 n
+ i=1
+
+ So the routines in Chapter E04 of the Library are relevant as
+ well as the special nonlinear equations routines.
+
+ There are two routines (C05NBF and C05PBF) for solving a set of
+ nonlinear equations. These routines require the f (and possibly
+ i
+ their derivatives) to be calculated in user-supplied routines.
+ These should be set up carefully so the Library routines can work
+ as efficiently as possible.
+
+ The main decision which has to be made by the user is whether to
+ ddf
+ i
+ supply the derivatives ----. It is advisable to do so if
+ ddx
+ j
+ possible, since the results obtained by algorithms which use
+ derivatives are generally more reliable than those obtained by
+ algorithms which do not use derivatives.
+
+ C05PBF requires the user to provide the derivatives, whilst
+ C05NBF does not. C05NBF and C05PBF are easy-to-use routines. A
+ routine, C05ZAF, is provided for use in conjunction with C05PBF
+ to check the user-provided derivatives for consistency with the
+ functions themselves. The user is strongly advised to make use of
+ this routine whenever C05PBF is used.
+
+ Firstly, the calculation of the functions and their derivatives
+ should be ordered so that cancellation errors are avoided. This
+ is particularly important in a routine that uses these quantities
+ to build up estimates of higher derivatives.
+
+ Secondly, scaling of the variables has a considerable effect on
+ the efficiency of a routine. The problem should be designed so
+ that the elements of x are of similar magnitude. The same comment
+ applies to the functions, all the f should be of comparable
+ i
+ size.
+
+ The accuracy is usually determined by the accuracy parameters of
+ the routines, but the following points may be useful:
+
+ (i) Greater accuracy in the solution may be requested by
+ choosing smaller input values for the accuracy parameters.
+ However, if unreasonable accuracy is demanded, rounding
+ errors may become important and cause a failure.
+
+ (ii) Some idea of the accuracies of the x may be obtained by
+ i
+ monitoring the progress of the routine to see how many
+ figures remain unchanged during the last few iterations.
+
+ (iii) An approximation to the error in the solution x, given by e
+ where e is the solution to the set of linear equations
+
+ J(x)e=-f(x)
+
+ T
+ where f(x)=(f (x),f (x),...,f (x)) (see Chapter F04).
+ 1 2 n
+
+ (iv) If the functions f (x) are changed by small amounts
+ i
+ (epsilon) , for i=1,2,...,n, then the corresponding change
+ i
+ in the solution x is given approximately by (sigma), where
+ (sigma) is the solution of the set of linear equations
+
+ J(x)(sigma)=-(epsilon), (see Chapter F04).
+
+ Thus one can estimate the sensitivity of x to any
+ uncertainties in the specification of f (x), for
+ i
+ i=1,2,...,n.
+
+ 3.3. Index
+
+ Zeros of functions of one variable:
+ Bus and Dekker algorithm C05ADF
+ Zeros of functions of several variables:
+ easy-to-use C05NBF
+ easy-to-use, derivatives required C05PBF
+ Checking Routine:
+ Checks user-supplied Jacobian C05ZAF
+
+
+ C05 -- Roots of One or More Transcendental Equations
+ Contents -- C05
+ Chapter C05
+
+ Roots of One or More Transcendental Equations
+
+ C05ADF Zero of continuous function in given interval, Bus and
+ Dekker algorithm
+
+ C05NBF Solution of system of nonlinear equations using function
+ values only
+
+ C05PBF Solution of system of nonlinear equations using 1st
+ derivatives
+
+ C05ZAF Check user's routine for calculating 1st derivatives
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C05ADF(3NAG) Foundation Library (12/10/92) C05ADF(3NAG)
+
+ C05 -- Roots of One or More Transcendental Equations C05ADF
+ C05ADF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C05ADF locates a zero of a continuous function in a given
+ interval by a combination of the methods of linear interpolation,
+ extrapolation and bisection.
+
+ 2. Specification
+
+ SUBROUTINE C05ADF (A, B, EPS, ETA, F, X, IFAIL)
+ INTEGER IFAIL
+ DOUBLE PRECISION A, B, EPS, ETA, F, X
+ EXTERNAL F
+
+ 3. Description
+
+ The routine attempts to obtain an approximation to a simple zero
+ of the function f(x) given an initial interval [a,b] such that
+ f(a)*f(b)<=0. The zero is found by calls to C05AZF(*) whose
+ specification should be consulted for details of the method used.
+
+ The approximation x to the zero (alpha) is determined so that one
+ or both of the following criteria are satisfied:
+
+ (i) |x-(alpha)| 0.0.
+
+ 4: ETA -- DOUBLE PRECISION Input
+ On entry: a value such that if |f(x)|0.0.
+
+ IFAIL= 2
+ Too much accuracy has been requested in the computation,
+ that is, EPS is too small for the computer being used. The
+ final value of X is an accurate approximation to the zero.
+
+ IFAIL= 3
+ A change in sign of f(x) has been determined as occurring
+ near the point defined by the final value of X. However,
+ there is some evidence that this sign-change corresponds to
+ a pole of f(x).
+
+ IFAIL= 4
+ Indicates that a serious error has occurred in C05AZF(*).
+ Check all routine calls. Seek expert help.
+
+ 7. Accuracy
+
+ This depends on the value of EPS and ETA. If full machine
+ accuracy is required, they may be set very small, resulting in an
+ error exit with IFAIL = 2, although this may involve more
+ iterations than a lesser accuracy. The user is recommended to set
+ ETA = 0.0 and to use EPS to control the accuracy, unless he has
+ considerable knowledge of the size of f(x) for values of x near
+ the zero.
+
+ 8. Further Comments
+
+ The time taken by the routine depends primarily on the time spent
+ evaluating F (see Section 5).
+
+ If it is important to determine an interval of length less than
+ EPS containing the zero, or if the function F is expensive to
+ evaluate and the number of calls to F is to be restricted, then
+ use of C05AZF(*) is recommended. Use of C05AZF(*) is also
+ recommended when the structure of the problem to be solved does
+ not permit a simple function F to be written: the reverse
+ communication facilities of C05AZF(*) are more flexible than the
+ direct communication of F required by C05ADF.
+
+ 9. Example
+
+ -x
+ The example program below calculates the zero of e -x within the
+ interval [0,1] to approximately 5 decimal places.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C05NBF(3NAG) Foundation Library (12/10/92) C05NBF(3NAG)
+
+ C05 -- Roots of One or More Transcendental Equations C05NBF
+ C05NBF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C05NBF is an easy-to-use routine to find a solution of a system
+ of nonlinear equations by a modification of the Powell hybrid
+ method.
+
+ 2. Specification
+
+ SUBROUTINE C05NBF (FCN, N, X, FVEC, XTOL, WA, LWA, IFAIL)
+ INTEGER N, LWA, IFAIL
+ DOUBLE PRECISION X(N), FVEC(N), XTOL, WA(LWA)
+ EXTERNAL FCN
+
+ 3. Description
+
+ The system of equations is defined as:
+
+ f (x ,x ,...,x )=0, for i=1,2,...,n.
+ i 1 2 n
+
+ C05NBF is based upon the MINPACK routine HYBRD1 (More et al [1]).
+ It chooses the correction at each step as a convex combination of
+ the Newton and scaled gradient directions. Under reasonable
+ conditions this guarantees global convergence for starting points
+ far from the solution and a fast rate of convergence. The
+ Jacobian is updated by the rank-1 method of Broyden. At the
+ starting point the Jacobian is approximated by forward
+ differences, but these are not used again until the rank-1 method
+ fails to produce satisfactory progress. For more details see
+ Powell [2].
+
+ 4. References
+
+ [1] More J J, Garbow B S and Hillstrom K E User Guide for
+ MINPACK-1. Technical Report ANL-80-74. Argonne National
+ Laboratory.
+
+ [2] Powell M J D (1970) A Hybrid Method for Nonlinear Algebraic
+ Equations. Numerical Methods for Nonlinear Algebraic
+ Equations. (ed P Rabinowitz) Gordon and Breach.
+
+ 5. Parameters
+
+ 1: FCN -- SUBROUTINE, supplied by the user.
+ External Procedure
+ FCN must return the values of the functions f at a point x.
+ i
+
+ Its specification is:
+
+ SUBROUTINE FCN (N, X, FVEC, IFLAG)
+ INTEGER N, IFLAG
+ DOUBLE PRECISION X(N), FVEC(N)
+
+ 1: N -- INTEGER Input
+ On entry: the number of equations, n.
+
+ 2: X(N) -- DOUBLE PRECISION array Input
+ On entry: the components of the point x at which the
+ functions must be evaluated.
+
+ 3: FVEC(N) -- DOUBLE PRECISION array Output
+ On exit: the function values f (x) (unless IFLAG is
+ i
+ set to a negative value by FCN).
+
+ 4: IFLAG -- INTEGER Input/Output
+ On entry: IFLAG > 0. On exit: in general, IFLAG should
+ not be reset by FCN. If, however, the user wishes to
+ terminate execution (perhaps because some illegal point
+ X has been reached), then IFLAG should be set to a
+ negative integer. This value will be returned through
+ IFAIL.
+ FCN must be declared as EXTERNAL in the (sub)program
+ from which C05NBF is called. Parameters denoted as
+ Input must not be changed by this procedure.
+
+ 2: N -- INTEGER Input
+ On entry: the number of equations, n. Constraint: N > 0.
+
+ 3: X(N) -- DOUBLE PRECISION array Input/Output
+ On entry: an initial guess at the solution vector. On
+ exit: the final estimate of the solution vector.
+
+ 4: FVEC(N) -- DOUBLE PRECISION array Output
+ On exit: the function values at the final point, X.
+
+ 5: XTOL -- DOUBLE PRECISION Input
+ On entry: the accuracy in X to which the solution is
+ required. Suggested value: the square root of the machine
+ precision. Constraint: XTOL >= 0.0.
+
+ 6: WA(LWA) -- DOUBLE PRECISION array Workspace
+
+ 7: LWA -- INTEGER Input
+ On entry: the dimension of the array WA. Constraint:
+ LWA>=N*(3*N+13)/2.
+
+ 8: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ If on entry IFAIL = 0 or -1, explanatory error messages are
+ output on the current error message unit (as defined by X04AAF).
+
+ IFAIL< 0
+ The user has set IFLAG negative in FCN. The value of IFAIL
+ will be the same as the user's setting of IFLAG.
+
+ IFAIL= 1
+ On entry N <= 0,
+
+ or XTOL < 0.0,
+
+ or LWA 0.
+
+ 3: X(N) -- DOUBLE PRECISION array Input/Output
+ On entry: an initial guess at the solution vector. On
+ exit: the final estimate of the solution vector.
+
+ 4: FVEC(N) -- DOUBLE PRECISION array Output
+ On exit: the function values at the final point, X.
+
+ 5: FJAC(LDFJAC,N) -- DOUBLE PRECISION array Output
+ On exit: the orthogonal matrix Q produced by the QR
+ factorization of the final approximate Jacobian.
+
+ 6: LDFJAC -- INTEGER Input
+ On entry:
+ the first dimension of the array FJAC as declared in the
+ (sub)program from which C05PBF is called.
+ Constraint: LDFJAC >= N.
+
+ 7: XTOL -- DOUBLE PRECISION Input
+ On entry: the accuracy in X to which the solution is
+ required. Suggested value: the square root of the machine
+ precision. Constraint: XTOL >= 0.0.
+
+ 8: WA(LWA) -- DOUBLE PRECISION array Workspace
+
+ 9: LWA -- INTEGER Input
+ On entry: the dimension of the array WA. Constraint:
+ LWA>=N*(N+13)/2.
+
+ 10: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ If on entry IFAIL = 0 or -1, explanatory error messages are
+ output on the current error message unit (as defined by X04AAF).
+
+ IFAIL< 0
+ A negative value of IFAIL indicates an exit from C05PBF
+ because the user has set IFLAG negative in FCN. The value of
+ IFAIL will be the same as the user's setting of IFLAG.
+
+ IFAIL= 1
+ On entry N <= 0,
+
+ or LDFJAC < N,
+
+ or XTOL < 0.0,
+
+ or LWA>=
+
+ C06(3NAG) Foundation Library (12/10/92) C06(3NAG)
+
+
+
+ C06 -- Summation of Series Introduction -- C06
+ Chapter C06
+ Summation of Series
+
+ 1. Scope of the Chapter
+
+ This chapter is concerned with calculating the discrete Fourier
+ transform of a sequence of real or complex data values, and
+ applying it to calculate convolutions and correlations.
+
+ 2. Background to the Problems
+
+ 2.1. Discrete Fourier Transforms
+
+ 2.1.1. Complex transforms
+
+ Most of the routines in this chapter calculate the finite
+ discrete Fourier transform (DFT) of a sequence of n complex
+ numbers z , for j=0,1,...,n-1. The transform is defined by:
+ j
+
+ n-1
+ ^ 1 -- ( 2(pi)jk)
+ z = --- > z exp(-i -------) (1)
+ k -- j ( n )
+ \/n j=0
+
+ for k=0,1,...,n-1. Note that equation (1) makes sense for all
+ ^
+ integral k and with this extension z is periodic with period n,
+ k
+ ^ ^ ^ ^
+ i.e. z =z , and in particular z =z .
+ k k+-n -k n-k
+
+ ^ ^
+ If we write z =x +iy and z =a +ib , then the definition of z
+ j j j k k k k
+ may be written in terms of sines and cosines as:
+
+ n-1
+ 1 -- ( ( 2(pi)jk) ( 2(pi)jk))
+ a = --- > (x cos( -------)+y sin( -------))
+ k -- ( j ( n ) j ( n ))
+ \/n j=0
+
+ n-1
+ 1 -- ( ( 2(pi)jk) ( 2(pi)jk))
+ b = --- > (y cos( -------)-x sin( -------)).
+ k -- ( j ( n ) j ( n ))
+ \/n j=0
+
+ The original data values z may conversely be recovered from the
+ j
+ ^
+ transform z by an inverse discrete Fourier transform:
+ k
+
+
+ n-1
+ 1 -- ^ ( 2(pi)jk)
+ z = --- > z exp(+i -------) (2)
+ j -- k ( n )
+ \/n k=0
+
+ for j=0,1,...,n-1. If we take the complex conjugate of (2), we
+
+
+ ^
+ find that the sequence z is the DFT of the sequence z . Hence
+ j k
+ ^
+ the inverse DFT of the sequence z may be obtained by: taking the
+ k
+ ^
+ complex conjugates of the z ; performing a DFT; and taking the
+ k
+ complex conjugates of the result.
+
+ Notes: definitions of the discrete Fourier transform vary.
+ Sometimes (2) is used as the definition of the DFT, and (1) as
+
+
+ the definition of the inverse. Also the scale-factor of 1/\/n may
+ be omitted in the definition of the DFT, and replaced by 1/n in
+ the definition of the inverse.
+
+ 2.1.2. Real transforms
+
+ If the original sequence is purely real valued, i.e. z =x , then
+ j j
+
+ n-1
+ ^ 1 -- ( 2(pi)jk)
+ z =a +ib = --- > x exp(-i -------)
+ k k k -- j ( n )
+ \/n j=0
+
+ ^ ^
+ and z is the complex conjugate of z . Thus the DFT of a real
+ n-k k
+ sequence is a particular type of complex sequence, called a
+ Hermitian sequence, or half-complex or conjugate symmetric with
+ the properties:
+ a =a b =-b b =0 and, if n is even, b =0.
+ n-k k n-k k 0 n/2
+
+ Thus a Hermitian sequence of n complex data values can be
+ represented by only n, rather than 2n, independent real values.
+ This can obviously lead to economies in storage, the following
+ scheme being used in this chapter: the real parts a for
+ k
+ 0<=k<=n/2 are stored in normal order in the first n/2+1 locations
+ of an array X of length n; the corresponding non-zero imaginary
+ parts are stored in reverse order in the remaining locations of
+ X. In other words, if X is declared with bounds (0:n-1) in the
+ ^
+ user's (sub)program, the real and imaginary parts of z are
+ k
+ stored as follows:
+
+ if n=2s if n=2s-1
+
+ X(0) a a
+ 0 0
+
+ X(1) a a
+ 1 1
+
+ X(2) a a
+ 2 2
+
+ . . .
+
+ . . .
+
+ . . .
+
+ X(s-1) a a
+ s-1 s-1
+
+ X(s) a b
+ s s-1
+
+ X(s+1) b b
+ s-1 s-2
+
+ . . .
+
+ . . .
+
+ . . .
+
+ X(n-2) b b
+ 2 2
+
+ X(n-1) b b
+ 1 1
+
+
+ ( n/2-1 )
+ 1 ( -- ( ( 2(pi)jk) ( 2(pi)jk)) )
+ Hence x = ---(a +2 > (a cos( -------)-b sin( -------))+a )
+ j ( 0 -- ( k ( n ) k ( n )) n/2)
+ \/n( k=0 )
+
+ where a = 0 if n is odd.
+ n/2
+
+ 2.1.3. Fourier integral transforms
+
+ The usual application of the discrete Fourier transform is that
+ of obtaining an approximation of the Fourier integral transform
+
+
+ +infty
+ /
+ F(s)= | f(t)exp(-i2(pi)st)dt
+ /
+ -infty
+
+ when f(t) is negligible outside some region (0,c). Dividing the
+ region into n equal intervals we have
+
+ n-1
+ c --
+ F(s)~= - > f exp(-i2(pi)sjc/n)
+ n -- j
+ j=0
+
+ and so
+
+ n-1
+ c --
+ F ~= - > f exp(-i2(pi)jk/n)
+ k n -- j
+ j=0
+
+ for k=0,1,...,n-1, where f =f(jc/n) and F =F(k/c).
+ j k
+
+ Hence the discrete Fourier transform gives an approximation to
+ the Fourier integral transform in the region s=0 to s=n/c.
+
+ If the function f(t) is defined over some more general interval
+ (a,b), then the integral transform can still be approximated by
+ the discrete transform provided a shift is applied to move the
+ point a to the origin.
+
+ 2.1.4. Convolutions and correlations
+
+ One of the most important applications of the discrete Fourier
+ transform is to the computation of the discrete convolution or
+ correlation of two vectors x and y defined (as in Brigham [1])
+ by:
+
+ n-1
+ --
+ convolution: z = > x y
+ k -- j k-j
+ j=0
+
+ n-1
+ --
+ correlation: w = > x y
+ k -- j k+j
+ j=0
+
+ (Here x and y are assumed to be periodic with period n.)
+
+ Under certain circumstances (see Brigham [1]) these can be used
+ as approximations to the convolution or correlation integrals
+ defined by:
+
+ +infty
+ /
+ z(s)= | x(t)y(s-t)dt
+ /
+ -infty
+
+ and
+
+ +infty
+ /
+ w(s)= | x(t)y(s+t)dt, -infty~~10 ).
+
+ Group 2 routines are designed to perform several transforms in a
+ single call, all with the same value of n. They do however
+ require more working storage. Even on scalar processors, they may
+ be somewhat faster than repeated calls to Group 1 routines
+ because of reduced overheads and because they pre-compute and
+ store the required values of trigonometric functions. Group 2
+ routines impose no practical restrictions on the value of n;
+ however the fast Fourier transform algorithm ceases to be 'fast'
+ if applied to values of n which cannot be expressed as a product
+ of small prime factors. All the above routines are particularly
+ efficient if the only prime factors of n are 2, 3 or 5.
+
+ If extensive use is to be made of these routines, users who are
+ concerned about efficiency are advised to conduct their own
+ timing tests.
+
+ To compute inverse discrete Fourier transforms the above routines
+ should be used in conjunction with the utility routines C06GBF,
+ C06GCF and C06GQF which form the complex conjugate of a Hermitian
+ or general sequence of complex data values.
+
+ 3.2. Multi-dimensional Fourier Transforms
+
+ C06FUF computes a 2-dimensional discrete Fourier transform of a
+ 2-dimensional sequence of complex data values. This is defined by
+
+ n -1 n -1
+ 1 2 ( 2(pi)j k ) ( 2(pi)j k )
+ ^ 1 -- -- ( 1 1) ( 2 2)
+ z = ------- > > z exp(-i ---------)exp(-i ---------).
+ -- -- ( n ) ( n )
+ k k /n n j =0 j =0 j j ( 1 ) ( 2 )
+ 1 2 \/ 1 2 1 2 1 2
+
+ 3.3. Convolution and Correlation
+
+ C06EKF computes either the discrete convolution or the discrete
+ correlation of two real vectors.
+
+ 3.4. Index
+
+ Complex conjugate,
+ complex sequence C06GCF
+ Hermitian sequence C06GBF
+ multiple Hermitian sequences C06GQF
+ Complex sequence from Hermitian sequences C06GSF
+ Convolution or Correlation
+ real vectors C06EKF
+ Discrete Fourier Transform
+ two-dimensional
+ complex sequence C06FUF
+ one-dimensional, multiple transforms
+ complex sequence C06FRF
+ Hermitian sequence C06FQF
+ real sequence C06FPF
+ one-dimensional, single transforms
+ complex sequence C06ECF
+ Hermitian sequence C06EBF
+ real sequence C06EAF
+
+
+
+ C06 -- Summation of Series Contents -- C06
+ Chapter C06
+
+ Summation of Series
+
+ C06EAF Single 1-D real discrete Fourier transform, no extra
+ workspace
+
+ C06EBF Single 1-D Hermitian discrete Fourier transform, no extra
+ workspace
+
+ C06ECF Single 1-D complex discrete Fourier transform, no extra
+ workspace
+
+ C06EKF Circular convolution or correlation of two real vectors,
+ no extra workspace
+
+ C06FPF Multiple 1-D real discrete Fourier transforms
+
+ C06FQF Multiple 1-D Hermitian discrete Fourier transforms
+
+ C06FRF Multiple 1-D complex discrete Fourier transforms
+
+ C06FUF 2-D complex discrete Fourier transform
+
+ C06GBF Complex conjugate of Hermitian sequence
+
+ C06GCF Complex conjugate of complex sequence
+
+ C06GQF Complex conjugate of multiple Hermitian sequences
+
+ C06GSF Convert Hermitian sequences to general complex sequences
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06EAF(3NAG) Foundation Library (12/10/92) C06EAF(3NAG)
+
+ C06 -- Summation of Series C06EAF
+ C06EAF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06EAF calculates the discrete Fourier transform of a sequence of
+ n real data values. (No extra workspace required.)
+
+ 2. Specification
+
+ SUBROUTINE C06EAF (X, N, IFAIL)
+ INTEGER N, IFAIL
+ DOUBLE PRECISION X(N)
+
+ 3. Description
+
+ Given a sequence of n real data values x , for j=0,1,...,n-1,
+ j
+ this routine calculates their discrete Fourier transform defined
+ by:
+
+ n-1
+ ^ 1 -- ( 2(pi)jk)
+ z = --- > x *exp(-i -------), k=0,1,...,n-1.
+ k -- j ( n )
+ \/n j=0
+
+ 1
+ (Note the scale factor of --- in this definition.) The
+
+
+ \/n
+ ^
+ transformed values z are complex, but they form a Hermitian
+ k
+ ^ ^
+ sequence (i.e., z is the complex conjugate of z ), so they are
+ n-k k
+ completely determined by n real numbers (see also the Chapter
+ Introduction).
+
+ To compute the inverse discrete Fourier transform defined by:
+
+ n-1
+ ^ 1 -- ( 2(pi)jk)
+ w = --- > x *exp(+i -------),
+ k -- j ( n )
+ \/n j=0
+
+ this routine should be followed by a call of C06GBF to form the
+ ^
+ complex conjugates of the z .
+ k
+
+ The routine uses the fast Fourier transform (FFT) algorithm
+ (Brigham [1]). There are some restrictions on the value of n (see
+ Section 5).
+
+ 4. References
+
+ [1] Brigham E O (1973) The Fast Fourier Transform. Prentice-
+ Hall.
+
+ 5. Parameters
+
+ 1: X(N) -- DOUBLE PRECISION array Input/Output
+ On entry: if X is declared with bounds (0:N-1) in the (sub)
+ program from which C06EAF is called, then X(j) must contain
+ x , for j=0,1,...,n-1. On exit: the discrete Fourier
+ j
+ transform stored in Hermitian form. If the components of the
+ ^
+ transform z are written as a +ib , and if X is declared
+ k k k
+ with bounds (0:N-1) in the (sub)program from which C06EAF is
+ called, then for 0<=k<=n/2, a is contained in X(k), and for
+ k
+ 1<=k<=(n-1)/2, b is contained in X(n-k). (See also Section
+ k
+ 2.1.2 of the Chapter Introduction, and the Example Program.)
+
+ 2: N -- INTEGER Input
+ On entry: the number of data values, n. The largest prime
+ factor of N must not exceed 19, and the total number of
+ prime factors of N, counting repetitions, must not exceed
+ 20. Constraint: N > 1.
+
+ 3: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ IFAIL= 1
+ At least one of the prime factors of N is greater than 19.
+
+ IFAIL= 2
+ N has more than 20 prime factors.
+
+ IFAIL= 3
+ N <= 1.
+
+ 7. Accuracy
+
+ Some indication of accuracy can be obtained by performing a
+ subsequent inverse transform and comparing the results with the
+ original sequence (in exact arithmetic they would be identical).
+
+ 8. Further Comments
+
+ The time taken by the routine is approximately proportional to
+ n*logn, but also depends on the factorization of n. The routine
+ is somewhat faster than average if the only prime factors of n
+ are 2, 3 or 5; and fastest of all if n is a power of 2.
+
+ On the other hand, the routine is particularly slow if n has
+ several unpaired prime factors, i.e., if the 'square-free' part
+ of n has several factors. For such values of n, routine C06FAF(*)
+ (which requires an additional n elements of workspace) is
+ considerably faster.
+
+ 9. Example
+
+ This program reads in a sequence of real data values, and prints
+ their discrete Fourier transform (as computed by C06EAF), after
+ expanding it from Hermitian form into a full complex sequence.
+
+ It then performs an inverse transform using C06GBF and C06EBF,
+ and prints the sequence so obtained alongside the original data
+ values.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06EBF(3NAG) Foundation Library (12/10/92) C06EBF(3NAG)
+
+ C06 -- Summation of Series C06EBF
+ C06EBF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06EBF calculates the discrete Fourier transform of a Hermitian
+ sequence of n complex data values. (No extra workspace required.)
+
+ 2. Specification
+
+ SUBROUTINE C06EBF (X, N, IFAIL)
+ INTEGER N, IFAIL
+ DOUBLE PRECISION X(N)
+
+ 3. Description
+
+ Given a Hermitian sequence of n complex data values z (i.e., a
+ j
+ sequence such that z is real and z is the complex conjugate
+ 0 n-j
+ of z , for j=1,2,...,n-1) this routine calculates their discrete
+ j
+ Fourier transform defined by:
+
+ n-1
+ ^ 1 -- ( 2(pi)jk)
+ x = --- > z *exp(-i -------), k=0,1,...,n-1.
+ k -- j ( n )
+ \/n j=0
+
+ 1
+ (Note the scale factor of --- in this definition.) The
+
+
+ \/n
+ ^
+ transformed values x are purely real (see also the the Chapter
+ k
+ Introduction).
+
+ To compute the inverse discrete Fourier transform defined by:
+
+ n-1
+ ^ 1 -- ( 2(pi)jk)
+ y = --- > z *exp(+i -------),
+ k -- j ( n )
+ \/n j=0
+
+ this routine should be preceded by a call of C06GBF to form the
+ complex conjugates of the z .
+ j
+
+ The routine uses the fast Fourier transform (FFT) algorithm
+ (Brigham [1]). There are some restrictions on the value of n (see
+ Section 5).
+
+ 4. References
+
+ [1] Brigham E O (1973) The Fast Fourier Transform. Prentice-
+ Hall.
+
+ 5. Parameters
+
+ 1: X(N) -- DOUBLE PRECISION array Input/Output
+ On entry: the sequence to be transformed stored in
+ Hermitian form. If the data values z are written as x +iy ,
+ j j j
+ and if X is declared with bounds (0:N-1) in the subroutine
+ from which C06EBF is called, then for 0<=j<=n/2, x is
+ j
+ contained in X(j), and for 1<=j<=(n-1)/2, y is contained in
+ j
+ X(n-j). (See also Section 2.1.2 of the Chapter Introduction
+ and the Example Program.) On exit: the components of the
+ ^
+ discrete Fourier transform x . If X is declared with bounds
+ k
+ (0:N-1) in the (sub)program from which C06EBF is called,
+ ^
+ then x is stored in X(k), for k=0,1,...,n-1.
+ k
+
+ 2: N -- INTEGER Input
+ On entry: the number of data values, n. The largest prime
+ factor of N must not exceed 19, and the total number of
+ prime factors of N, counting repetitions, must not exceed
+ 20. Constraint: N > 1.
+
+ 3: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ IFAIL= 1
+ At least one of the prime factors of N is greater than 19.
+
+ IFAIL= 2
+ N has more than 20 prime factors.
+
+ IFAIL= 3
+ N <= 1.
+
+ 7. Accuracy
+
+ Some indication of accuracy can be obtained by performing a
+ subsequent inverse transform and comparing the results with the
+ original sequence (in exact arithmetic they would be identical).
+
+ 8. Further Comments
+
+ The time taken by the routine is approximately proportional to
+ n*logn, but also depends on the factorization of n. The routine
+ is somewhat faster than average if the only prime factors of n
+ are 2, 3 or 5; and fastest of all if n is a power of 2.
+
+ On the other hand, the routine is particularly slow if n has
+ several unpaired prime factors, i.e., if the 'square-free' part
+ of n has several factors. For such values of n, routine C06FBF(*)
+ (which requires an additional n elements of workspace) is
+ considerably faster.
+
+ 9. Example
+
+ This program reads in a sequence of real data values which is
+ assumed to be a Hermitian sequence of complex data values stored
+ in Hermitian form. The input sequence is expanded into a full
+ complex sequence and printed alongside the original sequence. The
+ discrete Fourier transform (as computed by C06EBF) is printed
+ out.
+
+ The program then performs an inverse transform using C06EAF and
+ C06GBF, and prints the sequence so obtained alongside the
+ original data values.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06ECF(3NAG) Foundation Library (12/10/92) C06ECF(3NAG)
+
+ C06 -- Summation of Series C06ECF
+ C06ECF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06ECF calculates the discrete Fourier transform of a sequence of
+ n complex data values. (No extra workspace required.)
+
+ 2. Specification
+
+ SUBROUTINE C06ECF (X, Y, N, IFAIL)
+ INTEGER N, IFAIL
+ DOUBLE PRECISION X(N), Y(N)
+
+ 3. Description
+
+ Given a sequence of n complex data values z , for j=0,1,...,n-1,
+ j
+ this routine calculates their discrete Fourier transform defined
+ by:
+
+ n-1
+ ^ 1 -- ( 2(pi)jk)
+ z = --- > z *exp(-i -------), k=0,1,...,n-1.
+ k -- j ( n )
+ \/n j=0
+
+ 1
+ (Note the scale factor of --- in this definition.)
+
+
+ \/n
+
+ To compute the inverse discrete Fourier transform defined by:
+
+ n-1
+ ^ 1 -- ( 2(pi)jk)
+ w = --- > z *exp(+i -------),
+ k -- j ( n )
+ \/n j=0
+
+ this routine should be preceded and followed by calls of C06GCF
+ ^
+ to form the complex conjugates of the z and the z .
+ j k
+
+ The routine uses the fast Fourier transform (FFT) algorithm
+ (Brigham [1]). There are some restrictions on the value of n (see
+ Section 5).
+
+ 4. References
+
+ [1] Brigham E O (1973) The Fast Fourier Transform. Prentice-
+ Hall.
+
+ 5. Parameters
+
+ 1: X(N) -- DOUBLE PRECISION array Input/Output
+ On entry: if X is declared with bounds (0:N-1) in the (sub)
+ program from which C06ECF is called, then X(j) must contain
+ x , the real part of z , for j=0,1,...,n-1. On exit: the
+ j j
+ real parts a of the components of the discrete Fourier
+ k
+ transform. If X is declared with bounds (0:N-1) in the (sub)
+ program from which C06ECF is called, then a is contained in
+ k
+ X(k), for k=0,1,...,n-1.
+
+ 2: Y(N) -- DOUBLE PRECISION array Input/Output
+ On entry: if Y is declared with bounds (0:N-1) in the (sub)
+ program from which C06ECF is called, then Y(j) must contain
+ y , the imaginary part of z , for j=0,1,...,n-1. On exit:
+ j j
+ the imaginary parts b of the components of the discrete
+ k
+ Fourier transform. If Y is declared with bounds (0:N-1) in
+ the (sub)program from which C06ECF is called, then b is
+ k
+ contained in Y(k), for k=0,1,...,n-1.
+
+ 3: N -- INTEGER Input
+ On entry: the number of data values, n. The largest prime
+ factor of N must not exceed 19, and the total number of
+ prime factors of N, counting repetitions, must not exceed
+ 20. Constraint: N > 1.
+
+ 4: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ IFAIL= 1
+ At least one of the prime factors of N is greater than 19.
+
+ IFAIL= 2
+ N has more than 20 prime factors.
+
+ IFAIL= 3
+ N <= 1.
+
+ 7. Accuracy
+
+ Some indication of accuracy can be obtained by performing a
+ subsequent inverse transform and comparing the results with the
+ original sequence (in exact arithmetic they would be identical).
+
+ 8. Further Comments
+
+ The time taken by the routine is approximately proportional to
+ n*logn, but also depends on the factorization of n. The routine
+ is somewhat faster than average if the only prime factors of n
+ are 2, 3 or 5; and fastest of all if n is a power of 2.
+
+ On the other hand, the routine is particularly slow if n has
+ several unpaired prime factors, i.e., if the 'square-free' part
+ of n has several factors. For such values of n, routine C06FCF(*)
+ (which requires an additional n real elements of workspace) is
+ considerably faster.
+
+ 9. Example
+
+ This program reads in a sequence of complex data values and
+ prints their discrete Fourier transform.
+
+ It then performs an inverse transform using C06GCF and C06ECF,
+ and prints the sequence so obtained alongside the original data
+ values.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06EKF(3NAG) Foundation Library (12/10/92) C06EKF(3NAG)
+
+ C06 -- Summation of Series C06EKF
+ C06EKF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06EKF calculates the circular convolution or correlation of two
+ real vectors of period n. No extra workspace is required.
+
+ 2. Specification
+
+ SUBROUTINE C06EKF (JOB, X, Y, N, IFAIL)
+ INTEGER JOB, N, IFAIL
+ DOUBLE PRECISION X(N), Y(N)
+
+ 3. Description
+
+ This routine computes:
+
+ if JOB =1, the discrete convolution of x and y, defined by:
+
+ n-1 n-1
+ -- --
+ z = > x y = > x y ;
+ k -- j k-j -- k-j j
+ j=0 j=0
+
+ if JOB =2, the discrete correlation of x and y defined by:
+
+ n-1
+ --
+ w = > x y .
+ k -- j k+j
+ j=0
+
+ Here x and y are real vectors, assumed to be periodic, with
+ period n, i.e., x =x =x =...; z and w are then also
+ j j+-n j+-2n
+ periodic with period n.
+
+ Note: this usage of the terms 'convolution' and 'correlation' is
+ taken from Brigham [1]. The term 'convolution' is sometimes used
+ to denote both these computations.
+
+ ^ ^ ^ ^
+ If x, y, z and w are the discrete Fourier transforms of these
+ sequences,
+
+ n-1
+ ^ 1 -- ( 2(pi)jk)
+ i.e., x = --- > x *exp(-i -------), etc,
+ k -- j ( n )
+ \/n j=0
+
+ ^ ^ ^
+ then z =\/n.x y
+ k k k
+
+
+
+ ^ ^ ^
+ and w =\/n.x y
+ k k k
+
+ (the bar denoting complex conjugate).
+
+ This routine calls the same auxiliary routines as C06EAF and
+ C06EBF to compute discrete Fourier transforms, and there are some
+ restrictions on the value of n.
+
+ 4. References
+
+ [1] Brigham E O (1973) The Fast Fourier Transform. Prentice-
+ Hall.
+
+ 5. Parameters
+
+ 1: JOB -- INTEGER Input
+ On entry: the computation to be performed:
+ n-1
+ --
+ if JOB = 1, z = > x y (convolution);
+ k -- j k-j
+ j=0
+
+ n-1
+ --
+ if JOB = 2, w = > x y (correlation).
+ k -- j k+j
+ j=0
+ Constraint: JOB = 1 or 2.
+
+ 2: X(N) -- DOUBLE PRECISION array Input/Output
+ On entry: the elements of one period of the vector x. If X
+ is declared with bounds (0:N-1) in the (sub)program from
+ which C06EKF is called, then X(j) must contain x , for
+ j
+ j=0,1,...,n-1. On exit: the corresponding elements of the
+ discrete convolution or correlation.
+
+ 3: Y(N) -- DOUBLE PRECISION array Input/Output
+ On entry: the elements of one period of the vector y. If Y
+ is declared with bounds (0:N-1) in the (sub)program from
+ which C06EKF is called, then Y(j) must contain y , for
+ j
+ j=0,1,...,n-1. On exit: the discrete Fourier transform of
+ the convolution or correlation returned in the array X; the
+ transform is stored in Hermitian form, exactly as described
+ in the document C06EAF.
+
+ 4: N -- INTEGER Input
+ On entry: the number of values, n, in one period of the
+ vectors X and Y. The largest prime factor of N must not
+ exceed 19, and the total number of prime factors of N,
+ counting repetitions, must not exceed 20. Constraint: N > 1.
+
+ 5: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ IFAIL= 1
+ At least one of the prime factors of N is greater than 19.
+
+ IFAIL= 2
+ N has more than 20 prime factors.
+
+ IFAIL= 3
+ N <= 1.
+
+ IFAIL= 4
+ JOB /= 1 or 2.
+
+ 7. Accuracy
+
+ The results should be accurate to within a small multiple of the
+ machine precision.
+
+ 8. Further Comments
+
+ The time taken by the routine is approximately proportional to
+ n*logn, but also depends on the factorization of n. The routine
+ is faster than average if the only prime factors are 2, 3 or 5;
+ and fastest of all if n is a power of 2.
+
+ The routine is particularly slow if n has several unpaired prime
+ factors, i.e., if the 'square free' part of n has several
+ factors. For such values of n, routine C06FKF(*) is considerably
+ faster (but requires an additional workspace of n elements).
+
+ 9. Example
+
+ This program reads in the elements of one period of two real
+ vectors x and y and prints their discrete convolution and
+ correlation (as computed by C06EKF). In realistic computations
+ the number of data values would be much larger.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06FPF(3NAG) Foundation Library (12/10/92) C06FPF(3NAG)
+
+ C06 -- Summation of Series C06FPF
+ C06FPF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06FPF computes the discrete Fourier transforms of m sequences,
+ each containing n real data values. This routine is designed to
+ be particularly efficient on vector processors.
+
+ 2. Specification
+
+ SUBROUTINE C06FPF (M, N, X, INIT, TRIG, WORK, IFAIL)
+ INTEGER M, N, IFAIL
+ DOUBLE PRECISION X(M*N), TRIG(2*N), WORK(M*N)
+ CHARACTER*1 INIT
+
+ 3. Description
+
+ p
+ Given m sequences of n real data values x , for j=0,1,...,n-1;
+ j
+ p=1,2,...,m, this routine simultaneously calculates the Fourier
+ transforms of all the sequences defined by:
+
+ n-1
+ ^p 1 -- p ( 2(pi)jk)
+ z = --- > x *exp(-i -------), k=0,1,...,n-1; p=1,2,...,m.
+ k -- j ( n )
+ \/n j=0
+
+ 1
+ (Note the scale factor --- in this definition.)
+
+
+ \/n
+
+ ^p
+ The transformed values z are complex, but for each value of p
+ k
+ ^p ^p
+ the z form a Hermitian sequence (i.e.,z is the complex
+ k n-k
+ ^p
+ conjugate of z ), so they are completely determined by mn real
+ k
+ numbers (see also the Chapter Introduction).
+
+ The discrete Fourier transform is sometimes defined using a
+ positive sign in the exponential term:
+
+ n-1
+ ^p 1 -- p ( 2(pi)jk)
+ z = --- > x *exp(+i -------).
+ k -- j ( n )
+ \/n j=0
+
+ To compute this form, this routine should be followed by a call
+ ^p
+ to C06GQF to form the complex conjugates of the z .
+ k
+
+ The routine uses a variant of the fast Fourier transform (FFT)
+ algorithm (Brigham [1]) known as the Stockham self-sorting
+ algorithm, which is described in Temperton [2]. Special coding is
+ provided for the factors 2, 3, 4, 5 and 6. This routine is
+ designed to be particularly efficient on vector processors, and
+ it becomes especially fast as M, the number of transforms to be
+ computed in parallel, increases.
+
+ 4. References
+
+ [1] Brigham E O (1973) The Fast Fourier Transform. Prentice-
+ Hall.
+
+ [2] Temperton C (1983) Fast Mixed-Radix Real Fourier Transforms.
+ J. Comput. Phys. 52 340--350.
+
+ 5. Parameters
+
+ 1: M -- INTEGER Input
+ On entry: the number of sequences to be transformed, m.
+ Constraint: M >= 1.
+
+ 2: N -- INTEGER Input
+ On entry: the number of real values in each sequence, n.
+ Constraint: N >= 1.
+
+ 3: X(M,N) -- DOUBLE PRECISION array Input/Output
+ On entry: the data must be stored in X as if in a two-
+ dimensional array of dimension (1:M,0:N-1); each of the m
+ sequences is stored in a row of the array. In other words,
+ if the data values of the pth sequence to be transformed are
+ p
+ denoted by x , for j=0,1,...,n-1, then the mn elements of
+ j
+ the array X must contain the values
+ 1 2 m 1 2 m 1 2 m
+ x ,x ,...,x , x ,x ,..., x ,...,x ,x ,...,x .
+ 0 0 0 1 1 1 n-1 n-1 n-1
+ On exit: the m discrete Fourier transforms stored as if in
+ a two-dimensional array of dimension (1:M,0:N-1). Each of
+ the m transforms is stored in a row of the array in
+ Hermitian form, overwriting the corresponding original
+ sequence. If the n components of the discrete Fourier
+ ^p p p p
+ transform z are written as a +ib , then for 0<=k<=n/2, a
+ k k k k
+ p
+ is contained in X(p,k), and for 1<=k<=(n-1)/2, b is
+ k
+ contained in X(p,n-k). (See also Section 2.1.2 of the
+ Chapter Introduction.)
+
+ 4: INIT -- CHARACTER*1 Input
+ On entry: if the trigonometric coefficients required to
+ compute the transforms are to be calculated by the routine
+ and stored in the array TRIG, then INIT must be set equal to
+ 'I' (Initial call).
+
+ If INIT contains 'S' (Subsequent call), then the routine
+ assumes that trigonometric coefficients for the specified
+ value of n are supplied in the array TRIG, having been
+ calculated in a previous call to one of C06FPF, C06FQF or
+ C06FRF.
+
+ If INIT contains 'R' (Restart then the routine assumes that
+ trigonometric coefficients for the particular value of n are
+ supplied in the array TRIG, but does not check that C06FPF,
+ C06FQF or C06FRF have previously been called. This option
+ allows the TRIG array to be stored in an external file, read
+ in and re-used without the need for a call with INIT equal
+ to 'I'. The routine carries out a simple test to check that
+ the current value of n is consistent with the array TRIG.
+ Constraint: INIT = 'I', 'S' or 'R'.
+
+ 5: TRIG(2*N) -- DOUBLE PRECISION array Input/Output
+ On entry: if INIT = 'S' or 'R', TRIG must contain the
+ required coefficients calculated in a previous call of the
+ routine. Otherwise TRIG need not be set. On exit: TRIG
+ contains the required coefficients (computed by the routine
+ if INIT = 'I').
+
+ 6: WORK(M*N) -- DOUBLE PRECISION array Workspace
+
+ 7: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ If on entry IFAIL = 0 or -1, explanatory error messages are
+ output on the current error message unit (as defined by X04AAF).
+
+ IFAIL= 1
+ On entry M < 1.
+
+ IFAIL= 2
+ N < 1.
+
+ IFAIL= 3
+ INIT is not one of 'I', 'S' or 'R'.
+
+ IFAIL= 4
+ INIT = 'S', but none of C06FPF, C06FQF or C06FRF has
+ previously been called.
+
+ IFAIL= 5
+ INIT = 'S' or 'R', but the array TRIG and the current value
+ of N are inconsistent.
+
+ 7. Accuracy
+
+ Some indication of accuracy can be obtained by performing a
+ subsequent inverse transform and comparing the results with the
+ original sequence (in exact arithmetic they would be identical).
+
+ 8. Further Comments
+
+ The time taken by the routine is approximately proportional to
+ nm*logn, but also depends on the factors of n. The routine is
+ fastest if the only prime factors of n are 2, 3 and 5, and is
+ particularly slow if n is a large prime, or has large prime
+ factors.
+
+ 9. Example
+
+ This program reads in sequences of real data values and prints
+ their discrete Fourier transforms (as computed by C06FPF). The
+ Fourier transforms are expanded into full complex form using
+ C06GSF and printed. Inverse transforms are then calculated by
+ calling C06GQF followed by C06FQF showing that the original
+ sequences are restored.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06FQF(3NAG) Foundation Library (12/10/92) C06FQF(3NAG)
+
+ C06 -- Summation of Series C06FQF
+ C06FQF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06FQF computes the discrete Fourier transforms of m Hermitian
+ sequences, each containing n complex data values. This routine is
+ designed to be particularly efficient on vector processors.
+
+ 2. Specification
+
+ SUBROUTINE C06FQF (M, N, X, INIT, TRIG, WORK, IFAIL)
+ INTEGER M, N, IFAIL
+ DOUBLE PRECISION X(M*N), TRIG(2*N), WORK(M*N)
+ CHARACTER*1 INIT
+
+ 3. Description
+
+ p
+ Given m Hermitian sequences of n complex data values z , for
+ j
+ j=0,1,...,n-1; p=1,2,...,m, this routine simultaneously
+ calculates the Fourier transforms of all the sequences defined
+ by:
+
+ n-1
+ ^p 1 -- p ( 2(pi)jk)
+ x = --- > z *exp(-i -------), k=0,1,...,n-1; p=1,2,...,m.
+ k -- j ( n )
+ \/n j=0
+
+ 1
+ (Note the scale factor --- in this definition.)
+
+
+ \/n
+
+ The transformed values are purely real (see also the Chapter
+ Introduction).
+
+ The discrete Fourier transform is sometimes defined using a
+ positive sign in the exponential term
+
+ n-1
+ ^p 1 -- p ( 2(pi)jk)
+ x = --- > z *exp(+i -------).
+ k -- j ( n )
+ \/n j=0
+
+ To compute this form, this routine should be preceded by a call
+ ^p
+ to C06GQF to form the complex conjugates of the z .
+ j
+
+ The routine uses a variant of the fast Fourier transform (FFT)
+ algorithm (Brigham [1]) known as the Stockham self-sorting
+ algorithm, which is described in Temperton [2]. Special code is
+ included for the factors 2, 3, 4, 5 and 6. This routine is
+ designed to be particularly efficient on vector processors, and
+ it becomes especially fast as m, the number of transforms to be
+ computed in parallel, increases.
+
+ 4. References
+
+ [1] Brigham E O (1973) The Fast Fourier Transform. Prentice-
+ Hall.
+
+ [2] Temperton C (1983) Fast Mixed-Radix Real Fourier Transforms.
+ J. Comput. Phys. 52 340--350.
+
+ 5. Parameters
+
+ 1: M -- INTEGER Input
+ On entry: the number of sequences to be transformed, m.
+ Constraint: M >= 1.
+
+ 2: N -- INTEGER Input
+ On entry: the number of data values in each sequence, n.
+ Constraint: N >= 1.
+
+ 3: X(M,N) -- DOUBLE PRECISION array Input/Output
+ On entry: the data must be stored in X as if in a two-
+ dimensional array of dimension (1:M,0:N-1); each of the m
+ sequences is stored in a row of the array in Hermitian form.
+ p p p
+ If the n data values z are written as x +iy , then for
+ j j j
+ p
+ 0<=j<=n/2, x is contained in X(p,j), and for 1<=j<=(n-1)/2,
+ j
+ p
+ y is contained in X(p,n-j). (See also Section 2.1.2 of the
+ j
+ Chapter Introduction.) On exit: the components of the m
+ discrete Fourier transforms, stored as if in a two-
+ dimensional array of dimension (1:M,0:N-1). Each of the m
+ transforms is stored as a row of the array, overwriting the
+ corresponding original sequence. If the n components of the
+ ^p
+ discrete Fourier transform are denoted by x , for
+ k
+ k=0,1,...,n-1, then the mn elements of the array X contain
+ the values
+ ^1 ^2 ^m ^1 ^2 ^m ^1 ^2 ^m
+ x ,x ,...,x , x ,x ,..., x ,...,x ,x ,...,x .
+ 0 0 0 1 1 1 n-1 n-1 n-1
+
+ 4: INIT -- CHARACTER*1 Input
+ On entry: if the trigonometric coefficients required to
+ compute the transforms are to be calculated by the routine
+ and stored in the array TRIG, then INIT must be set equal to
+ 'I' (Initial call).
+
+ If INIT contains 'S' (Subsequent call), then the routine
+ assumes that trigonometric coefficients for the specified
+ value of n are supplied in the array TRIG, having been
+ calculated in a previous call to one of C06FPF, C06FQF or
+ C06FRF.
+
+ If INIT contains 'R' (Restart), then the routine assumes
+ that trigonometric coefficients for the particular value of
+ N are supplied in the array TRIG, but does not check that
+ C06FPF, C06FQF or C06FRF have previously been called. This
+ option allows the TRIG array to be stored in an external
+ file, read in and re-used without the need for a call with
+ INIT equal to 'I'. The routine carries out a simple test to
+ check that the current value of n is compatible with the
+ array TRIG. Constraint: INIT = 'I', 'S' or 'R'.
+
+ 5: TRIG(2*N) -- DOUBLE PRECISION array Input/Output
+ On entry: if INIT = 'S' or 'R', TRIG must contain the
+ required coefficients calculated in a previous call of the
+ routine. Otherwise TRIG need not be set. On exit: TRIG
+ contains the required coefficients (computed by the routine
+ if INIT = 'I').
+
+ 6: WORK(M*N) -- DOUBLE PRECISION array Workspace
+
+ 7: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ If on entry IFAIL = 0 or -1, explanatory error messages are
+ output on the current error message unit (as defined by X04AAF).
+
+ IFAIL= 1
+ On entry M < 1.
+
+ IFAIL= 2
+ On entry N < 1.
+
+ IFAIL= 3
+ On entry INIT is not one of 'I', 'S' or 'R'.
+
+ IFAIL= 4
+ On entry INIT = 'S', but none of C06FPF, C06FQF and C06FRF
+ has previously been called.
+
+ IFAIL= 5
+ On entry INIT = 'S' or 'R', but the array TRIG and the
+ current value of n are inconsistent.
+
+ 7. Accuracy
+
+ Some indication of accuracy can be obtained by performing a
+ subsequent inverse transform and comparing the results with the
+ original sequence (in exact arithmetic they would be identical).
+
+ 8. Further Comments
+
+ The time taken by the routine is approximately proportional to
+ nm*logn, but also depends on the factors of n. The routine is
+ fastest if the only prime factors of n are 2, 3 and 5, and is
+ particularly slow if n is a large prime, or has large prime
+ factors.
+
+ 9. Example
+
+ This program reads in sequences of real data values which are
+ assumed to be Hermitian sequences of complex data stored in
+ Hermitian form. The sequences are expanded into full complex form
+ using C06GSF and printed. The discrete Fourier transforms are
+ then computed (using C06FQF) and printed out. Inverse transforms
+ are then calculated by calling C06FPF followed by C06GQF showing
+ that the original sequences are restored.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06FRF(3NAG) Foundation Library (12/10/92) C06FRF(3NAG)
+
+ C06 -- Summation of Series C06FRF
+ C06FRF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06FRF computes the discrete Fourier transforms of m sequences,
+ each containing n complex data values. This routine is designed
+ to be particularly efficient on vector processors.
+
+ 2. Specification
+
+ SUBROUTINE C06FRF (M, N, X, Y, INIT, TRIG, WORK, IFAIL)
+ INTEGER M, N, IFAIL
+ DOUBLE PRECISION X(M*N), Y(M*N), TRIG(2*N), WORK(2*M*N)
+ CHARACTER*1 INIT
+
+ 3. Description
+
+ p
+ Given m sequences of n complex data values z , for j=0,1,...,n-1;
+ j
+ p=1,2,...,m, this routine simultaneously calculates the Fourier
+ transforms of all the sequences defined by:
+
+ n-1
+ ^p 1 -- p ( 2(pi)jk)
+ z = --- > z *exp(-i -------), k=0,1,...,n-1; p=1,2,...,m.
+ k -- j ( n )
+ \/n j=0
+
+ 1
+ (Note the scale factor --- in this definition.)
+
+
+ \/n
+
+ The discrete Fourier transform is sometimes defined using a
+ positive sign in the exponential term
+
+ n-1
+ ^p 1 -- p ( 2(pi)jk)
+ z = --- > z *exp(+i -------).
+ k -- j ( n )
+ \/n j=0
+
+ To compute this form, this routine should be preceded and
+ followed by a call of C06GCF to form the complex conjugates of
+ p ^p
+ the z and the z .
+ j k
+
+ The routine uses a variant of the fast Fourier transform (FFT)
+ algorithm (Brigham [1]) known as the Stockham self-sorting
+ algorithm, which is described in Temperton [2]. Special code is
+ provided for the factors 2, 3, 4, 5 and 6. This routine is
+ designed to be particularly efficient on vector processors, and
+ it becomes especially fast as m, the number of transforms to be
+ computed in parallel, increases.
+
+ 4. References
+
+ [1] Brigham E O (1973) The Fast Fourier Transform. Prentice-
+ Hall.
+
+ [2] Temperton C (1983) Self-sorting Mixed-radix Fast Fourier
+ Transforms. J. Comput. Phys. 52 1--23.
+
+ 5. Parameters
+
+ 1: M -- INTEGER Input
+ On entry: the number of sequences to be transformed, m.
+ Constraint: M >= 1.
+
+ 2: N -- INTEGER Input
+ On entry: the number of complex values in each sequence, n.
+ Constraint: N >= 1.
+
+ 3: X(M,N) -- DOUBLE PRECISION array Input/Output
+
+ 4: Y(M,N) -- DOUBLE PRECISION array Input/Output
+ On entry: the real and imaginary parts of the complex data
+ must be stored in X and Y respectively as if in a two-
+ dimensional array of dimension (1:M,0:N-1); each of the m
+ sequences is stored in a row of each array. In other words,
+ if the real parts of the pth sequence to be transformed are
+ p
+ denoted by x , for j=0,1,...,n-1, then the mn elements of
+ j
+ the array X must contain the values
+ 1 2 m 1 2 m 1 2 m
+ x ,x ,...,x , x ,x ,...,x ,..., x ,x ,...,x .
+ 0 0 0 1 1 1 n-1 n-1 n-1
+ On exit: X and Y are overwritten by the real and imaginary
+ parts of the complex transforms.
+
+ 5: INIT -- CHARACTER*1 Input
+ On entry: if the trigonometric coefficients required to
+ compute the transforms are to be calculated by the routine
+ and stored in the array TRIG, then INIT must be set equal to
+ 'I' (Initial call).
+
+ If INIT contains 'S' (Subsequent call), then the routine
+ assumes that trigonometric coefficients for the specified
+ value of n are supplied in the array TRIG, having been
+ calculated in a previous call to one of C06FPF, C06FQF or
+ C06FRF.
+
+ If INIT contains 'R' (Restart) then the routine assumes that
+ trigonometric coefficients for the particular value of n are
+ supplied in the array TRIG, but does not check that C06FPF,
+ C06FQF or C06FRF have previously been called. This option
+ allows the TRIG array to be stored in an external file, read
+ in and re-used without the need for a call with INIT equal
+ to 'I'. The routine carries out a simple test to check that
+ the current value of n is compatible with the array TRIG.
+ Constraint: INIT = 'I', 'S' or 'R'.
+
+ 6: TRIG(2*N) -- DOUBLE PRECISION array Input/Output
+ On entry: if INIT = 'S' or 'R', TRIG must contain the
+ required coefficients calculated in a previous call of the
+ routine. Otherwise TRIG need not be set. On exit: TRIG
+ contains the required coefficients (computed by the routine
+ if INIT = 'I').
+
+ 7: WORK(2*M*N) -- DOUBLE PRECISION array Workspace
+
+ 8: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ If on entry IFAIL = 0 or -1, explanatory error messages are
+ output on the current error message unit (as defined by X04AAF).
+
+ IFAIL= 1
+ On entry M < 1.
+
+ IFAIL= 2
+ On entry N < 1.
+
+ IFAIL= 3
+ On entry INIT is not one of 'I', 'S' or 'R'.
+
+ IFAIL= 4
+ On entry INIT = 'S', but none of C06FPF, C06FQF and C06FRF
+ has previously been called.
+
+ IFAIL= 5
+ On entry INIT = 'S' or 'R', but the array TRIG and the
+ current value of n are inconsistent.
+
+ 7. Accuracy
+
+ Some indication of accuracy can be obtained by performing a
+ subsequent inverse transform and comparing the results with the
+ original sequence (in exact arithmetic they would be identical).
+
+ 8. Further Comments
+
+ The time taken by the routine is approximately proportional to
+ nm*logn, but also depends on the factors of n. The routine is
+ fastest if the only prime factors of n are 2, 3 and 5, and is
+ particularly slow if n is a large prime, or has large prime
+ factors.
+
+ 9. Example
+
+ This program reads in sequences of complex data values and prints
+ their discrete Fourier transforms (as computed by C06FRF).
+ Inverse transforms are then calculated using C06GCF and C06FRF
+ and printed out, showing that the original sequences are
+ restored.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06FUF(3NAG) Foundation Library (12/10/92) C06FUF(3NAG)
+
+ C06 -- Summation of Series C06FUF
+ C06FUF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06FUF computes the two-dimensional discrete Fourier transform of
+ a bivariate sequence of complex data values. This routine is
+ designed to be particularly efficient on vector processors.
+
+ 2. Specification
+
+ SUBROUTINE C06FUF (M, N, X, Y, INIT, TRIGM, TRIGN, WORK,
+ 1 IFAIL)
+ INTEGER M, N, IFAIL
+ DOUBLE PRECISION X(M*N), Y(M*N), TRIGM(2*M), TRIGN(2*N),
+ 1 WORK(2*M*N)
+ CHARACTER*1 INIT
+
+ 3. Description
+
+ This routine computes the two-dimensional discrete Fourier
+ transform of a bivariate sequence of complex data values z ,
+ j j
+ 1 2
+ where j =0,1,...,m-1, j =0,1,...,n-1.
+ 1 2
+
+ The discrete Fourier transform is here defined by:
+
+ m-1 n-1 ( ( j k j k ))
+ ^ 1 -- -- ( ( 1 1 2 2))
+ z = ---- > > z *exp(-2(pi)i( ----+ ----)),
+ -- -- j j ( ( m n ))
+ k k \/mn j =0 j =0 1 2
+ 1 2 1 2
+
+ where k =0,1,...,m-1, k =0,1,...,n-1.
+ 1 2
+
+ 1
+ (Note the scale factor of ---- in this definition.)
+
+
+ \/mn
+
+ To compute the inverse discrete Fourier transform, defined with
+ exp(+2(pi)i(...)) in the above formula instead of exp(-
+ 2(pi)i(...)), this routine should be preceded and followed by
+ calls of C06GCF to form the complex conjugates of the data values
+ and the transform.
+
+ This routine calls C06FRF to perform multiple one-dimensional
+ discrete Fourier transforms by the fast Fourier transform (FFT)
+ algorithm in Brigham [1]. It is designed to be particularly
+ efficient on vector processors.
+
+ 4. References
+
+ [1] Brigham E O (1973) The Fast Fourier Transform. Prentice-
+ Hall.
+
+ [2] Temperton C (1983) Self-sorting Mixed-radix Fast Fourier
+ Transforms. J. Comput. Phys. 52 1--23.
+
+ 5. Parameters
+
+ 1: M -- INTEGER Input
+ On entry: the number of rows, m, of the arrays X and Y.
+ Constraint: M >= 1.
+
+ 2: N -- INTEGER Input
+ On entry: the number of columns, n, of the arrays X and Y.
+ Constraint: N >= 1.
+
+ 3: X(M,N) -- DOUBLE PRECISION array Input/Output
+
+ 4: Y(M,N) -- DOUBLE PRECISION array Input/Output
+ On entry: the real and imaginary parts of the complex data
+ values must be stored in arrays X and Y respectively. If X
+ and Y are regarded as two-dimensional arrays of dimension
+ (0:M-1,0:N-1), then X(j ,j ) and Y(j ,j ) must contain the
+ 1 2 1 2
+ real and imaginary parts of z . On exit: the real and
+ j j
+ 1 2
+ imaginary parts respectively of the corresponding elements
+ of the computed transform.
+
+ 5: INIT -- CHARACTER*1 Input
+ On entry: if the trigonometric coefficients required to
+ compute the transforms are to be calculated by the routine
+ and stored in the arrays TRIGM and TRIGN, then INIT must be
+ set equal to 'I', (Initial call).
+
+ If INIT contains 'S', (Subsequent call), then the routine
+ assumes that trigonometric coefficients for the specified
+ values of m and n are supplied in the arrays TRIGM and
+ TRIGN, having been calculated in a previous call to the
+ routine.
+
+ If INIT contains 'R', (Restart), then the routine assumes
+ that trigonometric coefficients for the particular values of
+ m and n are supplied in the arrays TRIGM and TRIGN, but does
+ not check that the routine has previously been called. This
+ option allows the TRIGM and TRIGN arrays to be stored in an
+ external file, read in and re-used without the need for a
+ call with INIT equal to 'I'. The routine carries out a
+ simple test to check that the current values of m and n are
+ compatible with the arrays TRIGM and TRIGN. Constraint: INIT
+ = 'I', 'S' or 'R'.
+
+ 6: TRIGM(2*M) -- DOUBLE PRECISION array Input/Output
+
+ 7: TRIGN(2*N) -- DOUBLE PRECISION array Input/Output
+ On entry: if INIT = 'S' or 'R',TRIGM and TRIGN must contain
+ the required coefficients calculated in a previous call of
+ the routine. Otherwise TRIGM and TRIGN need not be set.
+
+ If m=n the same array may be supplied for TRIGM and TRIGN.
+ On exit: TRIGM and TRIGN contain the required coefficients
+ (computed by the routine if INIT = 'I').
+
+ 8: WORK(2*M*N) -- DOUBLE PRECISION array Workspace
+
+ 9: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ If on entry IFAIL = 0 or -1, explanatory error messages are
+ output on the current error message unit (as defined by X04AAF).
+
+ IFAIL= 1
+ On entry M < 1.
+
+ IFAIL= 2
+ On entry N < 1.
+
+ IFAIL= 3
+ On entry INIT is not one of 'I', 'S' or 'R'.
+
+ IFAIL= 4
+ On entry INIT = 'S', but C06FUF has not previously been
+ called.
+
+ IFAIL= 5
+ On entry INIT = 'S' or 'R', but at least one of the arrays
+ TRIGM and TRIGN is inconsistent with the current value of M
+ or N.
+
+ 7. Accuracy
+
+ Some indication of accuracy can be obtained by performing a
+ subsequent inverse transform and comparing the results with the
+ original sequence (in exact arithmetic they would be identical).
+
+ 8. Further Comments
+
+ The time taken by the routine is approximately proportional to
+ mn*log(mn), but also depends on the factorization of the
+ individual dimensions m and n. The routine is somewhat faster
+ than average if their only prime factors are 2, 3 or 5; and
+ fastest of all if they are powers of 2.
+
+ 9. Example
+
+ This program reads in a bivariate sequence of complex data values
+ and prints the two-dimensional Fourier transform. It then
+ performs an inverse transform and prints the sequence so
+ obtained, which may be compared to the original data values.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06GBF(3NAG) Foundation Library (12/10/92) C06GBF(3NAG)
+
+ C06 -- Summation of Series C06GBF
+ C06GBF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06GBF forms the complex conjugate of a Hermitian sequence of n
+ data values.
+
+ 2. Specification
+
+ SUBROUTINE C06GBF (X, N, IFAIL)
+ INTEGER N, IFAIL
+ DOUBLE PRECISION X(N)
+
+ 3. Description
+
+ This is a utility routine for use in conjunction with C06EAF,
+ C06EBF, C06FAF(*) or C06FBF(*) to calculate inverse discrete
+ Fourier transforms (see the Chapter Introduction).
+
+ 4. References
+
+ None.
+
+ 5. Parameters
+
+ 1: X(N) -- DOUBLE PRECISION array Input/Output
+ On entry: if the data values z are written as x +iy and
+ j j j
+ if X is declared with bounds (0:N-1) in the (sub)program
+ from which C06GBF is called, then for 0<=j<=n/2, X(j) must
+ contain x (=x ), while for n/2= 1.
+
+ 3: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ IFAIL= 1
+ N < 1.
+
+ 7. Accuracy
+
+ Exact.
+
+ 8. Further Comments
+
+ The time taken by the routine is negligible.
+
+ 9. Example
+
+ This program reads in a sequence of real data values, calls
+ C06EAF followed by C06GBF to compute their inverse discrete
+ Fourier transform, and prints this after expanding it from
+ Hermitian form into a full complex sequence.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06GCF(3NAG) Foundation Library (12/10/92) C06GCF(3NAG)
+
+ C06 -- Summation of Series C06GCF
+ C06GCF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06GCF forms the complex conjugate of a sequence of n data
+ values.
+
+ 2. Specification
+
+ SUBROUTINE C06GCF (Y, N, IFAIL)
+ INTEGER N, IFAIL
+ DOUBLE PRECISION Y(N)
+
+ 3. Description
+
+ This is a utility routine for use in conjunction with C06ECF or
+ C06FCF(*) to calculate inverse discrete Fourier transforms (see
+ the Chapter Introduction).
+
+ 4. References
+
+ None.
+
+ 5. Parameters
+
+ 1: Y(N) -- DOUBLE PRECISION array Input/Output
+ On entry: if Y is declared with bounds (0:N-1) in the (sub)
+ program which C06GCF is called, then Y(j) must contain the
+ imaginary part of the jth data value, for 0<=j<=n-1. On
+ exit: these values are negated.
+
+ 2: N -- INTEGER Input
+ On entry: the number of data values, n. Constraint: N >= 1.
+
+ 3: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ IFAIL= 1
+ N < 1.
+
+ 7. Accuracy
+
+ Exact.
+
+ 8. Further Comments
+
+ The time taken by the routine is negligible.
+
+ 9. Example
+
+ This program reads in a sequence of complex data values and
+ prints their inverse discrete Fourier transform as computed by
+ calling C06GCF, followed by C06ECF and C06GCF again.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06GQF(3NAG) Foundation Library (12/10/92) C06GQF(3NAG)
+
+ C06 -- Summation of Series C06GQF
+ C06GQF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06GQF forms the complex conjugates of m Hermitian sequences,
+ each containing n data values.
+
+ 2. Specification
+
+ SUBROUTINE C06GQF (M, N, X, IFAIL)
+ INTEGER M, N, IFAIL
+ DOUBLE PRECISION X(M*N)
+
+ 3. Description
+
+ This is a utility routine for use in conjunction with C06FPF and
+ C06FQF to calculate inverse discrete Fourier transforms (see the
+ Chapter Introduction).
+
+ 4. References
+
+ None.
+
+ 5. Parameters
+
+ 1: M -- INTEGER Input
+ On entry: the number of Hermitian sequences to be
+ conjugated, m. Constraint: M >= 1.
+
+ 2: N -- INTEGER Input
+ On entry: the number of data values in each Hermitian
+ sequence, n. Constraint: N >= 1.
+
+ 3: X(M,N) -- DOUBLE PRECISION array Input/Output
+ On entry: the data must be stored in array X as if in a
+ two-dimensional array of dimension (1:M,0:N-1); each of the
+ m sequences is stored in a row of the array in Hermitian
+ p p p
+ form. If the n data values z are written as x +iy , then
+ j j j
+ p
+ for 0<=j<=n/2, x is contained in X(p,j), and for 1<=j<=(n-
+ j
+ p
+ 1)/2, y is contained in X(p,n-j). (See also Section 2.1.2
+ j
+ of the Chapter Introduction.) On exit: the imaginary parts
+ p p
+ y are negated. The real parts x are not referenced.
+ j j
+
+ 4: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ If on entry IFAIL = 0 or -1, explanatory error messages are
+ output on the current error message unit (as defined by X04AAF).
+
+ IFAIL= 1
+ On entry M < 1.
+
+ IFAIL= 2
+ On entry N < 1.
+
+ 7. Accuracy
+
+ Exact.
+
+ 8. Further Comments
+
+ None.
+
+ 9. Example
+
+ This program reads in sequences of real data values which are
+ assumed to be Hermitian sequences of complex data stored in
+ Hermitian form. The sequences are expanded into full complex form
+ using C06GSF and printed. The sequences are then conjugated
+ (using C06GQF) and the conjugated sequences are expanded into
+ complex form using C06GSF and printed out.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C06GSF(3NAG) Foundation Library (12/10/92) C06GSF(3NAG)
+
+ C06 -- Summation of Series C06GSF
+ C06GSF -- NAG Foundation Library Routine Document
+
+ Note: Before using this routine, please read the Users' Note for
+ your implementation to check implementation-dependent details.
+ The symbol (*) after a NAG routine name denotes a routine that is
+ not included in the Foundation Library.
+
+ 1. Purpose
+
+ C06GSF takes m Hermitian sequences, each containing n data
+ values, and forms the real and imaginary parts of the m
+ corresponding complex sequences.
+
+ 2. Specification
+
+ SUBROUTINE C06GSF (M, N, X, U, V, IFAIL)
+ INTEGER M, N, IFAIL
+ DOUBLE PRECISION X(M*N), U(M*N), V(M*N)
+
+ 3. Description
+
+ This is a utility routine for use in conjunction with C06FPF and
+ C06FQF (see the Chapter Introduction).
+
+ 4. References
+
+ None.
+
+ 5. Parameters
+
+ 1: M -- INTEGER Input
+ On entry: the number of Hermitian sequences, m, to be
+ converted into complex form. Constraint: M >= 1.
+
+ 2: N -- INTEGER Input
+ On entry: the number of data values, n, in each sequence.
+ Constraint: N >= 1.
+
+ 3: X(M,N) -- DOUBLE PRECISION array Input
+ On entry: the data must be stored in X as if in a two-
+ dimensional array of dimension (1:M,0:N-1); each of the m
+ sequences is stored in a row of the array in Hermitian form.
+ p p p
+ If the n data values z are written as x +iy , then for
+ j j j
+ p
+ 0<=j<=n/2, x is contained in X(p,j), and for 1<=j<=(n-1)/2,
+ j
+ p
+ y is contained in X(p,n-j). (See also Section 2.1.2 of the
+ j
+ Chapter Introduction.)
+
+ 4: U(M,N) -- DOUBLE PRECISION array Output
+
+ 5: V(M,N) -- DOUBLE PRECISION array Output
+ On exit: the real and imaginary parts of the m sequences of
+ length n, are stored in U and V respectively, as if in two-
+ dimensional arrays of dimension (1:M,0:N-1); each of the m
+ sequences is stored as if in a row of each array. In other
+ words, if the real parts of the pth sequence are denoted by
+ p
+ x , for j=0,1,...,n-1 then the mn elements of the array U
+ j
+ contain the values
+ 1 2 m 1 2 m 1 2 m
+ x ,x ,...,x , x ,x ,...,x ,..., x ,x ,...,x
+ 0 0 0 1 1 1 n-1 n-1 n-1
+
+ 6: IFAIL -- INTEGER Input/Output
+ On entry: IFAIL must be set to 0, -1 or 1. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ On exit: IFAIL = 0 unless the routine detects an error (see
+ Section 6).
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ If on entry IFAIL = 0 or -1, explanatory error messages are
+ output on the current error message unit (as defined by X04AAF).
+
+ IFAIL= 1
+ On entry M < 1.
+
+ IFAIL= 2
+ On entry N < 1.
+
+ 7. Accuracy
+
+ Exact.
+
+ 8. Further Comments
+
+ None.
+
+ 9. Example
+
+ This program reads in sequences of real data values which are
+ assumed to be Hermitian sequences of complex data stored in
+ Hermitian form. The sequences are then expanded into full complex
+ form using C06GSF and printed.
+
+ The example program is not reproduced here. The source code for
+ all example programs is distributed with the NAG Foundation
+ Library software and should be available on-line.
+
+@
\pagehead{NagSeriesSummationPackage}{NAGC06}
\pagepic{ps/v104nagseriessummationpackage.ps}{NAGC06}{1.00}
diff --git a/books/ps/v104applicationprograminterface.ps b/books/ps/v104applicationprograminterface.ps
new file mode 100644
index 0000000..fd9fbb7
--- /dev/null
+++ b/books/ps/v104applicationprograminterface.ps
@@ -0,0 +1,281 @@
+%!PS-Adobe-2.0
+%%Creator: Graphviz version 2.16.1 (Mon Jul 7 18:20:33 UTC 2008)
+%%For: (root) root
+%%Title: pic
+%%Pages: (atend)
+%%BoundingBox: (atend)
+%%EndComments
+save
+%%BeginProlog
+/DotDict 200 dict def
+DotDict begin
+
+/setupLatin1 {
+mark
+/EncodingVector 256 array def
+ EncodingVector 0
+
+ISOLatin1Encoding 0 255 getinterval putinterval
+EncodingVector 45 /hyphen put
+
+% Set up ISO Latin 1 character encoding
+/starnetISO {
+ dup dup findfont dup length dict begin
+ { 1 index /FID ne { def }{ pop pop } ifelse
+ } forall
+ /Encoding EncodingVector def
+ currentdict end definefont
+} def
+/Times-Roman starnetISO def
+/Times-Italic starnetISO def
+/Times-Bold starnetISO def
+/Times-BoldItalic starnetISO def
+/Helvetica starnetISO def
+/Helvetica-Oblique starnetISO def
+/Helvetica-Bold starnetISO def
+/Helvetica-BoldOblique starnetISO def
+/Courier starnetISO def
+/Courier-Oblique starnetISO def
+/Courier-Bold starnetISO def
+/Courier-BoldOblique starnetISO def
+cleartomark
+} bind def
+
+%%BeginResource: procset graphviz 0 0
+/coord-font-family /Times-Roman def
+/default-font-family /Times-Roman def
+/coordfont coord-font-family findfont 8 scalefont def
+
+/InvScaleFactor 1.0 def
+/set_scale {
+ dup 1 exch div /InvScaleFactor exch def
+ scale
+} bind def
+
+% styles
+/solid { [] 0 setdash } bind def
+/dashed { [9 InvScaleFactor mul dup ] 0 setdash } bind def
+/dotted { [1 InvScaleFactor mul 6 InvScaleFactor mul] 0 setdash } bind def
+/invis {/fill {newpath} def /stroke {newpath} def /show {pop newpath} def} bind def
+/bold { 2 setlinewidth } bind def
+/filled { } bind def
+/unfilled { } bind def
+/rounded { } bind def
+/diagonals { } bind def
+
+% hooks for setting color
+/nodecolor { sethsbcolor } bind def
+/edgecolor { sethsbcolor } bind def
+/graphcolor { sethsbcolor } bind def
+/nopcolor {pop pop pop} bind def
+
+/beginpage { % i j npages
+ /npages exch def
+ /j exch def
+ /i exch def
+ /str 10 string def
+ npages 1 gt {
+ gsave
+ coordfont setfont
+ 0 0 moveto
+ (\() show i str cvs show (,) show j str cvs show (\)) show
+ grestore
+ } if
+} bind def
+
+/set_font {
+ findfont exch
+ scalefont setfont
+} def
+
+% draw text fitted to its expected width
+/alignedtext { % width text
+ /text exch def
+ /width exch def
+ gsave
+ width 0 gt {
+ [] 0 setdash
+ text stringwidth pop width exch sub text length div 0 text ashow
+ } if
+ grestore
+} def
+
+/boxprim { % xcorner ycorner xsize ysize
+ 4 2 roll
+ moveto
+ 2 copy
+ exch 0 rlineto
+ 0 exch rlineto
+ pop neg 0 rlineto
+ closepath
+} bind def
+
+/ellipse_path {
+ /ry exch def
+ /rx exch def
+ /y exch def
+ /x exch def
+ matrix currentmatrix
+ newpath
+ x y translate
+ rx ry scale
+ 0 0 1 0 360 arc
+ setmatrix
+} bind def
+
+/endpage { showpage } bind def
+/showpage { } def
+
+/layercolorseq
+ [ % layer color sequence - darkest to lightest
+ [0 0 0]
+ [.2 .8 .8]
+ [.4 .8 .8]
+ [.6 .8 .8]
+ [.8 .8 .8]
+ ]
+def
+
+/layerlen layercolorseq length def
+
+/setlayer {/maxlayer exch def /curlayer exch def
+ layercolorseq curlayer 1 sub layerlen mod get
+ aload pop sethsbcolor
+ /nodecolor {nopcolor} def
+ /edgecolor {nopcolor} def
+ /graphcolor {nopcolor} def
+} bind def
+
+/onlayer { curlayer ne {invis} if } def
+
+/onlayers {
+ /myupper exch def
+ /mylower exch def
+ curlayer mylower lt
+ curlayer myupper gt
+ or
+ {invis} if
+} def
+
+/curlayer 0 def
+
+%%EndResource
+%%EndProlog
+%%BeginSetup
+14 default-font-family set_font
+1 setmiterlimit
+% /arrowlength 10 def
+% /arrowwidth 5 def
+
+% make sure pdfmark is harmless for PS-interpreters other than Distiller
+/pdfmark where {pop} {userdict /pdfmark /cleartomark load put} ifelse
+% make '<<' and '>>' safe on PS Level 1 devices
+/languagelevel where {pop languagelevel}{1} ifelse
+2 lt {
+ userdict (<<) cvn ([) cvn load put
+ userdict (>>) cvn ([) cvn load put
+} if
+
+%%EndSetup
+setupLatin1
+%%Page: 1 1
+%%PageBoundingBox: 36 36 114 152
+%%PageOrientation: Portrait
+0 0 1 beginpage
+gsave
+36 36 78 116 boxprim clip newpath
+1 1 set_scale 0 rotate 40 40 translate
+0.167 0.600 1.000 graphcolor
+newpath -4 -4 moveto
+-4 716 lineto
+536 716 lineto
+536 -4 lineto
+closepath fill
+1 setlinewidth
+0.167 0.600 1.000 graphcolor
+newpath -4 -4 moveto
+-4 716 lineto
+536 716 lineto
+536 -4 lineto
+closepath stroke
+% API
+gsave
+[ /Rect [ 8 72 62 108 ]
+ /Border [ 0 0 0 ]
+ /Action << /Subtype /URI /URI (bookvol10.4.pdf#nameddest=APPRULE) >>
+ /Subtype /Link
+/ANN pdfmark
+0.939 0.733 1.000 nodecolor
+newpath 62 108 moveto
+8 108 lineto
+8 72 lineto
+62 72 lineto
+closepath fill
+1 setlinewidth
+filled
+0.939 0.733 1.000 nodecolor
+newpath 62 108 moveto
+8 108 lineto
+8 72 lineto
+62 72 lineto
+closepath stroke
+0.000 0.000 0.000 nodecolor
+14.00 /Times-Roman set_font
+24 85.9 moveto 22 (API) alignedtext
+grestore
+% ORDSET
+gsave
+[ /Rect [ 0 0 70 36 ]
+ /Border [ 0 0 0 ]
+ /Action << /Subtype /URI /URI (bookvol10.2.pdf#nameddest=ORDSET) >>
+ /Subtype /Link
+/ANN pdfmark
+0.606 0.733 1.000 nodecolor
+newpath 70 36 moveto
+2.73566e-14 36 lineto
+6.33868e-15 1.06581e-14 lineto
+70 0 lineto
+closepath fill
+1 setlinewidth
+filled
+0.606 0.733 1.000 nodecolor
+newpath 70 36 moveto
+2.73566e-14 36 lineto
+6.33868e-15 1.06581e-14 lineto
+70 0 lineto
+closepath stroke
+0.000 0.000 0.000 nodecolor
+14.00 /Times-Roman set_font
+7.5 13.9 moveto 55 (ORDSET) alignedtext
+grestore
+% API->ORDSET
+gsave
+1 setlinewidth
+0.000 0.000 0.000 edgecolor
+newpath 35 72 moveto
+35 64 35 55 35 46 curveto
+stroke
+0.000 0.000 0.000 edgecolor
+newpath 38.5001 46 moveto
+35 36 lineto
+31.5001 46 lineto
+closepath fill
+1 setlinewidth
+solid
+0.000 0.000 0.000 edgecolor
+newpath 38.5001 46 moveto
+35 36 lineto
+31.5001 46 lineto
+closepath stroke
+grestore
+endpage
+showpage
+grestore
+%%PageTrailer
+%%EndPage: 1
+%%Trailer
+%%Pages: 1
+%%BoundingBox: 36 36 114 152
+end
+restore
+%%EOF
diff --git a/changelog b/changelog
index d1d652d..970fe18 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,9 @@
+20090302 tpd src/axiom-website/patches.html 20090302.02.mxr.patch
+20090302 tpd src/algebra/exposed.lsp add ApplicationProgramInterface
+20090302 tpd src/algebra/Makefile add API ApplicationProgramInterface
+20090302 tpd src/interp/util.lisp add domainsOf to browser autoload
+20090302 mxr books/bookvol10.4 add API ApplicationProgramInterface
+20090302 mxr books/ps/v104applicationprograminterface.ps
20090302 tpd src/axiom-website/patches.html 20090302.01.tpd.patch
20090302 tpd books/bookvol5 add user command documentation
20090302 tpd books/bookvol0 fix typo
diff --git a/src/algebra/Makefile.pamphlet b/src/algebra/Makefile.pamphlet
index bd760c9..08f7c1a 100644
--- a/src/algebra/Makefile.pamphlet
+++ b/src/algebra/Makefile.pamphlet
@@ -790,7 +790,7 @@ OASGP PDRING
<>=
LAYER2=\
- ${OUT}/ASP29.o \
+ ${OUT}/API.o ${OUT}/ASP29.o \
${OUT}/ATRIG.o ${OUT}/ATRIG-.o ${OUT}/BMODULE.o ${OUT}/CACHSET.o \
${OUT}/CHARNZ.o ${OUT}/CHARZ.o ${OUT}/DVARCAT.o ${OUT}/DVARCAT-.o \
${OUT}/ELEMFUN.o ${OUT}/ELEMFUN-.o ${OUT}/ESTOOLS2.o ${OUT}/EVALAB.o \
@@ -854,9 +854,10 @@ LAYER2=\
/*"ABELSG-" -> {"CABMON"; "ABELMON"; "SGROUP"; "MONOID"}*/
"ABELSG-" -> "LMODULE/SGROUP"
-"ASP29" [color="#88FF44",href="bookvol10.3.pdf#nameddest=ASP29"]
-"ASP29" -> "FORTCAT"
-/*"ASP29" -> {"TYPE"; "KOERCE"; "BOOLEAN"}*/
+"API" [color="#FF4488",href="bookvol10.4.pdf#nameddest=API"]
+/*"API" -> {"INT"; "LIST"; "ILIST"}*/
+"API" -> "ORDSET"
+/*"API" -> {"SETCAT"; "BASTYPE"; "KOERCE"}*/
"ATRIG" [color="#4488FF",href="bookvol10.2.pdf#nameddest=ATRIG"]
/*"ATRIG" -> {"RING"; "RNG"; "ABELGRP"; "CABMON"; "ABELMON"; "ABELSG"}*/
@@ -16431,6 +16432,7 @@ This keeps the regression test list in the algebra Makefile.
HELPFILE=${INT}/doc/help.helplist
SPADHELP=\
+ ${HELP}/ApplicationProgramInterface.help \
${HELP}/ArrayStack.help \
${HELP}/AssociationList.help ${HELP}/BalancedBinaryTree.help \
${HELP}/BasicOperator.help ${HELP}/BinaryExpansion.help \
@@ -16504,6 +16506,7 @@ is put into a int/Makefile.algebra and then executed by make.
TESTSYS= ${OBJ}/${SYS}/bin/interpsys
REGRESS=\
+ ApplicationProgramInterface.regress \
ArrayStack.regress \
AssociationList.regress BalancedBinaryTree.regress \
BasicOperator.regress BinaryExpansion.regress \
@@ -16589,6 +16592,18 @@ all: ${REGRESS}
@echo algebra test cases complete.
@
<>=
+${HELP}/ApplicationProgramInterface.help: ${BOOKS}/bookvol10.4.pamphlet
+ @echo 6999 create ApplicationProgramInterface.help from \
+ ${BOOKS}/bookvol10.4.pamphlet
+ @${TANGLE} -R"ApplicationProgramInterface.help" \
+ ${BOOKS}/bookvol10.4.pamphlet \
+ >${HELP}/ApplicationProgramInterface.help
+ @cp ${HELP}/ApplicationProgramInterface.help ${HELP}/API.help
+ @${TANGLE} -R"ApplicationProgramInterface.input" \
+ ${BOOKS}/bookvol10.4.pamphlet \
+ >${INPUT}/ApplicationProgramInterface.input
+ @echo "ApplicationProgramInterface (API)" >>${HELPFILE}
+
${HELP}/ArrayStack.help: ${BOOKS}/bookvol10.3.pamphlet
@echo 7000 create ArrayStack.help from ${BOOKS}/bookvol10.3.pamphlet
@${TANGLE} -R"ArrayStack.help" ${BOOKS}/bookvol10.3.pamphlet \
diff --git a/src/algebra/exposed.lsp.pamphlet b/src/algebra/exposed.lsp.pamphlet
index e211650..cc850fc 100644
--- a/src/algebra/exposed.lsp.pamphlet
+++ b/src/algebra/exposed.lsp.pamphlet
@@ -58,6 +58,7 @@
(|AlgebraGivenByStructuralConstants| . ALGSC)
(|Any| . ANY)
(|AnyFunctions1| . ANY1)
+ (|ApplicationProgramInterface| . API)
(|ArrayStack| . ASTACK)
(|AssociatedJordanAlgebra| . JORDAN)
(|AssociatedLieAlgebra| . LIE)
diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html
index c71f4fc..1838797 100644
--- a/src/axiom-website/patches.html
+++ b/src/axiom-website/patches.html
@@ -981,5 +981,7 @@ bookvol4 Hyperdoc tutorial on making new pages~~

gcl-2.6.8pre3.o.read.d.patch fix read-char-no-hang

20090302.01.tpd.patch
bookvol5 add user command documentation

+20090302.02.mxr.patch
+bookvol10.4 add API ApplicationProgrammingInterface

diff --git a/src/interp/util.lisp.pamphlet b/src/interp/util.lisp.pamphlet
index 85dac2f..a754b43 100644
--- a/src/interp/util.lisp.pamphlet
+++ b/src/interp/util.lisp.pamphlet
@@ -403,6 +403,7 @@ if you use the browse function of the {\bf hypertex} system.
|abSearch|
|detailedSearch|
|ancestorsOf|
+ |domainsOf|
|aPage|
|dbGetOrigin|
|dbGetParams|