Transcription
View metadata, citation and similar papers at core.ac.ukbrought to you byCOREprovided by Elsevier - Publisher ConnectorJ. Symbolic Computation (1998) 25, 367–382Generation and Verification of Algorithms forSymbolic-Numeric ProcessingLADISLAV KOCBACH† AND RICHARD LISKA‡†‡Department of Physics, University of Bergen, Allégaten 55, N-5007 Bergen, NorwayFaculty of Nuclear Sciences and Physical Engineering, Czech Technical University in PragueBřehová 7, 115 19 Prague, Czech RepublicSome large scale physical computations require algorithms performing symbolic computations with a particular class of algebraic formulas in a numerical code. Developingand implementing such algorithms in a numerical programming language is a tediousand error prone task. The algorithms can be developed in a computer algebra systemand their correctness can be checked by comparison with build-in facilities of the systemso that the system is used as an advanced debugging tool. After that a numerical codefor the algorithms is automatically generated from the same source code. The proposedmethodology is explained in detail on a simple example. Real applications to calculationof matrix elements of Coulomb interaction and two-centre exchange integrals needed inatomic collision codes, are described. The method makes the developing and debuggingof such algorithms easier and faster.c 1998 Academic Press Limited1. IntroductionIn certain large scale physical numerical calculations in e.g. quantum physics, oneneeds to include analytical operations like integration or differentiation of a particulartype of algebraic expressions, e.g. products of polynomials and exponentials. Most oftensuch operations are performed manually and the resulting expressions are coded in alanguage suitable for numerical computations, in some other cases large parts are performed entirely numerically. Recently it has been found advantageous to include some ofthe symbolic manipulation inside of the numerical codes, see e.g. Hansen (1990). Suchprocedures might yield higher precision, improve the efficiency and, in particular, resultin more general codes.There exist environments which allow symbolic and numeric algorithms to be usedtogether in compiled code, e.g. AXIOM-XL, the AXIOM Extension Language (Wattet al., 1994), in Axiom (Jenks and Sutor, 1992), but using such an environment wouldnot meet the aims of the applications as it would not produce portable numerical codeand its speed would be less than the speed of a purely numerical code, coded in e.g.,FORTRAN.Developing and implementing algorithms for symbolic manipulation in a language suchas FORTRAN is a tedious and error prone task, while the computer algebra systems0747–7171/98/030367 16 25.00/0sy970182c 1998 Academic Press Limited
368L. Kocbach and R. Liska(CASs) are designed for symbolic manipulation and it is quite straightforward to carryout both purely symbolic and mixed symbolic and numerical evaluations in many of theavailable CASs.To distinguish the algorithms which perform symbolic manipulation in a numericalprogramming language from other methods, we propose calling them symbolic-numericalgorithms. These typically use the fixed length integer and floating point representations,and thus lose the absolute precision of CASs, but they evaluate large quantities of similarexpressions orders of magnitude faster than CAS and are thus applicable to large scalecalculations. By symbolic algorithms we mean the algorithms dealing with formulas andusing the facilities for formula processing implemented in CASs.This paper presents a new application of CAS to the development of numerical codes,which might contain a certain amount of symbolic manipulation. In the method presented here, the algorithms suitable for implementation in a numerical code are designedby humans and coded in language of a CAS. The algorithms are tested by comparingtheir results on representative input data sets with the results of the symbolic algorithmsincluded in the CAS. The comparison is made at the algebraic level in the CAS. Strictlyspeaking by this comparision we are not attempting to prove the correctness of the algorithms but we are verifying the correctness of their implementation on finite inputdata sets. The numerical implementation is then automatically generated from the samesource, provided that the CAS in question contains a facility to convert both mathematical functions and its own control language into a code in programming languagesuitable for the numerical applications. The method uses a CAS for advanced debugging of symbolic-numeric algorithms and also allows comparisons of the algorithms indifferent floating-point arithmetics. The main advantage of our method is the improvement and speeding up of the development-debugging process. We assert that debuggingsymbolic-numeric algorithms can be done much more effectively in a CAS than in anumeric programming language.Using such an approach one can have strong confidence that the numerical code iscorrect. So here the knowledge from the computer algebra system is used to verify thecorrectness of the proposed algorithms. In the work reported here, the computer algebrasystem REDUCE (Hearn, 1995) with the standard code generation package GENTRAN(Gates, 1986) is used to develop codes in FORTRAN.Generally we deal here with the development of a particular symbolic processing algorithm which is usually used as part of a large numerical code. Typically the algorithmdeals only with a special domain of formulas. Many papers (e.g. Cook, 1990; Dewarand Richardson, 1990; Kant, 1993; Steinberg and Roache, 1985; Wang, 1986) have dealtwith the code generation of numerical algorithms. Some work has also been reported onprogram transformation techniques (e.g. Zippel, 1992) and automatic differentiation ofnumerical codes (Rostaing et al. , 1993); however, we are unaware of any work using codegeneration of symbolic algorithms.The paper is organized as follows. Section 2 provides a general description of themethod presented for the development of verified symbolic-numeric algorithms which isexplained in full detail on an elementary example in Section 3. Two particular applications of the method to practical problems in atomic physics are described in detail inSections 4 and 5. The first application deals with the calculation of the matrix elementsof a Coulomb interaction between two bound hydrogenic states and the second dealswith the calculation of three-dimensional two-centre exchange integrals with travellingorbitals. Both examples are of interest in the impact-parameter treatment of atom–atom
Generation and Verification of Algorithms for Symbolic-Numeric Processing369or ion–atom collisions. In these applications the calculations must be repeated for manycollision parameters and fast methods for numerical evaluation are essential. At the sametime, large amounts of symbolic evaluations are needed in order to set up the formulaeused in the numerical work. The test applications of the generated symbolic-numericcodes are discussed briefly in Section 6.2. Development of Verified AlgorithmsWe need to implement a particular symbolic processing algorithm A in a numericalprogramming language L. In general this implementation and mainly its debugging couldbe quite tedious, while the algorithm A can be usually implemented very simply in acomputer algebra system (CAS) as the CAS already includes many symbolic algorithmswhich are typically used as parts of the algorithm A. The algorithm A is dealing withformulas from the domain D. To enhance the debugging and verify the correctness ofthe developed algorithm we can use the approach described in general in this section.In this general description we use a CAS and the programming language L while in theapplications we have used the CAS REDUCE and the numerical programming languageFORTRAN. For a simple example illustrating this method see Section 3.2.1. symbolic implementationThe algorithm A is implemented in a CAS. The formulas from D are represented in theCAS as its standard formulas. The implementation and debugging is usually quite simple.For debugging and verification the CAS offers much better tools than the language L.The symbolic implementation is assumed to be error free. We believe that we can assumethis as the symbolic implementaion is really very simple, see Fig. 1 and Equations (4.1)and (5.4) showing formulas implemented in practical applications.2.2. data representationTo implement the algorithm A in the language L we have to choose the representationR of formulas from the domain D in the data structures of the language L. For thismethod we further need that the used data structures and the control commands (e.g.loops and conditions) are also supported by the CAS and that the CAS supports thecode generation of these structures and control commands in the language L. Typicallythese structures include only integers, floats and their arrays, e.g. a polynomial in onevariable can be represented by an array of its coefficients. To represent the multivariatepolynomials, that appear in parts of the processed formulas, we use either sparse representation, i.e. storing coefficients and degrees, or dense representation, i.e. storing allcoefficients. These representations have been chosen so that the implementation of thesymbolic-numeric algorithms using them might be quite simple.We should note here that the representation R does not usually need to be absolutelyprecise, i.e. including big integers, as the developed symbolic-numeric algorithm will befinally used in the numerical code which does not require absolute precision. Howeverwe have to be aware of possible rounding effects during the development of a symbolicnumeric algorithm, e.g. testing whether a number is zero when a rational number isreplaced by a float number.
370L. Kocbach and R. Liska2.3. symbolic-numeric implementationThe algorithm A is implemented again in the CAS, however now we use for the formulasfrom the domain D the representation R and use only the operations and semanticssupported by the language L. It may appear strange to represent in the CAS formulas bythe representation R, e.g. polynomials by arrays of their coefficients, but this is preciselywhat is needed for verification of the symbolic-numeric implementation. The symbolicnumeric implementation is necessarily much lengthier and complicated than the symbolicimplementation. This is because many subalgorithms of the algorithm A are known tothe CAS, e.g. the addition of two polynomials, and can be directly used in the symbolicimplementation while these subalgorithms have to be coded in the symbolic-numericimplementation. In other words, the two implementations are using different tools. Inthe symbolic implementation, the CAS is used with all the facilities it supports, whilein the symbolic-numeric implementation, though coded in the CAS language, only thefacilities (data structures, semantics, algorithms) supported by the language L can beused.2.4. verificationNow we have two implementations of the algorithm A, the symbolic implementationand the symbolic-numeric implementation, both implemented in the CAS. The symbolic implementation is assumed to be correct so we can verify the symbolic-numericimplementation by comparing the results of both on a representative set of input datato algorithm A. If the two implementations produce different results, a bug from thesymbolic-numeric implementation has to be removed. The comparision is done in theprecise arithmetics of the CAS. If needed in critical cases one can check the numericalquality of symbolic-numeric implementation by comparing in the CAS the results in tworounded arithmetics of different precisions. At the end of this step we have confidencethat the symbolic-numeric implementation is error free.It might be argued that the same type of error might somehow propagate into boththe symbolic and the symbolic-numeric representation. Since two entirely different approaches are used, the chance that an accidental error would lead to a fortutious agreement over a wide range of parameters seems to be very small indeed. On the other hand,the method cannot provide a rigorous proof of the correctness, only demonstrate thevalidity in the tested domain. It does, however, represent a major improvement over themethods usually used in verification of numerical codes. One further advantage is that ifthe algorithms are to be implemented in a new application with certain change of scope,the testing can be performed once more with stress on the new features. This is shortlydiscussed below in Section 6.Usually the verification is done over input data from a direct product of small integersets. However the results also depend on some other parameters, e.g. R and a in Section 3.If we are checking the symbolic-numeric implementation in the CAS these parameters are,as parts of formulas, treated symbolically while if we were to debug these algorithms in anumerical environment we would also need to check the results for many numerical valuesof these parameters. This is one of the points which makes the development-debuggingprocess easier by the presented method.
Generation and Verification of Algorithms for Symbolic-Numeric Processing3712.5. code generationBy using a code generation facility of the CAS the required implementation of the algorithm A in the language L is automatically generated from the verified symbolic-numericimplementation. The final result is the verified source code in language L implementingthe symbolic-numeric algorithm A.The numerical procedures obtained and tested by the described method might notbe optimal, neither from the point of speed nor accuracy (e.g. avoiding accumulation ofround-off errors), but the code is, in principle, error-free in our sense. It can thus serveas a base for further optimalization which can be performed by changing the symbolicnumeric implementation with regard to these aspects. The existence of an error-freecode when developing new versions will be recognized as a great advantage by anybodywho has worked on similar problems. It should also be mentioned that some specialevaluation techniques can be encoded in the symbolic-numeric implementation and thusbe automatically verified by the described method (as, e.g., the Horner scheme evaluationof polynomials).3. Simple ExampleFor better understanding of the development method described in the previous section,a very simple example, for which the method will be presented in detail, is included here.The problem considered is to implement in FORTRAN the program which calculates theintegralZ Rxn e ax dx,(3.1)I(R, n, a) 0for input parameters R, n, a where n 0 is integer and R 0, a 0 are floatingpoint numbers. The integral can be numerically integrated, however for any n it can beevaluated tonXI(R, n, a) e aRCi Rji A(3.2)i 0which would give a faster and more precise routine. We will develop this routine byapplying our method.3.1. symbolic implementationThe REDUCE program for calculation of integral (3.1), presented in Figure 1, is reallysimple and does not need any comments.Note that in this example we could proceed by calculating (3.1) for let say n 0, . . . , 20with R, a as parameters and then generate the FORTRAN routine including the results inthe form (3.2). However, such approach has disadvantages, e.g. later we need to calculateI(R, 25, a) and have to make another code generation, and it is impossible to apply suchan approach to more complicated cases where the number of necessary formulas can bevery large (e.g. of the order 104 as in Section 4). So we need to perform the manipulationwith formulas on a numerical level in FORTRAN. The actual evaluation of integral (3.1)suitable for the numerical work is done in the alternative way described below.
372L. Kocbach and R. Liskaprocedure integ(r,n,a);% Calculates the definite integral%int 0 r x n exp(-a x) d x% Input: r,n,a - parameters of the integral, n has to be non-negative integer% Output: value of the procedure - the definite integralbeginscalar y;y : int(x n*e ( - a*x),x);return (sub(x r,y) - sub(x 0,y));end;Figure 1. Symbolic implementation of (3.1), file integ.alg3.2. data representationAll formulas needed for calculating (3.1) have the form (3.2) which we need to takewith particular values of parameters n, a keeping R as variable. Such formula will berepresented by two floats a a, abs A, an array of integer exponents oexp(i) ji 1and an array of floating-point coefficients ocof(i) Ci 1 , where we have made the shiftby 1 in indices so that the FORTRAN arrays will begin from the standard index 1.3.3. symbolic-numeric implementationTo implement the calculation of integral (3.1) in terms of array representation describedin the previous section without the use of the REDUCE operator for integration int,we need to derive explicit formula (3.2) for calculation of (3.1). Applying several timesintegration per partes we getZ R1nn(n 1) n 2 axI(R, n, a) xn e ax dx Rn e ax 2 Rn 1 e ax Re aaa30n!n!· · · n 1 e ax n 1 ,aa(the last term n!/an 1 comes from the zero limit of the integral) from which we candeduce the recurrence relations for the degrees ji and coefficients Ci and a formula forthe absolute term A in (3.2):j0 n, ji n i,n 1 i1, i 1, . . . , n,C0 , Ci Ci 1aan!A n 1 .a(3.3)Recurrence relations appear regularly in symbolic-numeric implementations. Thesymbolic-numeric implementation based on (3.3) is shown on Figure 2. Note that theprocedure pinteg could be split into two procedures, one implementing (3.3) and theother (3.2), where the first procedure has to be called only after the change of n or a. Forcomments on declarations scalar, operator, literal and declare see Appendix A.
Generation and Verification of Algorithms for Symbolic-Numeric Processingprocedure fact(n);beginliteral"c Calculates Factorial of ndeclare fact:function;fact,f:real*8;n,i:integer ;f: if n 0 then 1else for i: 1:n product i;return fend;procedure pinteg(r,n,a);beginscalar abs,res;operator ocof,oexp;literal"c Calculates the definite integralliteral"cint 0 r x n exp(-a x) d xliteral"c Input: r,n,a - parameters of the integralliteral"cn has to be non-negative integerliteral"c Output: value of the procedure - the definite integraldeclare 0),r,a:real*8 ;literal"c Symbolic calculation of the integraloexp(1) : n;ocof(1) : -1/a;for i: 2:n 1 do oexp(i) : n 1 - i;ocof(i) : ocof(i-1)*(n 2-i)/a ;abs : fact(n)/a (n 1);literal"c Evaluation of the integralres : 0;for i: 1:n 1 do res : res ocof(i)*r oexp(i);res : e (-a*r)*res abs;return !*;",cr!*;",cr!*;Figure 2. Symbolic-numeric implementation of (3.1) based on (3.3), file integ.pro3.4. verificationThe verification of the symbolic-numeric implementation (see Figure 2) has been performed by comparing its results with symbolic implementation (see Figure 1). The verification code is presented on Figure 3. Its last line gives the result zero proving that thesymbolic-numeric implementation is correct for n 0, . . . , 50 from which we assume itto be correct for all non-negative integer n.3.5. code generationHaving verified the symbolic-numeric implementation in the file integ.pro we cangenerate the FORTRAN code directly from this file. The code generation commandsshown in Figure 4 are really very simple. Note that we generate the code directly from
374L. Kocbach and R. Liskain "genproc.red"; % to read and eliminate DECLARE and LITERALin "integ.pro";% definition of pinteg - in terms of array-operatorsin "integ.alg";% definition of integ - algebraic algorithm% testing of procedures in integ.pro by comparing with symbolic calculationfor n: 0:50 sum abs(integ(r,n,a) - pinteg(r,n,a));Figure 3. Verification of symbolic-numeric implementation, file integ.tstin "genproc.red";% loads gentran, defines switch genprocon genproc;gentranout "integ.f";in "integ.pro";gentranshut "integ.f";Figure 4. Code generation of FORTRAN symbolic-numeric implementation, file integ.genthe file integ.pro. To be able to use exactly the same file without any modificationswhich could introduce errors, we have introduced a new switch, genproc, described inAppendix A. This switch and the associated actions are responsible for code generationby interfacing the code generation package GENTRAN (Gates, 1986). The generatedFORTRAN code implementing the evaluation of the integral (3.1) is presented in Figure5. This code is guaranteed (if we have not made an error in the few lines shown in Figure1) to be error free for N 50. For evaluating of (3.1) it uses the developed symbolicnumeric implementation. Limitation on the maximum n given by the array bounds (here50) can be avoided by increasing the array bounds. Then the algorithm is also assumedto be error free also for N 50.4. Matrix Elements of Coulomb InteractionThe methodology outlined in Section 2 and demonstrated on the simple example inthe previous section has been applied for preparation of numerical FORTRAN codefor calculating the matrix elements of Coulomb interaction of two bound hydrogenicstates. These calculations can be separated into radial and angular parts. For the actualcalculations codes both parts are constructed with the assistance of the CAS REDUCE,however here we discuss only the radial parts where the described techniques are used.The radial matrix elements are given byZ RZ 11n1 n 2l 2lr Rn1 l1 (r)Rn2 l2 (r) dr RR(r)Rn2 l2 (r) dr, (4.1)Ml1 l2 (l) l 1l 1 n1 l1Rr0Rwhere n1 , n2 , l1 , l2 are quantum numbers, l1 l2 l l1 l2 is the transferred angular
Generation and Verification of Algorithms for Symbolic-Numeric Processing375REAL*8 FUNCTION FACT(N)REAL*8 FINTEGER N,Ic Calculates Factorial of nIF (N.EQ.0.0) THENF 1.0ELSEF 1DO 25001 I 1,NF F*I25001CONTINUEENDIFFACT FRETURNENDREAL*8 FUNCTION PINTEG(R,N,A)INTEGER N,I,OEXP(50)REAL*8 OCOF(50),R,Ac Calculates the definite integralcint 0 r x n exp(-a x) d xc Input: r,n,a - parameters of the integralcn has to be non-negative integerc Output: value of the procedure - the definite integralc Symbolic calculation of the integralOEXP(1) NOCOF(1) -(1.0/A)DO 25002 I 2,N 1OEXP(I) N (1-I)OCOF(I) OCOF(I-1)*((N (2.0-I))/A)25002 CONTINUEABS FACT(N)/A**(N 1)c Evaluation of the integralRES 0.0DO 25003 I 1,N 1RES RES OCOF(I)*R**OEXP(I)25003 CONTINUERES EXP(REAL(-(A*R)))*RES ABSPINTEG RESRETURNENDFigure 5. FORTRAN code (file integ.f) generated by GENTRAN from the file integ.pro (Figure 2)is the resulting implementation of the symbolic-numeric algorithm in FORTRAN.momentum quantum number and Rnl (r) are radial hydrogenic functionsRnl (r) qR R̄nl (r) 0,2 (r)r 2R̄nlR̄nl (r) wheredrLn l2l 12rn rl e r/n ,(4.2)where Lkj (r) are generalized Laguere polynomialsLkj (r)dj jdr dk rk e rerdrk .(4.3)
376L. Kocbach and R. Liska4.1. symbolic implementationThe symbolic implementation for calculating the matrix elements (4.1) is done by a fewlines of code implementing formulas (4.1), (4.2) and (4.3) using the operators performingdifferentiation and integration.4.2. data representationThe limitation on l and properties of the radial hydrogenic functions Rnl (r) guaranteethat any matrix element (4.1) can be expressed asZ RZ 1P1 (r)e ar dr RlP2 (r)e ar dr,(4.4)Mln11l2n2 (l) l 1R0Rwhere a 1/n1 1/n2 and Pi (r), i 1, 2 are polynomials in r. The matrix elements(4.4) result in the formula of the formMln11l2n2 (l) e aRnXCi Rji Ca Rja(4.5)i 0which is very similar to (3.2). So for all processing we need to represent polynomials inone variable which we represent as in Section 3.2 by integer array of exponents ji andfloating point array of coefficients Ci .4.3. symbolic-numeric implementationFor the symbolic-numeric implementation, calculating the radial hydrogenic functionsRnl (r), the formulas in l 12(n l)!22l 1 (n l 1)! r/n Xi 1rl ie( 1)Rnl (r) l 2n(n l)!3n (n l i 1)!(2l i 1)!i!i 0(4.6)has been derived. The integration of a polynomial multiplied by e ar which is neededin (4.4) is transformed into a linear combination of the integrals (3.1) which are evaluated by (3.3) (actually a generalization of (3.3) working with polynomials has beendeveloped). In addition the subalgorithms for polynomial addition, multiplication andcalculation of the absolute term of a polynomial, which are too long to be reproducedhere, have been implemented in the array representation. Finally the algorithm in thearray representation for calculation of the matrix elements (4.1) has been built from allthe foregoing subalgorithms.4.4. verification and code generationThe symbolic-numeric implementation has been verified by comparison of its resultswith the symbolic implementation. The matrix elements (4.1) for the quantum numbers0 n1 6, 0 n2 n1 (formulas are symmetric in n1 , n2 ) 0 l1 n1 , 0 l2 n2 , l1 l2 l l1 l2 (these restrictions are physical limitations on quantum numbers)have been calculated identically by both implementations. The typical quantum numbersused in the applications are small, usually the greatest one is around 4, so our verification
Generation and Verification of Algorithms for Symbolic-Numeric Processing377test has covered most of the relevant region of quantum numbers. Only the verificationof the whole algorithm has been done over this range of very small integers while theverification of all used subalgorithms which were mentioned above has been done bythe same method up to the degree 20. Again as in the case of the simple example inSection 3.5 the same source file which includes the symbolic-numeric implementationhas been used for the generation of a FORTRAN symbolic-numeric implementation. Insuch a way we have constructed the FORTRAN program for the analytical calculationof matrix elements (4.1). The algorithms used in the code have been verified.5. Exchange Integrals of Heavy-particle CollisionsThe three-dimensional overlap exchange integrals in the impact-parameter treatmentof heavy-particle collisions have the form (McDowell and Coleman, 1970)Z(5.1)I(n1 , l1 , m1 , n2 , l2 , m2 ) ψn 1 l1 m1 (r1 ) eia·r1 ib·r2 ψn2 l2 m2 (r2 ) dr1where the star denotes complex conjugation and the hydrogenlike wavefunctions ψnlmwith the quantum numbers n, l, m have the formψnlm (r) Rnl (r)Ylm (r),r r ,(5.2)where Rnl are radial hydrogen functions (4.2) and Ylm are spherical harmonics functions.The position vectors r1 ,r2 measured from the two centers are related by r2 r1 R,where R is the vector connecting the two centers.The wavefunctions can be expressed in cartesian coordinates as †Xψnlm (r) e αrCj rnj xlxj y lyj z lzj .(5.3)jSubstituting this into (5.1) we get the exchange integral (5.1) as a linear combination ofintegralsi(na , lax , lay , laz , nb , lbx , lby , lbz ) Zllr1na 1 xl1ax y1ay z1laz r2nb 1 xl2bx y2by z2lbz eia·r1 ib·r2 αr1 βr2 dr1 .(5.4)Using the method in Shakeshaft (1975) the integrals (5.4) can be transformed into onedimensional integralsi(na , la , nb , lb ) 2π( i)(la lb )·1 na nbZ 1 eix·R R y dw. laa lbb α βy0where we have used notationx bw a(1 w),y α2 (1 w) β 2 w (a b)2 w(1 w),l (lx , ly , lz ),† Here n , l , l , l are integer degrees, not quantum numbers.j xj yj zj(5.5)
378L. Kocbach and R. Liska ka (k ,k ,k ) (axx ,ayy ,azz ) ax k x ay k y az k z1 (1, 1, 1)The exchange matrix elements are obtained by numerical evaluation of the one-dimensionalintegrals (5.5). The symbolic-numeric methodology is applied to the evaluation of the integrands.5.1. symbolic implementationThe expression for the integrand in formula (5.5) can be directly implemented symbolically which allows us, having also implemented calculating the wavefunctions (5.2)in the form (5.3), to calculate the exchange integrals (5.1).5.2. symbolic-numeric implementation and data representationHere prior to proposing the data representation which will be used in the symbolicnumeric implementation we have developed a new method for calculating the derivativesin (5.5). In Kocbach and Liska (1994) it has been shown that (5.5) can be written in theclosed formi(na , la , nb , lb ) 2π( i)(la lb )·1 ( 1)na nb Z 1lblaXXlalb( iw1 R)ka (iwR)kb eix·Rkakb0(5.6)ka (0,0,0) kb (0,0,0)bna /2c bnb /2cXXbM/2cXnanbMDmDmDm(2αw1 )na 2ma (2w1 )ma (2βw)nb 2mbabma 0 mb 0 m (0,0,0)(2w)mb [2(a b)ww1 ]M 2m (2ww1 )m·1Ne tR X N jA t dw ,2N t2N 1 j 0 jwhere the notation w1 1 w,t y,N na nb (M m) · 1,M la lb ka kb ,m(mx ,my ,mz ) vxmx vymy vzmz ,v (vx , vy , vz )bMc (bMx c, bMy c, bMz c), kxk(kx , ky , kz )kykz ,mxmymzm(mx , my , mz )kkxkykz DmDmDm,DmxyznXk jXnx ,ny ,nz k (jx ,jy ,jz ) nynxnzXXXkx jx ky jy kz jzis used (bn/2c denotes truncated integer part of n/2, i.e. the greatest integer n/2) andwhere the coefficients Djn and Anj are defined by the recurrence relationsD01 1,
Generation and Verification of Algorithms for Symbolic-Numeric ProcessingD0n 1,n 1Djn (n 1 2j)Dj 1 Djn 1 ,A10 1, A11 R,( 2n 1),An0 An 10n 2, j 1, . . . , bn/2c,Anj An 1(j 2n 1) RAn 1jj 1 ,379(5.7), n 2, j 0, . . . , n.The closed form formula (5.6) contains a 12-tuple summation, however, in any particularcase most of the sums reduce to a single term. Keeping only the variables w, w1 , t whichdepend on the integration variable w the exchange integral (5.1) can be, by using (5.6),written in the closed form asZ 1jXmax kXmax lXmaxe tRCjkl wj w1k tl dw.(5.8)I(n1 , l1 , m1 , n2 , l2 , m2 ) 0j 0 k 0 l lminTo represent the wavefunctions (5.3) we use a floating-point array for coefficients Cj ,four integer arrays for degrees nj , lxj , lyj , lzj and of course a number of terms in the sum(α is represented by a special way by other physical quantities). The polynomial in w, w1 , tin the res
system REDUCE (Hearn, 1995) with the standard code generation package GENTRAN (Gates, 1986) is used to develop codes in FORTRAN. Generally we deal here with the development of a particular symbolic processing al-gorithm which is usually used as part of a large numerical code. Typically the algorithm deals only with a special domain of formulas.
