chris/js8call
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
5d64ac37c5
js8call/lib/qra/qracodes/pdmath.c

386 lines
10 KiB
C
Raw Normal View History

Initial Commit
2018-02-08 21:28:33 -05:00
// pdmath.c
// Elementary math on probability distributions
//
// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy
// ------------------------------------------------------------------------------
// This file is part of the qracodes project, a Forward Error Control
// encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes.
//
// qracodes 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.
// qracodes 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 qracodes source distribution.
// If not, see <http://www.gnu.org/licenses/>.
#include "pdmath.h"
typedef const float *ppd_uniform;
typedef void (*ppd_imul)(float*,const float*);
typedef float (*ppd_norm)(float*);
// define vector size in function of its logarithm in base 2
static const int pd_log2dim[7] = {
1,2,4,8,16,32,64
};
// define uniform distributions of given size
static const float pd_uniform1[1] = {
1.
};
static const float pd_uniform2[2] = {
1./2., 1./2.
};
static const float pd_uniform4[4] = {
1./4., 1./4.,1./4., 1./4.
};
static const float pd_uniform8[8] = {
1./8., 1./8.,1./8., 1./8.,1./8., 1./8.,1./8., 1./8.
};
static const float pd_uniform16[16] = {
1./16., 1./16., 1./16., 1./16.,1./16., 1./16.,1./16., 1./16.,
1./16., 1./16., 1./16., 1./16.,1./16., 1./16.,1./16., 1./16.
};
static const float pd_uniform32[32] = {
1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32.,
1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32.,
1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32.,
1./32., 1./32., 1./32., 1./32.,1./32., 1./32.,1./32., 1./32.
};
static const float pd_uniform64[64] = {
1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64.,
1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64.,
1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64.,
1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64.,
1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64.,
1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64.,
1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64.,
1./64., 1./64., 1./64., 1./64.,1./64., 1./64.,1./64., 1./64.
};
static const ppd_uniform pd_uniform_tab[7] = {
pd_uniform1,
pd_uniform2,
pd_uniform4,
pd_uniform8,
pd_uniform16,
pd_uniform32,
pd_uniform64
};
// returns a pointer to the uniform distribution of the given logsize
const float *pd_uniform(int nlogdim)
{
return pd_uniform_tab[nlogdim];
}
// in-place multiplication functions
// compute dst = dst*src for any element of the distrib
static void pd_imul1(float *dst, const float *src)
{
dst[0] *= src[0];
}
static void pd_imul2(float *dst, const float *src)
{
dst[0] *= src[0]; dst[1] *= src[1];
}
static void pd_imul4(float *dst, const float *src)
{
dst[0] *= src[0]; dst[1] *= src[1];
dst[2] *= src[2]; dst[3] *= src[3];
}
static void pd_imul8(float *dst, const float *src)
{
dst[0] *= src[0]; dst[1] *= src[1]; dst[2] *= src[2]; dst[3] *= src[3];
dst[4] *= src[4]; dst[5] *= src[5]; dst[6] *= src[6]; dst[7] *= src[7];
}
static void pd_imul16(float *dst, const float *src)
{
dst[0] *= src[0]; dst[1] *= src[1]; dst[2] *= src[2]; dst[3] *= src[3];
dst[4] *= src[4]; dst[5] *= src[5]; dst[6] *= src[6]; dst[7] *= src[7];
dst[8] *= src[8]; dst[9] *= src[9]; dst[10]*= src[10]; dst[11]*= src[11];
dst[12]*= src[12]; dst[13]*= src[13]; dst[14]*= src[14]; dst[15]*= src[15];
}
static void pd_imul32(float *dst, const float *src)
{
pd_imul16(dst,src);
pd_imul16(dst+16,src+16);
}
static void pd_imul64(float *dst, const float *src)
{
pd_imul16(dst, src);
pd_imul16(dst+16, src+16);
pd_imul16(dst+32, src+32);
pd_imul16(dst+48, src+48);
}
static const ppd_imul pd_imul_tab[7] = {
pd_imul1,
pd_imul2,
pd_imul4,
pd_imul8,
pd_imul16,
pd_imul32,
pd_imul64
};
// in place multiplication
// compute dst = dst*src for any element of the distrib give their log2 size
// arguments must be pointers to array of floats of the given size
void pd_imul(float *dst, const float *src, int nlogdim)
{
pd_imul_tab[nlogdim](dst,src);
}
static float pd_norm1(float *ppd)
{
float t = ppd[0];
ppd[0] = 1.f;
return t;
}
static float pd_norm2(float *ppd)
{
float t,to;
t =ppd[0]; t +=ppd[1];
if (t<=0) {
pd_init(ppd,pd_uniform(1),pd_log2dim[1]);
return t;
}
to = t;
t = 1.f/t;
ppd[0] *=t; ppd[1] *=t;
return to;
}
static float pd_norm4(float *ppd)
{
float t,to;
t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3];
if (t<=0) {
pd_init(ppd,pd_uniform(2),pd_log2dim[2]);
return t;
}
to = t;
t = 1.f/t;
ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t;
return to;
}
static float pd_norm8(float *ppd)
{
float t,to;
t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3];
t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7];
if (t<=0) {
pd_init(ppd,pd_uniform(3),pd_log2dim[3]);
return t;
}
to = t;
t = 1.f/t;
ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t;
ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t;
return to;
}
static float pd_norm16(float *ppd)
{
float t,to;
t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3];
t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7];
t +=ppd[8]; t +=ppd[9]; t +=ppd[10]; t +=ppd[11];
t +=ppd[12]; t +=ppd[13]; t +=ppd[14]; t +=ppd[15];
if (t<=0) {
pd_init(ppd,pd_uniform(4),pd_log2dim[4]);
return t;
}
to = t;
t = 1.f/t;
ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t;
ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t;
ppd[8] *=t; ppd[9] *=t; ppd[10] *=t; ppd[11] *=t;
ppd[12] *=t; ppd[13] *=t; ppd[14] *=t; ppd[15] *=t;
return to;
}
static float pd_norm32(float *ppd)
{
float t,to;
t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3];
t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7];
t +=ppd[8]; t +=ppd[9]; t +=ppd[10]; t +=ppd[11];
t +=ppd[12]; t +=ppd[13]; t +=ppd[14]; t +=ppd[15];
t +=ppd[16]; t +=ppd[17]; t +=ppd[18]; t +=ppd[19];
t +=ppd[20]; t +=ppd[21]; t +=ppd[22]; t +=ppd[23];
t +=ppd[24]; t +=ppd[25]; t +=ppd[26]; t +=ppd[27];
t +=ppd[28]; t +=ppd[29]; t +=ppd[30]; t +=ppd[31];
if (t<=0) {
pd_init(ppd,pd_uniform(5),pd_log2dim[5]);
return t;
}
to = t;
t = 1.f/t;
ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t;
ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t;
ppd[8] *=t; ppd[9] *=t; ppd[10] *=t; ppd[11] *=t;
ppd[12] *=t; ppd[13] *=t; ppd[14] *=t; ppd[15] *=t;
ppd[16] *=t; ppd[17] *=t; ppd[18] *=t; ppd[19] *=t;
ppd[20] *=t; ppd[21] *=t; ppd[22] *=t; ppd[23] *=t;
ppd[24] *=t; ppd[25] *=t; ppd[26] *=t; ppd[27] *=t;
ppd[28] *=t; ppd[29] *=t; ppd[30] *=t; ppd[31] *=t;
return to;
}
static float pd_norm64(float *ppd)
{
float t,to;
t =ppd[0]; t +=ppd[1]; t +=ppd[2]; t +=ppd[3];
t +=ppd[4]; t +=ppd[5]; t +=ppd[6]; t +=ppd[7];
t +=ppd[8]; t +=ppd[9]; t +=ppd[10]; t +=ppd[11];
t +=ppd[12]; t +=ppd[13]; t +=ppd[14]; t +=ppd[15];
t +=ppd[16]; t +=ppd[17]; t +=ppd[18]; t +=ppd[19];
t +=ppd[20]; t +=ppd[21]; t +=ppd[22]; t +=ppd[23];
t +=ppd[24]; t +=ppd[25]; t +=ppd[26]; t +=ppd[27];
t +=ppd[28]; t +=ppd[29]; t +=ppd[30]; t +=ppd[31];
t +=ppd[32]; t +=ppd[33]; t +=ppd[34]; t +=ppd[35];
t +=ppd[36]; t +=ppd[37]; t +=ppd[38]; t +=ppd[39];
t +=ppd[40]; t +=ppd[41]; t +=ppd[42]; t +=ppd[43];
t +=ppd[44]; t +=ppd[45]; t +=ppd[46]; t +=ppd[47];
t +=ppd[48]; t +=ppd[49]; t +=ppd[50]; t +=ppd[51];
t +=ppd[52]; t +=ppd[53]; t +=ppd[54]; t +=ppd[55];
t +=ppd[56]; t +=ppd[57]; t +=ppd[58]; t +=ppd[59];
t +=ppd[60]; t +=ppd[61]; t +=ppd[62]; t +=ppd[63];
if (t<=0) {
pd_init(ppd,pd_uniform(6),pd_log2dim[6]);
return t;
}
to = t;
t = 1.0f/t;
ppd[0] *=t; ppd[1] *=t; ppd[2] *=t; ppd[3] *=t;
ppd[4] *=t; ppd[5] *=t; ppd[6] *=t; ppd[7] *=t;
ppd[8] *=t; ppd[9] *=t; ppd[10] *=t; ppd[11] *=t;
ppd[12] *=t; ppd[13] *=t; ppd[14] *=t; ppd[15] *=t;
ppd[16] *=t; ppd[17] *=t; ppd[18] *=t; ppd[19] *=t;
ppd[20] *=t; ppd[21] *=t; ppd[22] *=t; ppd[23] *=t;
ppd[24] *=t; ppd[25] *=t; ppd[26] *=t; ppd[27] *=t;
ppd[28] *=t; ppd[29] *=t; ppd[30] *=t; ppd[31] *=t;
ppd[32] *=t; ppd[33] *=t; ppd[34] *=t; ppd[35] *=t;
ppd[36] *=t; ppd[37] *=t; ppd[38] *=t; ppd[39] *=t;
ppd[40] *=t; ppd[41] *=t; ppd[42] *=t; ppd[43] *=t;
ppd[44] *=t; ppd[45] *=t; ppd[46] *=t; ppd[47] *=t;
ppd[48] *=t; ppd[49] *=t; ppd[50] *=t; ppd[51] *=t;
ppd[52] *=t; ppd[53] *=t; ppd[54] *=t; ppd[55] *=t;
ppd[56] *=t; ppd[57] *=t; ppd[58] *=t; ppd[59] *=t;
ppd[60] *=t; ppd[61] *=t; ppd[62] *=t; ppd[63] *=t;
return to;
}
static const ppd_norm pd_norm_tab[7] = {
pd_norm1,
pd_norm2,
pd_norm4,
pd_norm8,
pd_norm16,
pd_norm32,
pd_norm64
};
float pd_norm(float *pd, int nlogdim)
{
return pd_norm_tab[nlogdim](pd);
}
void pd_memset(float *dst, const float *src, int ndim, int nitems)
{
int size = PD_SIZE(ndim);
while(nitems--) {
memcpy(dst,src,size);
dst +=ndim;
}
}
void pd_fwdperm(float *dst, float *src, const int *perm, int ndim)
{
// TODO: non-loop implementation
while (ndim--)
dst[ndim] = src[perm[ndim]];
}
void pd_bwdperm(float *dst, float *src, const int *perm, int ndim)
{
// TODO: non-loop implementation
while (ndim--)
dst[perm[ndim]] = src[ndim];
}
float pd_max(float *src, int ndim)
{
// TODO: faster implementation
float cmax=0; // we assume that prob distributions are always positive
float cval;
while (ndim--) {
cval = src[ndim];
if (cval>=cmax) {
cmax = cval;
}
}
return cmax;
}
int pd_argmax(float *pmax, float *src, int ndim)
{
// TODO: faster implementation
float cmax=0; // we assume that prob distributions are always positive
float cval;
int idxmax=-1; // indicates that all pd elements are <0
while (ndim--) {
cval = src[ndim];
if (cval>=cmax) {
cmax = cval;
idxmax = ndim;
}
}
if (pmax)
*pmax = cmax;
return idxmax;
}
Reference in New Issue Copy Permalink