[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Axiom-developer] 20090302.02.mxr.patch (bookvol10.4 ApplicationProgramm
|
From: |
daly |
|
Subject: |
[Axiom-developer] 20090302.02.mxr.patch (bookvol10.4 ApplicationProgrammingInterface) |
|
Date: |
Tue, 3 Mar 2009 18:26:28 -0600 |
A new package, ApplicationProgrammingInterface (API) was added to the system
with a single function called getDomains. This new package will be the
cover for exporting Axiom internal functions to the algebra level.
The getDomains function was originally proposed by Martin Rubey.
This is a slightly different implementation to keep the compiler happy.
Help documentation and a regression test were written and added.
Additionally, some Nag domains were documented.
No code was changed for the Nag routines.
Tim
======================================================================
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}
+<<ApplicationProgramInterface.input>>=
+)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.help>>=
+====================================================================
+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}
+
+<<package API ApplicationProgramInterface>>=
+)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.dotabb>>=
+"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}
+<<NagPolynomialRootsPackage.help>>=
+
+ 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}
+<<NagRootFindingPackage.help>>=
+
+ 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)|<EPS,
+
+ (ii) |f(x)|<ETA.
+
+ 4. References
+
+ None.
+
+ 5. Parameters
+
+ 1: A -- DOUBLE PRECISION Input
+ On entry: the lower bound of the interval, a.
+
+ 2: B -- DOUBLE PRECISION Input
+ On entry: the upper bound of the interval, b. Constraint: B
+ /= A.
+
+ 3: EPS -- DOUBLE PRECISION Input
+ On entry: the absolute tolerance to which the zero is
+ required (see Section 3). Constraint: EPS > 0.0.
+
+ 4: ETA -- DOUBLE PRECISION Input
+ On entry: a value such that if |f(x)|<ETA, x is accepted as
+ the zero. ETA may be specified as 0.0 (see Section 7).
+
+ 5: F -- DOUBLE PRECISION FUNCTION, supplied by the user.
+ External Procedure
+ F must evaluate the function f whose zero is to be
+ determined.
+
+ Its specification is:
+
+ DOUBLE PRECISION FUNCTION F (XX)
+ DOUBLE PRECISION XX
+
+ 1: XX -- DOUBLE PRECISION Input
+ On entry: the point at which the function must be
+ evaluated.
+ F must be declared as EXTERNAL in the (sub)program from
+ which C05ADF is called. Parameters denoted as Input
+ must not be changed by this procedure.
+
+ 6: X -- DOUBLE PRECISION Output
+ On exit: the approximation to the zero.
+
+ 7: IFAIL -- INTEGER Input/Output
+ Before entry, IFAIL must be assigned a value. For users not
+ familiar with this parameter (described in the Essential
+ Introduction) the recommended value is 0.
+
+ Unless the routine detects an error (see Section 6), IFAIL
+ contains 0 on exit.
+
+ 6. Error Indicators and Warnings
+
+ Errors detected by the routine:
+
+ IFAIL= 1
+ On entry EPS <= 0.0,
+
+ or A = B,
+
+ or F(A)*F(B)>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<N*(3*N+13)/2.
+
+ IFAIL= 2
+ There have been at least 200*(N+1) evaluations of FCN.
+ Consider restarting the calculation from the final point
+ held in X.
+
+ IFAIL= 3
+ No further improvement in the approximate solution X is
+ possible; XTOL is too small.
+
+ IFAIL= 4
+ The iteration is not making good progress. This failure exit
+ may indicate that the system does not have a zero, or that
+ the solution is very close to the origin (see Section 7).
+ Otherwise, rerunning C05NBF from a different starting point
+ may avoid the region of difficulty.
+
+ 7. Accuracy
+
+ ^
+ If x is the true solution, C05NBF tries to ensure that
+
+ ^ ^
+ ||x-x||<=XTOL*||x||.
+
+ -k
+ If this condition is satisfied with XTOL=10 , then the larger
+ components of x have k significant decimal digits. There is a
+ danger that the smaller components of x may have large relative
+ errors, but the fast rate of convergence of C05NBF usually avoids
+ this possibility.
+
+ If XTOL is less than machine precision, and the above test is
+ satisfied with the machine precision in place of XTOL, then the
+ routine exits with IFAIL = 3.
+
+ Note: this convergence test is based purely on relative error,
+ and may not indicate convergence if the solution is very close to
+ the origin.
+
+ The test assumes that the functions are reasonably well behaved.
+ If this condition is not satisfied, then C05NBF may incorrectly
+ indicate convergence. The validity of the answer can be checked,
+ for example, by rerunning C05NBF with a tighter tolerance.
+
+ 8. Further Comments
+
+ The time required by C05NBF to solve a given problem depends on n
+ , the behaviour of the functions, the accuracy requested and the
+ starting point. The number of arithmetic operations executed by
+ 2
+ C05NBF to process each call of FCN is about 11.5*n . Unless FCN
+ can be evaluated quickly, the timing of C05NBF will be strongly
+ influenced by the time spent in FCN.
+
+ Ideally the problem should be scaled so that at the solution the
+ function values are of comparable magnitude.
+
+ 9. Example
+
+ To determine the values x ,..., x which satisfy the tridiagonal
+ 1 9
+ equations:
+
+ (3-2x )x -2x =-1
+ 1 1 2
+
+ -x -1+(3-2x )x -2x =-1, i=2,3,...,8
+ i i i i+1
+
+ -x +(3-2x )x =-1.
+ 8 9 9
+
+ 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.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ C05PBF(3NAG) Foundation Library (12/10/92) C05PBF(3NAG)
+
+ C05 -- Roots of One or More Transcendental Equations C05PBF
+ C05PBF -- 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
+
+ C05PBF is an easy-to-use routine to find a solution of a system
+ of nonlinear equations by a modification of the Powell hybrid
+ method. The user must provide the Jacobian.
+
+ 2. Specification
+
+ SUBROUTINE C05PBF (FCN, N, X, FVEC, FJAC, LDFJAC, XTOL,
+ 1 WA, LWA, IFAIL)
+ INTEGER N, LDFJAC, LWA, IFAIL
+ DOUBLE PRECISION X(N), FVEC(N), FJAC(LDFJAC,N), XTOL, WA
+ 1 (LWA)
+ EXTERNAL FCN
+
+ 3. Description
+
+ The system of equations is defined as:
+
+ f (x ,x ,...,x )=0, i=1,2,...,n.
+ i 1 2 n
+
+ C05PBF is based upon the MINPACK routine HYBRJ1 (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 calculated, but it is not
+ recalculated 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
+ Depending upon the value of IFLAG, FCN must either return
+ the values of the functions f at a point x or return the
+ i
+ Jacobian at x.
+
+ Its specification is:
+
+ SUBROUTINE FCN (N, X, FVEC, FJAC, LDFJAC, IFLAG)
+ INTEGER N, LDFJAC, IFLAG
+ DOUBLE PRECISION X(N), FVEC(N), FJAC(LDFJAC,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 or the Jacobian must be evaluated.
+
+ 3: FVEC(N) -- DOUBLE PRECISION array Output
+ On exit: if IFLAG = 1 on entry, FVEC must contain the
+ function values f (x) (unless IFLAG is set to a
+ i
+ negative value by FCN). If IFLAG = 2 on entry, FVEC
+ must not be changed.
+
+ 4: FJAC(LDFJAC,N) -- DOUBLE PRECISION array Output
+ On exit: if IFLAG = 2 on entry, FJAC(i,j) must contain
+ ddf
+ i
+ the value of ---- at the point x, for i,j=1,2,...,n
+ ddx
+ j
+ (unless IFLAG is set to a negative value by FCN).
+
+ If IFLAG = 1 on entry, FJAC must not be changed.
+
+ 5: LDFJAC -- INTEGER Input
+ On entry: the first dimension of FJAC.
+
+ 6: IFLAG -- INTEGER Input/Output
+ On entry: IFLAG = 1 or 2:
+ if IFLAG = 1, FVEC is to be updated;
+
+ if IFLAG = 2, FJAC is to be updated.
+ 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 C05PBF 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: 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<N*(N+13)/2.
+
+ IFAIL= 2
+ There have been 100*(N+1) evaluations of the functions.
+ Consider restarting the calculation from the final point
+ held in X.
+
+ IFAIL= 3
+ No further improvement in the approximate solution X is
+ possible; XTOL is too small.
+
+ IFAIL= 4
+ The iteration is not making good progress. This failure exit
+ may indicate that the system does not have a zero or that
+ the solution is very close to the origin (see Section 7).
+ Otherwise, rerunning C05PBF from a different starting point
+ may avoid the region of difficulty.
+
+ 7. Accuracy
+
+ ^
+ If x is the true solution, C05PBF tries to ensure that
+
+ ^ ^
+ ||x-x|| <=XTOL*||x|| .
+ 2 2
+
+ -k
+ If this condition is satisfied with XTOL=10 , then the larger
+ components of x have k significant decimal digits. There is a
+ danger that the smaller components of x may have large relative
+ errors, but the fast rate of convergence of C05PBF usually avoids
+ the possibility.
+
+ If XTOL is less than machine precision and the above test is
+ satisfied with the machine precision in place of XTOL, then the
+ routine exits with IFAIL = 3.
+
+ Note: this convergence test is based purely on relative error,
+ and may not indicate convergence if the solution is very close to
+ the origin.
+
+ The test assumes that the functions and Jacobian are coded
+ consistently and that the functions are reasonably well behaved.
+ If these conditions are not satisfied then C05PBF may incorrectly
+ indicate convergence. The coding of the Jacobian can be checked
+ using C05ZAF. If the Jacobian is coded correctly, then the
+ validity of the answer can be checked by rerunning C05PBF with a
+ tighter tolerance.
+
+ 8. Further Comments
+
+ The time required by C05PBF to solve a given problem depends on n
+ , the behaviour of the functions, the accuracy requested and the
+ starting point. The number of arithmetic operations executed by
+ 2
+ C05PBF is about 11.5*n to process each evaluation of the
+ 3
+ functions and about 1.3*n to process each evaluation of the
+ Jacobian. Unless FCN can be evaluated quickly, the timing of
+ C05PBF will be strongly influenced by the time spent in FCN.
+
+ Ideally the problem should be scaled so that, at the solution,
+ the function values are of comparable magnitude.
+
+ 9. Example
+
+ To determine the values x ,..., x which satisfy the tridiagonal
+ 1 9
+ equations:
+
+ (3-2x )x -2x =-1
+ 1 1 2
+
+ -x +(3-2x )x -2x =-1, i=2,3,...,8.
+ i-1 i i i+1
+
+ -x +(3-2x )x =-1.
+ 8 9 9
+
+ 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{NagRootFindingPackage}{NAGC05}
\pagepic{ps/v104nagrootfindingpackage.ps}{NAGC05}{1.00}
@@ -63664,6 +65282,2213 @@ NagRootFindingPackage(): Exports == Implementation
where
@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{package NAGC06 NagSeriesSummationPackage}
+<<NagSeriesSummationPackage.help>>=
+
+ 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<s<+infty.
+ /
+ -infty
+
+ For more general advice on the use of Fourier transforms, see
+ Hamming [2]; more detailed information on the fast Fourier
+ transform algorithm can be found in Van Loan [3] and Brigham [1].
+
+ 2.2. References
+
+ [1] Brigham E O (1973) The Fast Fourier Transform. Prentice-
+ Hall.
+
+ [2] Hamming R W (1962) Numerical Methods for Scientists and
+ Engineers. McGraw-Hill.
+
+ [3] Van Loan C (1992) Computational Frameworks for the Fast
+ Fourier Transform. SIAM Philadelphia.
+
+ 3. Recommendations on Choice and Use of Routines
+
+ 3.1. One-dimensional Fourier Transforms
+
+ The choice of routine is determined first of all by whether the
+ data values constitute a real, Hermitian or general complex
+ sequence. It is wasteful of time and storage to use an
+ inappropriate routine.
+
+ Two groups, each of three routines, are provided
+
+ Group 1 Group 2
+
+ Real sequences C06EAF C06FPF
+
+ Hermitian sequences C06EBF C06FQF
+
+ General complex C06ECF C06FRF
+ sequences
+
+ Group 1 routines each compute a single transform of length n,
+ without requiring any extra working storage. The Group 1 routines
+ impose some restrictions on the value of n, namely that no prime
+ factor of n may exceed 19 and the total number of prime factors
+ (including repetitions) may not exceed 20 (though the latter
+ 6
+ restriction only becomes relevant when n>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<j<=n-1, X(j) must contain
+ j n-j
+ -y (=y ). In other words, X must contain the Hermitian
+ j n-j
+ sequence in Hermitian form. (See also Section 2.1.2 of the
+ Chapter Introduction). On exit: the imaginary parts y are
+ j
+ negated. The real parts x are not referenced.
+ j
+
+ 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 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>>=
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.
@
<<spadhelp>>=
+${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<br/>
gcl-2.6.8pre3.o.read.d.patch fix read-char-no-hang<br/>
<a href="patches/20090302.01.tpd.patch">20090302.01.tpd.patch</a>
bookvol5 add user command documentation<br/>
+<a href="patches/20090302.02.mxr.patch">20090302.02.mxr.patch</a>
+bookvol10.4 add API ApplicationProgrammingInterface<br/>
</body>
</html>
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|
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Axiom-developer] 20090302.02.mxr.patch (bookvol10.4 ApplicationProgrammingInterface),
daly <=