IDA/ida_gf64.cpp

   1
   2
   3 #include <iostream>
   4 #include <fstream>
   5 #include <chrono>
   6 #include <sstream>
   7 #include <ctime>
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <string.h>
  11 #include <algorithm>    // std::random_shuffle
  12 #include <vector>       // std::vector
  13 extern "C" {
  14   #include "jerasure.h"
  15 }
  16
  17 typedef unsigned long mylong;
  18 #define LLUI (long long unsigned int)
  19
  20
  21 using namespace std;
  22
  23
  24 void display(mylong *mat, int r, int c) {
  25  for(int i=0;i<r;i++) {
  26     for(int j=0;j<c;j++) {
  27       printf("%016llu ",LLUI mat[i*c+j]);
  28     }
  29     printf("\n");
  30   }
  31  printf("\n");
  32 }
  33
  34 mylong *matrix_multiply(gf_t *gf, mylong *m1, mylong *m2, int r1, int c1, int r2, int c2, int w)
  35 {
  36   mylong *product;
  37   int i, j, k;
  38
  39   product = (mylong *) malloc(sizeof(mylong)*r1*c2);
  40   for (i = 0; i < r1*c2; i++) product[i] = 0;
  41
  42   for (i = 0; i < r1; i++) {
  43     for (j = 0; j < c2; j++) {
  44       for (k = 0; k < r2; k++) {
  45         product[i*c2+j] ^= gf->multiply.w64(gf,m1[i*c1+k], m2[k*c2+j]);
  46       }
  47     }
  48   }
  49   return product;
  50 }
  51
  52
  53
  54 int invert_matrix(gf_t *gf, mylong *mat, mylong *inv, int rows)
  55 {
  56   int cols, i, j, k, x, rs2;
  57   int row_start;
  58   mylong tmp, inverse;
  59
  60   cols = rows;
  61
  62   k = 0;
  63   for (i = 0; i < rows; i++) {
  64     for (j = 0; j < cols; j++) {
  65       inv[k] = (i == j) ? 1 : 0;
  66       k++;
  67     }
  68   }
  69 //  display(inv, rows, rows);
  70 //  printf("\n");
  71
  72   /* First -- convert into upper triangular  */
  73   for (i = 0; i < cols; i++) {
  74     row_start = cols*i;
  75
  76     /* Swap rows if we ave a zero i,i element.  If we can't swap, then the
  77        matrix was not invertible  */
  78
  79     if (mat[row_start+i] == 0) {
  80       for (j = i+1; j < rows && mat[cols*j+i] == 0; j++) ;
  81       if (j == rows) return -1;
  82       rs2 = j*cols;
  83       for (k = 0; k < cols; k++) {
  84         tmp = mat[row_start+k];
  85         mat[row_start+k] = mat[rs2+k];
  86         mat[rs2+k] = tmp;
  87         tmp = inv[row_start+k];
  88         inv[row_start+k] = inv[rs2+k];
  89         inv[rs2+k] = tmp;
  90       }
  91     }
  92
  93     /* Multiply the row by 1/element i,i  */
  94     tmp = mat[row_start+i];
  95     if (tmp != 1) {
  96       inverse = gf->divide.w64(gf,1, tmp);
  97       for (j = 0; j < cols; j++) {
  98         mat[row_start+j] = gf->multiply.w64(gf,mat[row_start+j], inverse);
  99         inv[row_start+j] = gf->multiply.w64(gf,inv[row_start+j], inverse);
 100       }
 101     }
 102
 103     /* Now for each j>i, add A_ji*Ai to Aj  */
 104     k = row_start+i;
 105     for (j = i+1; j != cols; j++) {
 106       k += cols;
 107       if (mat[k] != 0) {
 108         if (mat[k] == 1) {
 109           rs2 = cols*j;
 110           for (x = 0; x < cols; x++) {
 111             mat[rs2+x] ^= mat[row_start+x];
 112             inv[rs2+x] ^= inv[row_start+x];
 113           }
 114         } else {
 115           tmp = mat[k];
 116           rs2 = cols*j;
 117           for (x = 0; x < cols; x++) {
 118             mat[rs2+x] ^= gf->multiply.w64(gf,tmp, mat[row_start+x]);
 119             inv[rs2+x] ^= gf->multiply.w64(gf,tmp, inv[row_start+x]);
 120           }
 121         }
 122       }
 123     }
 124   }
 125
 126   /* Now the matrix is upper triangular.  Start at the top and multiply down  */
 127
 128   for (i = rows-1; i >= 0; i--) {
 129     row_start = i*cols;
 130     for (j = 0; j < i; j++) {
 131       rs2 = j*cols;
 132       if (mat[rs2+i] != 0) {
 133         tmp = mat[rs2+i];
 134         mat[rs2+i] = 0;
 135         for (k = 0; k < cols; k++) {
 136           inv[rs2+k] ^= gf->multiply.w64(gf,tmp, inv[row_start+k]);
 137         }
 138       }
 139     }
 140   }
 141
 142 /*  printf("mat\n");
 143   display(mat, rows, rows);
 144   printf("\n");
 145    printf("inv\n");
 146   display(inv, rows, rows);
 147   printf("\n");
 148 */
 149   return 0;
 150 }
 151
 152
 153
 154
 155 int invertible_matrix(gf_t *gf, int *mat, int rows, int w)
 156 {
 157   int cols, i, j, k, x, rs2;
 158   int row_start;
 159   mylong tmp, inverse;
 160
 161   cols = rows;
 162
 163   /* First -- convert into upper triangular  */
 164   for (i = 0; i < cols; i++) {
 165     row_start = cols*i;
 166
 167     /* Swap rows if we ave a zero i,i element.  If we can't swap, then the
 168        matrix was not invertible  */
 169
 170     if (mat[row_start+i] == 0) {
 171       for (j = i+1; j < rows && mat[cols*j+i] == 0; j++) ;
 172       if (j == rows) return 0;
 173       rs2 = j*cols;
 174       for (k = 0; k < cols; k++) {
 175         tmp = mat[row_start+k];
 176         mat[row_start+k] = mat[rs2+k];
 177         mat[rs2+k] = tmp;
 178       }
 179     }
 180
 181     /* Multiply the row by 1/element i,i  */
 182     tmp = mat[row_start+i];
 183     if (tmp != 1) {
 184       inverse =  gf->divide.w64(gf,1, tmp);
 185       for (j = 0; j < cols; j++) {
 186         mat[row_start+j] =  gf->multiply.w64(gf,mat[row_start+j], inverse);
 187       }
 188     }
 189
 190     /* Now for each j>i, add A_ji*Ai to Aj  */
 191     k = row_start+i;
 192     for (j = i+1; j != cols; j++) {
 193       k += cols;
 194       if (mat[k] != 0) {
 195         if (mat[k] == 1) {
 196           rs2 = cols*j;
 197           for (x = 0; x < cols; x++) {
 198             mat[rs2+x] ^= mat[row_start+x];
 199           }
 200         } else {
 201           tmp = mat[k];
 202           rs2 = cols*j;
 203           for (x = 0; x < cols; x++) {
 204             mat[rs2+x] ^=  gf->multiply.w64(gf,tmp,mat[row_start+x]);
 205           }
 206         }
 207       }
 208     }
 209   }
 210   return 1;
 211 }
 212
 213
 214
 215
 216
 217 mylong*  readFullFile(int n, int t, mylong& sizeFile, mylong & padded_size) {
 218
 219 //  ifstream stream("lena.png", ios::in | ios::binary | ios::ate);
 220 //  ifstream stream("lena_small.png", ios::in | ios::binary | ios::ate);
 221   ifstream stream("/home/couturie/Downloads/CARCARIASS.zip", ios::in | ios::binary | ios::ate);
 222
 223   sizeFile=stream.tellg();
 224   std::cout << sizeFile << std::endl;
 225   stream.seekg(0, ios::beg);
 226
 227
 228
 229
 230
 231
 232   vector<uint8_t> contents((istreambuf_iterator<char>(stream)), istreambuf_iterator<char>());
 233
 234
 235
 236
 237
 238
 239
 240   //make padding, we need to pad to be divisible by 8*t, we
 241   if((sizeFile+8)%(8*t)!=0) {
 242     cout<<(int)(sizeFile/(8*t))<<endl;
 243
 244     int remainder=(8*t)*(1+(int)((sizeFile+8)/(8*t)))-sizeFile-8;
 245     cout << "remainder " << remainder << endl;
 246     uint8_t vec[remainder];
 247     for(int i=0;i<remainder;i++)
 248       vec[i]=0;
 249     //add remainder elements at the end
 250     contents.insert(contents.end(),vec,vec+remainder);
 251     //add 8 elements at the beginning for the size
 252     contents.insert(contents.begin(),vec,vec+8);
 253   }
 254
 255
 256
 257   cout << "res contents " << contents.size() << endl;
 258
 259 //  for(int i=0;i<contents.size();i++)
 260 //    cout << (int)contents[i] << " ";
 261   cout << endl;
 262
 263   uint8_t *p_contents=&contents[0];
 264 //  mylong *p_contents2=reinterpret_cast<mylong*>(p_contents);
 265
 266   padded_size=contents.size()/8;
 267
 268   mylong *p_contents2=new mylong[padded_size];
 269   memcpy(p_contents2,p_contents,sizeof(mylong)*padded_size);
 270          //mylong *p_contents2=(mylong*)p_contents;
 271
 272   p_contents2[0]=sizeFile;
 273
 274
 275
 276
 277 /*  for(int i=0;i<padded_size;i++)
 278     cout << p_contents2[i] << " ";
 279   cout << endl;
 280 */
 281   /* long res=0;
 282   for(int i=8-1;i>=0;i--) {
 283     res<<=8;
 284     res+=p_contents[i];
 285   }
 286
 287   cout << "convert val " << (long)res << endl;
 288
 289   res=0;
 290   for(int i=16-1;i>=8;i--) {
 291     res<<=8;
 292     res+=p_contents[i];
 293   }
 294
 295   cout << "convert val " << (long)res << endl;
 296   */
 297
 298   return p_contents2;
 299 }
 300
 301
 302 void saveFile(uint8_t *data, const char *fileName,long size_file) {
 303   cout<<"size file "<<size_file<<endl;
 304   FILE* pFile = fopen (fileName, "wb");
 305   fwrite (data , sizeof(uint8_t), size_file, pFile);
 306
 307   fclose (pFile);
 308 }
 309
 310
 311 void readFile(uint8_t *data, const char *fileName,long size_file) {
 312   cout<<"size file "<<size_file<<endl;
 313   FILE* pFile = fopen (fileName, "rb");
 314
 315   fseek(pFile, 0L, SEEK_END);
 316   int sz = ftell(pFile);
 317   rewind(pFile);
 318   if(sz!=size_file) {
 319     cout << "error : size of this chunk is not correct !! " << endl;
 320     exit(0);
 321   }
 322   fread (data , sizeof(uint8_t), size_file, pFile);
 323
 324   fclose (pFile);
 325 }
 326
 327
 328 void saveFile8(mylong *data, const char *fileName,long size_file) {
 329
 330   cout<<"size file 8 "<<size_file<<endl;
 331       display(data,1,1);
 332     FILE* pFile = fopen (fileName, "wb");
 333
 334
 335
 336     fwrite (data , sizeof(mylong), size_file, pFile);
 337     fclose (pFile);
 338
 339
 340
 341 }
 342
 343
 344 int main(int argc, char **argv)
 345 {
 346   unsigned int  w;
 347   int invert;
 348   mylong *matS;
 349   mylong *matG;
 350   mylong *matC;
 351   mylong *inverse;
 352   mylong *identity;
 353
 354
 355   int size=5000000;
 356   int t=4;
 357   int n=8;
 358
 359
 360   w=64;
 361   gf_t gf;
 362   gf_init_easy(&gf, w);
 363
 364   mylong padded_size;
 365   mylong sizeFile;
 366   matS=readFullFile(n,t,sizeFile,padded_size);
 367
 368   cout<<padded_size*8<<endl;
 369
 370 /* for(int i=0;i<padded_size;i++)
 371     cout << matS[i] << " ";
 372   cout << endl;
 373 */
 374
 375   int len=padded_size/t;
 376 //  int len=size/t;
 377   cout<<" len "<<len<<" "<<padded_size<<endl;
 378
 379   // matS = malloc(sizeof(mylong)*t*len);
 380   matG = malloc(sizeof(mylong)*t*n);
 381
 382
 383
 384
 385   srand48(time(0));
 386   std::srand ( unsigned ( std::time(0) ) );
 387
 388 /*
 389   for(int i=0;i<t;i++) {
 390     for(int j=0;j<len;j++) {
 391       matS[i*len+j]=lrand48()<<32|lrand48();
 392     }
 393   }
 394   */
 395
 396   for(int i=0;i<n;i++) {
 397     for(int j=0;j<t;j++) {
 398       matG[i*t+j]=lrand48()<<32|lrand48();
 399     }
 400   }
 401
 402
 403
 404   printf("Matrix S:\n");
 405 //  display(matS, len, t);
 406   display(matS, t, t);
 407
 408   printf("Matrix G:\n");
 409   display(matG, n, t);
 410
 411
 412   auto start = std::chrono::system_clock::now();
 413   matC=matrix_multiply(&gf, matG, matS, n, t, t, len, w);
 414   auto end = std::chrono::system_clock::now();
 415   std::chrono::duration<double> elapsed_seconds = end-start;
 416   std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";
 417
 418 //  display(matC,t,t);
 419
 420
 421   //Save trunks
 422   for(int i=0;i<n;i++) {
 423     stringstream ss;
 424     ss <<"lena_"<<i<<".png";
 425     string str = ss.str();
 426     saveFile((uint8_t*)&matC[i*len], str.c_str(),len*sizeof(mylong));
 427   }
 428
 429
 430
 431   mylong *matCs = malloc(sizeof(mylong)*t*len);
 432   mylong *matGs = malloc(sizeof(mylong)*t*t);
 433
 434
 435   std::vector<int> myvector;
 436
 437   // set some values:
 438   for (int i=0; i<n; ++i) myvector.push_back(i); // 1 2 3 4 5 6 7 8 9
 439
 440   // using built-in random generator:
 441   random_shuffle ( myvector.begin(), myvector.end() );
 442
 443
 444
 445   std::cout << "random chunk" << std::endl;
 446   int ind=0;
 447   for(int ii=0;ii<t;ii++) {
 448 //  for(int i=n-1;i>=t;i--) {
 449 //  for(int i=0;i<n;i+=2) {
 450
 451     auto i=myvector[ii];
 452     std::cout << myvector[i] << " ";
 453
 454
 455     stringstream ss;
 456     ss <<"lena_"<<i<<".png";
 457     string str = ss.str();
 458     readFile((uint8_t*)&matCs[ind*len], str.c_str(),len*sizeof(mylong));
 459
 460 //    display(&matCs[ind*len],1,1);
 461 //    display(&matC[i*len],1,1);
 462
 463     /*for(int j=0;j<len;j++) {
 464       matCs[ind*len+j]=matC[i*len+j];
 465       }*/
 466
 467
 468
 469     for(int j=0;j<t;j++) {
 470       matGs[ind*t+j]=matG[i*t+j];
 471     }
 472     ind++;
 473   }
 474   std::cout << std::endl;
 475
 476
 477   printf("Matrix Gs:\n");
 478   display(matGs, t, t);
 479
 480 //  printf("Matrix Cs:\n");
 481 //  display(matCs, t, len);
 482
 483   mylong* matGs_copy = malloc(sizeof(mylong)* t*t);
 484
 485   //WARNING invert changes the matrix
 486 //  invert = invertible_matrix(&gf,matGs, t, w);
 487 //  printf("\nInvertible Gs: %s\n", (invert == 1) ? "Yes" : "No");
 488
 489   memcpy(matGs_copy, matGs, sizeof(mylong)*t*t);
 490   mylong *matInvGs = malloc(sizeof(mylong)*t*t);
 491   invert_matrix(&gf, matGs_copy, matInvGs, t);
 492
 493   /*
 494   WARNING   this changes matGs
 495   identity=matrix_multiply(&gf, matInvGs, matGs, t, t, t, t, w);
 496   printf("Identity:\n");
 497   display(identity, t, t);
 498   */
 499
 500
 501   //printf("Matrix Gs:\n");
 502   //display(matGs, t, t);
 503
 504   mylong *matS2 = malloc(sizeof(mylong)*t*len);
 505
 506
 507   start = std::chrono::system_clock::now();
 508   matS2=matrix_multiply(&gf, matInvGs, matCs, t, t, t, len, w);
 509   end = std::chrono::system_clock::now();
 510   elapsed_seconds = end-start;
 511   std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";
 512
 513
 514   printf("Matrix S2:\n");
 515 //  display(matS2, len, t);
 516   display(matS2, t, t);
 517
 518   int equal=1;
 519   for(int i=0;i<padded_size && equal;i++) {
 520       equal=matS[i]==matS2[i];
 521 //      if(!equal)
 522 //      printf("ARGH\n");
 523   }
 524   if(equal)
 525     printf("EQUAL !!!\n");
 526
 527
 528   mylong new_size=matS2[0];
 529   cout << "size of data " << new_size << endl;
 530
 531
 532   //first elements that contains the size is removed
 533   uint8_t *reconstucted_data=reinterpret_cast<uint8_t*>(&matS2[1]);
 534   saveFile(reconstucted_data, "lena2.png",new_size);
 535   return 0;
 536 }
 537
 538