//////////////////////////////////////////////////////////////////////////////
//
//          Shamir's Secret Sharing demonstration (version 2008/11/10)
//
// Shamir A.,
// How to Share a Secret,
// Communications of the ACM, 22, 1979, pp. 612--613.
//
//                  Copyright (C) 2008 Christophe DAVID
//        http://www.christophedavid.org/w/c/w.php/Files/Security
//
//  This program is free software; you  can redistribute it and/or modify it
//  under the terms of the GNU General Public License version 2 as published
//  by the Free Software Foundation (http://www.gnu.org/licenses/gpl.html).
//
//////////////////////////////////////////////////////////////////////////////

//                This program was last tested with Scilab 4.1.2
//                            http://www.scilab.org

//////////////////////////////////////////////////////////////////////////////

// start reading at the very end of this file.

//////////////////////////////////////////////////////////////////////////////


errclear      // clear all errors
clear;        // kill variables
clearglobal   // kill global variables
clc;          // clear command window
tohome;       // move the cursor to the upper left corner of the Command Window
clf;          // clear or reset the current graphic figure (window) to default values
clear_pixmap; // erase the pixmap buffer


//////////////////////////////////////////////////////////////////////////////

function [ReturnValue] = FindMultipeBytesFmShares(SharesMatrix)
   // disp(SharesMatrix);
   ReturnValue = [];

   SizeOfSharesMatrix = size(SharesMatrix);

   X = SharesMatrix(:,1);
   // disp(X);

   for i = 2 : 1 : SizeOfSharesMatrix(1, 2);
      tmp(:, 1) =  X;  
      tmp(:, 2) =  SharesMatrix(:,i);
      //disp(tmp);
      ReturnValue(1, i - 1) = FindSingleByteFromShares(tmp);          
   end
endfunction

//////////////////////////////////////////////////////////////////////////////

function [ReturnValue] = FindSingleByteFromShares(SharesMatrix)
   //disp(SharesMatrix);

   SizeOfSharesMatrix = size(SharesMatrix);
   NumberOfShares     = SizeOfSharesMatrix(1, 1);

   CoefficientsMatrix = [];
   for i = 1 : 1 : NumberOfShares
      temp = 1;
      for j = 1 : 1 : NumberOfShares
         if (i <> j)
            temp = pmodulo(  (  (- temp)..
                             * SharesMatrix(j, 1)..
                             * ModuloPrimeInverse(  SharesMatrix(i, 1)..
                                                  - SharesMatrix(j, 1)..
                                                 )..
                            ),..
                          Prime..
                         );
         end
      end
      CoefficientsMatrix(1, i) = temp;
   end
   // disp(CoefficientsMatrix);

   ResultsMatrix = SharesMatrix(1:$, 2:2);
   //disp(ResultsMatrix);  

   ReturnValue = 0;
   for i = 1 : 1 : NumberOfShares
      ReturnValue = pmodulo(  (  ReturnValue..
                               + (   SharesMatrix(i, 2)..
                                   * CoefficientsMatrix(1, i)..
                                 )..
                             ),..
                           Prime..
                          );
   end
   //disp(ReturnValue);
endfunction

//////////////////////////////////////////////////////////////////////////////

function [ReturnValue] = ModuloPrimeInverse(InputValue)
   if InputValue < 0
      // ReturnValue = pmodulo((InversesMatrix(258 + InputValue)), Prime);      
      ReturnValue = pmodulo(InversesMatrix(258 + InputValue), Prime);
   else
      ReturnValue = InversesMatrix(InputValue + 1);
   end
   // printf("%d %d\n", InputValue, ReturnValue);
endfunction

//////////////////////////////////////////////////////////////////////////////

function [ReturnValue] = PolynomialFunction(CoefficientsMatrix, XMatrix)
   ReturnValue = 0;

   SizeOfCoefficientsMatrix = size(CoefficientsMatrix);
   SizeOfXMatrix            = size(XMatrix);
   NumberOfX                = SizeOfXMatrix(1,1);
   NumberOfCoefficients     = SizeOfCoefficientsMatrix(1, 2);

   for i = 1 : 1 : NumberOfX
      SubTotal = 0;
      for j = 1 : 1 : NumberOfCoefficients
        SubTotal = pmodulo( ( i  * XMatrix(i, 1)) +  CoefficientsMatrix(1, j),  Prime);
      end
   ReturnValue(i, 1) = SubTotal;
   end
endfunction

//////////////////////////////////////////////////////////////////////////////

function [ReturnValue] = PopulateInversesMatrix()
   ReturnValue = 0;

   x = 1;
   y = 1;
   ReturnValue(1) = 0;

   for i = 2 : 1 : Prime
      ReturnValue(x + 1) = y;
      x = pmodulo( 3 * x, Prime); //  3 is a generator of the ring Z/Z257
      y = pmodulo(86 * y, Prime); // 86 is the inverse of 3 modulo 257      
   end
   // disp(ReturnValue);
endfunction

//   0   1 129  86 193 103  43 147 225 200
// 180 187 150 178 202 120 241 121 100 230
//  90  49 222 190  75  72  89 238 101 195
//  60 199 249 148 189 235  50 132 115 145
//  45 163 153   6 111  40  95 175 166  21
//  36 126 173  97 119 243 179 248 226  61
//  30  59 228 102 253  87  74 234 223 149
// 246 181  25 169  66  24 186 247 201 244
// 151 165 210  96 205 127   3  65 184  26
//  20 209 176 152 216  46  83  53 139 135
//  18  28  63   5 215 164 177 245 188 224
// 250  44 218 116 124  38 113 134 159  54
//  15  17 158 140 114 220  51  85 255   2
// 172 206  37 143 117  99 240 242 203  98
// 123 144 219 133 141  39 213   7  33  69
//  12  80  93  42 252 194 229 239 122 118
// 204 174 211  41 105  81  48 237 231  73
// 192 254 130  52 161  47  92 106  13  56
//  10  71 233 191  88 232  76  11 108  34
//  23 183 170   4 155  29 198 227 196  31
//   9  78  14 138 160  84 131 221 236  91
//  82 162 217 146 251 104  94 212 112 142
// 125 207  22  68 109   8  58 197  62 156
//  19 168 185 182  67  35 208 167  27 157
// 136  16 137  55  79 107  70  77  57  32
// 110 214 154  64 171 128 256

//////////////////////////////////////////////////////////////////////////////

function [ReturnValue] = ShareMultipleBytes(MM, k, n)
   ReturnValue = (1:n);
   ReturnValue = ReturnValue';
   // disp(ReturnValue);  

   col = 2;
   for M = MM
      SingleByteShares = ShareSingleByte(M, k, n);
      //disp(SingleByteShares);
      for i = 1 : 1 : n
          ReturnValue(i, col) = SingleByteShares(i, 2);
      end
   col = col + 1;
   end
// disp(ReturnValue);  
endfunction


//////////////////////////////////////////////////////////////////////////////

function [ReturnValue] = ShareSingleByte(M, k, n)
   ReturnValue = 0;

   Coefficients       = rand(1, k - 1);
   Coefficients       = floor(Coefficients * Prime);

   Coefficients(1, k) = M;
   // disp(Coefficients);

   for i = 1 : 1 : n
      ReturnValue(i, 1) = i;
      ReturnValue(i, 2) = pmodulo(PolynomialFunction(Coefficients, i), Prime);
   end

   // disp(ReturnValue);

   if 1
      Polynomial         = "y = pmodulo((0";
      for i = 1 : 1 : k
         Polynomial      =   Polynomial + ' + ('  + string(Coefficients(1, i)) + '*(x**' + string(k-i) +  + '))';
      end
      Polynomial =   Polynomial + '), Prime)';
      // disp(Polynomial);
      deff("[y]=f(x)", Polynomial)
      fplot2d(0:Prime, f);
   end
endfunction

//////////////////////////////////////////////////////////////////////////////

function [ReturnValue] = WriteShares(SharesMatrix, k ,n)
   // disp(SharesMatrix);
   ReturnValue = [];

   SizeOfSharesMatrix = size(SharesMatrix);

   for row = 1 : 1 : SizeOfSharesMatrix(1, 1);
      ReturnValue(row) = sprintf("%02X%02X", SharesMatrix(row, 1), k);
      for col = 1 : 1 : SizeOfSharesMatrix(1, 2);
         if SharesMatrix(row, col) < 256
            PseudoHex = sprintf("%02X", SharesMatrix(row, col))
         else
            PseudoHex = 'G0'
         end
         ReturnValue(row) = ReturnValue(row) + PseudoHex;
      end
   end
endfunction


//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////

Prime          = 257;
InversesMatrix = PopulateInversesMatrix();


// Test PolynomialFunction
if 0
   X = [10; 20; 30];
   Coefficients = [2, 3, 4];
   test = PolynomialFunction(Coefficients, X);
end


// test ShareSingleByte
if 0
   M =  floor(rand() * 255);
   k =   3;
   n = 222;

   Shares = ShareSingleByte(M, k, n);
   //disp(Shares);

   if 0
      printf("\n");          
      for i = 1 : 1 : n
         printf("(:%d:%d:%d:)\n", k, Shares(i, 1), Shares(i, 2));
      end
   end  

   // solve first k shares
   // printf("\nM=%d, result=%d\n", M, FindSingleByteFromShares(Shares(1:k, 1:2)));          

   // solve last k shares
   printf("\nM=%d, result=%d\n", M, FindSingleByteFromShares(Shares(n-k+1:n, 1:2)));
end


// test ShareMultipleBytes
if 0
   k = 3;
   n = 5;
   SecretBytes = [3, 5, 7, 9, 11, 13, 15, 127, 256];
   Shares = ShareMultipleBytes(SecretBytes, k, n);
   // disp(Shares);
   recovered = FindMultipeBytesFmShares(Shares);
   disp(SecretBytes);
   disp(recovered);
end

// test ShareMultipleBytes - string
if 0
   k =   3;
   n = 255;
   SecretString = 'the quick brown fox jumps over the lazy dog';  
   disp(SecretString);

   SecretBytes = ascii(SecretString);
   Shares = ShareMultipleBytes(SecretBytes, k, n);
   // disp(Shares);

   if 0
      Txt = WriteShares(Shares);

      if 0
         Txt = [
                'E303E37CA2CE28DBF1';
                'F703F75027A5370B55';
                'D903D9B96FB75412B0';
               ];

         Shares = [];      
      end

      if 1
         FileName = '~/test/shares.txt';
         mdelete(FileName);
         f = file('open', FileName,'new');
         for i = 1 : 1 : n
            fprintf(f, "%s\n", Txt(i));
         end
         file('close', f);
      end


      SizeOfTxt = size(Txt);      
      for i = 1 : 1 : SizeOfTxt(1, 1);
         // printf("%s %d\n", txt(i), length(txt(1)));
         printf("%s\n", Txt(i));          
         A = strsplit(Txt(i), (2:2:length(Txt(i)) - 2));
         //  disp(A);

         SizeOfA = size(A);
         for j = 3 : 1 : SizeOfA(1, 1);          
            //disp(sscanf(A(j, 1), "%X"));        
            if A(j, 1) == 'G0'
               Shares(i, j - 2) = 256;
            else
               Shares(i, j - 2) = sscanf(A(j, 1), "%02X");
            end
         end
      end
   end

   //disp(Shares);


   recovered = ascii(FindMultipeBytesFmShares(Shares(1:k, 1:$)));
   disp(recovered);
end

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////


// global test

if 1
   k =  3;                                               // number of shares required to reconstruct the secret
   n =  5;                                               // total number of shares
   SecretString = 'My secret';              
   disp(SecretString);                                   // the secret
   SecretBytes  = ascii(SecretString);
   Shares       = ShareMultipleBytes(SecretBytes, k, n);
   disp(WriteShares(Shares));                            // the shares    
   reconstructed   = ascii(FindMultipeBytesFmShares(Shares(1:k, 1:$))); // the secret, reconstructed
   disp(reconstructed);
end


//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////