c*********************************
c     RLS algorithm
c*********************************
      parameter (ndata = 100)
      complex*16 rxxi(2,2)
      complex*16 xin(2,ndata),y(ndata),ww(2,0:ndata),r(ndata)
      complex*16 dwave, uwave, dsig, usig, zj, er
      complex*16 rtmp(2)
      real*8 pow1, pow2, ang1, ang2, pi, derad
      real*8 ee(0:ndata), alpha, mean_error, gamma
      real noise, st, thmr, thmi
      integer i, idst/584287/, iust/350107/, inst/455782/
      common /rdata/pow1,ang1,pow2,ang2,noise
      pi = acos(-1.0d0)
      derad = pi / 180.0d0
      zj = (0.0d0, 1.0d0)

c  signal parameters
      pow1 = 1.0d0
      pow2 = 1.0d0
      noise = 0.01d0
      st = sqrt(noise/2.0)
      ang1 = 30.0d0
      ang2 = -60.0d0

      open( 9,file='error-rls.dat')
      open(10,file='weight-rls.dat')
      write(10,*) ndata+1

c initial state
      rxxi(1,1) = 1.0d0
      rxxi(2,2) = 1.0d0
      rxxi(1,2) = 0.0d0
      rxxi(2,1) = 0.0d0
      ww(1,0) = (1.0d0, 0.0d0)
      ww(2,0) = (0.0d0, 0.0d0)
c
      alpha = 0.9d0
c
c  weight update loop
      do i = 1, ndata
         dsig = sqrt(pow1)*dwave(idst)
         usig = sqrt(pow2)*uwave(iust)

c  input at 1st element
         call normal(inst,0.0,st,thmr)
         call normal(inst,0.0,st,thmi)
         xin(1,i) = dsig + usig + dble(thmr)+zj*dble(thmi)

c  input at 2nd element
         call normal(inst,0.0,st,thmr)
         call normal(inst,0.0,st,thmi)
         xin(2,i) = dsig*cdexp(-zj*pi*sin(ang1*derad))
     &        + usig*cdexp(-zj*pi*sin(ang2*derad))
     &        + dble(thmr) + zj*dble(thmi)

c  array output
         y(i) = dconjg(ww(1,i-1))*xin(1,i)
     &        +dconjg(ww(2,i-1))*xin(2,i)

c  reference signal
         r(i) = dsig

c  mean squared error
         ee(i-1) = mean_error(ww(1,i-1))

c  error signal
         er = r(i) - y(i)

c  weight update
         rtmp(1)=rxxi(1,1)*xin(1,i)+rxxi(1,2)*xin(2,i)
         rtmp(2)=rxxi(2,1)*xin(1,i)+rxxi(2,2)*xin(2,i)

         gamma = dconjg(xin(1,i))*rtmp(1)
     &        + dconjg(xin(2,i))*rtmp(2)
     &        + alpha
         gamma = 1.0d0 / gamma

         ww(1,i) = ww(1,i-1)
     &        + gamma*dconjg(er)*rtmp(1)
         ww(2,i) = ww(2,i-1)
     &        + gamma*dconjg(er)*rtmp(2)

c  Rxx^{-1} update
         rxxi(1,1) = rxxi(1,1)/alpha
     &        - rtmp(1)*dconjg(rtmp(1))*gamma/alpha
         rxxi(1,2) = rxxi(1,2)/alpha
     &        - rtmp(1)*dconjg(rtmp(2))*gamma/alpha
         rxxi(2,1) = rxxi(2,1)/alpha
     &        - rtmp(2)*dconjg(rtmp(1))*gamma/alpha
         rxxi(2,2) = rxxi(2,2)/alpha
     &        - rtmp(2)*dconjg(rtmp(2))*gamma/alpha

      end do

      ee(ndata) = mean_error(ww(1,ndata))

c  results
      do i = 0, ndata
         write( 9,*) i, float(ee(i))
         write(10,*) ww(1,i),ww(2,i)
      end do

      close(9)
      close(10)

      stop
      end

c  desied signal
      complex*16 function dwave(idst)
      complex*16 zj/(0.0d0,1.0d0)/
      real*8 pi,phase
      real pp
      integer idst
      pi = acos(-1.0d0)
      call random(idst,pp)
      phase = 2.0d0*pi*pp
      dwave = cdexp(zj*phase)
      return
      end

c  undesired signal
      complex*16 function uwave(iust)
      complex*16 zj/(0.0d0,1.0d0)/
      real*8 pi,phase
      real pp
      integer iust
      pi = acos(-1.0d0)
      call random(iust,pp)
      phase = 2.0d0*pi*pp
      uwave = cdexp(zj*phase)
      return
      end

c  mean squared error
      real*8 function mean_error(weight)
      complex*16 weight(*)
      complex*16 rxx(2,2), rv(2), zj
      real*8 pow1, pow2, ang1, ang2, pi, derad
      real noise
      common /rdata/pow1,ang1,pow2,ang2,noise
      pi = acos(-1.0d0)
      derad = pi/180.0d0
      zj = (0.0d0, 1.0d0)
      rxx(1,1) = pow1+pow2+noise
      rxx(2,2) = rxx(1,1)
      rxx(1,2) = pow1*cdexp(zj*pi*sin(ang1*derad))
     &     + pow2*cdexp(zj*pi*sin(ang2*derad))
      rxx(2,1) = dconjg(rxx(1,2))
      rv(1) = pow1
      rv(2) = pow1*cdexp(-zj*pi*sin(ang1*derad))
      mean_error = pow1
     &  - 2.0d0*real(dconjg(weight(1))*rv(1)+dconjg(weight(2))*rv(2))
     &  + dconjg(weight(1))*(rxx(1,1)*weight(1)+rxx(1,2)*weight(2))
     &  + dconjg(weight(2))*(rxx(2,1)*weight(1)+rxx(2,2)*weight(2))
      return
      end

c  normal random
      subroutine normal(x,av,st,rn)
      integer x
      real av,st,rn,r1,r2,w1,w2
      call random(x,r1)
      call random(x,r2)
      w1=log(r1)
      w1=sqrt(-2.*w1)
      w2=sin(6.2831852*r2)
      rn=av+st*w1*w2
      return
      end

c  uniform random
      subroutine random(x,ru)
      integer x
      real ru
      x=x*65539
      if (x.lt.0) x=(x+2147483647)+1
      ru=x*0.4656613e-9
      return
      end