// Simulates an n-stage 4 x 4 arbiter PUF with one-hot encoded 24-valued inputs and 24-valued outputs.
// Each 24-digit input value is first mapped into a permutation (p1,p2,p3,p4), pi in {0,1,2,3}
// and then each pi is one-hot encoded into 4 bits. So one 24-valued digit is represented by 16 binary bits.
// Prints challenge-response pairs one per line, separated by space.
// e.g. if the numeber of stages n = 4 prints 
// c1,1...c1,16 c2,1... c2,16 c3,1... c3,16 c4,1... c4,16 r1
// where r1 is the 24-valued response to the one-hot encoded challeenge (c1,1...c1,16 ... c4,1... c4,16)


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <assert.h>

#include <cstdlib>
#include <cmath>
#include <limits>
#include <cstdio>
#include <iostream>

using namespace std;

#define SIZE                     4       // number of stages in APUF
#define M                        4
#define NUMBER_CHALLENGES        10
#define MEAN                     1.0
#define STD_DEV                  0.001


typedef struct sbox {
  float                 o1;                       /* 1st output*/
  float                 o2;                       /* 2nd output */
  float                 o3;                       /* 3rd output */
  float                 o4;                       /* 4th output */
  float                 d11;                      /* delay of wire connecting input 1 and ouptut 1 */
  float                 d12;                         
  float                 d13;                         
  float                 d14; 
  float                 d21;                         
  float                 d22;                         
  float                 d23;                         
  float                 d24; 
  float                 d31;                         
  float                 d32;                         
  float                 d33;                         
  float                 d34; 
  float                 d41;                         
  float                 d42;                         
  float                 d43;                         
  float                 d44;                         
} sbox;

static float map(unsigned a, unsigned b);
double generateGaussianNoise(double mu, double sigma);
static unsigned          perm[SIZE][4];   /* challenge bits as a permutation (a,b,c,d) */
static float             map(unsigned a, unsigned b);
static unsigned          ord(float a, float b, float c,  float d);

/* --------------------------------------------------------------------------- */
int 
main(int argc, char *argv[]) 
{
  unsigned          i, j, k, tmp, output, m, p;
  sbox              puf[SIZE+1];
  unsigned          x[SIZE][4][4]; /* one-hot encoded challenge bits {0,1} */
  unsigned          y[SIZE];     /* challenge bits {0,...,23} */
  unsigned          counter=0;


  int seed=(int)time(NULL);
  srand(seed);

/* Assign delays d1-d4 to each sbox */
  for(i=0; i < SIZE; i++) {
    puf[i].d11 = generateGaussianNoise(MEAN,STD_DEV);
    puf[i].d12 = generateGaussianNoise(MEAN,STD_DEV);   
    puf[i].d13 = generateGaussianNoise(MEAN,STD_DEV);   
    puf[i].d14 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d21 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d22 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d23 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d24 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d31 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d32 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d33 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d34 = generateGaussianNoise(MEAN,STD_DEV);  
    puf[i].d41 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d42 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d43 = generateGaussianNoise(MEAN,STD_DEV); 
    puf[i].d44 = generateGaussianNoise(MEAN,STD_DEV);

  }

  //for(i = 0; i < SIZE; i++) {
    //printf("i = %d:\nd11 = %f, d12 = %f, d13 = %f, d14 = %f\n", i, puf[i].d11, puf[i].d12, puf[i].d13, puf[i].d14);
    //printf("d21 = %f, d22 = %f, d23 = %f, d24 = %f\n", puf[i].d21, puf[i].d22, puf[i].d23, puf[i].d24);
    //printf("d31 = %f, d32 = %f, d33 = %f, d34 = %f\n", puf[i].d31, puf[i].d32, puf[i].d33, puf[i].d34);
    //printf("d41 = %f, d42 = %f, d43 = %f, d44 = %f\n\n", puf[i].d41, puf[i].d42, puf[i].d43, puf[i].d44);
  //}

  for(k = 0; k < NUMBER_CHALLENGES; k++) {
    for(j = 0; j < SIZE; j++) {
      y[j] = rand()%24; 
      //printf("y[%d] = %d\n", j, y[j]);
      map(y[j],j);
      for(m = 0; m < 4; m++) {
        //printf("perm[%d][%d] = %d\n", j, m, perm[j][m]);
        for(p = 0; p < 4; p++) {
          if(p == perm[j][m]) {
            x[j][m][p] = 1;
          }
          else {
             x[j][m][p] = 0;
          }
          printf("%d ", x[j][m][p]);
        }
        //printf(")\n");
      }
    }

    puf[0].o1 = 0;
    puf[0].o2 = 0;
    puf[0].o3 = 0;
    puf[0].o4 = 0;
    //printf("puf[0].o1 = %f, puf[0].o2 = %f, puf[0].o3 = %f, puf[0].o3 = %f\n", puf[0].o1, puf[0].o2, puf[0].o3, puf[0].o4);

    for(i = 0; i < SIZE; i++) {
      puf[i+1].o1 = (puf[i].o1 + puf[i].d11)*x[i][0][0] + (puf[i].o2 + puf[i].d21)*x[i][1][0] + (puf[i].o3 + puf[i].d31)*x[i][2][0] + (puf[i].o4 + puf[i].d41)*x[i][3][0];
      puf[i+1].o2 = (puf[i].o1 + puf[i].d12)*x[i][0][1] + (puf[i].o2 + puf[i].d22)*x[i][1][1] + (puf[i].o3 + puf[i].d32)*x[i][2][1] + (puf[i].o4 + puf[i].d42)*x[i][3][1];
      puf[i+1].o3 = (puf[i].o1 + puf[i].d13)*x[i][0][2] + (puf[i].o2 + puf[i].d23)*x[i][1][2] + (puf[i].o3 + puf[i].d33)*x[i][2][2] + (puf[i].o4 + puf[i].d43)*x[i][3][2];
      puf[i+1].o4 = (puf[i].o1 + puf[i].d14)*x[i][0][3] + (puf[i].o2 + puf[i].d24)*x[i][1][3] + (puf[i].o3 + puf[i].d34)*x[i][2][3] + (puf[i].o4 + puf[i].d44)*x[i][3][3];

      //printf("puf[%d].o1 = %f, puf[%d].o2 = %f, puf[%d].o3 = %f, puf[%d].o4 = %f\n", i+1, puf[i+1].o1, i+1, puf[i+1].o2, i+1, puf[i+1].o3, i+1, puf[i+1].o4);

    }
    //printf("\n");

    output = ord(puf[SIZE].o1, puf[SIZE].o2, puf[SIZE].o3, puf[SIZE].o4);

    printf("%d\n", output);

    //tmp = output;
    //for(j = 0; j < 2; j++) {
      //x[j] = tmp%2; 
      //tmp = tmp/2;
      //printf("%d", x[j]);
    //}
  }
  //printf("\n");

return(0);
}



static float 
map(unsigned a, unsigned b)
{

  switch(a) {
      case 0:
         perm[b][0] = 0;
         perm[b][1] = 1;
         perm[b][2] = 2;
         perm[b][3] = 3;
         break;
      case 1:
         perm[b][0] = 0;
         perm[b][1] = 1;
         perm[b][2] = 3;
         perm[b][3] = 2;
         break;
      case 2:
         perm[b][0] = 0;
         perm[b][1] = 2;
         perm[b][2] = 1;
         perm[b][3] = 3;
         break;
      case 3:
         perm[b][0] = 0;
         perm[b][1] = 2;
         perm[b][2] = 3;
         perm[b][3] = 1;
         break; 
      case 4:
         perm[b][0] = 0;
         perm[b][1] = 3;
         perm[b][2] = 1;
         perm[b][3] = 2;
         break;
      case 5:
         perm[b][0] = 0;
         perm[b][1] = 3;
         perm[b][2] = 2;
         perm[b][3] = 1;
         break;
      case 6:
         perm[b][0] = 1;
         perm[b][1] = 0;
         perm[b][2] = 2;
         perm[b][3] = 3;
         break;
      case 7:
         perm[b][0] = 1;
         perm[b][1] = 0;
         perm[b][2] = 3;
         perm[b][3] = 2;
         break;
      case 8:
         perm[b][0] = 1;
         perm[b][1] = 2;
         perm[b][2] = 0;
         perm[b][3] = 3;
         break;
      case 9:
         perm[b][0] = 1;
         perm[b][1] = 2;
         perm[b][2] = 3;
         perm[b][3] = 0;
         break; 
      case 10:
         perm[b][0] = 1;
         perm[b][1] = 3;
         perm[b][2] = 0;
         perm[b][3] = 2;
         break;
      case 11:
         perm[b][0] = 1;
         perm[b][1] = 3;
         perm[b][2] = 2;
         perm[b][3] = 0;
         break; 
      case 12:
         perm[b][0] = 2;
         perm[b][1] = 0;
         perm[b][2] = 1;
         perm[b][3] = 3;
         break;
      case 13:
         perm[b][0] = 2;
         perm[b][1] = 0;
         perm[b][2] = 3;
         perm[b][3] = 1;
         break;
      case 14:
         perm[b][0] = 2;
         perm[b][1] = 1;
         perm[b][2] = 0;
         perm[b][3] = 3;
         break;
      case 15:
         perm[b][0] = 2;
         perm[b][1] = 1;
         perm[b][2] = 3;
         perm[b][3] = 0;
         break; 
      case 16:
         perm[b][0] = 2;
         perm[b][1] = 3;
         perm[b][2] = 0;
         perm[b][3] = 1;
         break;
      case 17:
         perm[b][0] = 2;
         perm[b][1] = 3;
         perm[b][2] = 1;
         perm[b][3] = 0;
         break;
      case 18:
         perm[b][0] = 3;
         perm[b][1] = 0;
         perm[b][2] = 1;
         perm[b][3] = 2;
         break;
      case 19:
         perm[b][0] = 3;
         perm[b][1] = 0;
         perm[b][2] = 2;
         perm[b][3] = 1;
         break; 
      case 20:
         perm[b][0] = 3;
         perm[b][1] = 1;
         perm[b][2] = 0;
         perm[b][3] = 2;
         break;
      case 21:
         perm[b][0] = 3;
         perm[b][1] = 1;
         perm[b][2] = 2;
         perm[b][3] = 0;
         break;
      case 22:
         perm[b][0] = 3;
         perm[b][1] = 2;
         perm[b][2] = 0;
         perm[b][3] = 1;
         break;
      case 23:
         perm[b][0] = 3;
         perm[b][1] = 2;
         perm[b][2] = 1;
         perm[b][3] = 0;
         break; 
      default :
         printf("Error\n" );
  }
}



static unsigned 
ord(float a, float b, float c,  float d)
{

  if((a < b) && (b < c) && (c < d)) return(0);
  if((a < b) && (b < d) && (d < c)) return(1);
  if((a < c) && (c < b) && (b < d)) return(2);
  if((a < c) && (c < d) && (d < b)) return(3);
  if((a < d) && (d < b) && (b < c)) return(4);
  if((a < d) && (d < c) && (c < b)) return(5);
  if((b < a) && (a < c) && (c < d)) return(6);
  if((b < a) && (a < d) && (d < c)) return(7);
  if((b < c) && (c < a) && (a < d)) return(8);
  if((b < c) && (c < d) && (d < a)) return(9);
  if((b < d) && (d < a) && (a < c)) return(10);
  if((b < d) && (d < c) && (c < a)) return(11);
  if((c < a) && (a < b) && (b < d)) return(12);
  if((c < a) && (a < d) && (d < b)) return(13);
  if((c < b) && (b < a) && (a < d)) return(14);
  if((c < b) && (b < d) && (d < a)) return(15);
  if((c < d) && (d < a) && (a < b)) return(16);
  if((c < d) && (d < b) && (b < a)) return(17);
  if((d < a) && (a < b) && (b < c)) return(18);
  if((d < a) && (a < c) && (c < b)) return(19);
  if((d < b) && (b < a) && (a < c)) return(20);
  if((d < b) && (b < c) && (c < a)) return(21);
  if((d < c) && (c < a) && (a < b)) return(22);
  if((d < c) && (c < b) && (b < a)) return(23);
    
  return(0);
} 

 
// code from  https://en.wikipedia.org/wiki/Box-Muller_transform
double generateGaussianNoise(double mu, double sigma)
{
  static const double epsilon = std::numeric_limits<double>::min();
  static const double two_pi = 2.0*3.14159265358979323846;

  static double z1;
  static bool generate=false;
  generate = !generate;

  if (!generate)
     return z1 * sigma + mu;

  double u1, u2;
  do
   {
     u1 = rand() * (1.0 / RAND_MAX);
     u2 = rand() * (1.0 / RAND_MAX);
     //cout << "u1:" << u1 << endl;
     //cout << "u2:" << u2 << endl;
   }
  while ( u1 <= epsilon );

  double z0;
  z0 = sqrt(-2.0 * log(u1)) * cos(two_pi * u2);
  z1 = sqrt(-2.0 * log(u1)) * sin(two_pi * u2);
  return z0 * sigma + mu;
}