      program freud3
c                               2d        2d-2
c                    rho   -c0 x     -c1 x        ...
c recurrence  for |x|     e                            on (-inf.,+inf.)
c
c                                    d       d-1
c                   (rho-1)/2   -c0 x   -c1 x    ...
c      =           x           e                       on (0,inf.)
c
      integer*2 s(162374)
      real*8 f(2500),jac(2500,-9:0),x(-20:2500)
      real*8 coeff(0:10),cc,rho,err,term,fm
      integer degmax,deg,p,pp,p1,p2,ds(10)
      logical fex
      character*80 fname0,fname
      data ds/0,12,57,221,772,2492,7557,21807,60474,162374/
      do 100 i=-20,2500
  100 x(i)=0
      print *,' max degree ?  '
      read *,degmax
      if(degmax.gt.10)then
        print *, '  degree must be  <= 10  ... '
        stop
      endif
      is=ds(degmax)
      fname0='freud0xx.dat'
      fname='    '
      do i1=2,0,-1
       fname0(7:7)=char(48+i1)
c  '(assuming ascii)
       do i2=9,0,-1
        if(10*i1+i2.lt.degmax)goto 15
        fname0(8:8)=char(48+i2)
        inquire(file=fname0, exist=fex)
        if(fex)fname=fname0
        if(fex)print *,fname
       enddo
      enddo
   15 continue
      open(unit=1,file=fname)
      is=ds(degmax)
      read(1,'(20i4)')(s(i),i=1,is)
      print '(9i4,'' .. '',9i4)',(s(i),i=1,9),(s(i),i=is-8,is)
    6 print *,' degree (end if <0), number of recurrence coef.?  '
      read *,deg,nf
      if((deg.lt.0).or.(deg.gt.degmax))stop
      print *,' give rho and polynomial coef. (starting with largest degree) '
      read *,rho,(coeff(i),i=0,deg-1)
      cc=2*coeff(0)
      do 110 m=2,deg
       fm=m-1
  110 cc=cc*(4*m-2)/fm
      initx=initx+1
       do 12 n=1,nf+deg-1
        if(x(n).eq.0)x(n)=(n/cc)**(1.0d0/deg)
   12  continue
      print 11, rho,cc,( coeff(deg-i),i,i=deg,1,-1)
   11 format(' ',f5.2,f9.1,' ;',7(f10.4,' x',i3,1x))
      err=1.0d25
      do 1 iter=1,20
       p=1
       do 20 n=1,nf
        f(n)=coeff(deg-1)*x(n)*2
        jac(n,0)=f(n)
        do 20 k=1-deg,-1
   20  jac(n,k)=0
       do 2 m=2,deg
c  F
        do 21 i=0,s(p+1)-1
         pp=p+2+i*m
         p2=pp+m-1
         do 210 n=1,nf
          term=2*m*s(pp)*coeff(deg-m)*x(n)
          do 215 j=1,m-1
  215     term=term*x(n+s(pp+j))
         jac(n,0)=jac(n,0)+term
  210    f(n)=f(n)+term
   21   continue
c  dF/dx
        do 22 k=1-m,0
         p1=p2+1
         do 22 i=0,s(p1)-1
c monomes de dFn / dx(n+k)
          pp=p1+1+i*(m-1)
          p2=pp+m-2
          do 220 n=1,nf
           term=2*m*s(pp)*coeff(deg-m)*x(n)*x(n+k)
           do 2200 j=pp+1,p2
 2200      term=term*x(n+s(j))
  220     jac(n,k)=jac(n,k)+term
  22     continue
  2     p=p2+1
       print 291,(sqrt(x(i)),i=1,10)
       print 291,(f(i),i=1,10)
c      do 29 k=1-deg,0
c  29  print 291,(jac(i,k),i=1,10)
       print 291,(jac(i,0),i=1,10)
       if(iter.gt.1)print *,iter,' iter. err ',ierr,':',err
  291  format(10f8.4)
       if(err.lt.1.0d-13) goto 5
Cholesky
       do 3 n=1,nf
        mn=max(1-n,1-deg)
        f(n)=f(n)-n-rho*mod(n,2)
        do 3 k=mn,0
         do 30 l=mn,k-1
   30    jac(n,k)=jac(n,k)-jac(n,l)*jac(n+k,l-k)
         if(k.lt.0)then
           jac(n,k)=jac(n,k)/jac(n+k,0)
           f(n)=f(n)-jac(n,k)*f(n+k)
          elseif(jac(n,0).le.1.0d-15)then
                print *,' err Cholesky n= ',n,jac(n,0)
                print *,' sorry, try another intermediate polynomial'
                stop
           else
                jac(n,0)=sqrt(jac(n,0))
                f(n)=f(n)/jac(n,0)
         endif
 3     continue
Correction Newton
       err=0
       do 4 n=nf,1,-1
        do 40 i=1,min(deg-1,nf-n)
   40   f(n)=f(n)-f(n+i)*jac(n+i,-i)
        f(n)=f(n)/jac(n,0)
        if( abs(f(n)).gt.err )then
         ierr=n
         err=abs(f(n))
        endif
  4    x(n)=x(n)*exp(-f(n))
1     continue
 5    print *,' ----------- result ----------- '
      print *,deg
      print '(5f16.10)',rho,(coeff(i),i=0,deg-1)
      print *,nf
      do 55 i=0,nf
      print '(i5,f20.14)',i,x(i)
      print '(i5,f20.14)',i,sqrt(x(i))
      if(mod(i,2).ne.0)print '(i5,2f20.14)',i/2,x(i-1)+x(i),
     &  x(i-2)*x(i-1)
   55 continue
      goto 6
      end
