[MN-dev] [mndiff]: r77 - in trunk/noe: Makefile filerw.c filerw.h galois.c galois.h gfft.c gfft.h noe_internal.h rs.c rs.h rw.c rw.h test.c
michael
subversion at mplayerhq.hu
Tue Jul 10 00:15:17 CEST 2007
Author: michael
Date: Tue Jul 10 00:15:17 2007
New Revision: 77
Log:
version from 2002-11-26 13:12
Added:
trunk/noe/filerw.c
trunk/noe/filerw.h
trunk/noe/gfft.c
trunk/noe/gfft.h
trunk/noe/rw.c
trunk/noe/rw.h
Modified:
trunk/noe/Makefile
trunk/noe/galois.c
trunk/noe/galois.h
trunk/noe/noe_internal.h
trunk/noe/rs.c
trunk/noe/rs.h
trunk/noe/test.c
Modified: trunk/noe/Makefile
==============================================================================
--- trunk/noe/Makefile (original)
+++ trunk/noe/Makefile Tue Jul 10 00:15:17 2007
@@ -1,4 +1,4 @@
-OBJS = test.o galois.o rs.o
+OBJS = test.o galois.o rs.o gfft.o rw.o filerw.o
SRCS = $(OBJS:.o=.c) $(ASM_OBJS:.o=.s)
CFLAGS = -g -Wall -O4 $(OPTFLAGS) -I. $(EXTRA_INC)
Added: trunk/noe/filerw.c
==============================================================================
--- (empty file)
+++ trunk/noe/filerw.c Tue Jul 10 00:15:17 2007
@@ -0,0 +1,107 @@
+/*
+ * Copyright (C) 2002 Michael Niedermayer <michaelni at gmx.at>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
+
+#include <inttypes.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "rw.h"
+#include "filerw.h"
+
+static uint64_t fsize(FILE *f){
+ uint64_t cur= ftell(f);
+ uint64_t size;
+
+ fseek(f, 0, SEEK_END); //FIXME 64bit stuff?
+ size= ftell(f);
+ fseek(f, cur, SEEK_SET);
+
+ return size;
+}
+
+static int fileRead(noe_RwContext *c, uint8_t *data, uint64_t pos, int len){
+ int i;
+ noe_FileRwContext *p= (noe_FileRwContext *)c->priv;
+
+ for(i=0; i<p->fileCount; i++){
+ if(p->offset[i] > pos){
+ int e;
+ uint64_t currentLen= pos - p->offset[i];
+ if(currentLen > len) currentLen= len;
+ if(currentLen <= 0) return -1;
+
+ e= fread(data, currentLen, 1, p->file[i]);
+ if(e!=1) return -1;
+
+ len -= currentLen;
+ data+= currentLen;
+ pos += currentLen;
+ if(len==0) return 0;
+ }
+ }
+ return -1;
+}
+
+static int fileWrite(noe_RwContext *c, uint8_t *data, uint64_t pos, int len){
+ return -1;
+}
+
+static void fileUninit(noe_RwContext *c){
+ int i;
+ noe_FileRwContext *p;
+
+ if(c==NULL) return;
+
+ p = (noe_FileRwContext *)c->priv;
+
+ for(i=0; i<p->fileCount; i++){
+ fclose(p->file[i]);
+ }
+ free(p->file);
+ free(p->offset);
+ free(p);
+ free(c);
+}
+
+noe_RwContext *noe_getFileRwContext(char **fileNames, int fileCount, char *access){
+ int i;
+ uint64_t pos;
+ noe_RwContext *c= malloc(sizeof(noe_RwContext));
+ noe_FileRwContext *p= malloc(sizeof(noe_FileRwContext));
+
+ c->priv= p;
+ c->parent= NULL;
+
+ p->file= malloc(sizeof(FILE)*fileCount);
+ p->offset= malloc(sizeof(uint64_t)*fileCount);
+
+ c->read= fileRead;
+ c->write= fileWrite;
+ c->uninit= fileUninit;
+
+ pos=0;
+ for(i=0; i<fileCount; i++){
+ p->file[i]= fopen(fileNames[i], access);
+ if(p->file[i]==NULL) return NULL;
+
+ pos+= fsize(p->file[i]);
+ p->offset[i]= pos;
+ }
+
+ return c;
+}
Added: trunk/noe/filerw.h
==============================================================================
--- (empty file)
+++ trunk/noe/filerw.h Tue Jul 10 00:15:17 2007
@@ -0,0 +1,25 @@
+/*
+ * Copyright (C) 2002 Michael Niedermayer <michaelni at gmx.at>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
+
+typedef struct noe_FileRwContext{
+ FILE **file;
+ uint64_t *offset;
+ int fileCount;
+}noe_FileRwContext;
+
+//FIXME move to .c if this is the only thing here
\ No newline at end of file
Modified: trunk/noe/galois.c
==============================================================================
--- trunk/noe/galois.c (original)
+++ trunk/noe/galois.c Tue Jul 10 00:15:17 2007
@@ -16,58 +16,75 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
+//FIXME try to use 16bit table with check in noe_log/exp()
+
#include <inttypes.h>
#include <assert.h>
#include <stdio.h>
#include "noe_internal.h"
#include "galois.h"
+#include "gfft.h"
-GF2h16Element noe_exp[65536];
-GF2h16Element noe_log[65536];
-uint32_t noe_mod[256];
+GFF4Element noe_exp[SIZE];
+GFF4Element noe_log[SIZE];
void noe_init(){
- GF2h16Element ge= 1;
+ GFF4Element ge= 1;
int i;
- for(i=0; i<256; i++){
- noe_mod[i]= gf2poly_mod(i<<16, PRIMITIVE_POLYNOM);
- noe_mod[i]^= i<<16;
- }
-
for(i=0; i<SIZE; i++){
noe_exp[i]= ge;
noe_log[ge]= i;
- ge= gf2poly_mod(ge+ge, PRIMITIVE_POLYNOM);
+ ge= prod(ge, PRIMITIVE_ELEMENT);
}
noe_log[0]= 0xFFFF;
noe_log[1]= 0;
+
+ noe_gfft_init();
}
+/**
+ * ...
+ */
+int noe_prodPoly(GFF4Element *dst, GFF4Element *src1, GFF4Element *src2, int order1, int order2){
+ int order= order1 + order2;
+ int i;
-//int zech_log(GaloisField gf, GF2Polynom b){
-//}
+ if(order1*order2 < (order1+order2)*64){
+ for(i=order; i>=0; i--){
+ unsigned int acc=0;
+ int j, end;
-void noe_gfft(GF2h16Element *src, GF2h16Element *dst){
+ j = NOE_MAX(0, i - order2);
+ end= NOE_MIN(i, order1);
-}
+ for(; j<=end; j++){
+ acc+= prod(src1[j], src2[i-j]);
+ }
+ dst[i]= reduce(acc);
+ }
+ }else{
+ const int logSize= noe_log2(order2 + order1)+1;
+ const int size= 1<<logSize;
+ GFF4Element temp[2][size]; //FIXME is this fast?
+ const GFF4Element scale= inv(size);
-int noe_prodPoly(GF2h16Element *dst, GF2h16Element *src1, GF2h16Element *src2, int order1, int order2){
- int order= order1 + order2;
- int i;
-
- for(i=order; i>=0; i--){
- int acc=0;
- int j, end;
-
- j = NOE_MAX(0, i - order2);
- end= NOE_MIN(i, order1);
+ //FIXME avoid mem* (but note the memcpy/set only takes 2% of the prodPoly time)
+ memcpy(temp[0], src1, sizeof(GFF4Element)*(order1+1));
+ memset(temp[0] + order1 + 1, 0, sizeof(GFF4Element)*(size - order1 - 1));
+ memcpy(temp[1], src2, sizeof(GFF4Element)*(order2+1));
+ memset(temp[1] + order2 + 1, 0, sizeof(GFF4Element)*(size - order2 - 1));
- for(; j<=end; j++){
- acc^= prod(src1[j], src2[i-j]);
+ noe_gfft(temp[0], temp[0], logSize);
+ noe_gfft(temp[1], temp[1], logSize);
+
+ for(i=0; i<size; i++){
+ temp[0][i]= prod(prod(temp[0][i], temp[1][i]), scale);
}
- dst[i]= acc;
+
+ noe_igfft(temp[0], temp[0], logSize);
+ memcpy(dst, temp[0], sizeof(GFF4Element)*(order+1));
}
return order;
@@ -75,28 +92,51 @@ int noe_prodPoly(GF2h16Element *dst, GF2
/**
* returns the first order2 coefficients of the product of the 2 polynoms.
+ * Note dst must be at least 2^log2(order1+order2) big FIXME benchmark if its worth it
*/
-int noe_partialProdPoly(GF2h16Element *dst, GF2h16Element *src1, GF2h16Element *src2, int order1, int order2){
+int noe_partialProdPoly(GFF4Element *dst, GFF4Element *src1, GFF4Element *src2, int order1, int order2){
const int order= order2;
int i;
- for(i=order; i>=0; i--){
- int acc=0;
- int j, end;
-
- j = NOE_MAX(0, i - order2);
- end= NOE_MIN(i, order1);
+ if(order1*order2 < (order1+order2)*64){
+ for(i=order; i>=0; i--){
+ unsigned int acc=0;
+ int j, end;
- for(; j<=end; j++){
- acc^= prod(src1[j], src2[i-j]);
+ j = NOE_MAX(0, i - order2);
+ end= NOE_MIN(i, order1);
+
+ for(; j<=end; j++){
+ acc+= prod(src1[j], src2[i-j]);
+ }
+ dst[i]= reduce(acc);
}
- dst[i]= acc;
+ }else{
+ const int logSize= noe_log2(order2 + order1)+1;
+ const int size= 1<<logSize;
+ GFF4Element temp[size]; //FIXME is this fast?
+ const GFF4Element scale= inv(size);
+
+ //FIXME avoid mem* (but note the memcpy/set only takes 2% of the prodPoly time)
+ memcpy(temp, src2, sizeof(GFF4Element)*(order2+1));
+ memset(temp + order2 + 1, 0, sizeof(GFF4Element)*(size - order2 - 1));
+ memcpy(dst, src1, sizeof(GFF4Element)*(order1+1));
+ memset(dst + order1 + 1, 0, sizeof(GFF4Element)*(size - order1 - 1));
+
+ noe_gfft(dst, dst, logSize);
+ noe_gfft(temp, temp, logSize);
+
+ for(i=0; i<size; i++){
+ dst[i]= prod(prod(dst[i], temp[i]), scale);
+ }
+
+ noe_igfft(dst, dst, logSize);
}
return order;
}
-void noe_printPoly(GF2h16Element *src, int order){
+void noe_printPoly(GFF4Element *src, int order){
int i;
for(i=order; i>=0; i--){
@@ -108,59 +148,154 @@ void noe_printPoly(GF2h16Element *src, i
/**
* Evaluates the src polynom at x.
*/
-GF2h16Element noe_evalPoly(GF2h16Element *src, int order, GF2h16Element x){
- int acc=0;
+GFF4Element noe_evalPoly(GFF4Element *src, int order, GFF4Element x){
+ unsigned int acc=0;
int j;
for(j=order; j>=0; j--){
- acc = prod(acc, x) ^ src[j];
+ acc = sum(prod(acc, x), src[j]);
}
return acc;
}
-void noe_getDerivative(GF2h16Element *dst, GF2h16Element *src, int order){
+int noe_getDerivative(GFF4Element *dst, GFF4Element *src, int order){
int i;
- for(i=0; i<order-1; i+=2){
- dst[i ]= src[i+1];
- dst[i+1]= 0;
+ for(i=0; i<order; i++){
+ dst[i ]= prod(src[i+1], i+1);
}
- if(i<order)
- dst[i ]= src[i+1];
+
+ return order-1;
}
-static inline int rev(int x){
- x = (((x & 0xaaaa) >> 1) | ((x & 0x5555) << 1));
- x = (((x & 0xcccc) >> 2) | ((x & 0x3333) << 2));
- x = (((x & 0xf0f0) >> 4) | ((x & 0x0f0f) << 4));
- x = (( x >> 8) | ((x & 0x00ff) << 8));
- return(x);
+/**
+ * ...
+ */
+void noe_synthPoly(GFF4Element *dst, GFF4Element *src, int count){
+ if(count<20){
+ int i;
+
+ dst[0]= 1;
+
+ for(i=0; i<count; i++){
+ noe_prodPoly(dst, dst, src+2*i, i, 1);
+ }
+ }else{
+ int pass, i, countLeft;
+ const int passCount= noe_log2(count-1)+1;
+ GFF4Element temp[2][4<<noe_log2(count-1)];
+
+ assert(passCount>0);
+
+ countLeft= count;
+
+ // noe_printPoly(src, 1);
+ // noe_printPoly(src+2, 1);
+
+ for(i=0; i<(count>>1); i++){
+ noe_prodPoly(temp[0]+4*i, src+4*i, src+4*i+2, 1, 1);
+ }
+ if(count&1){
+ memcpy(temp[0]+4*i, src+4*i, 2*sizeof(GFF4Element));
+ temp[0][4*i+2]= 0;
+ countLeft++;
+ }
+ countLeft>>=1;
+
+ // noe_printPoly(temp[0], 2);
+
+ for(pass=1; countLeft>1; pass++){
+ const int step= 2<<pass;
+ GFF4Element *temp1= temp[ pass&1 ];
+ GFF4Element *temp2= temp[(pass&1)^1];
+
+ for(i=0; i<(countLeft>>1); i++){
+ noe_prodPoly(temp1+2*step*i, temp2+2*step*i, temp2+2*step*i+step, step>>1, step>>1);
+ }
+
+ if(countLeft&1){
+ int len= (step>>1) + 1;
+ memcpy(temp1+2*step*i , temp2+2*step*i, len*sizeof(GFF4Element));
+ memset(temp1+2*step*i + len , 0, step*sizeof(GFF4Element));
+ countLeft++;
+ }
+ countLeft>>=1;
+ }
+
+ assert(pass==passCount);
+
+ // noe_printPoly(temp[0], count);
+
+ memcpy(dst, temp[(pass&1)^1], (count+1)*sizeof(GFF4Element));
+ // noe_printPoly(dst, count);
+ }
}
-void noe_gfft_dit2(GF2h16Element *dst, GF2h16Element *src){
- int i, j, pass;
-
- for(i=0; i<SIZE; i++){
- dst[i]= src[rev(i)];
+int noe_getOrder(GFF4Element *src, int maxOrder){
+ while(maxOrder && src[maxOrder]==0)
+ maxOrder--;
+
+ return maxOrder;
+}
+
+/**
+ * ...
+ * Note: rem can be identical to nom
+ */
+int noe_divPoly(GFF4Element *quot, GFF4Element *rem, GFF4Element *nom, GFF4Element *denom, int nomOrder, int denomOrder){
+ int i;
+ int remOrder= nomOrder;
+ int quotOrder= 0;
+
+ if(rem != nom)
+ memcpy(rem, nom, (nomOrder+1)*sizeof(GFF4Element));
+
+ denomOrder= noe_getOrder(denom, denomOrder);
+
+ quot[0]=0;
+
+ for(; remOrder >= denomOrder; remOrder--){
+ GFF4Element scale;
+
+ if(rem[remOrder]==0) continue;
+
+ scale= prod(rem[remOrder], inv(denom[denomOrder]));
+
+ for(i=0; i<=denomOrder; i++){ //FIXME optimize
+ const int remIndex= i + remOrder - denomOrder;
+ rem[remIndex]= sum(rem[remIndex], SIZE - prod(scale, denom[i]));
+ }
+
+ quot[remOrder - denomOrder]= scale;
+ if(quotOrder==0)
+ quotOrder= remOrder - denomOrder;
}
+
+ return noe_getOrder(rem, remOrder);
+}
- for(pass=0; pass<EXPONENT; pass++){
- int block;
- const int np= 1<<pass;
- const int blockCount= 1<<(EXPONENT-pass-1);
-
- for(block=0; block<blockCount; block++){
- int n;
- const int top= block*np<<1;
- const int bot= top + np;
-
- for(n=0; n<np; n++){
- int a= dst[top + n];
- int b= dst[bot + n];
-
- dst[top + n]= a^b;
+int noe_diffPoly(GFF4Element *dst, GFF4Element *src1, GFF4Element *src2, int order1, int order2){
+ int i;
+ int minOrder= NOE_MIN(order1, order2);
+
+ for(i=0; i<=minOrder; i++){
+ dst[i]= sum(src1[i], SIZE - src2[i]);
+ }
+
+ if(order1==order2){
+ return noe_getOrder(dst, order1);
+ }else{
+ if(dst!=src1){
+ for(; i<=order1; i++){
+ dst[i]= src1[i];
}
}
+
+ for(; i<=order2; i++){
+ dst[i]= SIZE - src2[i];
+ }
+
+ return NOE_MAX(order1, order2);;
}
}
Modified: trunk/noe/galois.h
==============================================================================
--- trunk/noe/galois.h (original)
+++ trunk/noe/galois.h Tue Jul 10 00:15:17 2007
@@ -16,113 +16,48 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
-#define PRIMITIVE_ELEMENT 2
-#define PRIMITIVE_POLYNOM 0x0001002D
-#define EXPONENT 16
-#define SIZE (1<<EXPONENT)
-
-typedef uint16_t GF2h16Element;
+#define PRIMITIVE_ELEMENT 3
+#define SIZE 0x10001
-typedef uint32_t GF2Polynom;
+typedef unsigned int GFF4Element;
-extern GF2h16Element noe_exp[65536];
-extern GF2h16Element noe_log[65536];
-extern uint32_t noe_mod[256];
+extern GFF4Element noe_exp[SIZE];
+extern GFF4Element noe_log[SIZE];
-int noe_prodPoly(GF2h16Element *dst, GF2h16Element *src1, GF2h16Element *src2, int order1, int order2);
-int noe_partialProdPoly(GF2h16Element *dst, GF2h16Element *src1, GF2h16Element *src2, int order1, int order2);
+int noe_prodPoly(GFF4Element *dst, GFF4Element *src1, GFF4Element *src2, int order1, int order2);
+int noe_partialProdPoly(GFF4Element *dst, GFF4Element *src1, GFF4Element *src2, int order1, int order2);
void noe_init();
-void noe_printPoly(GF2h16Element *src, int order);
-void noe_getDerivative(GF2h16Element *dst, GF2h16Element *src, int order);
-GF2h16Element noe_evalPoly(GF2h16Element *src, int order, GF2h16Element x);
-
+void noe_printPoly(GFF4Element *src, int order);
+int noe_getDerivative(GFF4Element *dst, GFF4Element *src, int order);
+GFF4Element noe_evalPoly(GFF4Element *src, int order, GFF4Element x);
+void noe_synthPoly(GFF4Element *dst, GFF4Element *src, int count);
+int noe_divPoly(GFF4Element *quot, GFF4Element *rem, GFF4Element *nom, GFF4Element *denom, int nomOrder, int denomOrder);
+int noe_diffPoly(GFF4Element *dst, GFF4Element *src1, GFF4Element *src2, int order1, int order2);
+int noe_getOrder(GFF4Element *src, int maxOrder);
-static inline GF2Polynom gf2poly_prod(GF2Polynom a, GF2Polynom b){
- int i;
- unsigned int ret=0;
-#if 0
- if(a&1) ret^=b;
- b+=b;
- if(a&2) ret^=b;
- b+=b;
- if(a&4) ret^=b;
- b+=b;
- if(a&8) ret^=b;
- b+=b;
- if(a&16) ret^=b;
- b+=b;
- if(a&32) ret^=b;
- b+=b;
- if(a&64) ret^=b;
- b+=b;
- if(a&128) ret^=b;
- b+=b;
- if(a&256) ret^=b;
- b+=b;
- if(a&512) ret^=b;
- b+=b;
- if(a&1024) ret^=b;
- b+=b;
- if(a&2048) ret^=b;
- b+=b;
- if(a&4096) ret^=b;
- b+=b;
- if(a&8192) ret^=b;
- b+=b;
- if(a&16384) ret^=b;
- b+=b;
- if(a&32768) ret^=b;
- return ret;
-#elif 1
- unsigned int tmp[4];
- tmp[0]= 0;
- tmp[1]= b;
- tmp[2]= b + b;
- tmp[3]= (b + b) ^ b;
+static inline GFF4Element reduce(GFF4Element a){
+ a = (a&0xFFFF) - (a>>16);
+ a += (((int)a)>>31)&0x10001;
- ret ^= tmp[a&3] ; a>>=2;
- ret ^= tmp[a&3]<<2 ; a>>=2;
- ret ^= tmp[a&3]<<4 ; a>>=2;
- ret ^= tmp[a&3]<<6 ; a>>=2;
- ret ^= tmp[a&3]<<8 ; a>>=2;
- ret ^= tmp[a&3]<<10; a>>=2;
- ret ^= tmp[a&3]<<12; a>>=2;
- ret ^= tmp[a&3]<<14;
- return ret;
-#else
- if(a==0 || b==0) return 0;
- ret= noe_log[a] + noe_log[b];
- if(ret >= 0xFFFF) ret-= 0xFFFF;
- return noe_exp[ret];
-#endif
+ return a;
}
-static inline GF2Polynom gf2poly_mod(GF2Polynom a, GF2Polynom b){
- int i;
+static inline GFF4Element prod(GFF4Element a, GFF4Element b){
+//printf("%d %d\n", a,b);
+ if((a&b)==0x10000) return 1;
+ else return reduce(a*b);
+}
- b<<= EXPONENT - 1;
+static inline GFF4Element sum(GFF4Element a, GFF4Element b){
- for(i=0; i<EXPONENT; i++){
- if((a^b) < a) a^=b;
- b>>=1;
- }
+ a += b - 0x10001;
+ a += (((int)a)>>31)&0x10001;
return a;
}
-static inline GF2Polynom gf2poly_mod_m(GF2Polynom a){
- a ^= noe_mod[a>>24]<<8;
- a ^= noe_mod[a>>16];
-
- return a;
-}
-
-static inline GF2h16Element prod(GF2h16Element a, GF2h16Element b){
- return gf2poly_mod_m(gf2poly_prod(a, b));
-}
-
-static inline GF2h16Element inv(GF2h16Element a){
+static inline GFF4Element inv(GFF4Element a){
assert(a!=0);
- return noe_exp[(1<<EXPONENT) - 1 - noe_log[a]];
+ return noe_exp[SIZE - 1 - noe_log[a]];
}
Added: trunk/noe/gfft.c
==============================================================================
--- (empty file)
+++ trunk/noe/gfft.c Tue Jul 10 00:15:17 2007
@@ -0,0 +1,288 @@
+/*
+ * Copyright (C) 2002 Michael Niedermayer <michaelni at gmx.at>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
+
+#include <inttypes.h>
+#include <assert.h>
+#include <stdio.h>
+
+#include "noe_internal.h"
+#include "galois.h"
+#include "gfft.h"
+
+uint8_t noe_revTable[256];
+
+void noe_gfft_init(){
+ int i;
+
+ for(i=0; i<256; i++){
+ int x= i;
+ x = (((x & 0xaa) >> 1) | ((x & 0x55) << 1));
+ x = (((x & 0xcc) >> 2) | ((x & 0x33) << 2));
+ x = (( x >> 4) | ((x & 0x0f) << 4));
+
+ noe_revTable[i]= x;
+ }
+}
+
+int noe_log2(unsigned int v) // from ffmpeg
+{
+ int n;
+
+ n = 0;
+ if (v & 0xffff0000) {
+ v >>= 16;
+ n += 16;
+ }
+ if (v & 0xff00) {
+ v >>= 8;
+ n += 8;
+ }
+ if (v & 0xf0) {
+ v >>= 4;
+ n += 4;
+ }
+ if (v & 0xc) {
+ v >>= 2;
+ n += 2;
+ }
+ if (v & 0x2) {
+ n++;
+ }
+ return n;
+}
+
+void noe_permute(GFF4Element *dst, GFF4Element *src, int outCount){
+ int i;
+
+ if(dst==src){
+ for(i=1; i<outCount; i++){
+ const int r= noe_bitReverse(i);
+ unsigned int t;
+ if(r<=i) continue;
+ t= dst[i];
+ dst[i]= dst[r];
+ dst[r]= t;
+ }
+ }else{
+ for(i=0; i<outCount; i++){
+ dst[i]= src[noe_bitReverse(i)];
+ }
+ }
+}
+
+static inline void fft2(GFF4Element *p){
+ const unsigned int a= p[0];
+ const unsigned int b= p[1];
+
+ p[0]= sum(a, b);
+ p[1]= sum(a, SIZE - b);
+}
+
+static inline void fft4(GFF4Element *p){
+ unsigned int a,b,c,d;
+
+ d= p[0] + SIZE*0x8000ULL;
+
+ a= d +p[2];
+ b= p[1]+p[3];
+ c= d -p[2];
+ d= (p[3] - p[1])<<8 /*(p[1] - p[3])*noe_exp[16384]*/;
+
+ p[0]= reduce(a+b);
+ p[1]= reduce(a-b);
+ p[2]= reduce(c+d);
+ p[3]= reduce(c-d);
+}
+
+/*
+aWi + bWiWn
+(aWi - bWiWn)W2i
+
+(a + bWn)Wi
+(a - bWn)WiW2i
+*/
+static inline void fft8(GFF4Element *p){
+ unsigned int a,b;
+
+ a= p[0];
+ b= p[4];
+ p[0]= a+b;
+ p[4]= a-b;
+
+ a= p[1];
+ b= p[5];
+ p[1]= a+b;
+ p[5]= reduce(SIZE*0x8000ULL + ((a-b)<<12))/*(a-b)*noe_exp[8192]*/;
+
+ a= p[2];
+ b= p[6];
+ p[2]= a+b;
+ p[6]= (b-a)<<8 /*(a-b)*noe_exp[16384]*/;
+
+ a= p[3];
+ b= p[7];
+ p[3]= a+b;
+ p[7]= (a-b)<<4/*(a-b)*noe_exp[3*8192]*/;
+
+ fft4(p);
+ fft4(p+4);
+}
+
+static void fft16(GFF4Element *p){
+ int n;
+
+ for(n=0; n<8; n++){
+ const unsigned int w= noe_exp[ n<<12 ];
+ const unsigned int a= p[n ];
+ const unsigned int b= p[n + 8];
+
+ p[n ]= sum(a, b);
+ p[n + 8]= reduce(sum(a, SIZE - b)*w);
+ }
+
+ fft8(p);
+ fft8(p+8);
+}
+
+static void fftn(GFF4Element *p, int logSize){
+ int n;
+ const int size= 1<<(logSize-1);
+
+ for(n=0; n<size; n++){
+ const unsigned int w= noe_exp[ n<<(16-logSize) ];
+ const unsigned int a= p[n ];
+ const unsigned int b= p[n + size];
+
+ p[n ]= sum(a, b);
+ p[n + size]= reduce(sum(a, SIZE - b)*w);
+ }
+
+ if(logSize==5){
+ fft16(p);
+ fft16(p+16);
+ }else{
+ fftn(p, logSize-1);
+ fftn(p+size, logSize-1);
+ }
+}
+
+static void fftn2(GFF4Element *p, GFF4Element *src, int logSize){
+ int n;
+ const int size= 1<<(logSize-1);
+
+ for(n=0; n<size; n++){
+ const unsigned int w= noe_exp[ n<<(16-logSize) ];
+ const unsigned int a= src[n ];
+ const unsigned int b= src[n + size];
+
+ p[n ]= sum(a, b);
+ p[n + size]= reduce(sum(a, SIZE - b)*w);
+ }
+
+ if(logSize==5){
+ fft16(p);
+ fft16(p+16);
+ }else{
+ fftn(p, logSize-1);
+ fftn(p+size, logSize-1);
+ }
+}
+
+static inline void ifft2(GFF4Element *p){
+ const unsigned int a= p[0];
+ const unsigned int b= p[1];
+
+ p[0]= sum(a, b);
+ p[1]= sum(a, SIZE - b);
+}
+
+static void ifftn(GFF4Element *p, int logSize){
+ int n;
+ const int size= 1<<(logSize-1);
+
+ if(logSize==2){
+ ifft2(p);
+ ifft2(p+2);
+ }else{
+ ifftn(p, logSize-1);
+ ifftn(p+size, logSize-1);
+ }
+
+ for(n=0; n<size; n++){
+ const unsigned int w= noe_exp[ (1<<16) - (n<<(16-logSize)) ];
+ const unsigned int a= p[n ];
+ const unsigned int b= prod(p[n + size], w);
+
+ p[n ]= sum(a, b);
+ p[n + size]= sum(a, SIZE - b);
+ }
+
+}
+
+void noe_gfft(GFF4Element *dst, GFF4Element *src, int logSize){
+// int i, j, pass;
+//printf("%X %X\n", noe_exp[4096], noe_exp[3*4096]);
+
+ assert(logSize>=5);
+
+ if(src==dst)
+ fftn(dst, logSize);
+ else
+ fftn2(dst, src, logSize);
+#if 0
+ memcpy(dst, src, sizeof(GFF4Element)<<logSize);
+ for(pass=0; pass<12; pass++){
+ int block;
+ const int np= 1<<(16-pass-1);
+ const int blockCount= 1<<pass;
+
+ for(block=0; block<blockCount; block++){
+ int n;
+ const int top= block*np<<1;
+ const int bot= top + np;
+
+ for(n=0; n<np; n++){
+ const unsigned int w= noe_exp[ n<<pass ];
+ const unsigned int a= dst[top + n];
+ const unsigned int b= dst[bot + n];
+
+ dst[top + n]= sum(a, b);
+ dst[bot + n]= prod(sum(a, SIZE - b), w);
+ }
+ }
+ }
+
+ {
+ int block;
+ for(block=0; block<4096; block++){
+{START_TIMER
+ fft16(dst + block*16);
+STOP_TIMER}
+ }
+ }
+#endif
+}
+
+void noe_igfft(GFF4Element *dst, GFF4Element *src, int logSize){
+ assert(logSize>=5);
+
+ if(src!=dst)
+ memcpy(dst, src, sizeof(GFF4Element)<<logSize);
+
+ ifftn(dst, logSize);
+}
Added: trunk/noe/gfft.h
==============================================================================
--- (empty file)
+++ trunk/noe/gfft.h Tue Jul 10 00:15:17 2007
@@ -0,0 +1,30 @@
+/*
+ * Copyright (C) 2002 Michael Niedermayer <michaelni at gmx.at>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
+
+void noe_gfft_init();
+void noe_gfft(GFF4Element *dst, GFF4Element *src, int size);
+void noe_igfft(GFF4Element *dst, GFF4Element *src, int size);
+void noe_permute(GFF4Element *dst, GFF4Element *src, int outCount);
+int noe_log2(unsigned int v);
+
+extern uint8_t noe_revTable[256];
+
+static inline unsigned int noe_bitReverse(unsigned int x){
+ return (noe_revTable[x&255]<<8) | (noe_revTable[x>>8]);
+}
+
Modified: trunk/noe/noe_internal.h
==============================================================================
--- trunk/noe/noe_internal.h (original)
+++ trunk/noe/noe_internal.h Tue Jul 10 00:15:17 2007
@@ -18,3 +18,33 @@
#define NOE_MAX(a,b) ((a) > (b) ? (a) : (b))
#define NOE_MIN(a,b) ((a) > (b) ? (b) : (a))
+
+static inline long long rdtsc()
+{
+ long long l;
+ asm volatile( "rdtsc\n\t"
+ : "=A" (l)
+ );
+// printf("%d\n", int(l/1000));
+ return l;
+}
+
+#define START_TIMER \
+static uint64_t tsum=0;\
+static int tcount=0;\
+static int tskip_count=0;\
+uint64_t tend;\
+uint64_t tstart= rdtsc();\
+
+
+#define STOP_TIMER \
+tend= rdtsc();\
+if(tcount<2 || tend - tstart < 4*tsum/tcount){\
+ tsum+= tend - tstart;\
+ tcount++;\
+}else\
+ tskip_count++;\
+if(256*256*256*64%tcount==0){\
+ fprintf(stderr, "%Ld dezicycles, %d runs, %d skips\n", tsum*10/tcount, tcount, tskip_count);\
+ fflush(stderr);\
+}
\ No newline at end of file
Modified: trunk/noe/rs.c
==============================================================================
--- trunk/noe/rs.c (original)
+++ trunk/noe/rs.c Tue Jul 10 00:15:17 2007
@@ -19,61 +19,63 @@
#include <inttypes.h>
#include <assert.h>
#include <malloc.h>
+#include <stdio.h>
#include "noe_internal.h"
#include "galois.h"
#include "rs.h"
+#include "gfft.h"
+
/**
*
*/
-void noe_getRsGenerator(GF2h16Element *dst, int first, int order){
- GF2h16Element factor[2];
+void noe_getRsGenerator(GFF4Element *dst, int first, int order){
int i;
-
- dst[0]= 1;
-
- factor[0]= noe_exp[first];
- factor[1]= 1;
+ GFF4Element factor[2*order];
- for(i=0; i<order; i++){
- int j;
+ factor[0]= SIZE - noe_exp[first];
+ factor[1]= 1;
- noe_prodPoly(dst, dst, factor, i, 1);
- factor[0]= prod(factor[0], PRIMITIVE_ELEMENT);
+ for(i=1; i<order; i++){
+ factor[2*i+0]= prod(factor[2*i-2], PRIMITIVE_ELEMENT);
+ factor[2*i+1]= 1;
}
+
+ noe_synthPoly(dst, factor, order);
}
-void noe_getSyndrom(GF2h16Element *syn, GF2h16Element *src, int first, int order){
- int i;
- //FIXME use GFFT if order is large
+void noe_getSyndrom(GFF4Element *syn, GFF4Element *src, int first, int order){
+ noe_gfft(syn, src, 16);
+ noe_permute(syn, syn, order+1);
syn[0]= 1;
-
+#if 0
for(i=0; i<order; i++){
- syn[i+1]= noe_evalPoly(src, (1<<EXPONENT)-2, noe_exp[first + i]);
+ syn[i+1]= noe_evalPoly(src, SIZE-2, noe_exp[first + i]);
}
+#endif
}
#if 1
noe_RsEncoder *noe_getRsEncoder(int parityCount){
- const int dataCount= (1<<EXPONENT) - parityCount - 1;
+ const int dataCount= SIZE - parityCount - 1;
int i;
- GF2h16Element *locator, *factor;
- GF2h16Element locatorDerivative[parityCount];
+ GFF4Element *locator, *factor;
+ GFF4Element locatorDerivative[parityCount];
noe_RsEncoder *encoder= memalign(16, sizeof(noe_RsEncoder));
- locator= encoder->parityLocator= memalign(16, sizeof(GF2h16Element)*(parityCount+1));
- factor= encoder->parityFactor = memalign(16, sizeof(GF2h16Element)* parityCount );
+ locator= encoder->parityLocator= memalign(16, sizeof(GFF4Element)*(parityCount+1));
+ factor= encoder->parityFactor = memalign(16, sizeof(GFF4Element)* parityCount );
encoder->parityCount= parityCount;
locator[0]= 1;
for(i=0; i<parityCount; i++){
- GF2h16Element locationFactor[2];
+ GFF4Element locationFactor[2];
locationFactor[0] = 1;
- locationFactor[1] = noe_exp[dataCount + i];
- noe_prodPoly(locator, locator, locationFactor, i, 1);
+ locationFactor[1] = SIZE - noe_exp[dataCount + i];
+ noe_prodPoly(locator, locator, locationFactor, i, 1); //FIXME use synthPoly
}
#if 0
for(i=0; i<parityCount; i++){
@@ -84,17 +86,18 @@ for(i=0; i<parityCount; i++){
if(!noe_evalPoly(locator, parityCount, inv(noe_exp[dataCount - 1])))
printf("Internal Error 2\n");
#endif
- noe_getDerivative(locatorDerivative, locator, parityCount);
+ i= noe_getDerivative(locatorDerivative, locator, parityCount);
+ assert(i==parityCount-1);
for(i=0; i<parityCount; i++){
- GF2h16Element X= noe_exp[dataCount + i];
- GF2h16Element invX= noe_exp[(1<<EXPONENT) - dataCount - i - 1];
- GF2h16Element denom;
+ GFF4Element X= noe_exp[dataCount + i];
+ GFF4Element invX= noe_exp[SIZE - dataCount - i - 1];
+ GFF4Element denom;
assert(X == inv(invX));
- denom= noe_evalPoly(locatorDerivative, parityCount-1, invX);
- factor[i]= prod(X, inv(denom));
+ denom= noe_evalPoly(locatorDerivative, parityCount-1, invX); //FIXME do in freq domain if parityCount>1000
+ factor[i]= prod(SIZE - X, inv(denom));
}
return encoder;
@@ -115,35 +118,222 @@ void noe_freeRsEncoder(noe_RsEncoder *en
/**
*
*/
-void noe_rsEncode(noe_RsEncoder *encoder, GF2h16Element *data){
+void noe_rsEncode(noe_RsEncoder *encoder, GFF4Element *data){
int i;
const int parityCount= encoder->parityCount;
- const int dataCount= (1<<EXPONENT) - parityCount - 1;
- GF2h16Element syn[parityCount+1];
- GF2h16Element omega[2*parityCount+1];
- GF2h16Element *locator= encoder->parityLocator;
- GF2h16Element *factor= encoder->parityFactor;
+ const int dataCount= SIZE - parityCount - 1;
+ GFF4Element syn[SIZE - 1];
+ GFF4Element omega[SIZE - 1];
+ GFF4Element *locator= encoder->parityLocator;
+ GFF4Element *factor= encoder->parityFactor;
- memset(data + dataCount, 0, parityCount*sizeof(GF2h16Element));
+ memset(data + dataCount, 0, parityCount*sizeof(GFF4Element));
noe_getSyndrom(syn, data, 1, parityCount);
- if(parityCount < 1000 || 1){
- noe_partialProdPoly(omega, locator, syn, parityCount, parityCount);
-
+ noe_partialProdPoly(omega, locator, syn, parityCount, parityCount);
+
+ if(parityCount < 1000){
for(i=0; i<parityCount; i++){
- GF2h16Element invX= noe_exp[(1<<EXPONENT) - dataCount - i - 1];
+ GFF4Element invX= noe_exp[SIZE - dataCount - i - 1];
- GF2h16Element parity= prod(noe_evalPoly(omega, parityCount, invX), factor[i]);
+ GFF4Element parity= prod(noe_evalPoly(omega, parityCount, invX), factor[i]);
- data[dataCount + i]= parity;
+ data[dataCount + i]= SIZE - parity;
}
}else{
- //FIXME freq domain encoding
+//START_TIMER
+#if 0
+ for(i=0; i<parityCount; i++){
+ GFF4Element invX= noe_exp[SIZE - dataCount - i - 1];
+
+ GFF4Element parity= prod(noe_evalPoly(omega, parityCount, invX), factor[i]);
+
+ data[dataCount + i]= SIZE - parity;
+ }
+#else
+ memset(omega + parityCount + 1, 0, (SIZE - 2 - parityCount)*sizeof(GFF4Element));
+ noe_gfft(omega, omega, 16);
+ for(i=0; i<parityCount; i++){
+ int index= SIZE - dataCount - i - 1;
+
+ GFF4Element parity= prod(omega[noe_bitReverse(index)], factor[i]);
+
+ data[dataCount + i]= SIZE - parity;
+ }
+#endif
+//STOP_TIMER
+ }
+}
+#endif
+
+/**
+ * ...
+ * @param omega must be parityCount+2 big
+ */
+int noe_rsEuclid(GFF4Element *locator[2], GFF4Element *omega[2], int parityCount, int erasureCount, int *locatorOrderP, int *omegaOrderP){
+ int i, quotOrder;
+ int locatorOrder[2]= {0,0};
+ int omegaOrder[2]= {parityCount+1, -1};
+ GFF4Element quot[parityCount+2]; //FIXME size ok?
+
+ memset(omega[0], 0, (parityCount+1)*sizeof(GFF4Element));
+ omega[0][parityCount+1]= 1;
+ locator[0][0]=0;
+ locator[1][0]=1;
+ omegaOrder[1]= noe_getOrder(omega[1], parityCount);
+
+ for(i=0; 1; i++){
+ const int di= i&1;
+ const int si= 1-di;
+#ifdef DEBUG_EUCLID
+printf("Locator0:");noe_printPoly(locator[di], locatorOrder[di]);
+printf("Omega0:");noe_printPoly(omega[di], omegaOrder[di]);
+printf("Locator1:");noe_printPoly(locator[si], locatorOrder[si]);
+printf("Omega1:");noe_printPoly(omega[si], omegaOrder[si]);
+#endif
+ if(2*omegaOrder[si] <= parityCount + erasureCount)
+ break;
+
+ quotOrder= omegaOrder[di] - omegaOrder[si];
+
+ omegaOrder[di]= noe_divPoly(quot, omega[di], omega[di], omega[si], omegaOrder[di], omegaOrder[si]);
+
+ quotOrder= noe_prodPoly(quot, quot, locator[si], quotOrder, locatorOrder[si]);
+ locatorOrder[di]= noe_diffPoly(locator[di], locator[di], quot, locatorOrder[di], quotOrder);
+#ifdef DEBUG_EUCLID
+printf("quot:");noe_printPoly(quot, quotOrder);
+printf("Locator:");noe_printPoly(locator[di], locatorOrder[di]);
+printf("Omega:");noe_printPoly(omega[di], omegaOrder[di]);
+#endif
}
+
+ i++;
+
+ *locatorOrderP= locatorOrder[i&1];
+ *omegaOrderP= omegaOrder[i&1];
+
+ return i&1;
+}
+
+int noe_rsDecode(GFF4Element *data, int *erasure, int erasureCount, int parityCount){
+ int i;
+ int errorCount, psiOrder, gfftEval, omegaOrder, phantomErrorCount;
+ const int dataCount= SIZE - 1 - parityCount;
+ GFF4Element erasureLocator[erasureCount + 1];
+ GFF4Element erasureSynthSrc[erasureCount*2 + 1];
+ GFF4Element locator0[parityCount + 1];
+ GFF4Element locator1[parityCount + 1];
+ GFF4Element *locators[2]={locator0, locator1};
+ GFF4Element omega0[SIZE - 1];
+ GFF4Element omega1[SIZE - 1];
+ GFF4Element *omegas[2]={omega0, omega1}; //FIXME clean this shit up
+ GFF4Element *fErrorLocator, *errorLocator, *omega, *syn, *psi;
+
+ syn= omegas[1];
+ /* kill erased symbols */
+ for(i=0; i<erasureCount; i++){
+//printf("%d\n", erasure[i]);
+ data[ erasure[i] ]=0;
+ }
+
+ noe_getSyndrom(syn, data, 1, parityCount);
+ for(i=1; i<=parityCount; i++){
+ if(syn[i]) break;
+ }
+ if(i>parityCount)
+ return 0;
+
+ phantomErrorCount=0;
+ //FIXME check truncated symbols syndrom
+
+ for(i=0; i<erasureCount; i++){
+ erasureSynthSrc[2*i + 0] = 1;
+ erasureSynthSrc[2*i + 1] = SIZE - noe_exp[erasure[i]];
+ }
+ noe_synthPoly(erasureLocator, erasureSynthSrc, erasureCount);
+// noe_printPoly(erasureLocator, erasureCount);
+#if 0
+for(i=0; i<erasureCount; i++){
+ int eval= noe_evalPoly(erasureLocator, erasureCount, inv(noe_exp[erasure[i]]));
+ assert( eval==0);
}
#endif
+ noe_partialProdPoly(syn, erasureLocator, syn, erasureCount, parityCount);
+ i= noe_rsEuclid(locators, omegas, parityCount, erasureCount, &errorCount, &omegaOrder);
+//printf("errorCount %d \n", errorCount);
+
+ if(i<0) return -1;
+
+ omega= omegas[i];
+ errorLocator= locators[i];
+ fErrorLocator= omegas[1-i]; //reuse some unused space
+ psi= locators[1-i]; //reuse some unused space
+
+ gfftEval= errorCount>20;
+ if(gfftEval && errorCount){
+ memcpy(fErrorLocator, errorLocator, (errorCount+1)*sizeof(GFF4Element));
+ memset(fErrorLocator + errorCount+1, 0, (SIZE - errorCount-2)*sizeof(GFF4Element));
+ noe_gfft(fErrorLocator, fErrorLocator, 16);
+ }
+
+ psiOrder= noe_prodPoly(psi, errorLocator, erasureLocator, errorCount, erasureCount);
+ psiOrder= noe_getDerivative(psi, psi, psiOrder);
+
+ if(errorCount){
+ for(i=0; i<SIZE-1; i++){
+ int r, e;
+ GFF4Element X, invX, nom, denom;
+
+ if(gfftEval){ //FIXME optimize
+ if(fErrorLocator[i])
+ continue;
+
+ r= (- noe_bitReverse(i))&0xFFFF;
+ }else{
+ if(noe_evalPoly(errorLocator, errorCount, noe_exp[i]))
+ continue;
+
+ r= (-i)&0xFFFF;
+ }
+//printf("E%6d %6d\n", i, r);
+
+ X = noe_exp[r];
+ invX= noe_exp[SIZE - r - 1];
+ nom= prod(X, noe_evalPoly(omega, omegaOrder, invX));
+ denom= noe_evalPoly(psi, psiOrder, invX);
+
+ assert(r>=0 && r<=SIZE-2);
+ e= prod(nom, inv(denom));
+ if(data[r]==0 && e==0x10000){
+ if(r < dataCount)
+ return -1;
+ else
+ phantomErrorCount++;
+ }
+ data[r]= sum(data[r], e);
+ }
+ }
+
+ for(i=0; i<erasureCount; i++){
+ const int r= erasure[i];
+ const GFF4Element X = noe_exp[r];
+ const GFF4Element invX= noe_exp[SIZE - r - 1];
+ GFF4Element nom, denom;
+
+ assert(r>=0 && r<0x10000);
+ assert(prod(X, invX)==1);
+
+ nom= prod(X, noe_evalPoly(omega, omegaOrder, invX));
+ denom= noe_evalPoly(psi, psiOrder, invX);
+
+ assert(data[r]==0);
+ data[r]= prod(nom, inv(denom));
+ }
+
+ return erasureCount + errorCount - phantomErrorCount;
+}
Modified: trunk/noe/rs.h
==============================================================================
--- trunk/noe/rs.h (original)
+++ trunk/noe/rs.h Tue Jul 10 00:15:17 2007
@@ -15,25 +15,27 @@
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
-
+
typedef struct noe_RsEncoder{
int parityCount;
- GF2h16Element *parityLocator;
- GF2h16Element *parityFactor;
+ GFF4Element *parityLocator;
+ GFF4Element *parityFactor;
}noe_RsEncoder;
-void noe_getRsGenerator(GF2h16Element *dst, int first, int order);
+void noe_getRsGenerator(GFF4Element *dst, int first, int order);
/**
* gets the syndroms.
* @param src a 2^16 entry long polynom
* @param syn the syndrom polynom will be stored here
*/
-void noe_getSyndrom(GF2h16Element *syn, GF2h16Element *src, int first, int order);
+void noe_getSyndrom(GFF4Element *syn, GFF4Element *src, int first, int order);
noe_RsEncoder *noe_getRsEncoder(int parityOrder);
void noe_freeRsEncoder(noe_RsEncoder *encoder);
/**
*
*/
-void noe_rsEncode(noe_RsEncoder *encoder, GF2h16Element *data);
+void noe_rsEncode(noe_RsEncoder *encoder, GFF4Element *data);
+
+int noe_rsDecode(GFF4Element *data, int *erased, int erasedCount, int parityCount);
Added: trunk/noe/rw.c
==============================================================================
--- (empty file)
+++ trunk/noe/rw.c Tue Jul 10 00:15:17 2007
@@ -0,0 +1,93 @@
+/*
+ * Copyright (C) 2002 Michael Niedermayer <michaelni at gmx.at>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
+
+#include <inttypes.h>
+
+#include "rw.h"
+
+typedef struct noe_Map{
+ uint64_t src;
+ uint64_t dst;
+ uint64_t len;
+}noe_Map;
+
+typedef struct noe_RemapRwContext{
+ Map *map;
+ int count;
+}noe_RemapRwContext;
+
+
+static int noe_remapRead(noe_RwContext *c, uint8_t *data, uint64_t pos, int len){
+ int i;
+ noe_RemapRwContext *p= (noe_RemapRwContext *)c->priv;
+ noe_RwContext *parent= c->parent;
+
+ //FIXME do binary search or cache last
+
+ for(i=0; i<p->count; i++){
+ if(pos < p->map[i].src + p->map[i].len){
+ const uint64_t diff= p->map[i].dst - p->map[i].src;
+ uint64_t currentLen= p->map[i].src + p->map[i].len - pos;
+ if(currentLen > len) currentLen= len;
+
+ parent->read(parent, data, pos + diff, currentLen);
+
+ data+= currentLen;
+ pos+= currentLen;
+ len-= currentLen;
+ if(len == 0) return 0;
+ }
+ }
+
+ return -1;
+}
+
+static void remapUninit(noe_RwContext *c){
+ int i;
+ noe_RemapRwContext *p;
+
+ if(c==NULL) return;
+
+ p = (noe_RemapRwContext *)c->priv;
+
+ free(p->map);
+ free(p);
+ free(c);
+}
+
+noe_RwContext *noe_getRemapRwContext(noe_RwContext *parent, uint64_t *map, int count){
+ int i;
+ uint64_t pos;
+ noe_RwContext *c= malloc(sizeof(noe_RwContext));
+ noe_RemapRwContext *p= malloc(sizeof(noe_RemapRwContext));
+
+ c->priv= p;
+ c->parent= parent;
+
+ p->map= malloc(sizeof(noe_Map)*count);
+ memcpy(p->map, map, sizeof(noe_Map)*count);
+ p->count= count;
+
+ //FIXME sort map ?
+
+ c->read= remapRead;
+ c->write= remapWrite;
+ c->uninit= remapUninit;
+
+ return c;
+}
Added: trunk/noe/rw.h
==============================================================================
--- (empty file)
+++ trunk/noe/rw.h Tue Jul 10 00:15:17 2007
@@ -0,0 +1,25 @@
+/*
+ * Copyright (C) 2002 Michael Niedermayer <michaelni at gmx.at>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
+
+typedef struct noe_RwContext{
+ int (*read )(struct noe_RwContext *context, uint8_t *data, uint64_t pos, int len);
+ int (*write)(struct noe_RwContext *context, uint8_t *data, uint64_t pos, int len);
+ void (*uninit)(struct noe_RwContext *context);
+ struct noe_RwContext *parent;
+ void *priv;
+}noe_RwContext;
\ No newline at end of file
Modified: trunk/noe/test.c
==============================================================================
--- trunk/noe/test.c (original)
+++ trunk/noe/test.c Tue Jul 10 00:15:17 2007
@@ -24,18 +24,41 @@
#include "noe_internal.h"
#include "galois.h"
+#include "gfft.h"
#include "rs.h"
-main(){
+static const unsigned int gfftChecksums[20]={
+0xFCB3E495,
+0x5704AC39,
+0x07F19DC4,
+0xE8E3A16F,
+0xD3FDE16D,
+0xB2F03FC5,
+0xB42E9BB6,
+0x1281C15E,
+0x4F1C1BA2,
+0x0EA4927A,
+0x2E0EBEC7,
+0x2D80F238,
+0x6D40B8AB,
+0x3EF87629,
+0xD1649F79,
+0x24A3D9E4,
+0x4EDC2FE4,
+0xF24C2F4F,
+0x1BFE26C9,
+};
+
+int main(){
int i, j;
- GF2Polynom bp=1;
+ GFF4Element bp=1;
noe_init();
- printf("Testing GF(2^%d)\n", EXPONENT);
- for(i=0; i<(1<<EXPONENT)-1; i++){
+ printf("Testing GF(65537)\n");
+ for(i=0; i<SIZE-1; i++){
if((bp==1 || bp==0) && i>0){
- printf("FAIL\n");
+ printf("FAIL %d %d\n", i, bp);
break;
}
bp= prod(bp, PRIMITIVE_ELEMENT);
@@ -44,24 +67,68 @@ main(){
printf("Testing inverse\n");
bp=1;
- for(i=1; i<(1<<EXPONENT); i++){
+ for(i=1; i<SIZE; i++){
if(prod(bp, inv(bp)) != 1){
printf("FAIL %X %X\n", (int)bp, (int)inv(bp));
break;
}
bp= prod(bp, PRIMITIVE_ELEMENT);
}
+#if 0
+ printf("Experiment");
+ for(i=0; i<(1<<30); i++){
+ unsigned int a= reduce(i*0xFF01);
+ unsigned int b= reduce(reduce(i)*0xFF01);
+
+ if(a!=b) printf("%X %X %X %X\n", i, reduce(i), i*0xFF01, reduce(i)*0xFF01);
+ }
+#endif
+ srand (31415);
- srand (time (0));
-
- printf("Testing generator polynoms");
+ printf("Testing GFFT");
for(i=1; i<20; i++){
- int n= (1<<EXPONENT) - 1;
- int k= (1<<EXPONENT) - 1 - i;
- GF2h16Element syn[i+1];
- GF2h16Element gen[i+1];
- GF2h16Element data[k];
- GF2h16Element code[n];
+ int n= SIZE - 1;
+ GFF4Element syn[n];
+ GFF4Element code[n];
+ int acc=0;
+
+ for(j=0; j<n; j++)
+ code[j]= rand() % 0x10001;
+
+ noe_gfft(syn, code, 16);
+
+ for(j=0; j<n; j++){
+ acc+= syn[j];
+ acc*=2*j+1;
+ }
+ if(gfftChecksums[i-1] != acc)
+ printf("FAIL: 0x%08X != 0x%08X\n", gfftChecksums[i-1], acc);
+
+ noe_igfft(syn, syn, 16);
+ for(j=0; j<n; j++){
+ syn[j]= prod(syn[j], inv(1<<16));
+ }
+
+ for(j=0; j<n; j++){
+ if(syn[j] != code[j]) printf("IGFFT missmatch at %d %X!=%X\n", j, syn[j], code[j]);
+ }
+
+// noe_printPoly(gen, i);
+// noe_printPoly(syn, i-1);
+
+ printf("."); fflush(stdout);
+ }
+
+ srand (31415);
+
+ printf("\nTesting generator polynoms");
+ for(i=1; i<100; i+=5){
+ int n= SIZE - 1;
+ int k= SIZE - 1 - i;
+ GFF4Element syn[n];
+ GFF4Element gen[i+1];
+ GFF4Element data[k];
+ GFF4Element code[n];
for(j=0; j<k; j++)
data[j]= rand() & 0xFFFF;
@@ -82,13 +149,13 @@ main(){
}
printf("\nTesting encoder");
- for(i=1; i<50; i+=10){
- int n= (1<<EXPONENT) - 1;
- int k= (1<<EXPONENT) - 1 - i;
- GF2h16Element syn[i+1];
- GF2h16Element gen[i+1];
- GF2h16Element data[k];
- GF2h16Element code[n];
+ for(i=1; i<1000; i+=200){
+ int n= SIZE - 1;
+ int k= SIZE - 1 - i;
+ GFF4Element syn[n];
+// GFF4Element gen[i+1];
+ GFF4Element data[k];
+ GFF4Element code[n];
noe_RsEncoder *encoder;
int pass=5;
@@ -98,7 +165,7 @@ main(){
for(j=0; j<k; j++)
data[j]= rand() & 0xFFFF;
- memcpy(code, data, n*sizeof(GF2h16Element));
+ memcpy(code, data, n*sizeof(GFF4Element));
noe_rsEncode(encoder, code);
@@ -115,19 +182,87 @@ main(){
}
noe_freeRsEncoder(encoder);
}
+
+ printf("\nTesting decoder");
+ for(i=1; i<2000; i+=100){
+ int n= SIZE - 1;
+ int k= SIZE - 1 - i;
+ GFF4Element syn[n];
+ GFF4Element erased[i];
+ GFF4Element data[k];
+ GFF4Element code[n];
+ noe_RsEncoder *encoder;
+ int pass=5;
+
+ encoder= noe_getRsEncoder(i);
+
+ for(pass=5; pass; pass--){
+ int erasedCount, errorCount;
+ for(j=0; j<k; j++)
+ data[j]= rand() & 0xFFFF;
+
+ memcpy(code, data, n*sizeof(GFF4Element));
+
+ noe_rsEncode(encoder, code);
+
+ erasedCount=0;
+ if(rand()&0x40){
+ erasedCount= (rand()%i)+1;
+ for(j=0; j<erasedCount; j++){
+ int e=-1;
+
+ while(e==-1 || code[e]==0x10000) e=rand()%n;
+
+ erased[j]= e;
+ code[ e ]= 0x10000;
+ }
+ for(j=0; j<erasedCount; j++){
+ int e= erased[j];
+
+ code[ e ]= rand()%0xFFFF;
+ }
+ }
+ errorCount=0;
+ if(erasedCount==0 || (rand()&0x40)){
+ int max= (i - erasedCount)/2;
+ if(max){
+ errorCount= (rand() % max)+1;
+ for(j=0; j<errorCount; j++){
+ int pos= rand()%n;
+ code[pos] = rand()%0xFFFF;
+// printf("P%6d\n", pos);
+ }
+ }
+ }
+//printf("erased:%d errors:%d parity:%d\n", erasedCount, errorCount, i);
+ noe_rsDecode(code, erased, erasedCount, i);
+
+ // noe_printPoly(gen, i);
+ // noe_printPoly(syn, i-1);
+
+ for(j=0; j<k; j++){
+ if(data[j]!=code[j]){
+ printf("FAIL at:%d is %X %X\n", j, data[j], code[j]);
+ break;
+ }
+ }
+ printf("."); fflush(stdout);
+ }
+ noe_freeRsEncoder(encoder);
+ }
return 0;
printf("Testing erasure decoding\n");
for(i=1; i<20; i++){
- int n= (1<<EXPONENT) - 1;
- int k= (1<<EXPONENT) - 1 - i;
- GF2h16Element syn[i+1];
- GF2h16Element gen[i+1];
- GF2h16Element data[k];
- GF2h16Element code[n];
- GF2h16Element locator[i+1];
- GF2h16Element omega[2*i];
+ int n= SIZE - 1;
+ int k= SIZE - 1 - i;
+ GFF4Element syn[n];
+ GFF4Element gen[i+1];
+ GFF4Element data[k];
+ GFF4Element code[n];
+ GFF4Element locator[i+1];
+ GFF4Element omega[2*i];
int locatorOrder=0;
@@ -141,7 +276,7 @@ return 0;
locator[0] = 1;
for(j=0; j<i; j++){
int index;
- GF2h16Element locationFactor[2];
+ GFF4Element locationFactor[2];
if((rand()&7)==0 && j) continue;
@@ -153,8 +288,6 @@ return 0;
code[index] += rand()&0xFFFF;
}
- syn[0] ^= 1;
-
noe_getSyndrom(syn, code, 1, i);
noe_prodPoly(omega, syn, locator, i-1, locatorOrder);
// noe_getDeriative(locator, locator, locatorOrder);
More information about the Mndiff-dev
mailing list