Schoonschip manual examples. *end C Example 1 Z XX = (a + b)^2 *begin B a Z XX= (a + b)^2 *end C Example 2. P output !Only the first two characters of 'output' are significant. Digits 3 Rationalization 3 ! 0 leads to no rationalization attempt. Z XX=3.14 + 1/3*a1 + 7.1357689E20*a2 !E20 or E+20 implies 10^20. *yep Digits !This restores to the default option, 5 digits. Rationalization !Default: 22 digits, with check. *end C Example 3. BLOCK text{xx,yy} C This argument is substituted in a call: 'xx'. This not: xx. Arguments can be glued together: 'xx''yy'. ENDBLOCK BLOCK lines{A,B} 'A'C line 1 'A' line 2 'B' line 3 ENDBLOCK C Calling block text: text{yes,sir} DELETE text C Calling lines twice: lines{,~} lines{~,} DO var=1,3,2 P input !Normally only the first DO round is printed. C value of var: 'var' ENDDO DO var=-1,0 !The third number is 1 by default. P input C value of var: 'var'. C This may be part of a name, for example XX'var', or an index, for example YY(3+'var'). In the last case Schoonschip will interpret the % sign as minus if not on a C line and correctly work out the number. ENDDO *end C Example 4. A a1,a2,a3 !Algebraic symbol list. F f1,f2,f3 !Functions. I mu,nu !Indices. V p,q,k !Vectors. Z XX = p(mu)*{ a3^-20*a1*q(mu) + a2*a1^7*D(mu,nu)*k(nu) + a3*f2(mu,k,q)*f1} !Note the use of {}. + q(mu)*{a3*a2*a1 + f3*f2*f1} *end C Example 5. A a1,a2=i,a3=c,a4=c,a5=i=3 V p,q=z,k Z XX = Conjg(a1 + a2 + a3 + a4 + a5) Z YY = Conjg(a3C + a4C) Z ZZ = (a1 + a5)^5 Z AA = p(mu)*{p(mu) + q(mu)} + D(mu,3)*{p(mu) + q(mu)} + pDk Oldnew a4C=b4,k=K *end C Example 6. I mu,nu,la,ka Z XX = D(mu,nu)*f1(a1,a2,mu)*f2(a3,a4,nu) + f1(a1,a2,mu)*f2(a3,a4,nu)*D(mu,nu) Z YY = D(la,ka)*f1(a1,a2,la)*f2(a3,a4,ka) + f1(a1,a2,la)*f2(a3,a4,ka)*D(la,ka) Sum,la,ka *end C Example 7. A a1,a2=i,a3=c F f1,f2=i,f3=c,f4=u Z xx=Conjg{ f1(a1,a2,a3)*f2(a1,a2,a3,f3)*f4(a1,a2,a3) } *end C Example 8. Digits 20 Ratio 0 Z xx= DS(J,0,20,(1.),(1/J)) *begin Digits Ratio Z yy = DS(J,0,5,(X^J),(1/J)) *end C Example 9. X EXP(y) = DS(J,0,6,(y^J),(1/J)) Z COS = { EXP{ (i*x) } + EXP{ -(i*x) } }/2 *end C Example 10. B a11,a21,a31 D rc(n)=a11,a12,a13, a21,a22,a23, a31,a32,a33 X matrix(n,m)=rc(3*n+m-3) Z xxx = DS{J1,1,3,(DS{J2,1,3,(DS{J3,1,3, { DP(J1,J2,J3)*matrix(J1,1)*matrix(J2,2)*matrix(J3,3) } } ) } ) } *begin B a11,a21,a31,b11,b21,b31 D ra(n)=a11,a12,a13, a21,a22,a23, a31,a32,a33 D rb(n)=b11,b12,b13, b21,b22,b23, b31,b32,b33 X matrix(n,m,ra)=ra(3*n+m-3) X DET(ra) = DS{J1,1,3,(DS{J2,1,3,(DS{J3,1,3, { DP(J1,J2,J3)*matrix(J1,1,ra)*matrix(J2,2,ra)*matrix(J3,3,ra) } } ) } ) } Z detb=DET(rb) Z deta=DET(ra) *begin F f1 B b11,b21,b31 D ra(n)=(a1-a2),0, 0,(a1+a2) D rb(n)=b11,b12,b13, b21,b22,b23, b31,b32,b33 D rc(n)=c11 X matrix(k,n,m,f1)=f1(k*n+m-k) D DET(n,f1) = f1(1), DS{J1,1,2,(DS{J2,1,2, { DP(J1,J2)*matrix(n,J1,1,f1)*matrix(n,J2,2,f1) } } ) } , DS{K1,1,3,(DS{K2,1,3,(DS{K3,1,3, { DP(K1,K2,K3)*matrix(n,K1,1,f1)*matrix(n,K2,2,f1)*matrix(n,K3,3,f1) } } ) } ) } Z deta=DET(2,ra) Z detb=DET(3,rb) Z detc=DET(1,rc) *end C Example 11. T List(n)=a,b,c F f Z compl = DS(J,1,3,{ f(List(J)) } ) *begin A x=c T r1(n)=a11,a12,a13 T r2(n)=a21,a22,a23 T r3(n)=a31,a32,a33 T matr(n,m)=r1,r2,r3 T weird(n,a1,a2)=Conjg(a1+a2),Integ(3*a2) X XX(a1,a2)=a2*a1 Z sqa13 = DS(J,1,3,{f1(matr(1,J))*f1(matr(J,3))} ) Z weirdo=XX(weird(1,x,7),weird(2,x,7)) *end C Example 12. P lists Common yyy,ccc(0) Z xxx(a,b)=(a+b)^2 Z ccc(3,a,b)=(a+b)^3 Z ccc(4,a,b)=(a+b)^4 Keep xxx *next F f1,f2 V p Z yyy(p)=xxx(c,d) + p(nu)*f2(mu) Sum mu,nu *begin Write fileC *begin V q Z zzz=a1*yyy(q) Z Abc=xxx(a,b) Delete ccc,yyy !The actual delete when this section is done. *begin Z xyz=ccc(3,e,f) *end !The reader may want to delete fileC at this point. C Example 13. Z integr = a*x^2 +b*x Id,x^n~ = x2^(n+1)/(n+1) - x1^(n+1)/(n+1) *begin Z integr = (a*x^2 + b*x + c)*dx Id,x^n~*dx = x2^(n+1)/(n+1) - x1^(n+1)/(n+1) Al,dx = x2 - x1 *begin C Method to integrate expressions of the form x^n*sin(x) or x^n*cos(x). C The method is based on the equation (n even): C Integral( x^n*sin(x)) = - n! * cos(x) + Integral( sin(x) * n * { x^(n-1) - (n-2)*(n-1)*x^(n-3) + ...} ) + Integral( -cos(x) * { x^n - (n-1)*n*x^(n-2) + ...} ) and similar equations for odd n and also cos(x). The X-expressions Coefa and Coefb corresponds to the above sums enclosed in curly brackets. Below the function DK is used to separate even and odd n cases. Remember that Integ converts and truncates its argument to a short integer in the range -128, 127. X Coefa(n,m,x)=DS(j,0,m,{ x^(n-2*j-1) } , { -(n-2*j)*(n-2*j+1) } ) X Coefb(n,m,x)=DS(k,0,m,{ x^(n-2*k) } , { -(n-2*k+2)*(n-2*k+1) } ) X Icos(n,x) = DK(2*Integ(n/2),n)*{ - DB(n)*sin(x) ! n even. + cos(x)*n*Coefa(n,(n-1)/2,x) + sin(x)*Coefb(n,(n-1)/2-1,x) } + DK(2*Integ(n/2),n+1)*{ DB(n)*cos(x) ! n odd. + cos(x)*n*Coefa(n,(n-1)/2-1,x) + sin(x)*Coefb(n,(n-1)/2,x) } X Isin(n,x) = DK(2*Integ(n/2),n)*{ DB(n)*cos(x) ! n even. + sin(x)*n*Coefa(n,(n-1)/2,x) - cos(x)*Coefb(n,(n-1)/2-1,x) } + DK(2*Integ(n/2),n+1)*{ + DB(n)*sin(x) ! n odd. + sin(x)*n*Coefa(n,(n-1)/2-1,x) - cos(x)*Coefb(n,(n-1)/2,x) } Z xx = x^8*cos(x) Id,x^n~*cos(x) = Icos(n,x2) - Icos(n,x1) Al,x^n~*sin(x) = Isin(n,x2) - Isin(n,x1) *end C Example 14. V p,q A a,b,c,d,e,f,g,h I mu,nu F f1,f2,f3 X expr=f1(a,b,p)*f2(a,c,q)*f3(d,e)*f1(g,h)* { a^7 + a^-7 + a^2 + pDq^2 + pDp + p(mu) + p(nu) } *fix C Class 0, no keyword. Z xxx=expr Id,f1(a1~,a2~)=a1^10*a2^20 *begin C Class 0, keyword Always. Z xxx=expr Id,Always,f1(a1~,a2~)=a1^10*a2^20 *begin C Class 1, keyword Multi. Z xxx=expr Id,Multi,a^3 = xyz + hij *begin C Class 2, exponent 1, no keyword. Z xxx=expr Id,pDq = XYZ + HIJ *begin C Class 3, no keyword. Z xxx=expr Id,a^2 = a1^7/15 *begin C Class 5. Z xxx=expr Id,p(mu) = - q(mu) *begin C Class 6. Z xxx=expr Id,p(mu~) = - q(mu) *begin C Class 7, keyword Funct. Z xxx=expr Id,Funct,a = a27 Id,f1(a1~,a2~,a3~) = 200*a1*a2*a3 Al,f2(a1~,a2~,a3~) = a1^10*a2^11*a3^12 *begin C Class 8, keyword Once. Z xxx=expr Id,Once,a^2 = XXX *begin C Class 9. Z xxx=expr Id,a^2*f2(a1~,c,p~) = F2(a1,c,p) *begin C Class 10. Z xxx=expr Id,f1~(a,b~,p~) = F(a,b,p) *begin C Class 11, keyword Adiso. Z xxx = expr Id,Adiso,f1(g,h)*f1(a~,b~,c~) = F1(a,b,c,g,h) *begin C Class 12, keyword Ainbe. Z xxx=expr Id,Ainbe,f1(g,h)*f1(a~,b~,c~) = F1(a,b,c,g,h) Id,Ainbe,f1(a~,b~,c~)*f1(g,h) = F2(a,b,c,g,h) *begin C Class 13. Z xxx=expr Id,f1(x~,b,p)*f2(x~,c,q) = F(x,b,p,c,q) *begin C Class 15, keyword Dotpr. Z xxx=expr Id,Dotpr,p(mu~) = - q(mu) *begin C Class 16, keyword Funct. Z xxx=expr Id,Funct,p(mu~) = - F1(a,b,mu) *end C Example 15. Integration of polynomium. A a,b,c,d,e,x Z xxx = a*x^2 + b*x + c + d/x + e/x^2 IF NOT x^-1=[Log(x)] AND NOT x^n~=x^(n+1)/(n+1) Al,Addfa,x ENDIF *end C Example 16. V p,q F F I mu=N,nu=N Z xx= pDp^2 * pDq^3 * F1(p,p,q,p) * p(mu) * q(nu) Id,All,p,N,F P output *yep C Showing the dimensionality of the created indices: Id,F(i1~,i2~,i3~,i4~,i5~,i6~,i7~,i8~,i9~,i10~,i11~)= F(i1,i2,i3,i6,i7,i8,i9,i10,i11)*D(i4,i5) *begin C A more realistic example. I al,be,mu,nu A N,N_ V p,k F F,F20,F22 Z xx = G(1,1,al,be,p,al,be,p) Id,All,p,N,F P output *yep Id,F(i1~,i2~) = D(i1,i2)*F20(k) + k(i1)*k(i2)*F22(k) *end C Example 17. V p,q,k I mu A a,b,c,d,e,f,g,h F F1,F2,F3 Z xx = F1(e,d,c,b,a) + F2(e,d,c,b,a) + F3(e,d,c,b,a) + F2(-125,-30,-1,0,30,125) + F2(pDq,pDk,kDq,p(3),q(2),F2,7,mu,p,a) + f1(e,f,g,d,a,b,c,h,k) Id,Asymm,F1,2,3,4,F2,F3,4,5 Al,Asymm,f1:1-3:5-7: *end C Example 18. A a1,a2,a3,a4,alt1,agt2,aeq15,ai12 Z xxx= 0.5*a1 + 1.5*a2 + 2.5*a3 + 1.E-20*a4 IF Coef,lt,1. Al,Addfa,alt1 ENDIF IF Coef,gt,2. Al,Addfa,agt2 ENDIF IF Coef,eq,1.5 Al,Addfa,aeq15 ENDIF IF Coef,ib,1,2 Al,Addfa,ai12 ENDIF IF Coef,lt,1E-15 !Delete the term if coefficient < 1E-15. Al,Addfa,0 ENDIF *end C Example 19. A a1,a2,a3,a4,a5 F F1,F2,F3 Z xx= F3(a1,a2)*F1(a1,a2,a3)*F1(a4,a5) + F1(a1,a2,a4)*F2(a1,a2,a3)*F1(a1,a2,a3) + F3(a1,a4)*F2(a1,a2,a3)*F1(a1,a2,a4)*F2(a3,a4,a5) Id,Commu,F1,F3 *end C Example 20. F F1,F2 A a,b,c,d Z xx = F1("c,"b,"a,/,"e,*,a1) + F1("c,"b,"a,/,"e,*,a1,*,xy) + F1("c,"b,"a,/,"e,*,a1,*,xy,*,a2) + F1("c,"b,"a,/,"e,*,a1,*,xy,*,a2,*,xz) + F2("a,"b,/,"F,*,x,*,y,*,z1) + F2("b,"a,/,"F,*,x,*,y,*,z1) + F2("F,"c,"b,"a,*,*,z,*,y,*,x) + F2("F,"c,"b,"a,*,*,z,*,yp,*,x) Id,Compo,,F1,,F2 Id,Always,F2(F1~,a~,b~,c~,d~,e~) = F1(a,b,c) Id,F1~(a~,b~,z1) = F1(a,b) *begin C j1 i1 ---0-------0--- i3 j4 | |j2 | | i2 ---0-------0--- i4 j3 T Ch(n)="A,"B,"C T Cg(n)="A,"C,"B F BC,ABC A i1,i2,i3,i4,j1,j2,j3,j4 Z solu = square("A,"A,"A,"A) Id,square(c1~,c2~,c3~,c4~) = DS{L1,1,3,(DS{L2,1,3, (DS{L3,1,3,(DS{L4,1,3,( v3(c1,Ch(L1),Cg(L4),*,i1,*,-j1,*,j4)* con(Ch(L1),Cg(L1),*,j1,*,-j1)* v3(c3,Cg(L1),Ch(L2),*,-i3,*,j1,*,-j2)* con(Ch(L2),Cg(L2),*,j2,*,-j2)* v3(c4,Cg(L2),Ch(L3),*,-i4,*,j2,*,-j3)* con(Ch(L3),Cg(L3),*,j3,*,-j3)* v3(c2,Cg(L3),Ch(L4),*,i2,*,j3,*,-j4)* con(Ch(L4),Cg(L4),*,j4,*,-j4) ) } ) } ) } ) } Id,Compo,,v3,con Id,v3(f1~,a1~,a2~,a3~) = f1(a1,a2,a3) Al,con(f1~,a1~,a2~) = f1(a1) *end C Example 20. A x,y T tt(n) = "1,"2,"3,"4,"5 C Here a complicated way to make C pow(x) = a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5. Z pow(x) = DS(J,1,5,{f1(/,"a,tt(J))*x^J}) Id,Compo,f1 Id,f1(y~) = y Keep pow *next Z xx = pow(x) Id,Count,3,x,1 Keep pow *next Z xx = pow(y) Id,Count,f1,y,2,a3,10 Keep pow *next V p,q A AA Z xx = pow(y)*f1(a2,a3)*pDq^2 Id,Count,AA,y,2,f1,-4,p,1,q,3 *end C Example 22. A a,b,c,d,e Z xxx = F1(e,d,c,b,a) + F2(e,d,c,b,a) +F3(e,d,c,b,a) Id,Cyclic,F1,2,5,4 Id,Symme,F2,F3,2,3,4 *end C Example 23. V p,q,k,pp,qp,kp Z xxx = Epf(k,p,q)*Epf(kp,pp,qp) + Epf(i1,i2,i3,i4)*Epf(i1,i2,j3,j4) Id,Epfred *end C Example 24. F f1,f2,f3,f4 Z xxx = f1(-a1,-a2,a3,-a4) + f2(-a1,-a2,a3,-a4) + f3(-a1,-a2,a3,-a4) + f4(-a1,-a2,a3,-a4) Id,Even,f1,2,3,f2 Id,Odd,f3,2,3,f4 *end C Example 25. V p,q Z xxx = a1^2*pDq^3/p(4) Id,Numer,a1,2,pDq,1.E10,p(4),1.E5 *end C Example 26. A a1,a2,a2ma1 F f1 B b2,b3,b4,b5 Z xx=f1(8,4) Id,f1(n~,m~)= { b2*a1^n*a2^m + b3*a1^-n*a2^m + b4*a1^n*a2^-m + b5*a1^-n*a2^-m } Id,Ratio,a1,a2,a2ma1 *begin C Here the use of Ratio with [] names. B [b-a] Z xxx = 1/[x+a]^3 * 1/[x+b]^2 Id,Ratio,[x+a],[x+b],[b-a] P output *yep C Just checking... Id,[x+a]^n~ = (x+a)^(3+n)*(x+b)^2/[x+a]^3/[x+b]^2 Al,[x+b]^n~ = (x+a)^3*(x+b)^(2+n)/[x+a]^3/[x+b]^2 Id,b = [b-a] + a *end C Example 27. Muon decay. Masses and momenta: muon Mm,k; electron Me,p; anti-e-neutrino qp; mu-neutrino q. The length of the 3-dimensional parts of p, q and qp are denoted by pl, ql and qpl. The variable z is the cosine of the angle between p and q (3-dimensional parts). Also k0, p0 etc are the energies of the particles. Note that k4=i*k0 etc. Note also that q0=ql, qp0=qpl since the neutrino masses are taken to be zero. C Further quantities: Pi=3.14... and Mpr is the proton mass. C The evaluation is in the muon rest-system, thus kl=0 and k0=Mm. A Me,Mm,ql,qpl,pl,p0,Pi,b,z,bp,Mpr,Al V p,q,qp,k I mu,mup,i1,i2 F Dia=u Z Rate = Dia(mu,1)*Conjg(Dia(mup,10))/2 C Factor 1/2 for averaging over muon polarizations. Id,Dia(mu~,n~) = Ubg(n,Me,p)*G(n,n+1,mu,G6)*Ug(n+1,0,qp) *Ubg(n+2,0,q)*G(n+2,n+3,mu,G6)*Ug(n+3,Mm,k) Id,Spin,p,q,k,qp C There are two separate traces here. They may be unified to advantage as they have two indices in common (mu and mup). Id,Gammas,"U C Apply conservation of four-momentum qp=k-p-q and the mass-shell relations kDk=-Mm^2, pDp=-Me^2 and qDq=qpDqp=0. Id,qpDq=kDq-pDq Al,qpDk=-Mm^2-pDk-qDk Al,qpDp=pDk+Me^2-qDp Id,Funct,qp(mu~) = k(mu) - p(mu) - q(mu) *yep Id,pDq=pl*ql*z-p0*ql Al,pDk=-Mm*p0 Al,qDk=-ql*Mm Al,kDk=-Mm^2 Al,pDp=-Me^2 C Integration over all angles of the vector p gives a factor 4*Pi: Id,Addfa,4*Pi C The integration over the azimuthal angle of q gives a factor 2*Pi. Including the various factors 1/(2k0)..., a factor (2*Pi)^4/{(2*Pi)^3)}^3 and the factors from going to polar coordinates for q and p gives: Id,Addfa,2*Pi/16/Mm/p0/ql/qpl*pl^2*ql^2/32/Pi^5 C The remaining delta function is delta(k0-p0-q0-qp0). For given p0=Sqrt(pl^2+Me^2) and q0=ql this may be solved for z using qp0=qpl=length(p+q)= Sqrt(pl^2+ql^2+2*pl*ql*z) where z is the cosine of the angle of q with the third axis, taken along p. In this last comment p and q are three-vectors. C Note that integration over the delta-function gives a factor 1/Abs(F), where F is the derivative of the argument of the delta function with respect to z. This F is equal to ql*pl/qpl. Id,z=(0.5*Me^2+0.5*Mm^2-Mm*p0-ql*Mm+ql*p0)/pl/ql Al,Addfa,qpl/ql/pl *yep B Pi C The next integration is over the length ql. The endvalues are denoted by qma and qmi, with qma=(Mm-p0+pl)/2 and qmi=(Mm-p0-pl)/2. These values obtain for configurations where the electron and both neutrinos are aligned. The maximum value obtains if the momentum qp is in the direction of the electron momentum, while q is pointed in the opposite direction. Then pl+qpl=ql and q0=Mm-p0-qp0, and one solves ql(=q0)=(Mm-p0+pl)/2. C Actually, ql^-1 is not occurring here, but for completeness we include it in the integration. IF NOT ql^-1=[Log(qa/qb)] AND NOT ql^n~=qma^(n+1)/(n+1)-qmi^(n+1)/(n+1) Al,Addfa,qma-qmi ENDIF Id,qma=0.5*Mm-0.5*p0+0.5*pl Al,qmi=0.5*Mm-0.5*p0-0.5*pl *yep Id,Multi,p0^2 = pl^2+Me^2 Id,p0 = (pl^2+Me^2)/p0 C The answer contains essentially one kinematic variable, namely pl (or p0). The spectrum that follows from this equation is typical for a V-A theory. Using a variable x (shown below) this spectrum follows an equation derived by Michel provided a parameter ro in that equation, called Michel parameter, is set to 3/4. P output *yep C The final integration is over pl. A new variable x defined by x=pl+p0, with p0=Sqrt(pl^2+Me^2) is used in the following. One has dx/d(pl)=1+p/p0 and this may be rewritten as d(pl)/p0 = dx/x. Note that pl=(x+Me^2/x)/2. There are terms with and without p0^-1, and they are treated separately. The variable plmx is the maximum value of pl, i.e. plmx=(Mm-Me^2/Mm)/2. This is achieved if both neutrinos are aligned and opposite to the electron. The energy of the neutrinos is then simply pl, since the sum of the momenta must be equal to the electron momentum. In that case therefore Mm=p0+pl, which gives the stated value for plmx. Note that Mm is the endpoint value for x, corresponding to maximum pl. The minimum value for pl is 0, for x it is Me. IF p0^-1=1/x Al,pl^n~=(x/2-Me^2/x/2)^n Id,x^-1=[Log(Mm/Me)] Al,x^n~=Mm^(n+1)/(n+1)-Me^(n+1)/(n+1) ELSE ..IF NOT pl^n~=plmx^(n+1)/(n+1) ..Al,Addfa,plmx ..ENDIF ENDIF Id,plmx=Mm/2-Me^2/Mm/2 P output *yep Digits 7 C Numerical evaluation. One usually introduces the coupling constant g/Sqrt(2)/Mpr^2, and then g is to be taken such that the muon lifetime comes out to be 2.197134E-6. Note also the factor (h/(2*Pi))^-1 to convert Mev to 1/seconds. Id,Addfa,g^2/6.582173E-22/2/Mpr^4 Id,Numer,Mm,105.65946,Me,0.5110034,[Log(Mm/Me)],5.33160,Mpr,938.2796 Al,Numer,Pi,3.141592653589793238 C Multiply the Rate (= inverse lifetime) with the experimental lifetime. Then g can be solved from the requirement that the result must be one. This gives g = 1.024E-5, which corresponds to the well-known result for the fermi-coupling constant: Gf = 1.02 * 10^-5 /Sqrt(2) /Mpr^2 . Id,Addfa,2.197134E-6 P output *yep Id,g = 1.0246275E-5 *end C Example 28. C The problem is to make a series expansion of a complicated expression in terms of quantities a and b. The expression was given in terms of the quantities called X1,X2,Omx1a,X2b,Bmx2a,Fac1 and Fac2 here below. C The problem was solved in two separate runs. Here the first run is shown. The object is to give the series expansion for the quantities mentioned. Up to order 7 in a and/or b is required. The quantity ep is an expansion parameter, essentially counting the order in a and b. Terms ep^8 and higher are to be ignored, the sooner the better. For this the assignment ep=8 in the A-list is the best tool. C The calculation proceeds by building up one quantity after the next, step by step. Common X1,X2,Omx1a,X2b,Bmx2a,Fac1,Fac2 A ep=8,a,b X X(a)=ep^2*(a^2+b^2-2*a*b)-2*ep*(a+b) *fix C First certain roots called WO and WOI are needed. These are of the form Sqrt(1+XX) and 1/Sqrt(1+XX). The quantity XX itself is a function of a and b, given above as X(a). In the following DO-loop the powers of X(a) are computed: XX(1)=X(a), XX(2)=X(a)^2 etc. Z XX(0)=1. Keep XX *next DO K1=1,7 B ep Z XX('K1')=X(a)*XX('K1'-1) Nprint XX('K1') Keep XX *next ENDDO C Now the roots WO and WOI can be computed. The expansions used here are simply the expansions for Sqrt(1+x) and 1/Sqrt(1+x). B ep Z WO = 1 + XX(1)/2 - XX(2)/8 + XX(3)/16 - 5*XX(4)/128 + 7*XX(5)/256 - 21*XX(6)/1024 + 33/2048*XX(7) Z WOI = 1 - XX(1)/2 + 3*XX(2)/8 - 5*XX(3)/16 + 35*XX(4)/128 - 63*XX(5)/256 + 231*XX(6)/1024 - 429/2048*XX(7) Keep WO,WOI *next B ep Z X1=0.5+0.5*b*ep-0.5*a*ep+0.5*WO Z X2=0.5+0.5*b*ep-0.5*a*ep-0.5*WO Keep WOI *next B ep Z Omx1a=1/a/ep-X1/a/ep Z X2b=X2/b/ep Z Bmx2a=b/a-X2/a/ep Keep WOI *next B ep A ep=7 C For Fac1 and Fac2 only up to ep^6 is needed. Z Fac1=(b*ep-X1)*WOI Z Fac2=(X2-b*ep)*WOI *begin C Wrinting of common files. The printed list shows which common files were in existence here; a common file has =C after its name. The X-expression X uses two levels and has =X2 after its name. It is not written. P lists Write expP *end C Example 29. C In this second part (Example 28 is first part) we will not show the whole calculation of all expressions as they were needed, but just the evaluation of the first one, which is enough to clarify the procedure. First the output of example 27 must be entered. Enter expP A ep=7,a,b *fix B ep Z Arg(1) = X2 Z Arg(3) = (b-a)*X2/b Z Arg(4) = 1 - X2b Z Arg(5) = Omx1a*X2b - 1 Z Arg(6) = (b-a)*(1-X1)/a Z Arg(7) = (a-b)*ep*Omx1a Z Arg(8) = Bmx2a Keep Arg *next B ep C The Sp are certain functions of a variable x. Sp(1) = Log(1-x)/x, higher Sp are obtained by iteration, involving differentiation. In the end x must be set to (a-b)/a. The XX2 are powers of -Arg(1). Z Sp(1) = - Lomx/x Z XX2(1) = - Arg(1) Nprint Sp,XX2 Keep Sp,XX2,Arg *next B ep DO K1=2,6 Z Sp('K1') = Difx*Sp('K1'-1)/'K1' Z XX2('K1') = - XX2('K1'-1)*Arg(1) Id,Difx*Lomx = - 1/Omx + Lomx*Difx Id,Difx*Omx^n~ = - n*Omx^(n-1) + Omx^n*Difx Id,Difx*x^n~ = n*x^(n-1) + x^n*Difx Id,Difx = 0 C Note that only negative exponents of x occur here. Id,x^n~*Omx^-1 = DS(J1,1,-n,(x^-J1)) + 1/Omx Nprint Sp,XX2 Keep Sp,XX2,Arg *next ENDDO B Lga,Lgb DO K2=1,6 Z Sp('K2') = Sp('K2')*x^'K2' ENDDO C Omx = 1-x appears only with negative exponent. It is rewritten in terms of x/(1-x) = Xomx, which is possible here since the exponent of x is here always larger than minus the exponent of Omx. Id,x^n~*Omx^m~ = Xomx^-m*x^(n+m) Id,x = (a-b)/a Al,Xomx = (a-b)/b Id,Lomx = Lgb - Lga Nprint Sp Keep Sp,XX2,Arg *next C Now finally the first of the desired expressions is worked out. B ep,Lga,Lgb Z Exp(1) = DS(J3,1,6,{Sp(J3)*XX2(J3)}) *end