//**Turbo碼譯碼，Log-MAP算法**//
//**假設是AWGN信道上的BPSK傳輸,信噪比Eb/N0=1dB,Lc=4Eb/N0**//
#include<stdio.h>
#include<math.h>
#include <stdlib.h>
#define SIZE 9
#define L 3
#define Lc 40
#define Es 1
#define Pi 3.14159265358979   
#define Epsilon exp(1)        


double NB_RND()    //隨機數列發生器[0，1）//
{double nA, nB, nC, nD;
static double nClock=0;
static double nSeed=0;

if(nClock==0)
     nClock=nSeed=rand();
if(nClock==rand())
     nSeed=nSeed*Pi/Epsilon;
else
     nClock=nSeed=rand();
if(nSeed==0)                 
     nSeed=Pi+Epsilon;    
nA=(nSeed*100)/(Pi-int(Pi)+Epsilon-int(Epsilon));
nB=(nA -int(nA))/nSeed;
if(nB==0)
nC=nB+1/Pi+Epsilon;
else
nC=nB+1/nB;
nD=exp(nB+log(nC))-int(exp(nB+log(nC)));
return(nD);
}
double mk(double a,double s,double p,int c,int u)  //Mk(e)的計算//
{double mk;
 if(u==1)
     mk=a-log(1+exp(a))+1/2*s+1/2*p*(2*c-1);
 else
	 mk=-log(1+exp(a))-1/2*s+1/2*p*(2*c-1);
 return(mk);
}
double abk(double t1,double k1,double t2,double k2)  //Ak(e)和Bk(e)的計算//
{double s1,s2,s;
 s1=exp(t1+k1);s2=exp(t2+k2);
 s=log(s1+s2);
 return(s);
}
//分量譯碼器//
void DEC(double a[SIZE+1],double ys[SIZE+1],double yp[SIZE+1],double e[SIZE])
{double me1[SIZE+1],me2[SIZE+1],me3[SIZE+1],me4[SIZE+1],me5[SIZE+1],me6[SIZE+1],me7[SIZE+1],me8[SIZE+1];
 double a0[SIZE],a1[SIZE],a2[SIZE],a3[SIZE],b0[SIZE+1],b1[SIZE+1],b2[SIZE+1],b3[SIZE+1];
 int i,u,c;
 for(i=1;i<=SIZE;i++)
 {c=0;u=0;
  me1[i]=mk(a[i],ys[i],yp[i],c,u);
  c=1;u=1;
  me2[i]=mk(a[i],ys[i],yp[i],c,u);
  c=0;u=1;
  me3[i]=mk(a[i],ys[i],yp[i],c,u);
  c=1;u=0;
  me4[i]=mk(a[i],ys[i],yp[i],c,u);
  c=0;u=0;
  me5[i]=mk(a[i],ys[i],yp[i],c,u);
  c=1;u=1;
  me6[i]=mk(a[i],ys[i],yp[i],c,u);
  c=0;u=1;
  me7[i]=mk(a[i],ys[i],yp[i],c,u);
  c=1;u=0;
  me8[i]=mk(a[i],ys[i],yp[i],c,u);
 }
 a0[0]=1;a1[0]=0;a2[0]=0;a3[0]=0;
 b0[SIZE]=0;b1[SIZE]=0;b2[SIZE]=1;b3[SIZE]=0;
 for(i=1;i<SIZE;i++)
 {a0[i]=abk(a0[i-1],me1[i],a2[i-1],me6[i]);
  a1[i]=abk(a0[i-1],me2[i],a2[i-1],me5[i]);
  a2[i]=abk(a1[i-1],me4[i],a3[i-1],me7[i]);
  a3[i]=abk(a1[i-1],me3[i],a3[i-1],me8[i]);
 }
 for(i=SIZE-1;i>=1;i--)
 {b0[i]=abk(b0[i+1],me1[i+1],b1[i+1],me2[i+1]);
  b1[i]=abk(b3[i+1],me3[i+1],b2[i+1],me4[i+1]);
  b2[i]=abk(b1[i+1],me5[i+1],b0[i+1],me6[i+1]);
  b3[i]=abk(b2[i+1],me7[i+1],b3[i+1],me8[i+1]);
 }
 for(i=1;i<SIZE;i++)
    e[i]=log(exp(a0[i-1]+1/2*yp[i]+b1[i])+exp(a1[i-1]-1/2*yp[i]+b2[i])+exp(a2[i-1]+1/2*yp[i]+b0[i])+exp(a3[i-1]-1/2*yp[i]+b3[i]))-log(exp(a0[i-1]-1/2*yp[i]+b0[i])+exp(a1[i-1]+1/2*yp[i]+b3[i])+exp(a2[i-1]-1/2*yp[i]+b1[i])+exp(a3[i-1]+1/2*yp[i]+b2[i]));
}
void interlace(double a[SIZE+1],double b[SIZE+1])  //交織器//
{double interlace[L][L];
 int i,j,k; 
 k=1;
 for(i=0;i<L;i++)
    for(j=0;j<L;j++)
	{interlace[i][j]=a[k];
     k++;
	}
 k=1;
for(j=0;j<L;j++)
   for(i=0;i<L;i++)
      {b[k]= interlace[i][j]; 
       k++;
      }
}
void uninterlace(double a[SIZE+1],double b[SIZE+1])  //解交織器//
{double interlace[L][L];
 int i,j,k; 
 k=1;
 for(j=0;j<L;j++)
    for(i=0;i<L;i++)
	{interlace[i][j]=a[k];
     k++;
	}
 k=1;
for(i=0;i<L;i++)
   for(j=0;j<L;j++)
      {b[k]= interlace[i][j]; 
       k++;
      }
}
void main()
{int x[SIZE][SIZE],y[SIZE+1][3],y0[SIZE+1],y1[SIZE+1],y2[SIZE+1],u[SIZE+1];
 double y0_in[SIZE+1],y00_in[SIZE+1],y1_in[SIZE+1],y2_in[SIZE+1];
 double a[SIZE+1],e[SIZE];
 double out1[SIZE+1],out2[SIZE+1];
 int i,j,k;
 printf("Please input the stream:\n");
 for(i=1;i<=SIZE;i++)
	 for(j=1;j<3;j++)
		 scanf("%d",&x[i][j]);
 printf("\n");
 for(i=1;i<=SIZE;i++)
	 for(j=1;j<3;j++)
         y[i][j]=sqrt(Es)*(2*x[i][j]-1);
 for(i=1;i<=SIZE;i++)   //串并轉換，信道置信度加權//
	 for(j=1;j<3;j++)
	 {if(j==1)
		 {y0[i]=y[i][j];
		  y0_in[i]=Lc*(y0[i]+1/10*NB_RND());
		 }
		 else
		 {if(i%2==1)
			 {y1[i]=y[i][j];
			 y1_in[i]=Lc*(y1[i]+1/10*NB_RND());
			 y2[i]=0;
			 y2_in[i]=Lc*(y2[i]+1/10*NB_RND());
			 }
			 else
			 {y1[i]=0;
			 y1_in[i]=Lc*(y1[i]+1/10*NB_RND());
			 y2[i]=y[i][j];
			 y2_in[i]=Lc*(y2[i]+1/10*NB_RND());
			 }
		 }
	 }
 for(i=1;i<=SIZE;i++)
	 a[i]=0;
 interlace(y0_in,y00_in);
 for(k=1;k<6;k++)    //迭代6次//
 {DEC(a,y0_in,y1_in,e);
  interlace(e,a);
  DEC(a,y00_in,y2_in,e);
  uninterlace(e,a);
 }
 DEC(a,y0_in,y1_in,e);
 interlace(e,a);
 DEC(a,y00_in,y2_in,e);
 for(i=1;i<=SIZE;i++)
    out1[i]=a[i]+e[i]+y00_in[i];
 uninterlace(out1,out2);
 for(i=1;i<=SIZE;i++)  //硬判決，輸出譯碼后的碼字//
 {if(out2[i]>=0)
	u[i]=1;
  else
	u[i]=0;
  printf("%d",u[i]);
 }
}