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