C     Last change:  LS   25 Jan 100    6:06 pm
      program boulanger

      implicit none
c 16 est le nombre de decimales binaires, 10000 est le nb d'iterations
      REAL A(16,10000)
      real x1(16),x2(16),xinit(10016)
c taille xinit nb de decimales binaires  plus nb iterations
      real y1,y2,fin
      real beta,xsol,xdepart,ailleurs,temp,tempp,temppp,xdec
      INTEGER n,i,k,itime,iseed

c
      open (2,FILE='resuXY.xls')

c n est le nombre d iterations
      n=10000
c beta est le coefficielnt de contraction
      beta=0.7
c y1 est l ordonnee initiale
      y1=0.2
      write (*,*) y1

c ecriture aleatoire du nombre x1 initial en binaire

      do i=1,16
        call RANDOM_NUMBER(x1(i))
      	if (x1(i).le.0.5) then
      	x1(i)=0
        else
      	x1(i)=1
      	end if
      xinit(i)=x1(i)
      end do

c debut de la boucle sur les iterations
      do k=1,n

c boucle sur les decimales de x1
      do i=1,16

c decalage du vecteur d un cran vers la gauche (partie fractionnaire de 2x1)
      x2(i)=x1(i+1)
      end do
c determination aleatoire du 0 ou 1 rajoute a la fin du vecteur apres le decalage
      call RANDOM_NUMBER(fin)

      	if (fin.le.0.5) then
      	fin=0
        else
      	fin=1
        end if

      x2(16)=fin


c ecriture des termes rajoutes a xinit (ecriture de sa vraie valeur)
      xinit(16+k)=x2(16)


c xdec est l'ecriture en decimale de x au cours des iterations
      xdec=0
      do i=1,16
      x1(i)=x2(i)

c ecriture du vecteur x1dans la matrice A a chaque iteration
      A(i,k)=x1(i)
c passage en decimal
      tempp=2**(i)
      tempp=1/tempp
      xdec=xdec+x1(i)*tempp

      end do


c calcul des y avec contraction beta
      if (A(1,k)==0.) then
      y2=beta*y1/2
      else
      y2=beta*(y1+1)/2
      end if
      y1=y2

c on ecrit les xdec au cours des iterations mais pas le xinit dont on ne dispose pas encore
      write (2,*) xdec,y1
      end do
c fin de la boucle sur les iterations

c ecriture du vecteur initial dans le xinit.xls
c la boucle va de 1 a nb iterations plus nb de decimales binaires

c le calcul s'arrete et ecrit des etoiles a partir de 32 calculs a cause de la puissance de 2 donc on approxime a 30
      do i=1,30
      temppp=2**(i)
      temppp=1/temppp
      xdepart=xdepart+xinit(i)*temppp
      end do
      write (*,*) xdepart

c passage de binaire en decimal pour le vecteur de sortie x1 (image de xinit au bout des n iterations)

      do i=1,16
      temp=2**(i)
      temp=1/temp
      xsol=xsol+x1(i)*temp
      end do
      write (*,*) xsol
      write (*,*) y1


      close (2)

      end program