program polynomial
      implicit none
      integer*4 n0
      parameter (n0=1000)
      integer*4 i,k,m,kmax,lmax,
     1n ,ni( n0),nl( n0),ki( n0),mi( n0),ml( n0),kl( n0),
     1np,nip(n0),nlp(n0),kip(n0),mip(n0),mlp(n0),klp(n0),
     1nc,nic(n0),nlc(n0),kic(n0),mic(n0),mlc(n0),klc(n0)
      real*8 p(n0),pi,pp(n0),c(n0)
      real*8 ti,fi,ri,tl,fl,rl,xi,yi,zi,xl,yl,zl,dl,
     1sr,si,slr,sli,sir,sii

      call Initfac
      pi=4.0d0*atan(1.0d0)
c     sum(k=0..7,m=-k,k)4pi(2k+1) ri^ni/rl(^nl+1)Yki,miYklml
c     Ykm - spherical functions
      kmax=7
      lmax=7

      i=0
      do 1 k=0,kmax
      do 1 m=-k,k
      i=i+1
      if(i.gt.n0)call report('too many terms')
      p(i)=4.0d0*pi/dble(2*k+1)
      if(mod(m,2).ne.0)p(i)=-p(i)
      ni(i)=k
      nl(i)=k+1
      ki(i)=k
      mi(i)=m
      kl(i)=k
1     ml(i)=-m
      n=i
      write(6,*)n,' terms'

      call pq(      n0,n ,p ,ni ,nl, ki, mi, kl, ml,
     1                 n ,p ,ni ,nl, ki, mi, kl, ml,
     1                 np,pp,nip,nlp,kip,mip,klp,mlp,lmax,kmax)
      write(6,*)'squared -',np,' terms'
      call pq(      n0,np,pp,nip,nlp,kip,mip,klp,mlp,
     1                 np,pp,nip,nlp,kip,mip,klp,mlp,
     1                 nc,c ,nic,nlc,kic,mic,klc,mlc,lmax,kmax)
      write(6,*)'fourth power -',nc,' terms'
      do 5 i=1,nc
 5    write(6,609)i,c(i),nic(i),nlc(i),kic(i),mic(i),klc(i),mlc(i)
 609  format(i4,f12.6,2i3,2x,2i3,2x,2i3)
      write(6,*)
      write(6,*)'1/r  1/r (expansion):'
      do m=1,5
      if(m.eq.1)ti=45.0d0*pi/180.0d0
      if(m.eq.1)fi=45.0d0*pi/180.0d0
      if(m.eq.1)ri=1.0d0
      if(m.eq.1)tl=145.0d0*pi/180.0d0
      if(m.eq.1)fl=-45.0d0*pi/180.0d0
      if(m.eq.1)rl=4.0d0
      if(m.eq.2)ti=-45.0d0*pi/180.0d0
      if(m.eq.2)fi=45.0d0*pi/180.0d0
      if(m.eq.2)ri=2.0d0
      if(m.eq.2)tl=145.0d0*pi/180.0d0
      if(m.eq.2)fl=-35.0d0*pi/180.0d0
      if(m.eq.2)rl=4.0d0
      if(m.eq.3)ti=45.0d0*pi/180.0d0
      if(m.eq.3)fi=245.0d0*pi/180.0d0
      if(m.eq.3)ri=1.50d0
      if(m.eq.3)tl=0.0d0*pi/180.0d0
      if(m.eq.3)fl=-15.0d0*pi/180.0d0
      if(m.eq.3)rl=4.0d0
      if(m.eq.4)ti=5.0d0*pi/180.0d0
      if(m.eq.4)fi=5.0d0*pi/180.0d0
      if(m.eq.4)ri=1.0d0
      if(m.eq.4)tl=-45.0d0*pi/180.0d0
      if(m.eq.4)fl=45.0d0*pi/180.0d0
      if(m.eq.4)rl=4.0d0
      if(m.eq.5)ti=40.0d0*pi/180.0d0
      if(m.eq.5)fi=45.0d0*pi/180.0d0
      if(m.eq.5)ri=3.0d0
      if(m.eq.5)tl=25.0d0*pi/180.0d0
      if(m.eq.5)fl=-45.0d0*pi/180.0d0
      if(m.eq.5)rl=4.0d0
      xi=ri*sin(ti)*cos(fi)
      yi=ri*sin(ti)*sin(fi)
      zi=ri*cos(ti)
      xl=rl*sin(tl)*cos(fl)
      yl=rl*sin(tl)*sin(fl)
      zl=rl*cos(tl)
      dl=dsqrt((xi-xl)**2+(yi-yl)**2+(zi-zl)**2)
c     reproduce 1/r:
      sr=0.0d0
      si=0.0d0
      do 2 i=1,n
      call sph(kl(i),ml(i),tl,fl,slr,sli)
      call sph(ki(i),mi(i),ti,fi,sir,sii)
      sr=sr+ri**ni(i)/rl**nl(i)*p(i)*(sir*slr-sii*sli)
 2    si=si+ri**ni(i)/rl**nl(i)*p(i)*(sii*slr+sir*sli)
      write(6,600)1.0d0/dl,sr,si
600   format(3f10.5)
      write(6,*)
      write(6,*)'1/r^2  1/r^2 (expansion):'
c     reproduce 1/r^2:
      sr=0.0d0
      si=0.0d0
      do 3 i=1,np
      call sph(klp(i),mlp(i),tl,fl,slr,sli)
      call sph(kip(i),mip(i),ti,fi,sir,sii)
      sr=sr+ri**nip(i)/rl**nlp(i)*pp(i)*(sir*slr-sii*sli)
 3    si=si+ri**nip(i)/rl**nlp(i)*pp(i)*(sii*slr+sir*sli)
      write(6,600)1.0d0/dl**2,sr,si
      write(6,*)
      write(6,*)'1/r^4  1/r^4 (expansion):'
c     reproduce 1/r^4:
      sr=0.0d0
      si=0.0d0
      do 4 i=1,nc
      call sph(klc(i),mlc(i),tl,fl,slr,sli)
      call sph(kic(i),mic(i),ti,fi,sir,sii)
      sr=sr+ri**nic(i)/rl**nlc(i)*c(i)*(sir*slr-sii*sli)
 4    si=si+ri**nic(i)/rl**nlc(i)*c(i)*(sii*slr+sir*sli)
      write(6,600)1.0d0/dl**4,sr,si
      enddo
      end

      subroutine sph(l,m,t,f,re,im)
c     spherical function Yl,m(t,f) ... r, im, real and imaginaty parts
      implicit none
      integer*4 l,m
      real*8 t,f,re,im,n,norm,pol,p
      norm=n(l,m)
      pol=p(m,l,t)
      re=norm*pol*cos(dble(m)*f)
      im=norm*pol*sin(dble(m)*f)
      return
      end

      function p(m,l,t)
c     associated Legendre polynomial              
      implicit none
      integer*4 nfac0
      parameter (nfac0=100)
      integer*4 l,m,ll,i,mp
      real*8 t,z,p,sgn,pmm,pmmp,po
      real*8 dm,d0,df
      common/dpar/dm,d0,df(nfac0)
      real*8 f0,ff
      common/npar/f0,ff(nfac0)
      real*8,allocatable::pp(:)
      if(l.eq.0.and.m.eq.0)then
       p=1.0d0
       return
      endif
      z=cos(t)
c     pp is p for positive m:
      mp=abs(m)
c     Pm,m
      pmm=df(2*mp-1)*(1.0d0-z**2)**(dble(mp)/2.0d0)
      if(mp.eq.l)then
       po=pmm
      else
c      Pm,m+1
       pmmp=z*dble(2*mp+1)*pmm
       if(mp+1.eq.l)then
        po=pmmp
       else
        allocate(pp(l-mp+1))
        pp(1)=pmm
        pp(2)=pmmp
        do 1 ll=mp+2,l
        i=ll-mp+1
1       pp(i)=
     1  (z*dble(2*ll-1)*pp(i-1)-dble(ll+mp-1)*pp(i-2))/dble(ll-mp)
        po=pp(l-mp+1)
       endif
      endif
      if(m.gt.0)then
       p=po
      else
       if(mod(mp,2).eq.0)then
        sgn= 1.0d0
       else
        sgn=-1.0d0
       endif
       p=sgn*ff(l-mp)/ff(l+mp)*po
      endif
      return
      end

      function n(l,m)
      implicit none
      real*8 n,fp
      integer*4 l,m
      integer*4 nfac0
      parameter (nfac0=100)
      real*8 f0,factorial
      common/npar/f0,factorial(nfac0)
      fp=4.0d0*4.0d0*datan(1.0d0)
c     n=dsqrt(dble(2*l+1)*factorial(l-abs(m))/fp/factorial(l+abs(m)))
      n=dsqrt(dble(2*l+1)*factorial(l-m)/fp/factorial(l+m))
      return
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      subroutine pq(n0,np,p ,nip,nlp,kip,mip,klp,mlp,
     1                 nq,q ,niq,nlq,kiq,miq,klq,mlq,
     1                  k,s ,nis,nls,kis,mis,kls,mls,lmax,kmax)
c     s = p x q ;special polynoms:
c     p = sum (j) ri^nij/rL^nLj Ykjmj(i)YkLjmL(j(L)
      
      integer*4 n0,lmax,i,j,lp,mp,lq,mq,kk,kmax,
     1np,nip(*),nlp(*),kip(*),mip(*),klp(*),mlp(*),
     1nq,niq(*),nlq(*),kiq(*),miq(*),klq(*),mlq(*),
     1 k,nis(*),nls(*),kis(*),mis(*),kls(*),mls(*)
      real*8 p(*),q(*),s(*),c,yyy,cpq,tol
c
c     generate all possible Yi x YL terms and digest
      k=0
      do 7 i=1,np
      do 7 j=1,nq
c     Y(i)p x Y(i)q = sum YYY Ylm(i), etc
      cpq=p(i)*q(j)
      if(cpq.ne.0.0d0.and.nip(i)+niq(j).le.lmax)then
       do 71 lp=0,kmax
       do 71 mp=-lp,lp
       do 71 lq=0,kmax
       do 71 mq=-lq,lq
       c=yyy(lp,mp,kip(i),mip(i),kiq(j),miq(j))
     1 * yyy(lq,mq,klp(i),mlp(i),klq(j),mlq(j))
       if(c.ne.0.0d0)then
c       try to find this term:
        do 8 kk=1,k
        if(nis(kk).eq.nip(i)+niq(j).and.
     1     nls(kk).eq.nlp(i)+nlq(j).and.
     1     kis(kk).eq.lp.and.
     1     mis(kk).eq.mp.and.
     1     kls(kk).eq.lq.and.
     1     mls(kk).eq.mq)then
c        this term exists, add it and go to next:
         s(kk)=s(kk)+cpq*c
         goto 71
        endif
8       continue
c       this term was not found, add it:
        k=k+1
        if(k.gt.n0)call report('too many terms - Y stage')
        nis(k)=nip(i)+niq(j)
        nls(k)=nlp(i)+nlq(j)
        kis(k)=lp
        mis(k)=mp
        kls(k)=lq
        mls(k)=mq
        s(kk)=cpq*c
       endif
71     continue
      endif
7     continue
      write(6,*)k,' final unique terms'
c     delete small ones
      tol=1.0d-5
99    do 9 i=1,k
      if(dabs(s(i)).lt.tol)then
       do 91 kk=i,k-1
       s(kk)=s(kk+1)
       nis(kk)=nis(kk+1)
       nls(kk)=nls(kk+1)
       kis(kk)=kis(kk+1)
       mis(kk)=mis(kk+1)
       kls(kk)=kls(kk+1)
91     mls(kk)=mls(kk+1)
       k=k-1
       goto 99
      endif
9     continue
      write(6,*)k,' large terms left'

      return
      end
c     ============================================================
      function yyy(l1,m1,l2,m2,l3,m3)
c     <l1 m1 |Y_l2,m2|l3 m3>
      implicit none
      real*8 yyy,pi,la,ma,lb,mb,lc,mc,rthrj,z0,one,two,four
      integer*4 l1,m1,l2,m2,l3,m3
      data z0,one,two,pi,four/0.0d0,1.0d0,2.0d0,
     13.14159265358979d0,4.0d0/
      la=DBLE(l1)
      ma=DBLE(m1)
      lb=DBLE(l2)
      mb=DBLE(m2)
      lc=DBLE(l3)
      mc=DBLE(m3)
c
c    int Y(*)l1m1 Yl2m2 Yl3m3 = dsqrt((2l1+1)(2l2+1)(2l3+1)/(4pi))
c
c    x    (l1 l2 l3) (l1 l2 l3) x (-1)^m1
c         (m1 m2 m3) (0  0  0 )
c
      yyy=dsqrt((two*la+one)*(two*lb+one)*(two*lc+one)/(four*pi))
     1*rthrj(la,lb,lc,ma,mb,mc)*rthrj(la,lb,lc,z0,z0,z0)
      if(mod(m1,2).ne.0)yyy=-yyy
      return
      end
      
      Function rthrj(x1,x2,x3,y1,y2,y3)
c     {real three - J symbols}
      implicit none
      integer*4 nfac0
      parameter (nfac0=100)
      real*8 f0,factorial
      common/npar/f0,factorial(nfac0)
      real*8  rthrj,x1,x2,x3,y1,y2,y3,ar,u,sum,msign
      integer*4 j1d,j2d,j3d,m1d,m2d,m3d,s,ni,i1,i2,i3,i4,i5,i6
      j1d=NINT(x1*2.0d0)
      j2d=NINT(x2*2.0d0)    
      j3d=NINT(x3*2.0d0)   
      m1d=NINT(y1*2.0d0)
      m2d=NINT(y2*2.0d0)    
      m3d=NINT(y3*2.0d0)   
      if(abs(m1d).gt.j1d.or.abs(m2d).gt.j2d.or.abs(m3d).gt.j3d.or.
     1m1d+m2d .ne.-m3d.or.j1d+j2d.lt.j3d.or.j2d+j3d.lt.j1d.or.
     2j3d+j1d.lt.j2d   .or.mod(j1d-m1d,2).ne.0 .or.
     3mod(j2d-m2d,2)   .ne.0.or.mod(j3d-m3d,2).ne. 0.or.
     4mod(j1d+j2d+j3d,2).ne.0) then 
       ar=0.0d0
      else 
       if (mod((j1d-j2d-m3d)/2,2).eq.0) then
        msign=1.0d0 
       else
        msign=-1.0d0   
       endif
       s=NINT(x1+x2+x3)   
       ar=factorial(s-j3d)*factorial(s-j2d)*factorial(s-j1d)
     1 /factorial(s+1)
       ar=ar*factorial((j1d+m1d) / 2)*factorial((j1d-m1d) / 2)   
       ar=ar*factorial((j2d+m2d) / 2)*factorial((j2d-m2d) / 2)   
       ar=ar*factorial((j3d+m3d) / 2)*factorial((j3d-m3d) / 2)   
       ar=msign*dsqrt(ar)   
       sum=0.0d0   
       do 1 ni=0 , s 
       i1=ni                       
       i2=(j1d+j2d-j3d) / 2-ni   
       i3=(j1d-m1d) / 2-ni       
       i4=(j2d+m2d) / 2-ni   
       i5=(j3d-j2d+m1d) / 2+ni   
       i6=(j3d-j1d-m2d) / 2+ni   
       if( mod(ni, 2) .ne. 0) then
        msign=-1.0d0
       else
        msign=1.0d0   
       endif
       if(((((i2.ge.0).and.(i3.ge.0)).and.(i4.ge.0)).and.(i5.ge.0))
     1 .and.(i6.ge.0)) then
        u=factorial(i1)*factorial(i2)   
        u=u*factorial(i3)*factorial(i4)   
        u=u*factorial(i5)*factorial(i6)   
        sum=sum+msign/u   
       endif   
1      continue
       ar=ar*sum   
       if(abs(ar).lt.1.0d-38) ar=0.0d0
      endif
      rthrj=ar
      return
      end

      subroutine Initfac
      implicit none
      integer*4 nfac0,i
      parameter (nfac0=100)
      real*8 f0,factorial,fac
      common/npar/f0,factorial(nfac0)
      real*8 dm,d0,df,dfac
      common/dpar/dm,d0,df(nfac0)
      do  1 i=0 , nfac0
1     factorial(i)=fac(i)
      do  2 i=-1, nfac0
2     df(i)=dfac(i)
      return
      end

      Function DFac(n)
c     double factorial = n!!
      implicit none
      integer*4 n,ii
      real*8 DFac,bb
      if(n.eq.-1.or.n.eq.0)then
       bb=1.0d0
      else
       bb=dble(n)
       if(mod(n,2).ne.0)then
        do 1 ii=n-2,1,-2
1       bb=bb*dble(ii)
       else
        do 2 ii=n-2,2,-2
2       bb=bb*dble(ii)
       endif
      endif
      dfac=bb
      return
      end

      Function Fac(n)
      implicit none
      integer*4 n,ii
      real*8 Fac,bb
      bb=1.0d0
      do 1 ii=2 , n
1     bb=bb*dble(ii)
      fac=bb
      return
      end