/*******************************************************************************
*
* subroutines for LDPC codes study, version 3.2.2
* Copyright (C) 2003-2010 Misha Stepanov
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
*
******************************************************************************/
#ifdef LDPCC_RANDOM
#include
#include
#endif
#ifdef LDPCC_LP
#include
#endif
#define LDPCC_BIG_VALUE 1.e+10
#define ERROR(M) \
{ fprintf(stderr, M); exit(1); }
#define ALLOCATE(A, N, T, P) \
(A) = (P)malloc((N) * sizeof(T));\
if ((A) == NULL) ERROR("ERR: malloc\n")
typedef double real;
/*===type_structures==========================================================*/
typedef struct {
int bits; /* length of the codeword (in bits) */
int checks; /* number of parity checks */
int **Hi, **Ha; /* parity check matrix arrays */
int *punctured; /* 0 / 1 --- bit (is not) / (is) punctured */
int *message; /* encoded message, the values are +/- 1 */
real *h; /* log-likelihoods (magnetic fields) h_i */
real *m; /* decoding output (magnetization) m_i */
real **eta, **mu, **lambda, **old_eta, *pre_eta; /* messages */
int HiHa_grown, ID_grown;
/* 0 / 1 --- the corresponding arrays (are not) / (are) generated */
#ifdef LDPCC_LP
LPX *lp; /* linear programming structure */
int LP_grown; /* 0 / 1 --- is (is not) generated */
#endif
} code;
/*----------------------------------------------------------------------------*/
#define Gaussian_channel 0
#define AWGN_channel 0
#define Laplacian_channel 1
#define BSC_channel 8
#define BEC_channel 9
#define isotropic_erasure_channel 19
typedef struct {
real SNR; /* the ratio of *amplitudes* of signal and noise,
for AWGN differs from standard one by SNR -> SNR^2 */
real bit_damage_probability;
#ifdef LDPCC_RANDOM
gsl_rng *RNG;
int RNG_grown;
#endif
int type; /* variants are Gaussian_channel, ..., see definitions above */
} channel;
/*----------------------------------------------------------------------------*/
#define sum_product 0
#define min_sum 1
#define bit_flipping 2
typedef struct {
int n_iter; /* number of iterations in iterative decoding */
/* real (*node)(); */
int WCC; /* != 0 --- with codeword check on each iteration */
int relaxed; real tau; /* relaxation, slowing of iterations */
} iterative_decoder;
/*===channel_noise============================================================*/
#ifdef LDPCC_RANDOM
void kill_RNG(channel *pC)
{
if (pC->RNG_grown)
{
gsl_rng_free(pC->RNG);
pC->RNG_grown = 0;
}
}
void grow_RNG(channel *pC)
{
if (pC->RNG_grown) kill_RNG(pC);
gsl_rng_env_setup();
pC->RNG = gsl_rng_alloc(gsl_rng_ranlxd2);
gsl_rng_set(pC->RNG, (unsigned long int)time(NULL));
pC->RNG_grown = 1;
}
void toss_noise(channel *pC, real *xi, const int bits)
{
int i, sign;
for (i = 0; i < bits; i++)
switch (pC->type)
{
case Gaussian_channel:
xi[i] = gsl_ran_gaussian(pC->RNG, 1. / (pC->SNR));
break;
case Laplacian_channel:
sign = 2 * gsl_rng_uniform_int(pC->RNG, 1) - 1;
xi[i] = sign * gsl_ran_exponential(pC->RNG, 1. / (pC->SNR));
break;
case BSC_channel:
xi[i] = 0.;
if (gsl_rng_uniform(pC->RNG) < pC->bit_damage_probability)
xi[i] = 2.;
break;
case BEC_channel:
xi[i] = 0.;
if (gsl_rng_uniform(pC->RNG) < pC->bit_damage_probability)
xi[i] = 1.;
break;
case isotropic_erasure_channel:
xi[i] = 1. + gsl_ran_gaussian(pC->RNG, 1. / (pC->SNR));
break;
default: /* Gaussian channel */
xi[i] = gsl_ran_gaussian(pC->RNG, 1. / (pC->SNR));
break;
}
}
#endif
/*===general_decoding_subroutines=============================================*/
void send_msg_w2b_calc_h(code *pH, channel *pC, real *xi)
/* initializing h_i --- sending messages from world to bits */
/* h[i] = (1/2) log(p(x_i=+|x_out) / p(x_i=-|x_out)) */
{
int i;
real SNR, SNR2, *h, amp_h;
h = pH->h;
SNR = (pC->SNR); SNR2 = SNR * SNR;
if (pC->type == 8)
{
amp_h = pC->bit_damage_probability;
amp_h = (0.5 * log((1. - amp_h) / amp_h));
}
for (i = 0; i < (pH->bits); i++)
if ((pH->punctured)[i] == 0)
switch (pC->type)
{
case Gaussian_channel:
h[i] = ((real)((pH->message)[i]) - xi[i]) * SNR2; break;
case Laplacian_channel:
h[i] = ((real)((pH->message)[i]) - xi[i]) * SNR; break;
case BSC_channel:
if (xi[i] < 1.) h[i] = amp_h; else h[i] = -amp_h;
h[i] *= ((real)((pH->message)[i]));
break;
case BEC_channel:
if (xi[i] > 0.5) h[i] = 0.;
else h[i] = (real)((pH->message)[i]) * LDPCC_BIG_VALUE;
break;
default: /* Gaussian channel */
h[i] = ((real)((pH->message)[i]) - xi[i]) * SNR2; break;
}
else h[i] = 0.;
}
int check_codeword(code *pH)
/* checks if the decoding output m_i gives a codeword */
{
int a, I, flag, itmp;
for (flag = 1, a = 0; a < (pH->checks); a++)
{
for (itmp = 1, I = 1; I <= (pH->Hi)[a][0]; I++)
if ((pH->m)[(pH->Hi)[a][I]] < 0.) itmp = -itmp;
if (itmp < 0) { flag = 0; break; }
}
return flag;
}
int check_error(code *pH)
{
int i, flag;
for (flag = 0, i = 0; i < (pH->bits); i++)
if ((((real)(pH->message)[i]) * (pH->m)[i]) <= 0.)
{ flag = 1; break; }
return flag;
}
/*===iterative_decoding=======================================================*/
void kill_iterative_decoder(code *pH)
{
int i;
if (pH->ID_grown)
{
if (pH->HiHa_grown == 0)
ERROR("Iterative decoder is grown but not the code!\n")
free(pH->h); free(pH->pre_eta); free(pH->m);
for (i = 0; i < (pH->bits); i++)
{ free((pH->eta)[i]); free((pH->mu)[i]); free((pH->lambda)[i]); }
free(pH->eta); free(pH->mu); free(pH->lambda);
pH->ID_grown = 0;
}
}
void grow_iterative_decoder(code *pH)
{
int i;
if (pH->HiHa_grown == 0)
{
if (pH->ID_grown)
ERROR("Iterative decoder is grown but not the code!\n")
else return;
}
kill_iterative_decoder(pH);
ALLOCATE(pH->eta, pH->bits, real *, real **)
for (i = 0; i < (pH->bits); i++)
{ ALLOCATE((pH->eta)[i], (pH->Ha)[i][0] + 1, real, real *) }
ALLOCATE(pH->mu, pH->bits, real *, real **)
for (i = 0; i < (pH->bits); i++)
{ ALLOCATE((pH->mu)[i], (pH->Ha)[i][0] + 1, real, real *) }
ALLOCATE(pH->lambda, pH->bits, real *, real **)
for (i = 0; i < (pH->bits); i++)
{ ALLOCATE((pH->lambda)[i], (pH->Ha)[i][0] + 1, real, real *) }
ALLOCATE(pH->old_eta, pH->bits, real *, real **)
for (i = 0; i < (pH->bits); i++)
{ ALLOCATE((pH->old_eta)[i], (pH->Ha)[i][0] + 1, real, real *) }
ALLOCATE(pH->pre_eta, pH->bits, real, real *)
pH->ID_grown = 1;
}
int iterative_decoding(code *pH, channel *pC, iterative_decoder *pD, real *xi)
{
int it, bits, checks, i, j, a, J, A, B, sign;
real *h, **eta, **mu, **lambda, **old_eta, *pre_eta, *m;
real abs, min_abs;
if (pH->ID_grown == 0)
ERROR("Iterative decoder is not grown!\n")
/* making variables more local */
bits = pH->bits; checks = pH->checks;
h = pH->h; m = pH->m; eta = pH->eta; mu = pH->mu;
lambda = pH->lambda; old_eta = pH->old_eta; pre_eta = pH->pre_eta;
send_msg_w2b_calc_h(pH, pC, xi);
for (i = 0; i < bits; i++) m[i] = h[i];
/* setting eta_ia = mu_ia = lambda_ia = 0 */
for (i = 0; i < bits; i++)
for (A = 1; A <= (pH->Ha)[i][0]; A++) mu[i][A] = 0.;
if (pD->relaxed)
for (i = 0; i < bits; i++)
for (A = 1; A <= (pH->Ha)[i][0]; A++)
{ eta[i][A] = 0.; lambda[i][A] = 0.; }
/* iterations, main cycle of the decoder */
for (it = 0; it < (pD->n_iter); it++)
{
if (pD->WCC) if (check_codeword(pH)) break;
/* making backup of eta_ia */
if (pD->relaxed)
for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
old_eta[i][A] = eta[i][A];
/* sending messages from bits to checks */
for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
for (eta[i][A] = h[i], B = 1; B <= (pH->Ha)[i][0]; B++)
if (B != A) eta[i][A] += mu[i][B];
if (pD->relaxed)
{
/* calculating lambda_ia from eta_ia, old_eta_ia */
for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
lambda[i][A] += ((pD->tau) * (eta[i][A] - old_eta[i][A]));
/* calculating eta_ia from lambda_ia */
for (i = 0; i < bits; i++)
{
for (pre_eta[i] = 0., A = 1; A <= (pH->Ha)[i][0]; A++)
pre_eta[i] += lambda[i][A];
pre_eta[i] /= (real)((pH->Ha)[i][0] - 1.);
}
for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
eta[i][A] = pre_eta[i] - lambda[i][A];
}
/* sending messages from checks to bits */
for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
/* a = Ha[i][A] */
{
a = (pH->Ha)[i][A];
sign = 1; min_abs = -1.;
for (J = 1; J <= (pH->Hi)[a][0]; J++) if ((pH->Hi)[a][J] != i)
/* j = Hi[a][J] */
{
j = (pH->Hi)[a][J];
for (B = 1; B <= (pH->Ha)[j][0]; B++) if ((pH->Ha)[j][B] == a) break;
abs = eta[j][B];
if (abs < 0.) { sign = -sign; abs = -abs; }
if ((abs < min_abs) || (min_abs < 0.)) min_abs = abs;
}
mu[i][A] = sign * min_abs;
}
/* calculating m_i */
for (i = 0; i < bits; i++)
for (m[i] = h[i], A = 1; A <= (pH->Ha)[i][0]; A++) m[i] += mu[i][A];
}
return it;
}
/*===LP_decoding==============================================================*/
#ifdef LDPCC_LP
void kill_lp_decoder(code *pH)
{
if (pH->LP_grown)
{
if (pH->HiHa_grown == 0)
ERROR("Linear programming decoder is grown but not the code!\n")
glp_delete_prob(pH->lp);
pH->LP_grown = 0;
}
}
void grow_lp_decoder(code *pH)
{
int i, a, flag, rows, row, *ind;
real rflag, *val;
if (pH->HiHa_grown == 0)
{
if (pH->LP_grown)
ERROR("Linear programming decoder is grown but not the code!\n")
else return;
}
kill_lp_decoder(pH);
pH->lp = glp_create_prob();
glp_set_obj_dir(pH->lp, GLP_MAX);
glp_add_cols(pH->lp, pH->bits);
for (i = 0; i < (pH->bits); i++)
glp_set_col_bnds(pH->lp, i + 1, GLP_DB, -1., 1.);
ALLOCATE(ind, (pH->bits) + 1, int, int *)
ALLOCATE(val, (pH->bits) + 1, real, real *)
for (rows = 0, a = 0; a < (pH->checks); a++) if ((pH->Hi)[a][0] > 0)
rows += (1 << ((pH->Hi)[a][0] - 1));
glp_add_rows(pH->lp, rows);
for (row = 1, a = 0; a < (pH->checks); a++)
switch ((pH->Hi)[a][0])
{
case 0: break;
case 1: ind[1] = (pH->Hi)[a][1] + 1; val[1] = 1.;
glp_set_row_bnds(pH->lp, row, GLP_FX, 1., 1.);
glp_set_mat_row(pH->lp, row, 1, ind, val);
row++; break;
case 2: ind[1] = (pH->Hi)[a][1] + 1; val[1] = 1.;
ind[2] = (pH->Hi)[a][2] + 1; val[2] = -1.;
glp_set_row_bnds(pH->lp, row, GLP_FX, 0., 0.);
glp_set_mat_row(pH->lp, row, 2, ind, val);
row ++; break;
default: for (i = 1; i <= (pH->Hi)[a][0]; i++)
{ ind[i] = (pH->Hi)[a][i] + 1; val[i] = 1.; }
for (flag = 0; flag <= (pH->Hi)[a][0]; )
{
for (rflag = 1., i = 1; i <= (pH->Hi)[a][0]; i++)
rflag *= val[i];
if (rflag < 0.) /* not a local codeword */
{
glp_set_row_bnds(pH->lp, row, GLP_UP,
0., ((real)(pH->Hi)[a][0]) - 2.);
glp_set_mat_row(pH->lp, row, (pH->Hi)[a][0], ind, val);
row++;
}
for (flag = 1; flag <= (pH->Hi)[a][0]; flag++)
{
val[flag] *= -1.;
if (val[flag] < 0.) break;
}
}
break;
}
free(ind); free(val);
pH->LP_grown = 1;
}
void lp_decoding(code *pH, channel *pC, real *xi)
{
int i;
if (pH->LP_grown == 0)
{
printf("Linear programming decoder is not grown!\n");
return;
}
send_msg_w2b_calc_h(pH, pC, xi);
for (i = 0; i < (pH->bits); i++)
glp_set_obj_coef(pH->lp, i + 1, (pH->h)[i]);
glp_simplex(pH->lp, NULL);
for (i = 0; i < (pH->bits); i++)
(pH->m)[i] = glp_get_col_prim(pH->lp, i + 1);
}
#ifdef LDPCC_RANDOM
/* Finding low weight pseudocodewords using LP decoding and the
* procedure described in
* M. Chertkov, M.G. Stepanov, An efficient pseudocodeword search
* algorithm for linear programming decoding of LDPC codes, IEEE
* Transactions on Information Theory 54 (4) 1514-1520 (2008).
* [arxiv: cs.IT/0601113] */
real lp_search(code *pH, channel *pC, real *xi)
{
int i, flag;
real sxi, sxi2, old_sxi2;
for (flag = 0, old_sxi2 = 1.e+10, sxi2 = 1.e+9; sxi2 != old_sxi2; )
{
if (flag == 0) toss_noise(pC, xi, pH->bits);
lp_decoding(pH, pC, xi);
for (flag = 0, i = 0; i < (pH->bits); i++)
if ((((real)(pH->message)[i]) * (pH->m)[i]) <= 0.)
{ flag = 1; break; }
if (flag)
{
old_sxi2 = sxi2;
for (sxi = 0., sxi2 = 0., i = 0; i < (pH->bits); i++)
{
xi[i] = 1. - (pH->m)[i];
sxi += fabs(xi[i]); sxi2 += (xi[i] * xi[i]);
}
for (i = 0; i < (pH->bits); i++) xi[i] *= (sxi / sxi2);
for (sxi2 = 0., i = 0; i < (pH->bits); i++) sxi2 += (xi[i] * xi[i]);
}
}
return sxi2;
}
#endif
#endif
/*===parity_check_matrix======================================================*/
void kill_H_matrix(code *pH)
{
int i, a;
if (pH->HiHa_grown)
{
if (pH->ID_grown) kill_iterative_decoder(pH);
#ifdef LDPCC_LP
if (pH->LP_grown) kill_lp_decoder(pH);
#endif
free(pH->h); free(pH->m); free(pH->punctured); free(pH->message);
for (i = 0; i < pH->bits; i++) free((pH->Ha)[i]); free(pH->Ha);
for (a = 0; a < pH->checks; a++) free((pH->Hi)[a]); free(pH->Hi);
pH->HiHa_grown = 0;
}
}
void read_H_matrix(const char *filename, code *pH)
/*
* The format of the file is the following:
* --- begin ----------- | 3 2 |
* #bits #checks | | <- example of the H-matrix file
* matrix from 0s and 1s | 1 1 0 | for the (3, 1) repetition code
* --- end ------------- | 0 1 1 |
* matrix form: #checks rows, #bits columns
*/
{
int i, a, c, *ni, *na;
FILE *in;
if (pH->ID_grown) kill_iterative_decoder(pH);
#ifdef LDPCC_LP
if (pH->LP_grown) kill_lp_decoder(pH);
#endif
if (pH->HiHa_grown) kill_H_matrix(pH);
in = fopen(filename, "r");
fscanf(in, "%d %d", &i, &a); pH->bits = i; pH->checks = a;
ALLOCATE(na, pH->bits, int, int *)
ALLOCATE(ni, pH->checks, int, int *)
for (a = 0; a < (pH->checks); a++)
for (i = 0;i < (pH->bits); i++)
{ fscanf(in, "%d", &c); if (c) { na[i]++; ni[a]++; } }
fclose(in);
ALLOCATE(pH->Hi, pH->checks, int *, int **)
for (a = 0; a < pH->checks; a++)
{
ALLOCATE((pH->Hi)[a], ni[a] + 1, int, int *)
(pH->Hi)[a][0] = ni[a];
}
ALLOCATE(pH->Ha, pH->bits, int *, int **)
for (i = 0; i < pH->bits; i++)
{
ALLOCATE((pH->Ha)[i], na[i] + 1, int, int *)
(pH->Ha)[i][0] = na[i];
}
for (i = 0; i< (pH->bits); i++) na[i] = 0;
for (a = 0; a< (pH->checks); a++) ni[a] = 0;
in = fopen(filename, "r");
fscanf(in, "%d %d", &i, &a);
for (a = 0; a < (pH->checks); a++) for (i = 0; i < (pH->bits); i++)
{
fscanf(in, "%d", &c);
if (c)
{ na[i]++; (pH->Ha)[i][na[i]] = a; ni[a]++; (pH->Hi)[a][ni[a]] = i; }
}
fclose(in); free(na); free(ni);
ALLOCATE(pH->punctured, pH->bits, int, int *)
for (i = 0; i < (pH->bits); i++) (pH->punctured)[i] = 0;
ALLOCATE(pH->message, pH->bits, int, int *)
for (i = 0; i < (pH->bits); i++) (pH->message)[i] = 1;
ALLOCATE(pH->h, pH->bits, real, real *)
ALLOCATE(pH->m, pH->bits, real, real *)
pH->HiHa_grown = 1;
}
void read_short_H_matrix(const char *filename, code *pH)
/*
* The format of the file is the following:
* --- begin ----------- | 3 2 |
* #bits #checks | | <- example of the H-matrix file
* #checks rows with the | 2 0 1 | for the (3, 1) repetition code
* #bits connected and | 2 1 2 |
* their list |
* --- end ------------- |
*/
{
int i, a, I, *ni, *na;
FILE *in;
if (pH->ID_grown) kill_iterative_decoder(pH);
#ifdef LDPCC_LP
if (pH->LP_grown) kill_lp_decoder(pH);
#endif
if (pH->HiHa_grown) kill_H_matrix(pH);
in = fopen(filename, "r");
fscanf(in, "%d %d", &i, &a); pH->bits = i; pH->checks = a;
ALLOCATE(na, pH->bits, int, int *)
ALLOCATE(ni, pH->checks, int, int *)
ALLOCATE(pH->Hi, pH->checks, int *, int **)
for (i = 0; i < (pH->bits); i++) na[i] = 0;
for (a = 0; a < (pH->checks); a++)
{
fscanf(in, "%d", &(ni[a]));
ALLOCATE((pH->Hi)[a], ni[a] + 1, int, int *)
(pH->Hi)[a][0] = ni[a];
for (I = 1; I <= (pH->Hi)[a][0]; I++)
{ fscanf(in, "%d", &i); (pH->Hi)[a][I] = i; na[i]++; }
}
fclose(in); free(ni);
ALLOCATE(pH->Ha, pH->bits, int *, int **)
for (i = 0; i < pH->bits; i++)
{
ALLOCATE((pH->Ha)[i], na[i] + 1, int, int *)
(pH->Ha)[i][0] = na[i]; na[i] = 1;
for (a = 0; a < (pH->checks); a++)
for (I = 1; I <= (pH->Hi)[a][0]; I++)
if ((pH->Hi)[a][I] == i) { (pH->Ha)[i][na[i]] = a; na[i]++; break; }
}
free(na);
ALLOCATE(pH->punctured, pH->bits, int, int *)
for (i = 0; i < (pH->bits); i++) (pH->punctured)[i] = 0;
ALLOCATE(pH->message, pH->bits, int, int *)
for (i = 0; i < (pH->bits); i++) (pH->message)[i] = 1;
ALLOCATE(pH->h, pH->bits, real, real *)
ALLOCATE(pH->m, pH->bits, real, real *)
pH->HiHa_grown = 1;
}
void write_H_matrix(const char *filename, code *pH)
{
int i, a, A, flag;
FILE *out;
out = fopen(filename, "w");
fprintf(out, "%d %d\n\n", pH->bits, pH->checks);
for (a = 0; a < (pH->checks); a++)
{
for (i = 0; i < (pH->bits); i++)
{
for (flag = 0, A = 1; A <= (pH->Ha)[i][0]; A++)
if ((pH->Ha)[i][A] == a) { flag = 1; break; }
fprintf(out, "%d ", flag);
}
fprintf(out, "\n");
}
fclose(out);
}
void write_short_H_matrix(const char *filename, code *pH)
{
int a, I;
FILE *out;
out = fopen(filename, "w");
fprintf(out, "%d %d\n\n", pH->bits, pH->checks);
for (a = 0; a < (pH->checks); a++)
{
fprintf(out, "%d ", (pH->Hi)[a][0]);
for (I = 1; I <= (pH->Hi)[a][0]; I++)
fprintf(out, "%d ", (pH->Hi)[a][I]);
fprintf(out, "\n");
}
fclose(out);
}
/*===diluting_parity_check_matrix=============================================*/
/* "Dendro trick" for lowering connectivity of parity checks described in
* M. Chertkov, M. Stepanov, Pseudo-codeword Landscape,
* IEEE International Symposium on Information Theory 2007
* (24-29 June 2007, Nice, France), pp. 1546-1550.
* DOI: 10.1109/ISIT.2007.4557442
* [arxiv: cs.IT/0701084] */
int dilute_addition(int degree)
{
if (degree < 4) return 0;
else
{
if (degree == 4) return 1;
else
return ((degree >> 1) + dilute_addition((degree + 1) >> 1));
}
}
void dilute_H_matrix(code *pH1, code *pH2, int number, int ind[])
{
int i, j, I, J, a, b, A, B, *ni, *na, *map, flag;
kill_H_matrix(pH2);
pH2->bits = pH1->bits; pH2->checks = pH1->checks;
for (A = 0; A < number; A++)
{
j = dilute_addition((pH1->Hi)[ind[A]][0]);
pH2->bits += j; pH2->checks += j;
}
ALLOCATE(na, pH2->bits, int, int *)
ALLOCATE(ni, pH2->checks, int, int *)
ALLOCATE(pH2->Hi, pH2->checks, int *, int **)
ALLOCATE(pH2->Ha, pH2->bits, int *, int **)
for (i = 0; i < (pH1->bits); i++) na[i] = (pH1->Ha)[i][0];
for (; i < (pH2->bits); i++) na[i] = 2;
/* moving all to be diluted checks to the back, the a'th (not being
* diluted) check goes to map[a]'th check in a new H matrix */
ALLOCATE(map, pH1->bits, int, int *)
for (a = 0, b = 0; a < (pH1->checks); a++)
{
map[a] = -1;
for (flag = 1, A = 0; A < number; A++)
if (a == ind[A]) { flag = 0; break; }
if (flag)
{ map[a] = b; ni[b] = (pH1->Hi)[a][0]; b++; }
}
for (; b < (pH2->checks); b++) ni[b] = 3;
for (a = 0; a < pH2->checks; a++)
{
ALLOCATE((pH2->Hi)[a], ni[a] + 1, int, int *)
(pH2->Hi)[a][0] = ni[a];
}
for (i = 0; i < pH2->bits; i++)
{
ALLOCATE((pH2->Ha)[i], na[i] + 1, int, int *)
(pH2->Ha)[i][0] = na[i];
}
/* untouched part of parity check matrix */
for (i = 0; i < (pH1->bits); i++)
for (A = 1; A <= na[i]; A++) (pH2->Ha)[i][A] = map[(pH1->Ha)[i][A]];
for (; i < (pH2->bits); i++) { (pH2->Ha)[i][1] = (pH2->Ha)[i][2] = -1; }
for (a = 0; a < (pH1->checks); a++) if (map[a] != -1)
for (I = 1; I <= ni[map[a]]; I++) (pH2->Hi)[map[a]][I] = (pH1->Hi)[a][I];
/* diluted checks */
/* I / na is now number / list of bits we are to process */
/* J / B --- first available bit / check */
for (J = (pH1->bits), B = (pH1->checks) - number, A = 0; A < number; A++)
{
I = (pH1->Hi)[ind[A]][0];
for (i = 0; i < I; i++) na[i] = (pH1->Hi)[ind[A]][i + 1];
while (I > 3)
{
if (I == 4)
{
for (j = 0; j < 4; j++)
{
for (b = 1; b <= (pH2->Ha)[na[j]][0]; b++)
if ((pH2->Ha)[na[j]][b] == -1) break;
(pH2->Ha)[na[j]][b] = B + (j >> 1);
(pH2->Hi)[B + (j >> 1)][(j % 2) + 1] = na[j];
}
(pH2->Hi)[B][3] = J; (pH2->Hi)[B + 1][3] = J;
(pH2->Ha)[J][1] = B; (pH2->Ha)[J][2] = B + 1;
J++; B += 2; I = 0;
}
else
{
for (j = 0; j < (I >> 1); j++, J++)
{
for (b = 1; b <= (pH2->Ha)[na[2 * j]][0]; b++)
if ((pH2->Ha)[na[2 * j]][b] == -1) break;
(pH2->Ha)[na[2 * j]][b] = B + j;
for (b = 1; b <= (pH2->Ha)[na[2 * j + 1]][0]; b++)
if ((pH2->Ha)[na[2 * j + 1]][b] == -1) break;
(pH2->Ha)[na[2 * j + 1]][b] = B + j;
(pH2->Hi)[B + j][1] = na[2 * j];
(pH2->Hi)[B + j][2] = na[2 * j + 1];
(pH2->Hi)[B + j][3] = J; (pH2->Ha)[J][2] = B + j;
na[j] = J;
}
if ((I % 2) == 1) na[j] = na[I - 1];
B += (I >> 1); I = ((I + 1) >> 1);
if ((I % 2) == 1)
{
/* flipping the list of bits to make the diluting more symmetric */
for (j = 0; j < (I >> 1); j++)
{ i = na[j]; na[j] = na[I - j - 1]; na[I - j - 1] = i; }
}
}
}
if (I == 3)
{
for (j = 0; j < 3; j++)
{
for (b = 1; b <= (pH2->Ha)[na[j]][0]; b++)
if ((pH2->Ha)[na[j]][b] == -1) break;
(pH2->Ha)[na[j]][b] = B;
(pH2->Hi)[B][j + 1] = na[j];
}
B++;
}
}
free(na); free(ni);
ALLOCATE(pH2->punctured, pH2->bits, int, int *)
for (i = 0; i < (pH1->bits); i++) (pH2->punctured)[i] = (pH1->punctured)[i];
for (; i < (pH2->bits); i++) (pH2->punctured)[i] = 1;
ALLOCATE(pH2->message, pH2->bits, int, int *)
for (i = 0; i < (pH2->bits); i++) (pH2->message)[i] = 1;
ALLOCATE(pH2->h, pH2->bits, real, real *)
ALLOCATE(pH2->m, pH2->bits, real, real *)
pH2->HiHa_grown = 1;
}
/******************************************************************************/