gflib.c

density evolution for LDPC codes design and some related documents

/*
 * GFlib.c
 * Description: Library for arithmetics over GF(2^x) (1<x <=23)
 * Author: Tadashi WADAYAMA
 * Created:1997/11/22
 * Copyright (C) 1997 Tadashi WADAYAMA
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include "GFlib.h"

/*
The library "GFlib" offers arithmetics functions over GF(2^x) and
functions for polynomials over GF(2^x).

0. Initialization

Before using the functions  in the following, you need to 
make tables on Zech logarithm. The function make_table(x)
prepares the tables for GF(2^m).

1. Arithmetics on GF(2^x)

We must declare an variable which stores a value in GF(2^x) 
type "GF". 

    Addition        : add(a,b) or  a^b
    Multiplication  : mult(a,b)
    Inverse         : inv(a)
    x-th power of 
    the primitive 
    element         : alpha(x)
    Power a^b       : power(a,b)
    Division        : gfdiv(a,b)

2. Polynomial 

We must declare an variable which stores a polynomial over GF(2^x) 
type "GF_poly". The type GF_poly consists of two members:

  GF_poly A:
       A.deg  : the degree of the polynomial A
       A.x[i] : the coefficient (in GF(2^x)) corresponding to x^i

The followings are the functions you can use:

  void print_poly(FILE* file,GF_poly *a) 
  Display the polynomial.

  void copy_poly(GF_poly *A, GF_poly *B)
  Copy from A to B.

  void clear_poly(GF_poly *A,int pdeg)
  Clear the polynomial A, where  pdeg is the degree of A.

  void add_poly(GF_poly *A, GF_poly *B, GF_poly *C)
  The result of A + B  is stored in C
  
  void mult_poly(GF_poly *A, GF_poly *B, GF_poly *C)
  The result of A * B is stored in C
  
  void div_poly(GF_poly *X, GF_poly *Y, GF_poly *Q, GF_poly *R)
  The result of X/Y is stored in Q and the remainder is stored in R

  void print_roots(GF_poly *poly)
  The roots of the poly is displayed in standard output.


*/


int* etov; /* etov[exponent] : table for exponent to vector conversion  */
int* vtoe; /* vtoe[vector]   : table for vector to exponent conversion */


int MODULUS;
int MASK;

GF gfdiv(GF a, GF b)
{

  if (a == 0) return 0;

  if (b == 0) {
    fprintf(stderr, "gfdiv : div by zero\n");
    exit(-1);
  }

  return alpha(gflog(a) - gflog(b));
}

GF inv(GF a)
{
  if (a == 1) return  1;
  else return etov[MODULUS-vtoe[a]];
}

GF mult(GF a, GF b)
{
  if (a == 0 || b == 0) return 0;
  else return etov[(vtoe[a] + vtoe[b]) % MODULUS];
}

GF alpha(int x)	
{
  if (x >= 0) return etov[x % MODULUS ];
  else return etov[(MODULUS + x) % MODULUS ];
}

GF power(int a, int b)
{
  if (b == 0) return 1;
  if (a == 0) return 0;
  if (b < 0) {
    b = -b;
    a = inv(a);
  }
  return etov[(vtoe[a] * b) % MODULUS];
}

int make_table(int x)
{
  int i, poly;
  
  switch (x) { 

  case 2: 
    poly =   07; MASK=03;break; /* x^2 + x + 1 */
  case 3: 
    poly =  013; MASK=07;break; /* x^3 + x + 1 */
  case 4: 
    poly =  023; MASK=017;break; /* x^4 + x + 1 */
  case 5: 
    poly =  045; MASK=037;break; /* x^5 + x^2 + 1 */
  case 6: 
    poly = 0103; MASK=077;break; /* x^6 + x + 1 */
  case 7: 
    poly = 0203; MASK=0177;break; /* x^7 + x + 1 */
  case 8: 
    poly = 0435; MASK=0377;break; /* x^8+x^4+x^3+x^2+1 */
				/* 11/22 */
  case 9: 
    poly = 01021; MASK=0777;break; /*  1 000 010 001  */

  case 10: 
    poly = 02011; MASK=01777;break; /* 10 000 001 001 */

  case 11: 
    poly = 04005; MASK=03777;break; /* 100 000 000 101 */

  case 12: 
    poly = 010123; MASK=07777;break; /* 1 000 001 010 011 */

  case 13: 
    poly = 020033; MASK=017777;break; /* 10 000 000 011 011 */

  case 14: 
    poly = 042103; MASK=037777;break; /* 100 010 001 000 011 */

  case 15: 
    poly = 0100003; MASK=077777;break; /* 1 000 000 000 000 011 */

  case 16: 
    poly = 0210013; MASK=0177777;break; /* 10 001 000 000 001 011 */

  case 17: 
    poly = 0400011; MASK=0377777;break; /* 100 000 000 000 001 001 */

  case 18: 
    poly = 01000201; MASK=0777777;break; /* 1 000 000 000 010 000 001 */

  case 19: 
    poly = 02000047; MASK=01777777;break; /* 10 000 000 000 000 100 111 */

  case 20: 
    poly = 04000011; MASK=03777777;break; /* 100 000 000 000 000 001 001 */

  case 21: 
    poly = 010000005; MASK=07777777;break; /* 1 000 000 000 000 000 000 101 */

  case 22: 
    poly = 020000003; MASK=017777777;break; /* 10 000 000 000 000 000 000 011*/

  case 23: 
    poly = 040000041; MASK=037777777;break;/* 100 000 000 000 000 000 100 001*/

  default: fprintf(stderr,"maketable: invalid argument.\n");
    exit(-1);
  } 

  MODULUS = (1<<x)-1;
  if ((etov = (int*)malloc(sizeof(int)*(MODULUS+1))) == NULL) {
    fprintf(stderr,"Can't allocate memory\n");
    exit(-1);
  }
  if ((vtoe = (int*)malloc(sizeof(int)*(MODULUS+1))) == NULL) {
    fprintf(stderr,"Can't allocate memory\n");
    exit(-1);
  }

  etov[0] = 1;
  for (i = 1; i <= MODULUS - 1; i++) {
    etov[i] = etov[i-1] << 1;
    if (((etov[i]>>x)&1) == 1) etov[i] ^= poly;
    vtoe[etov[i]] = i;
  }
  vtoe[0] = -1;
  return MODULUS;
}

void print_poly(FILE* fp, GF_poly *a)
{
  int i;
  for (i = 0; i <= a->deg; i++) {
    if (a->x[i] != 0)  {
      fprintf(fp,"a[%d] x^%d",vtoe[a->x[i]],i);
      if (i != a->deg) fprintf(fp," + ");
    }
  }
  fprintf(fp,"\n");
}

void copy_poly(GF_poly *A, GF_poly *B)
{
  int i;

  B->deg = A->deg;
  for (i = 0; i <= A->deg; i++) {
    B->x[i] = A->x[i];
  }
}

void add_poly(GF_poly *A, GF_poly *B, GF_poly *C)
{
  int i;

    C->deg = A->deg; 

  if (A->deg > B->deg) {
    for (i = B->deg+1; i<=A->deg; i++) B->x[i] = 0;
  }

  if (A->deg < B->deg) {
    C->deg = B->deg; 
    for (i = A->deg+1; i<=B->deg; i++) A->x[i] = 0;
  }

  for (i = 0; i <= C->deg; i++) 
    C->x[i] = add(A->x[i], B->x[i]);

  i = C->deg;
  while ( (C->x[i] == 0) && (i > 0)) i--;
  C->deg = i;

}

void mult_poly(GF_poly *A, GF_poly *B, GF_poly *C)
{
  int i, ia, ib;

  if (((A->deg == 0) && (A->x[0] == 0)) || ((B->deg == 0) && (B->x[0] == 0))) {
    C->deg = 0;
    C->x[0] = 0;
    return;
  }

  C->deg = A->deg + B->deg;

  for (i = 0; i <= C->deg; i++) {
    C->x[i] = 0;
  }

  for (ib = 0; ib <= B->deg; ib++) {
    for(ia = 0; ia <= A->deg; ia++) {
      C->x[ia + ib] ^= mult(A->x[ia],B->x[ib]);
    }
  }
}

void div_poly(GF_poly *X, GF_poly *Y, GF_poly *Q, GF_poly *R)
{
  int i, qi, yi;
  GF inverse_y_head;
  GF quotient;
  if ((Y->deg == 0) && (Y->x[0] == 0)) {
    fprintf(stderr,"div_poly: div by zero\n");
    exit(-1);
  }

  while (X->x[X->deg] == 0)
    X->deg--;

  while (Y->x[Y->deg] == 0)
    Y->deg--;

  inverse_y_head = inv(Y->x[Y->deg]);
  Q->deg = X->deg - Y->deg;
  R->deg = X->deg;

  for (i = 0; i <= R->deg; i++) R->x[i] = X->x[i];

  if (X->deg < Y->deg) {
    Q->deg = 0;
    Q->x[0] = 0;
    R->deg = X->deg;
    return;
  }

  for (qi = 0; qi <= Q->deg; qi++) {
    quotient = mult( R->x[R->deg - qi], inverse_y_head);
    Q->x[Q->deg - qi] = quotient;
    for (yi = 0; yi <= Y->deg; yi++) {
      R->x[R->deg - yi - qi] ^= 
	    mult(quotient,Y->x[Y->deg - yi]);
    }
  }

  i = R->deg;
  while((R->x[i] == 0)&&(i>0)) i--;
  R->deg = i;
}

void print_roots(GF_poly *poly)
{
  GF i;
  printf("roots: ");
  for (i = 0; i <= MODULUS; i++) {
    if (subst(poly,i) == 0) {
      printf("a[%d]  ",gflog(i));
    }
  }
  printf("\n");
}

GF subst(GF_poly *a, GF in)
{
  int i;
  GF sum;

  sum = 0;
  
  if (in == 0) return a->x[0];

  for (i = 0; i <= a->deg; i++) {
    if (a->x[i] != 0) 
      sum ^= mult(a->x[i],etov[(vtoe[in]*i) % MODULUS]);
  }
  return sum;
}

void clear_poly(GF_poly *a,int pdeg)
{
  int i; 
  a->deg = 0;
  for (i = 0; i <= pdeg; i++) {
    a->x[i] = 0;
  }
  a->deg = pdeg;
}