'''
plplot.py Legendre polys plot
l Pl(z)=(2l-1) z P(l-1) - (l-1)P(l-2)
Pl(z):= (1/2^l l!) d^l (zz-1)^l
Plm(z,x) = (-)^m x^m d^m Pl(z) where x= sqrt(1-zz)~ sin(t); xx+zz=1
P0= 1
P1=z P11=-x
P2=(3zz-1)/2; P21= -3xz P22=3xx
P3=(5zzz-3z)/2 P31=-3/2x(5zz-1) P32=15zxx; P33=-15xxx
P4=(35zzzz-30zz+3)/8
'''
from math import exp
def div(p,q) : return int(p/q)
def divide(p,q) : # improve integer fractions
def hasfact(x,z) : return (x>0) and ((x % z)==0)
pn=[2,3,5,7,11,13,17,19,23,29,31,37]; ln= len(pn) # small prime nbrs
p,q,s = (abs(p),abs(q),-1) if ((p*q)<0) else (abs(p),abs(q),1)
if p==0 : q=1
for i in range(ln) :
z=pn[i]
while hasfact(p,z) and hasfact(q,z) : p=div(p,z); q=div(q,z)
return [s*p,q]
def polynscale(u,p,q) : # scale factor on polynomial
n=len(u); z=[[]]*n
for i in range(n) : z[i]= divide(p*u[i][0],q*u[i][1])
return z
def mult(a,b) : # multiply rationals or integers, result rational
ta=str(type(a)); inta= (ta.find('int')>0)
tb=str(type(b)); intb= (tb.find('int')>0)
pa,qa= (a,1) if inta else (a[0],a[1])
pb,qb= (b,1) if intb else (b[0],b[1])
return [pa*pb,qa*qb]
def add(a,b) : # add rational or int
ta=str(type(a)); inta= (ta.find('int')>0)
tb=str(type(b)); intb= (tb.find('int')>0)
pa,qa= (a,1) if inta else (a[0],a[1])
pb,qb= (b,1) if intb else (b[0],b[1])
return divide(pa*qb+pb*qa,qa*qb)
def polynadd(u,v,a,b) : # rational-coeff polynomial added, with int? factors
n=len(u); m=len(v); k=max(n,m); w=[[]]*k
for i in range(k):
uu= mult(a,u[i]) if (i<n) else [0,1]
vv= mult(b,v[i]) if (i<m) else [0,1]
q=uu[1]*vv[1]; p=uu[0]*vv[1]+vv[0]*uu[1]; w[i]= divide(p,q)
return w
def polynmul(u,v) : # multiply polynomials
n=len(u); m=len(v)
if (n==0)or(m==0) : w=[[0,1]]
else :
w=[[]]*(n+m-1)
for i in range(n) :
for k in range(m) :
a=u[i][0];b=u[i][1]; c=v[k][0];d=v[k][1]; z= divide(a*c,b*d); j=i+k
if w[j]==[] : w[j]=z # first use of j
else : p=z[0]*w[j][1]+z[1]*w[j][0]; q=z[1]*w[j][1]; w[j]=divide(p,q)
return w
def polyndiff(u) : # derivative of polynomial
n=len(u); w=[[]]*(n-1)
for i in range(1,n) : w[i-1]= [i*u[i][0],u[i][1]]
return w
def plzero(n) : # rational polys of z
# l Pl(z)=(2l-1) z P(l-1) - (l-1)P(l-2)
# l Pl(z,x)= (2l-1) z P(l-1)(z,x) - (l-1)(zz+[1~=xx])P(l-2)
pl=[[]]*20; pl[0]= [[1,1]]; z= [[0,1],[1,1]]; un=[[1,1],[0,1],[1,1]]
for i in range(1,n) :
y= polynmul(pl[i-1],z); y= polynscale(y, 2*i-1,1)
if (i-2)>= 0 : w= polynscale(pl[i-2],i-1,1); y=polynadd(y,w,1,-1)
pl[i]= polynscale(y,1,i)
return pl
def pleval(p,zmin,zmax,dz) :
nz= int((zmax-zmin)/(dz+0.0)) + 1; m=len(p); xy=[0]*(2*nz)
for i in range(nz) :
x=zmin+i*dz; y=0; zz=1
for k in range(m) : y += zz* p[k][0]/(p[k][1]+0.0); zz *=x
xy[2*i]=x; xy[2*i+1]=y
return xy
class svgdump :
def plotbegin(sd, w,h,f) : # args: width height zoomfactor
W=str(w); H=str(h); u=''; fw=str(f*w); fh=str(f*h)
s=( '<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n'+
'<!DOCTYPE svg>\n'+
'<svg xmlns="http://www.w3.org/2000/svg"\n'+
' xmlns:xlink="http://www.w3.org/1999/xlink"\n'+
' width="'+W+'" height="'+H+'" viewBox="0 0 '+fw+' '+fh+'">\n')
return s
def plotend(sd) : return '</svg>\n'
def move(sd,x,y) : return 'M'+str(x)+','+str(y)
def line(sd,x,y) : return 'L'+str(x)+','+str(y)
def labl(sd,s,x,y, size=14,fill='#000') :
f=' font-family="Arial"'
t='<tspan x="'+str(x)+'" y="'+str(y)+'">'+s+'</tspan>'; fsz= str(size)
return ('<text fill="'+fill+'"'+f+' font-size="'+fsz+'px">'+t+'</text>\n')
def path(sd,w,s,lc='#000') : # path w=with lc=linecolor
t='<path fill="none" stroke="'+lc+'" stroke-width="'+str(w)+'" '
t+= ' stroke-linecap="round" '
return t + 'd="\n'+s+'" />\n'
# end class svgdump
def aprx(x) : # 3-digit approximation
y= (-x) if(x<0) else x; z= int(1e3*y+1e-3); z= (-z) if (x<0) else z
return z/1000.0
def curve(sd,scale,xy, fill='#000', width=1) :
(xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= div(len(xy),2)
for i in range(n) :
x= aprx(xmin+fx*(xy[2*i]-xs)); y= aprx(ymin+fy*(xy[2*i+1]-ys))
s+= sd.move(x,y) if(i==0) else sd.line(x,y)
return sd.path(width,s,lc=fill)
def miniplot(fn, xys, text) :
colo=['#000','#000','#00a','#077','#0a0','#770','#a00']
xs,fx,xmin= -1, 750/2.0 ,25 # xscale start in user units, min in pixels
ys,fy,ymin= -1, -550/2.0 ,575 ; scale=(xs,fx,xmin, ys,fy,ymin)
sd=svgdump(); buf=sd.plotbegin(800,600,1)
buf += sd.labl(text,50,50); i=0
for xy in xys : buf += curve(sd,scale,xy, fill=colo[i%7]); i +=1
g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()
def legendretest() :
xyf=[1,-1, -1,-1, -1,1, 1,1, 1,-1, 0,-1, 0,1, 1,1, 1,0, -1,0]; xys=[xyf]
pl=plzero(6)
for i in range(6) : xy=pleval(pl[i], -1,1, 0.01); xys += [xy]
miniplot('plplot', xys, 'Legendre Pl_1 to Pl_5')
# radial funcs u(x)/x= x^l exp(-x) radpol(n,l,x) Laguerre-type polyns, n>l
def division(p,q) : # integer fractions
def hasfact(x,z) : return (x>0) and ((x % z)==0)
def div(p,q) : return int(p/q)
pn=[2,3,5,7,11,13,17,19,23]; ln= len(pn) # small prime nbrs
p,q,s = (abs(p),abs(q),-1) if ((p*q)<0) else (abs(p),abs(q),1)
for i in range(ln) :
z=pn[i]
while hasfact(p,z) and hasfact(q,z) : p=div(p,z); q=div(q,z)
return s*p,q
def radialpol(nmax=5) : # coeffs of radial hydrogen atom polyns
print('Radialpolynome.'); pols=[]; # n=1,2,3,4, l=0..n-1, a0=1.
for n in range(1,nmax+1) :
for l in range(n) :
j=2*l+2; k=0; p=1; q=1; t=''; pol=[]
while p != 0 :
p,q= division(p,q); pol += [[p,q]]
t += ' '+ (str(p) if(q==1) else (str(p)+'/'+str(q)))
q *= (k+1)*(k+j); p *= (-2*n+2*k+j); k+=1
print('n='+str(n)+' l='+str(l)+' '+t)
pols += [[n,l,pol]]
return pols
def radiallist(nlist,mode=0) : # laguerre-type radial functions
def pol(q,x) : # polynomial q
y=0; zz=1; m=len(q)
for k in range(m) : y += zz* q[k][0]/(q[k][1]+0.0); zz *=x
return y
testrun= (mode==0)
p=radialpol(); m=len(p);
dx=0.05; xmax=15.0;xmin=0.0; ymin=-1;ymax=1; xys=[]; minmax=[]; k=0
for i in range(m) :
n=p[i][0]; l=p[i][1]; q=p[i][2];
if n==nlist : # scan x^l expm(x) pol(q,x), 200 points xy[l]
x=0; y1=0;y=0; j=0; ymi=0;yma=0
while (x<1)or(abs(y)>1e-3)or(y1>1e-3) :
z= exp(-x)*pol(q,x)
for k in range(l) : z *=x
y1=abs(y); y=z; x +=dx; j+=1; ymi=min(ymi,y); yma=max(yma,y)
#print('n,l,j,min,max='+str([n,l,j,ymi,yma])); minmax += [(ymi,yma)]
if not testrun :
x=0; xy=[]; a,b= ymi,yma; b= (-a) if((-a)>b) else b; b=1.01*b
while x <= xmax :
z= exp(-x)*pol(q,x)
for k in range(l) : z *=x
xy += [x,z/b]; x += dx
xys += [xy]
k+=1 # after defining xys[k] and/or minmax[k]
xyf=[ xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin,
xmin,0, xmax,0] # frame
xs,fx,xpix= xmin, 750/(xmax-xmin),25 # scale start user units, min in pixels
ys,fy,ypix= ymin,-550/(ymax-ymin),575 ; scale=(xs,fx,xpix, ys,fy,ypix)
colo=['#000','#000','#00a','#077','#0a0','#770','#a00']
fn='radialf'; text='Radial n=5 l=0,1,2,3,4'; xys= [xyf] + xys
if not testrun : # strange plot, case l=0 tiny, l=4 huge.
sd=svgdump(); buf=sd.plotbegin(800,600,1)
buf += sd.labl(text,50,500); i=0
for t in xys : buf += curve(sd,scale, t, fill=colo[i%7]); i +=1
g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()
legendretest()
radiallist(5,1)
# Legendre polynomials, rational-number arith
def div(p,q) : return int(p/q)
def getprimes(n) : # want n prime numbers
p=[2]; i=1; plast=2
while i<=n :
q=plast+1; found=False
while not found:
j=0
while ((q%p[j])!=0) and ((p[j]*p[j])<=q) : j +=1
if (q%p[j] != 0) : p += [q]; plast=q; found=True; i+=1
else : q+=1
return p
def greatestdiv(x) : # greatest divisor, skip 0s in list
n=len(x); y=[0]*n; w=[0]*n; h=0; ymin=-1; gcd=1
pn=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]; ln= len(pn) # small primes
for i in range(n) :
z=abs(x[i])
if z>0 : ymin= z if(ymin<0) else min(ymin,z); w[h]=z; h+=1
w=w[:h]; ip=0
while (ip<ln)and(pn[ip]<=ymin) :
z=pn[ip]; ok=True
while ok :
for i in range(h) : ok= ok and ((w[i]%z)==0)
if ok :
for i in range(h) : w[i]= div(w[i],z)
ymin= div(ymin,z); gcd *= z
ip +=1
return gcd
def smallestmult(x) : # smallest multiple of nonzero numbers
n=len(x); y=abs(x[0])
for i in range(1,n) :
z=abs(x[i]); y = y if (z==0) else div(y*z,greatestdiv([y,z]))
return y
def divide(p,q) : # improve integer fractions
def hasfact(x,z) : return (x>0) and ((x % z)==0)
pn=[2,3,5,7,11,13,17,19,23,29,31,37]; ln= len(pn) # small prime nbrs
p,q,s = (abs(p),abs(q),-1) if ((p*q)<0) else (abs(p),abs(q),1)
if p==0 : q=1
for i in range(ln) :
z=pn[i]
while hasfact(p,z) and hasfact(q,z) : p=div(p,z); q=div(q,z)
return [s*p,q]
def intsolve(ai, bi) : # solve integer equation, return fractions as pairs
# ai is p-by-p matrix, bi is p-vector. solve Ax=b
def fracadd(x,y,f) : # x + f*y . x,y fractions, f integer.
a,b= tuple(x); c,d= tuple(y); c *=f; return divide(a*d+c*b,b*d)
p=len(bi); xg=[[]]*p; bg=[0]*p; ok= True; ag=[[]]*p; j=0
for i in range(p) : # local copy
ag[i]=[0]*p; bg[i]= int(bi[i])
for k in range(p) : ag[i][k] = int(ai[i][k])
while ok and (j<p) : # kill lower triangle columns, row by row
i=j # get pivot. swap first non-zero ag[j][j] row up
while (ag[i][j]==0) and (i<p) : i +=1
ok= (i<p)
if ok and (i>j) :
for h in range(p) : z=ag[i][h]; ag[i][h]=ag[j][h]; ag[j][h]=z
z=bg[i]; bg[i]=bg[j]; bg[j]=z
for i in range(j+1,p) :
rp= ag[i][j]; rq= ag[j][j] # line i = rq*(line i) - rp* (line j)
bg[i]= rq*bg[i]-rp*bg[j]; ag[i][j]= 0; ok= ok and (rq!=0)
for l in range(j+1,p) : ag[i][l]= rq*ag[i][l]-rp*ag[j][l]
j +=1
if ok : # matrix ag is upper triangular
i=p-1
while i>=0 : # make vector xg: (ag) * (xg) = (bg)
ax= [0,1]
for j in range(i+1,p) : ax= fracadd(ax,xg[j],ag[i][j])
z= fracadd([bg[i],1],ax,-1); xg[i]= divide(z[0],z[1]*ag[i][i]); i -=1
return ok, xg # junk if ok=False
def polynscale(u,p,q) : # scale factor on polynomial
n=len(u); z=[[]]*n
for i in range(n) : z[i]= divide(p*u[i][0],q*u[i][1])
return z
def mult(a,b) : # multiply rationals or integers, result rational
ta=str(type(a)); inta= (ta.find('int')>0)
tb=str(type(b)); intb= (tb.find('int')>0)
pa,qa= (a,1) if inta else (a[0],a[1])
pb,qb= (b,1) if intb else (b[0],b[1])
return [pa*pb,qa*qb]
def add(a,b) : # add rational or int
ta=str(type(a)); inta= (ta.find('int')>0)
tb=str(type(b)); intb= (tb.find('int')>0)
pa,qa= (a,1) if inta else (a[0],a[1])
pb,qb= (b,1) if intb else (b[0],b[1])
return divide(pa*qb+pb*qa,qa*qb)
def polynadd(u,v,a,b) : # rational-coeff polynomial added, with int? factors
n=len(u); m=len(v); k=max(n,m); w=[[]]*k
for i in range(k):
uu= mult(a,u[i]) if (i<n) else [0,1]
vv= mult(b,v[i]) if (i<m) else [0,1]
q=uu[1]*vv[1]; p=uu[0]*vv[1]+vv[0]*uu[1]; w[i]= divide(p,q)
return w
def polynmul(u,v) : # multiply polynomials
n=len(u); m=len(v); w=[[]]*(n+m-1)
for i in range(n) :
for k in range(m) :
a=u[i][0];b=u[i][1]; c=v[k][0];d=v[k][1]; z= divide(a*c,b*d); j=i+k
if w[j]==[] : w[j]=z # first use of j
else : p=z[0]*w[j][1]+z[1]*w[j][0]; q=z[1]*w[j][1]; w[j]=divide(p,q)
return w
def polyndiff(u) : # derivative of polynomial
n=len(u); w=[[]]*(n-1)
for i in range(1,n) : w[i-1]= [i*u[i][0],u[i][1]]
return w
def wholenb(x) : # x is a list rationals, scale to get integers only. bad.
n=len(x); q=1; y=[0]*n
for i in range(n) : r=x[i][1]; z=divide(q,r); q = z[0]*r
for i in range(n) : z=divide(q*x[i][0],x[i][1]); y[i]= z[0]
if y[0]<0 :
for i in range(n) : y[i]= -y[i]
return y
def polynpurge(u) : # omit leading zero terms and scale to integers
# return new polyn and scalefactor applied
n=len(u); j=n-1; k=0
while (u[j][0]==0)and (j>=0) : j -=1
while (u[k][0]==0)and (k<n) : k+=1
j +=1; v= u[:j]; p=[0]*j; q=[0]*j; sign= (-1) if (u[k][0]<0) else 1
for i in range(j) : p[i]= v[i][0]; q[i]=v[i][1]
z=smallestmult(q)
for i in range(j) : r= divide(p[i]*z,q[i]); p[i]=r[0]
y=greatestdiv(p);
for i in range(j) : r=divide(p[i],y); v[i]=[sign*r[0],1]
return v, divide(z,y)
def present(s,q,mode,noisy) :
# mode 0 no x, mode 1 complementary x, mode -n const x^n, mode>0 l=mode
n=len(q); dmax=1; j= (n-1) if (mode==1) else 0; s+='('
if mode>0 : j=mode # new
v= '' if (mode >=0) else ('x'+str(-mode))
for i in range(n) : dmax= q[i][1] if (q[i][1]>dmax) else dmax
for i in range(n) :
a=q[i][0]; b= q[i][1]; bad= (dmax % b) != 0; sig=''
if bad : print('present bug exit'); exit()
t=('x'+str(j)) if (j>0) else ''; u=('z'+str(i)) if (i>0) else ''
if a!=0 : cf= div(dmax,b)*a; sig='+' if (cf>0) else ''
if a!=0 : s+= ' '+sig+str(cf)+v+t+u
j -=1
s+=')/'+str(dmax)
if noisy : print(s)
return (s+'\n')
def plhomogenize(q,l,i) : # degree l, implicit x-power i.
# multiply with (1+zz) truncated poly up to k, again up to k-2 etc.
# powers x are implied, the complementary ones to z
# inverse plreduce on sphere: drive out all xx as (1-zz) leaving only x.
n=len(q); f=[[1,1],[0,1],[1,1]] # factor polyn 1+zz to apply
t= []; ok=True; dg=[0]*n # t target polyn
for k in range(n) :
if q[k][0] != 0 : j=l-(k+i); dg[k]=div(j,2); ok= ok and ((j%2)==0)
if not ok : print('plhomogenize bad input. exit.'); exit()
#print('dg='+str(dg)) # dbg
for k in range(n) :
u=[[0,1]]*(k+1); u[k]=[]+q[k]; h=dg[k]
# if h>0 : multiply coeff with f^h, then add to t
for j in range(h) : u=polynmul(u,f)
t= polynadd(t,u,1,1)
t,f= polynpurge(t); return t # divide out common factors
def plmore(pl,l,noisy) :
#Plm(z,x) = (-)^m x^m d^m Pl(z) where x= sqrt(1-zz)~ sin(t); xx+zz=1
# option: homogenize with products (xx+zz) when deficient
ql=pl;pwr=1;sign=-1; hm='';pi=''; rl=[[]]*(l+1);rl[0]=[]; qlist=[[l,0,[]+pl]]
for i in range(1,l+1) :
ql= polyndiff(ql); pwr += 1; ql=polynscale(ql,sign,1); qlist+=[[l,i,ql]]
pi += present('Plm['+str(l)+str(i)+']=', ql, -i, noisy)
# power i of x included, z^k must get (l-(k+i))/2 factors (1+zz).
rl[i]= plhomogenize(ql,l,i)
hm += present('Hlm['+str(l)+str(i)+']=', rl[i], l, noisy)
return hm,pi,rl,qlist
def plzero(n,noisy) : # rational polys of z
# l Pl(z)=(2l-1) z P(l-1) - (l-1)P(l-2)
# l Pl(z,x)= (2l-1) z P(l-1)(z,x) - (l-1)(zz+[1~=xx])P(l-2)
# returns hom,hpl: homog versions in (x,z); pl,pli: official Plm
# have either odd-only or even-only powers. let p be the maxpower
# multiply term p-2 with (xx+zz), (p-2k) with (xx+zz)^k.
# store in same array type, powers of x are implicit.
pl=[[]]*20; pl[0]= [[1,1]]; z= [[0,1],[1,1]]; un=[[1,1],[0,1],[1,1]]
ql=[[]]*20; ql[0]= [[1,1]]; pli='Plm[0,0]=(1)\n'; pllist=[]
hom=''; hpl=[[[1,1]]] # homogeneous variant, stringlist and poly list
for i in range(1,n) :
y= polynmul(pl[i-1],z); y= polynscale(y, 2*i-1,1)
if (i-2)>= 0 : w= polynscale(pl[i-2],i-1,1); y=polynadd(y,w,1,-1)
pl[i]= polynscale(y,1,i)
pli += present('Plm['+str(i)+'0]=',pl[i], 0, noisy)
y= polynmul(ql[i-1],z); y= polynscale(y, 2*i-1,1) # homog version
if (i-2)>= 0 : # wrong for i>2
w= polynscale(ql[i-2],i-1,1); w= polynmul(w,un); y=polynadd(y,w,1,-1)
ql[i]= polynscale(y,1,i); qq,fq = polynpurge(ql[i]) # must keep ql[i] as-is
hom += present('Hlm['+str(i)+'0]=',qq, i,noisy)
hm,pi,rl,qlist = plmore(pl[i],i, noisy)
pllist += qlist; rl[0]=qq; hpl += rl
hom+=hm;pli+=pi # expand powers x^(2n) as (1-zz)^n
return pl,hom,pli,hpl,pllist
def listpos(lst,i,j,k) :
n=len(lst); ix=-1
for h in range(n) : ix=h if([i,j,k]==lst[h]) else ix
if ix<0 : lst += [[i,j,k]]; ix=n
return lst,ix
def combine(l,m,base) :
n=len(base); cf=[[]]*(2*n); cfs=[]; icf=0
for h in range(n) :
i,j,k=tuple(base[h]); i1,j1,k1= i-1,j-1,k; i2,j2,k2= i,j,k-2; cfok=True
ok1= (i1<0)or(j1<0); ok2= (k2<0) # vanishing terms
if not ok1 : cf[icf]= [4*i*j,i1,j1,k1,h]; icf+=1
if not ok2 : cf[icf]= [k*(k-1),i2,j2,k2,h]; icf +=1
if n==1 : # check that delta deriv terms vanish, ie have a negative ix
i,j,k=tuple(base[0]); cfs=[1]
i1,j1,k1= i-1,j-1,k; i2,j2,k2= i,j,k-2
ok= ((i1<0)or(j1<0)) and (k2<0)
if not ok : print('combine bug: '+str([i1,j1,k1])+str([i2,j2,k2]))
if (n==2)and False : # must have 1 match and 2 negatives, old code
for h in range(1,4): cfok= cfok and (cf[0][h]==cf[1][h])
# in general extract all cf with same 123 triplets, return matrix
cfs=[cf[1][0],-cf[0][0]] # winning zero lincomb
print('combine cfs:'+str(cfok)+' '+str(cf)+' '+str(cfs))
if n>=2 : # make equation system, index different ijk with position p
mx=[[]]*n; lst=[]; b=[0]*n # mx has 2 nontrivial lines, common center elem?
for h in range(n) : mx[h]=[0]*n
for h in range(n) : mx[n-1][h]=1; b[n-1]=1 # dummy row
for h in range(icf) :
i,j,k= tuple(cf[h][1:4]);
lst,p= listpos(lst,i,j,k); q=cf[h][4] # q= basis ix
# make n*n matrix so that coeff sums for same p become zero
# add a line with maxvals of cols? or sum=1?
mx[p][q]= cf[h][0]
# print('matrix '+str(mx)); # showmatrix(m)
ok,xg=intsolve(mx,b); # boost to smallest ints and return as cfs
cfs= wholenb(xg)
if False : print('intsolve ok='+str(ok)+' '+str(xg)+' cfs='+str(cfs))
return cfs
def listterms(lmax,mode,noisy) :
# loop over (l,m), monomials of Ylm to polynomials
# in terms of Plm exp(im f), Plm(x,z) is obtained by x^(i+j) z^k
# when x= sin(teta); z= cos(teta) this written in c s smbols?
# note Pl(-m)=Plm: because i-j=0 iff j-i=-m.
full=(mode==1); data=[]; text=''; dm=-1
for l in range(lmax+1) :
m=l; lastm= (-l) if full else 0
if mode==2 : m=0; dm=1; lastm=l # growing order
while (dm*(lastm-m)) >= 0 :
i=l; s='l='+str(l)+' m='+str(m); base=[]; t=''+s; u=''+s; block=[l,m]
while i>=0 :
j=i-m; k=l-i-j
if (j>=0)and(k>=0) : s += ' '+str([i,j,k]); base += [[i,j,k]]
i-=1
cfs= combine(l,m,base); lb=len(base)
for h in range(lb) : # new neutral: cc becomes x, ss becomes z
z=str(cfs[h]); z= ('+'+z) if (cfs[h]>0) else z; i,j,k=tuple(base[h])
g=i+j; cc='' if (g==0) else ('x' if (g==1) else 'x^'+str(g))
ss='' if (k==0) else ('z' if (k==1) else 'z^'+str(k)); void=(cc+ss)==''
if not void : z= '+' if (cfs[h]==1) else ('-' if (cfs[h]==-1) else z)
t += ' '+z+str(base[h]); u+= ' '+z+cc+ss; block+=[base[h]+[cfs[h]]]
data += [block]
if noisy : print(u)
m += dm; text += u +'\n'
return data,text
'''
I(m,n)= integ 0 to pi of sin^m cos^n [ n even ]
- m += 2 implies I(m+2,n) = (m+1)/(m+n+2) *I(m,n)
- n += 2 implies I(m,n+2) = (n+1)/(m+n+2) *I(m,n)
- start functions : 1, sin, cos, sincos: 0,0 1,0 0,1, 1,1
0,0: pi. 1,0: 2 0,1 : 0 1,1 : 0
'''
def isico0pi(m,n) : # integral 0 to pi of sin^m(x) cos^n(x)
if (n%2 != 0) : return [0,1] # ugly trick: negative val encode pi factor
nn=0; mm= 0 if (m%2 == 0) else 1; p=-1 if (mm==0) else 2; q=1
while nn<n : p *= (nn+1); q *= (mm+nn+2); nn += 2
p,q= tuple(divide(p,q))
while mm<m : p *= (mm+1); q *= (mm+nn+2); mm += 2
return divide(p,q)
def csproduct(n,a,b,c, m,d,e,f) : # product of 2-var polynomials p(x,z)
# length n and m. a,d= lists of coeffs, b,e powers of x, c,f powers of z.
# returns list of coeffs and list of power pairs they target
# isico0pi uses these power pairs, results to be summed up then
# variant: a,d are pairs=rationals.
def ratpr(a,b) : return [a[0]*b[0],a[1]*b[1]]
def ratsum(a,b) :
if a==[] : a=[0,1]
p=a[0]*b[1]+a[1]*b[0];q=a[1]*b[1]; return divide(p,q)
def getindex(ix,p,q) : # where to store x^p y^q
i=0; n=len(ix); found=False
while (not found) and (i<n) and (ix[i]!=[]) :
found= (ix[i][0]==p)and(ix[i][1]==q); i= (i+1) if (not found) else i
if not found : ix[i]=[p,q] # new exponent combination
return i
t=str(type(a[0])); ratnum= (t.find('list')>0)
x= ([[]]*(n+m)) if ratnum else ([0]*(n+m))
ix=[[]]*(n+m); jmax=0
for i in range(n) :
for k in range(m) :
yy=b[i]+e[k]; zz=c[i]+f[k]; j=getindex(ix,yy,zz) # side effect on ix
if ratnum : x[j]= ratsum(x[j],ratpr(a[i],d[k]))
else : x[j] += a[i]*d[k]
if j>jmax : jmax=j
return x[:(jmax+1)],ix[:(jmax+1)]
def intprod(pa,ia,la, pb,ib,lb) : # integral of product of 2 plm(x,z) polyn
# pa(x,z)= sum(k,pa[k]*x^ia*z^k, pb same. nonhomog, const x power=i+j
# status ok for plzero
na=len(pa); cfa=[0]*na; pxa=[0]*na; pza=[0]*na
for i in range(na) : cfa[i]=pa[i]; pxa[i]=ia; pza[i]=i
nb=len(pb); cfb=[0]*nb; pxb=[0]*nb; pzb=[0]*nb
for i in range(nb) : cfb[i]=pb[i]; pxb[i]=ib; pzb[i]=i
cf,pw= csproduct(na,cfa,pxa,pza, nb,cfb,pxb,pzb); nt=len(cf); sum=[]
# print('intprod cf,pw='+str([cf,pw]))
for j in range(nt) :
pw[j][0] +=1 # sin factor integral measure
z= isico0pi(pw[j][0],pw[j][1]); a=z[0]; b=z[1]
z= [[0,1],[-a,b]] if (a<0) else [[a,b]]
sum= polynadd(sum,z, 1, cf[j])
# print('pa,sum='+str([pa,sum]))
return sum
def iproducts(dat,all,noisy) :
# got listterms data for Plm, eval int(Plm^2 sin(a)da)
# check, int(Pl1m1 Pl1ml2m sin(a) da) must vanish if l1 != l2.
# scalar prods either have a factor pi or not, in all terms.
mode=2 # mode2: cs are sinuses, sn cosinuses. Plm(x,z): z~cos, x~sin.
n=len(dat); csh=[0]*32; snh=[0]*32; cfh=[0]*32; bugs=0
csk=[0]*32; snk=[0]*32; cfk=[0]*32; tbl='' # tbl result table of squares
for h in range(n) :
x=dat[h]; l=x[0];m=x[1]; jh=len(x)-2; kmax= n if all else (h+1)
for i in range(jh) : # cf=coeff, cs=cosinus and sn=sinus powers of Plm
y=x[2+i]; cfh[i]=y[3]; csh[i]=y[0]+y[1]; snh[i]=y[2]
for k in range(h,kmax) : # (h,h+1) or (h,n)
x=dat[k]; lk=x[0];mk=x[1]; jk=len(x)-2; sum=[]
for i in range(jk) :
y=x[2+i]; cfk[i]=y[3]; csk[i]=y[0]+y[1]; snk[i]=y[2]
cf2,pw2= csproduct(jh,cfh,csh,snh,jk,cfk,csk,snk); nt=len(cf2) # nt terms
for j in range(nt) : pw2[j][0] += 1 # increment sin-power of integrand
for j in range(nt) : # arg order: sin-power, then cos-power
if mode==1 : # bad
z=isincos(pw2[j][1],pw2[j][0],0); sum= polynadd(sum,z, 1, cf2[j])
elif mode==2 : # better integrals, sin/cos swapped, 0 to pi. ok!
z= isico0pi(pw2[j][0],pw2[j][1]); a=z[0]; b=z[1]
z= [[0,1],[-a,b]] if (a<0) else [[a,b]]
sum= polynadd(sum,z, 1, cf2[j])
#print('l,m='+str([l,m])+' cf2='+str(cf2)
# +' pw2='+str(pw2)+' sum='+str(sum))
if noisy : print('l,m,lk,mk,nt='+str([l,m,lk,mk,nt])+' sum='+str(sum))
z= sum[0][0] if (len(sum)==1) else (abs(sum[0][0])+abs(sum[1][0]))
if (k==h) and (z==0) : print('BUG1'); bugs+=1
if k==h : tbl += 'l='+str(l)+' m='+str(m)+' '+str(sum)+'\n'
if (l != lk) and (m==mk) and (z != 0) : print('BUG2'); bugs +=1
if noisy : print('Bugs:'+str(bugs))
return tbl
def getdecimal(s,n,trig) :
d='0123456789'; i=0; m=len(s); p=[0]*n; t=''
if trig != '' : i=s.find(trig)
for k in range(n) :
while (i<m) and (d.find(s[i])<0) : i+=1
while (i<m) and (d.find(s[i])>=0) : t+=s[i]; i+=1
p[k]= int(t); t=''
return p
def stripblanks(s) :
n=len(s); t=''
for i in range(n) : t += '' if (s[i]==' ') else s[i]
return t
def plmtable(txt) : # latex-style table
u= txt.split('\n'); lu=len(u); v0='---'; v=[[]]*lu
if u[lu-1]=='' : lu-=1; u= u[:lu]
for i in range(lu) : v[i]= u[i].split(' ')
t='\\begin{array}{|l|l|l|l|l|}\n'
t+= ('\\ell & m & P_{\\ell m}(x,z) & P^h_{\\ell m}(x,z) homogen & '+
'\\int_0^{\\pi} d\\theta\\; \\sin(\\theta)\\; (P^h_{\\ell m}'+
'(\\sin(\\theta),\\cos(\\theta))^2 \\\\\n')
for i in range(lu) :
if (v[i][0] != v0) : t+= '\\hline\n'
t+=v[i][0]+' & '+v[i][1]+' & '+v[i][2]+' & '+v[i][3]+' & '+v[i][4]+'\\\\\n'
v0=v[i][0]
t+= '\\end{array}\n'
print(t); return(t)
def integlm(l,m) : # theoretical integral Plm^2, defining non-homog versions
p=2; q=2*l+1; j=l-m+1
while j<=(l+m) : p *= j; j +=1
# print('integlm l,m,p,q'+str([l,m,p,q]))
return divide(p,q)
def pleval(p,l,m) : # values at z=+-1
if m>0 : return [] # any power of x=0
zv=[1,-1]; nz=len(zv); m=len(p); yv=[[]]*nz
for i in range(nz) :
z=zv[i]; y=0; zz=1
for k in range(m) : y=add(y,mult(zz,p[k])); zz *=z
yv[i]=y
return yv
def testplm1(lmax=5) : # Plm build and table listing
noisy=False
plz,r,ri,hpl,pllist= plzero(lmax+1, noisy)
if noisy : print(r);print('') # recursion
d,txt = listterms(lmax,2, noisy) # lin-equation method
tbl=iproducts(d,True,noisy); # integrals, odd trigon powers no pi factor.
tx2=txt.split('\n'); tb2= tbl.split('\n'); n=len(tx2); m=len(tb2); dat=''
for i in range(n) :
u=tx2[i]; s=tb2[i]
if len(s) > 2 :
j=u.find('+'); v=u[j+1:]; w=u[:j]; v=stripblanks(v)
p= getdecimal(s,2,'[['); tb2[i]=(w+v+' '+str(p[0])+'/'+str(p[1]))
ri2=ri.split('\n'); m=len(ri2)
for i in range(m):
u=ri2[i]
if len(u)>0 :
a=u.find('('); v=stripblanks(u[a:]); b=v.find(')'); w='';c='?'
if v[b+1:]=='/1' : v=v[:b+1]
if v[1]=='+' : v=v[0]+v[2:]
for k in range(len(v)) :
block= (c=='x')or(c=='z');
if block : w+= c if(v[k]=='1') else (c+'^'+ v[k])
c=v[k]; w+= '' if ((c=='x')or(c=='z')or block) else c
ri2[i]=w; # print(tb2[i]+' '+ri2[i]) # remix to final string
for i in range(m) :
if len(ri2[i])>0 :
q= tb2[i].split(' '); # print('['+q[0]+'|'+q[1]+']')
x= getdecimal(q[0],1,''); y= getdecimal(q[1],1,'')
s=str(x[0])+' '+str(y[0])+' '+ri2[i]+' '+q[2]+' '+q[3]
dat= s if (dat=='') else (dat+'\n'+s)
plmtable(dat)
testplm1()
