IHAT(x):= AHAT(e):= MA(m,n):=SUM(m SUB j_,j_,1,n) rt:=[[0,r12,r13,r14,r15,r16,r17,r18,r19],[-r12,0,r23,r24,r25,r26,r27,r28,r29]~ ,[-r13,-r23,0,r34,r35,r36,r37,r38,r39],[-r14,-r24,-r34,0,r45,r46,r47,r48,r49]~ ,[-r15,-r25,-r35,-r45,0,r56,r57,r58,r59],[-r16,-r26,-r36,-r46,-r56,0,r67,r68,~ r69],[-r17,-r27,-r37,-r47,-r57,-r67,0,r78,r79],[-r18,-r28,-r38,-r48,-r58,-r68~ ,-r78,0,r89],[-r19,-r29,-r39,-r49,-r59,-r69,-r79,-r89,0]] e1:= et:=[e1,e2,e3,e4,e5,e6,e7,e8,e9] IAM(m):=IHAT(SUM(m SUB i*et SUB i,i,1,DIMENSION(m))) R_ij(i,j):=rt SUB i SUB j delta_i(i,m,n):=PRODUCT(IF(i/=j_,R_ij(j_,i)^m SUB j_,1),j_,1,n) eta_s(s,i,m,n):=-1/ABS(s)*SUM(IF(i/=j_,m SUB j_*R_ij(i,j_)^s,0),j_,1,n) "CS10(n) function is from Johann Wiesenbauer :" CS10(n,k_,n_:=2,s_:=[[[1]]],t_,u_):=PROG(LOOP(IF(n=1,exit),k_:=n_,t_:=[[n_]],~ LOOP(k_:-1,IF(k_=0,exit),u_:=s_ SUB (k_-n_),LOOP(IF(FIRST(FIRST(u_))<=k_,exit~ ),u_:=REST(u_)),t_:=APPEND(t_,VECTOR(APPEND([k_],v_),v_,u_))),s_:=ADJOIN(t_,s~ _),n_:+1,n:-1),u_:=FIRST(s_),s_:=0,LOOP(IF(u_=[],RETURN(s_)),t_:=PRODUCT(x SU~ B l_,l_,FIRST(u_)),t_:=t_/PRODUCT(l_!,l_,(FACTORS(t_)) COL 2),s_:+t_,u_:=REST~ (u_))) RE1710(s,i,m,n):=SUBST(CS10(ABS(s)),x,VECTOR(eta_s(r_*SIGN(s),i,m,n),r_,1,ABS~ (s))) C_s(s,i,m,n):=IF(s=0,1,RE1710(s,i,m,n)) E_p(p,j,s,h):=IF(p=0,1,SUM(R_ij(i_,j)*E_p(p-1,j,i_+1,h),i_,s,h)) RE1712A(q,j,h):=-2*AHAT((q+1)*et SUB j+SUM(et SUB i_,i_,1,h))+SUM((q+r_/2+1)*~ E_p(h-r_,j,1,h)*IHAT((q+r_)*et SUB j),r_,0,h) RE1712B(p,j,h):=IF(p>0,RE1712A(p-h,j,h),IF(j>h,RE1712A(p,j,h),RE1712A(p-1,j,h~ ))) RE1712F(i,p,j,h):= RE1712C(i,q,j,h):=(AHAT((q+1)*et SUB j+SUM(et SUB i_,i_,1,h))-SUM((q+r_/2+1)*~ E_p(h-r_,j,1,h)*RE1712F(i,q+r_,j,h),r_,0,h-1))/((q+h/2+1)*E_p(0,j,1,h)) RE1712D(i,q,j,h):=(AHAT((q+1)*et SUB j+SUM(et SUB i_,i_,1,h))-SUM((q+r_/2+1)*~ E_p(h-r_,j,1,h)*RE1712F(i,q+r_,j,h),r_,1,h))/((q+1)*E_p(h,j,1,h)) RE1712E(i,q,j,h):=(AHAT((q+1)*et SUB j+SUM(et SUB i_,i_,1,h))-SUM((q+r_/2+1)*~ E_p(h-r_,j,1,h)*RE1712F(i,q+r_,j,h),r_,2,h))/((q+1/2+1)*E_p(h-1,j,1,h)) RE1713A(i,q,j,h):=SUM(COMB(q,r_)*R_ij(j,i)^(q-r_)*RE1712F(i,r_,i,h),r_,0,q) RE1712F(i,p,j,h):=IF(p>0,IF(i/=j,RE1713A(i,p,j,h),IF(p=1 AND h=4,IHAT(et SUB ~ j),RE1712C(i,p-h,j,h))),IF(p>-2,IHAT(p*et SUB j),IF(j>h,RE1712D(i,p,j,h),RE17~ 12E(i,p-1,j,h)))) RE179(i,m,n,h,ma_):=PROG(ma_:=MA(m,n),IF(ma_>-1,SUM(C_s(ma_-q_,i,m,n)*RE1712F~ (i,q_,i,h),q_,0,ma_),0)+SUM(IF(-m SUB i_>0,delta_i(i_,m,n)*SUM(C_s(m SUB i_+q~ _,i_,m,n)*RE1712F(i,-q_,i_,h),q_,1,-m SUB i_),0),i_,1,n)) RE175(v,a,b,a_):=PROG(a_:=VECTOR(a SUB i_/b SUB i_,i_,1,DIM(a)),SUBST(v,APPEN~ D(VECTOR(VECTOR(R_ij(i_,j_),j_,i_+1,DIM(a)),i_,1,DIM(a)-1)),APPEND(VECTOR(VEC~ TOR(a_ SUB k_-a_ SUB l_,l_,k_+1,DIM(a)),k_,1,DIM(a)-1)))) S_t(h,a,b,t):=PRODUCT(SQRT(a SUB j_+b SUB j_*t),j_,1,h) P_t(m,n,a,b,t):=PRODUCT((a SUB j_+b SUB j_*t)^m SUB j_,j_,1,n) PS_t(m,h,n,a,b,t):=P_t(m,n,a,b,t)/S_t(h,a,b,t) BA(m,n,b):=PRODUCT(b SUB j_^m SUB j_,j_,1,n) AHAT_xy(m,h,n,a,b,x,y):=(PS_t(m,h,n,a,b,x)-PS_t(m,h,n,a,b,y))/BA(m,n,b) IHAT_t(m,h,n,a,b,t):=PS_t(m,h,n,a,b,t)/BA(m,n,b)