kva.cpp

包括RS码的编码,硬(BM)/软(KV)译码,AWGN信道调制解调仿真. 具体采用何种编译码方案和调制解调方式可在Profile.txt文件中指定(内有详细说明). 且扩展性极好,容易向其中加入新的调

#include <stdlib.h>
#include <stdafx.h>
#include <stdio.h>
#include <math.h>
#include "crtdbg.h"
#include "KVA.h"

//#define ListSize 3
//#define _DEBUGGING

DecAlgInterface InitKVA(const RSCodeParam& rsp, MPAStopCondit MPAStopCdt, float MPAThres)
{
	KVAMember& kv=*((KVAMember*)malloc(sizeof(KVAMember)));
	kv.gf=rsp.gf;
	kv.Q=rsp.N+1;
	kv.N=rsp.N;
	kv.K=rsp.K;
	kv.MPAStopCdt=MPAStopCdt;
	kv.MPAThres=(MPAStopCdt==GivenLambda)? MPAThres : (float)floor(MPAThres);

	switch(MPAStopCdt) {
		case GivenL:
			kv.MaxC=(kv.K-1)*((int)kv.MPAThres+1)*((int)kv.MPAThres+2)/2-1;
			break;
		case GivenC:
			kv.MaxC=(int)kv.MPAThres;
			break;
		case GivenS:
			if ((int)kv.MPAThres>rsp.M) kv.MaxC=(int)kv.MPAThres-rsp.M+(rsp.M+1)*rsp.M/2;
			else kv.MaxC=((int)kv.MPAThres+1)*(int)kv.MPAThres/2;
			break;
		case GivenM:
			kv.MaxC=kv.N*(int)kv.MPAThres*((int)kv.MPAThres+1)/2;
			break;
		case GivenLambda:
			kv.MaxC=kv.N*((int)MPAThres)*((int)MPAThres+1)/2;
			break;
		default:
			_ASSERT(1==2);
			break;
	}
	int qn=kv.Q*kv.N, i;
	kv.M_1D=(int*)malloc(sizeof(int)*qn);	
	kv.indx=(int*)malloc(sizeof(int)*qn);
	kv.P=(double*)malloc(sizeof(double)*qn);
	//快速排序法所用到的堆栈(s:左边界, t:右边界, r:分割点)
	kv.s=(int*)malloc(sizeof(int)*qn);
	kv.t=(int*)malloc(sizeof(int)*qn);
	kv.r=(int*)malloc(sizeof(int)*qn);
	//r[]同时还用作每个递归步执行状态的标志: -1表示初次进入; -2表示已对左半段排序完

	kv.MaxDY=(int)(1+sqrt(1+8*kv.MaxC/(kv.K-1)))/2-1;
    //kv.MaxDX=(int)(kv.MaxC/(kv.MaxDY+1)+kv.MaxDY*(kv.K-1)/2);//恐怕不够用
	kv.MaxDX=kv.MaxC;
	kv.g=(BiPolynom*)malloc(sizeof(BiPolynom)*(kv.MaxDY+1));
	//for (i=0; i<=kv.MaxDY; i++) InitBiPolynom(&kv.g[i], &kv.gf, kv.MaxDX+2, kv.MaxDY);//+2 ??
	for (i=0; i<=kv.MaxDY; i++) InitBiPolynom(&kv.g[i], &kv.gf, kv.MaxDX, kv.MaxDY);//+2 ??
	kv.WD=(int*)malloc(sizeof(int)*(kv.MaxDY+1));
	kv.delta=(GFE*)malloc(sizeof(GFE)*(kv.MaxDY+1));

	kv.YRoots=(Polynom*)malloc(sizeof(Polynom)*(kv.MaxDY+1));
	for (i=0; i<=kv.MaxDY; i++) InitPolynom(&kv.YRoots[i], rsp.gf, kv.K-1);
	kv.roots=(Polynom*)malloc(sizeof(Polynom)*(kv.K+1));
	for (i=0; i<kv.K+1; i++) InitPolynom(&kv.roots[i], rsp.gf, kv.MaxDY);
    kv.Q0=(BiPolynom*)malloc(sizeof(BiPolynom)*(kv.K+1));
	for (i=0; i<kv.K+1; i++) InitBiPolynom(&kv.Q0[i], &kv.gf, kv.MaxDX+(kv.MaxDY-1)*kv.K, kv.MaxDY);
	kv.f=NewPolynom(kv.gf, kv.MaxDY);
	kv.equalroot=(bool*)malloc(sizeof(bool)*(kv.K+1));

	kv.MLCode=(GFE*)malloc(sizeof(GFE)*kv.N);

	DecAlgInterface decalg;
	decalg.AlgStructSize=sizeof(KVAMember);
	decalg.CloseDecAlg=CloseKVA;
	decalg.DecAlgRegs=&kv;
	decalg.DecodeOneWord_Hard=NULL;
	decalg.DecodeOneWord_Soft=DecodeOneWord_KVA;
	decalg.isHardDecoder=false;

	return decalg;
}

void CloseKVA(void* A) {
	KVAMember& kv=*((KVAMember*)A);
	free(kv.M_1D);
	free(kv.indx);
	free(kv.P);
	free(kv.s); free(kv.t); free(kv.r);
	free(kv.WD);
	free(kv.delta);
	for (int i=0; i<=kv.MaxDY; i++) FreeBiPolynom(&kv.g[i]);
	free(kv.g);
	for (i=0; i<=kv.MaxDY; i++) FreePolynom(&kv.YRoots[i]);
	free(kv.YRoots);
	for (i=0; i<kv.K+1; i++) FreePolynom(&kv.roots[i]);
    free(kv.roots);
	for (i=0; i<kv.K+1; i++) FreeBiPolynom(&kv.Q0[i]);
	free(kv.Q0);	
	FreePolynom(&kv.f);
	free(kv.equalroot);
	free(kv.MLCode);	
}


int* MultiplicityMat_Gross(KVAMember& kv, const double* Pie_1D, float lambda) {
	int QN=kv.Q*kv.N;
	int *M_1D=kv.M_1D;
	for (int i=0; i<QN; i++) M_1D[i]=(int)(Pie_1D[i]*lambda);
	return M_1D;
}

int* MultiplicityMat_GS(KVAMember& kv, const double* Pie_1D, int M) {
	int *M_1D=kv.M_1D;
	int iMaxPr, QN=kv.Q*kv.N;
	double maxPr;
	for (int i=0; i<kv.N; i++)
	{
		maxPr=Pie_1D[i];
		M_1D[i]=M;
		iMaxPr=0;
		for (int j=kv.N; j<QN; j+=kv.N)
		{
			if (Pie_1D[j+i]>=maxPr)
			{
				M_1D[iMaxPr+i]=0;
				iMaxPr=j;
				maxPr=Pie_1D[j+i];
				M_1D[j+i]=M;
			}
			else
			{
				M_1D[j+i]=0;
			}
		}
	}
	return M_1D;
}

/*以Algorithm A进行重数分配.返回重数矩阵的首元地址*/
int* MultiplicityMat_AA(KVAMember& kv, const double* Pie_1D, int S, MPAStopCondit stopcdt) {
	int Q=kv.Q, N=kv.N, QN=Q*N;
	int K=kv.K;
	int *M_1D=kv.M_1D;
	int *indx=kv.indx, *s=kv.s, *t=kv.t, *r=kv.r;
	double *P=kv.P;
	//先用快速排序法对P的元素进行排序(实际是移动下标数组indx[]的元素,不改变P阵)
	
	for (int i=0; i<QN; i++)
	{
		M_1D[i]=0;
		indx[i]=i;
		P[i]=Pie_1D[i];
	}
	
	for (i=0; i<QN; i++) r[i]=-1;
	int curDepth=0, j;
	s[0]=0;  t[0]=QN-1;
	double v;//快速排序每一步分割的基准值
    while (curDepth>=0)
	{
		if (r[curDepth]==-1)
		{
		   if (t[curDepth]==s[curDepth])
		   {
			  curDepth--;
			  continue;
		   }

		   v=P[indx[t[curDepth]]];
		   i=s[curDepth]-1; j=t[curDepth]+1;
		   while (i<j)
		   {
			   while (P[indx[++i]]<v);
			   while (P[indx[--j]]>v);
			   if (i<j)
			   {
				  indx[i]^=indx[j];
				  indx[j]^=indx[i];
				  indx[i]^=indx[j];
			   }
		   }
		   r[curDepth]=i;
		   s[curDepth+1]=s[curDepth];
		   t[curDepth+1]=i-1;
		   curDepth++;
		}
		else if (r[curDepth]!=-2)
		{
		   s[curDepth+1]=r[curDepth];
		   t[curDepth+1]=t[curDepth];
		   r[curDepth]=-2;
		   curDepth++;
		}
		else
		{
		   r[curDepth]=-1;
		   curDepth--;
		}
	}

#ifdef _DEBUGGING
	//OnlyForDebug:测试排序是否正确
	for (i=0; i<QN-1; i++) _ASSERT(P[indx[i]]<=P[indx[i+1]]);
#endif
/*
	if (stopcdt==GivenM)
	{//如果是指定重数M则将最大的N个点的重数置为给定值
		for (i=QN-1; i>=QN-kv.N; i--) M_1D[indx[i]]=S;
		return M_1D;
	}
*/
	//否则执行Algorithm A
	//因P已有序,所以每次更新一个Pie*值后的重排序可用二分查找完成,提高了效率
	int C=0;
	//int delta=0, V, tmp, i1=1, i2=0, i3=0;//用于由C计算delta的方法中的变量
	//delta为使线性约束条件至少有C个的最小(1,k-1)加权阶数即delta([1,k-1], C)
	int k, mid;//用于对被更新的pie值进行插入排序的变量
	for (i=0; (stopcdt==GivenS && i<S) || (stopcdt!=GivenS); i++)
	{
		P[indx[QN-1]]/=(double)((++M_1D[indx[QN-1]])+1);//Algorithm A关键步
		C+=M_1D[indx[QN-1]];
		//if (stopcdt==GivenM && M_1D[indx[QN-1]]==S) {QN--; continue;}
		if (stopcdt!=GivenS)
		{	//限定了不超过列表容量L的条件
			/*MaxC已事先计算好,则以下部分可去
			//求出delta([1,k-1], C)
			if (V>=i1-i3)
			{
				V-=(i1-i3);
				i3=0;
				delta++;
				if (++i2==K-1) {i2=0; i1++;}

				tmp=(K-1-i2)*i1;
				if (V>=tmp)
				{
					V-=tmp;
					delta+=(K-1-i2);
					i2=0;
					i1++;
					tmp=i1*(K-1);
					while (V>=tmp)
					{
						V-=tmp;
						i1++;
						delta+=(K-1);
					}
					i2=V/i1;
					delta+=i2;
					i3=V-i1*i2;
				}
				else
				{
					delta+=V/i1;
					i2+=V/i1;
					V%=i1;
					i3=V;
				}
			}
			else i3+=V;
			*/
			//若超限,则恢复到上一步结果并结束迭代过程
			//if (delta/(K-1)>L) {C-=(M_1D[indx[QN-1]]--); break;}
			if (C>kv.MaxC) {C-=(M_1D[indx[QN-1]]--); break;}
		}

		//对改变后的最大pie值进行插入排序
		if (QN>1)
		{   
			v=P[indx[QN-1]];
			if (v<P[indx[QN-2]])
			{
			   j=0; k=QN-2;
			   while (j<k)
			   {
				   mid=(j+k)/2;
				   if (P[indx[mid]]>=v) k=mid; else j=mid+1;
			   }
			   mid=indx[QN-1];
			   for (k=QN-2; k>=j; k--) indx[k+1]=indx[k];
			   indx[j]=mid;
			}
		}
#ifdef _DEBUGGING
	//OnlyForDebug:测试排序是否正确
	for (j=0; j<QN-1; j++) _ASSERT(P[indx[j]]<=P[indx[j+1]]);
#endif

	}
	
#ifdef _DEBUGGING
	if (stopcdt==GivenS) _ASSERT(i==S);
	else _ASSERT(C<=kv.MaxC && (C+M_1D[indx[QN-1]]+1>kv.MaxC));
#endif

	return M_1D;
}


int max(int n1, int n2) {return (n1>n2? n1:n2);}

BiPolynom* Interpolation(KVAMember& kv, int* M_1D) {
  Using_Pow2Vec_Vec2Pow_NG(kv.gf);    
  
  int r, s;
  for (int j=0; j<=kv.MaxDY; j++)
  {
      for (r=0; r<=kv.MaxDX; r++)
	  {
		  for (s=0; s<=kv.MaxDY; s++) kv.g[j].coef[s][r]=0;
	  }
	  kv.g[j].coef[j][0]=1;
      kv.g[j].DegX=0;  kv.g[j].DegY=j;
	  if (j==0) kv.WD[j]=0; else kv.WD[j]=kv.WD[j-1]+kv.K-1;
  }

  int js, minWDeg, curWDeg;
  GFE expa, expb, alpha, beta;
  int globalMaxDegX=0, x, y;
  int maxDegX, maxDegY;
  //迭代C次
  for (int i=0; i<kv.Q*kv.N; i++)
  {
	  if (M_1D[i]==0) continue;
	  else
	  {
		 alpha=Pow2Vec[i%kv.N];
		 beta=(i<kv.N)? 0 : Pow2Vec[i/kv.N-1];
	  }
  	  for (r=0; r<M_1D[i]; r++)
	  for (s=0; s<M_1D[i]-r; s++)
	  {   //计算Delta[j]
		  for (j=0; j<=kv.MaxDY; j++) kv.delta[j]=0;
		  expa=1;
		  for (x=r; x<=globalMaxDegX; x++)
		  {
			  expb=1;
		      for (y=s; y<=kv.MaxDY; y++)
			  {
				  if (CombMod2(x,r) && CombMod2(y,s))
				  {
					  for (j=0; j<=kv.MaxDY; j++)
					  {
						  kv.delta[j]^=GFMPY_V2V(GFMPY_V2V(expa, expb, NG), 
							                  kv.g[j].coef[y][x], NG);
					  }