/*** Boomerang attack on COCONUT98 and 16-round Khufu implemented by Roman Dovgard ***/

#include<stdio.h>
#include<string.h>
#include<stdlib.h>

#define UINT                          unsigned short      // 16-bit blocks.  
#define WORD                          unsigned long       // 32-bit blocks. 
#define DWORD                         unsigned long long  // 64-bit blocks.
#define POWER32MINUS1                 4294967295U         // 2^32-1.
#define POWER31                       2147483648U         // 2^31.
#define NEWLINE                       printf("\n")        // New line printing.
#define SPACE                         printf(" ")         // Space printing.
#define RANGE                         500000              // Radius in which we search for the keys.
#define KHUFU_RIGHT_PAIRS             1024                // Number of right pairs we need to attack Khufu's S box.
#define KHUFU_POOL_SIZE               80000               // Pool of plaintexts/ciphertexts we use in Khufu ~> 2^16.
#define KHUFU_ITERS                   8                   // Number of a's we use in Khufu.
#define MAX_POOL_SIZE                 200000              // Maximal pool size in birthday paradox.
#define MAX_POLY_SIZE                 200                 // Size of polynomials used in decorrelation module.
#define DEBUG_                                            // Whether we debug or run the program.
#define SEARCH_RANGE                                      // Whether we have debugging or real version of key recovery.
#define DECORRELATION_MODULE                              // Whether we have decorrelation module in COCONUT98.
#define PRINT_                                            // Whether we print debugging information.

#ifdef DECORRELATION_MODULE
#define COCONUT98_CHARS               8                   // Number of characteristics used for k1 and k8 recovery.
#define COCONUT98_CHARS_SECOND        6                   // Number of characteristics used for k2 recovery.
#define COCONUT98_CHARS_THIRD         6                   // Number of characteristics used for k3 recovery.
#define COCONUT98_RIGHT_PAIRS         20                  // Number of right pairs used for k1 and k8 recovery.
#define COCONUT98_RIGHT_PAIRS_SECOND  8                   // Number of right pairs used for k2 recovery.
#define COCONUT98_RIGHT_PAIRS_THIRD   8                   // Number of right pairs used for k3 recovery.
#else
#define COCONUT98_CHARS               7                   
#define COCONUT98_CHARS_SECOND        7                   
#define COCONUT98_CHARS_THIRD         6
#define COCONUT98_RIGHT_PAIRS         10 
#define COCONUT98_RIGHT_PAIRS_SECOND  8  
#define COCONUT98_RIGHT_PAIRS_THIRD   8  
#endif


// Two auxiliary functions.

WORD min(WORD a, WORD b) { return (a<=b)?a:b; }
WORD max(WORD a, WORD b) { return (a>=b)?a:b; }

// Structure that keeps 64-bit block in comfortable way for encryption/decryption.

struct PAIR
{
  WORD x, y;
 
  PAIR() {}
  PAIR(WORD _x, WORD _y) : x(_x), y(_y) {} 

  PAIR operator+ (PAIR& pair) { return PAIR(x^pair.x, y^pair.y); }
};

// COCONUT98's constant c.

WORD coconut98_c;

// COCONUT98's S box.

WORD coconut98_S[256];

// COCONUT98's keys.

WORD coconut98_k1, coconut98_k2, coconut98_k3, coconut98_k4,
  coconut98_k5, coconut98_k6, coconut98_k7, coconut98_k8, 
  coconut98_K1, coconut98_K2, coconut98_K3, coconut98_K4,
  coconut98_K5, coconut98_K6, coconut98_K7, coconut98_K8;

// COCONUT98's fast decorrealtion module support.

PAIR coconut98_K5_K6_pair, coconut98_K7_K8_pair, coconut98_K5_K6_pair_inv, coconut98_K7_K8_pair_inv;

// Khufu's S boxes.

WORD khufu_S1[256], khufu_S2[256];

// Random 32-bit block generation.

inline void random_word(WORD& text);

// Random 64-bit block generation.

inline void random_dword(DWORD& text);

// Two random 32-bit blocks generation.

inline void random_pair(PAIR& pair);

// 32-bit block printing.

inline void print_word(WORD text);

// 64-bit block printing.

inline void print_dword(DWORD text);

// Two 32-bit blocks printing.

inline void print_pair(PAIR pair);

// COCONUT98's initialization.

void coconut98_random_init();

// COCONUT98's printing configuration.

void coconut98_print_config();

// COCONUT98's F function.

inline WORD coconut98_F(WORD text);

// COCONUT98's one round of encryption.

inline PAIR coconut98_encryption_round(PAIR plaintext, WORD key);

// COCONUT98's one round of decryption. 

inline PAIR coconut98_decryption_round(PAIR ciphertext, WORD key);

// COCONUT98's probability of characterstic for F.

double coconut98_calc_approx_prob_of_characteristic_for_F(WORD delta, WORD nabla, UINT trials = 10000);

// COCONUT98's Psi function.  

inline PAIR coconut98_Psi(PAIR& plaintext, WORD& k1, WORD& k2, WORD& k3, WORD& k4);

// Polynomial null.

inline void poly_null(UINT a[]);

// Polynomial printing.

inline void poly_print(UINT a[]);

// Polynomial one.

inline void poly_one(UINT a[]);

// Polynomial copy.

inline void poly_copy(UINT src[], UINT dst[]);

// Polynomial degree.

inline long poly_deg(UINT a[]);

// Polynomial addition.

inline void poly_add(UINT a[], UINT b[], UINT c[]);

// Polynomial multiplication.

inline void poly_mul(UINT a[], UINT b[], UINT c[]);

// Polynomial shift.

inline void poly_shift(UINT a[], WORD shift);

// Polynomial division with remainder.

inline long poly_div(UINT a[], UINT b[], UINT q[], UINT r[]);

// Polynomial extended euclidean algorithm.

inline void poly_ext_euclidean(UINT a[], UINT b[], UINT d[], UINT x[], UINT y[]);

// Convertion from DWORD to POLY.

inline void dword_to_poly(DWORD& dword, UINT poly[]);

// Convertion from POLY to DWORD.

inline void poly_to_dword(UINT poly[], DWORD& dword);

// Convertion from DWORD to PAIR.

inline void dword_to_pair(DWORD& dword, PAIR& pair);

// Convertion from PAIR to DWORD.

inline void pair_to_dword(PAIR& pair, DWORD& dword);

// Convertion from POLY to PAIR.

inline void poly_to_pair(UINT poly[], PAIR& pair);

// Convertion from PAIR to POLY.

inline void pair_to_poly(PAIR& pair, UINT poly[]);

// Pair multiplication using polynomials.

inline PAIR pair_mul(PAIR pair_a, PAIR pair_b);

// Pair fast multiplication using xtime.

inline PAIR pair_fast_mul(PAIR pair_a, PAIR pair_b);

// Pair inverse using polynomials

inline PAIR pair_inv(PAIR pair_a);

// Checking polynomial interface.

void poly_check();

// COCONUT98's encryption function.

inline PAIR coconut98_encrypt(PAIR& plaintext, WORD& k1, WORD& k2, WORD& k3, WORD& k4, WORD& k5,
			      WORD& k6, WORD& k7, WORD& k8, WORD& K5, WORD& K6, WORD& K7, WORD& K8);
// COCONUT98's decryption function.

inline PAIR coconut98_decrypt(PAIR& ciphertext, WORD& k1, WORD& k2, WORD& k3, WORD& k4, WORD& k5,
			      WORD& k6, WORD& k7, WORD& k8, WORD& K5, WORD& K6, WORD& K7, WORD& K8);

// COCONUT98's attacker encryption trial.

inline PAIR coconut98_encryption_trial(PAIR& plaintext);

// COCONUT98's attacker decryption trial.

inline PAIR coconut98_decryption_trial(PAIR& ciphertext);

// COCONUT98's attacker encryption trial on rounds 2 3 4 5 6 7 8,
// while k1 of the cipher is given via argument key to this function. 

PAIR coconut98_encryption_trial_2345678(PAIR& plaintext, WORD key);

// COCONUT98's attacker decryption trial on rounds 2 3 4 5 6 7 8,
// while k1 of the cipher is given via argument key to this function. 

PAIR coconut98_decryption_trial_2345678(PAIR& ciphertext, WORD key);

// COCONUT98's attacker encryption trial on rounds 3 4 5 6 7 8, while k1 and k2 of
// the cipher are given via arguments key1 and key2 respectively to this function. 

PAIR coconut98_encryption_trial_345678(PAIR& plaintext, WORD key1, WORD key2);

// COCONUT98's attacker decryption trial on rounds 3 4 5 6 7 8, while k1 and k2 of
// the cipher are given via arguments key1 and key2 respectively to this function. 

PAIR coconut98_decryption_trial_345678(PAIR& ciphertext, WORD key1, WORD key2);

// COCONUT98's decryption checking.

void coconut98_check_inverse();

// COCONUT98's boomerang attack.

void coconut98_boomerang_attack();

// Khufu's initialization.

void khufu_random_init();

// Khufu's printing configuration.

void khufu_print_config();

// Khufu's one round of encryption.

inline PAIR khufu_encryption_round(PAIR plaintext, WORD S[]);

// Khufu's one round of decryption. 

inline PAIR khufu_decryption_round(PAIR ciphertext, WORD S[]);

// Khufu's encryption function.

inline PAIR khufu_encrypt(PAIR& plaintext, WORD S1[], WORD S2[], WORD& k1, WORD& k2, WORD& k3, WORD& k4);

// Khufu's decryption function.

inline PAIR khufu_decrypt(PAIR& ciphertext, WORD S1[], WORD S2[], WORD& k1, WORD& k2, WORD& k3, WORD& k4);

// Khufu's attacker encryption trial.

inline PAIR khufu_encryption_trial(PAIR& plaintext);

// Khufu's attacker decryption trial.

inline PAIR khufu_decryption_trial(PAIR& ciphertext);

// Khufu's decryption checking.

void khufu_check_inverse();

// Partition procedure for Quicksort algorihtm.

WORD Partition(PAIR arr[], PAIR snd_arr[], WORD p, WORD r); 

// Quicksort for boomerang attack on Khufu.

void Quicksort(PAIR arr[], PAIR snd_arr[], WORD p, WORD r);

// Binary search for boomerang attack on Khufu.

long BinarySearch(PAIR arr[], WORD val, WORD p, WORD r);

// Second partition procedure for second Quicksort algorihtm.

WORD Partition(PAIR arr[], PAIR snd_arr[], WORD p, WORD r); 

// Second Quicksort for boomerang attack on Khufu.

void Quicksort(PAIR arr[], PAIR snd_arr[], WORD p, WORD r);

// Second binary search for boomerang attack on Khufu.

long BinarySearch(PAIR arr[], WORD val, WORD p, WORD r);

// Simple partition procedure for simple Quicksort algorihtm.

WORD PartitionSimple(WORD arr[], WORD p, WORD r);

// Simple Quicksort for birthday paradox.

void QuicksortSimple(WORD arr[], WORD p, WORD r);

// Simple binary search for birthday paradox. 

long BinarySearchSimple(WORD arr[], WORD val, WORD p, WORD r);

// Calculation of probability of birthday paradox.

double birthday_paradox_prob_calc(WORD pool, WORD trials);

// Khufu's boomerang attack.

void khufu_boomerang_attack();

// Khufu's right quartet probability calculation. 

double khufu_right_quartet_prob_calc(WORD trials);

int main(int argc, char **argv)
{
#ifdef DEBUG

  srand48(atoi(argv[1]));
  srandom(atoi(argv[1]));
  
  coconut98_random_init();

  //khufu_random_init();

  //coconut98_check_inverse();

  //poly_check();

  //double res = khufu_right_quartet_prob_calc(100);
  //printf("Probability of right quartet in Khufu is %f\n\n", res);

  int i;
  double arr[32];
  WORD delta, nabla;

  for(i = 0; i<32; i++)
    {
      delta = 1<<i;
      if(i + 11 <= 31) nabla = 1<<(i+11);
      else nabla = 1<<(i-21);      
      double res = coconut98_calc_approx_prob_of_characteristic_for_F(delta, nabla, 10000); //100
      arr[i] = res;
    }

  NEWLINE;
  for(i = 0; i<32; i++)
    {
      printf("%f ", arr[i]);
      if(i == 9 || i == 19 || i == 29) NEWLINE; 
    }
  NEWLINE;

  double average = 0, without_average = 0;
  
  for(i = 0; i<32; i++)
    {
      int t, j = (i + 11 <= 31)? i+11 : i-21;
      long double c = 1;
      for(t = 0; t < 8; t++) c *= arr[i];
      for(t = 0; t < 4; t++) c *= arr[j];
      if(c == 0) continue;
      c = 1;
      for(t = 0; t < 8; t++) c *= arr[i];
      //for(t = 0; t < 4; t++) c *= arr[j];
      c = 1 / c;
      if(i != 19 && i != 30) without_average += c;
      average += c;
      long long_c = (long)(c);
      //printf("c[%d] is %ld\n", i, long_c);
      printf("%ld & ", long_c);
    }

  printf("%ld\n", (long)(average/8.0));
  //printf("average without 19 and 30 is %ld\n", (long)(without_average/6.0));
  
  NEWLINE;
  NEWLINE;
  
  
#else
  
  NEWLINE;
  
  if(argc != 3)
    {
      printf("Usage: boomerang COCONUT98/Khufu/Birthday seed\n\n");
      return 0;
    }
  
  int seed = atoi(argv[2]);
  
  srand48(seed);
  srandom(seed);

  if(!strcmp(argv[1],"COCONUT98"))
    {
      coconut98_random_init();
      coconut98_print_config();
      NEWLINE;
      coconut98_boomerang_attack();
      NEWLINE;
      NEWLINE;
      return 0;
    }


  if(!strcmp(argv[1],"Khufu"))
    {
      khufu_random_init();
      khufu_print_config();
      NEWLINE;
      khufu_boomerang_attack();
      NEWLINE;
      NEWLINE;
      return 0;
    }

  if(!strcmp(argv[1],"Birthday"))
  {
    double res = birthday_paradox_prob_calc(KHUFU_POOL_SIZE, seed); // 65536
    printf("We calculate here probability of birthday paradox as follows.\n"
           "In each trial we generate two pools of size %d in domain of size 2^32.\n"
           "Counting in how many out of %d trials we have no intersection "
	   "between the\ntwo pools we obtained estimate which says "
	   "that probability of no collision\nis %f.\n\n", KHUFU_POOL_SIZE, seed, res);
    // with 65536 Probability of no collision is 0.349000 (1000 trials).
    // with KHUFU_POOL_SIZE Probability of no collision is 0.219000 (1000 trials).
    return 0;
  }

  printf("Usage: boomerang COCONUT98/Khufu/Birthday seed\n\n");
  
#endif
  
  return 0;
}

/* The following function preforms random word generation

   Input/Output: text
*/ 

void random_word(WORD& text)
{
  // Function lrand48 generates uniformly ditributed random number in the range [0,2^31-1].
  // Adding 0 or 2^31 randomly (and uniformly) to lrand48() we get uniformly ditributed number in
  // the range [0,2^32-1].

  text = (WORD)(lrand48()) + (random()&01) * POWER31; 
}

/* The following function performs random double word generation

   Input/Output: text
*/

void random_dword(DWORD& text)
{
  WORD high, low;
  random_word(low);
  random_word(high);
  text = high;
  text <<= 32;
  text += low;
}

/* The following function performs random pair generation

   Input/Output: pair
*/

void random_pair(PAIR& pair)
{
  random_word(pair.x);
  random_word(pair.y);
}

/* The following function prints a word

   Input: text
*/

void print_word(WORD text)
{
  // We print 32-bit numbers in hexadecimal representation.

  printf("%04X%04X", (UINT)(text>>16), (UINT)(text&0x0000FFFF));
}

/* The following function prints a double word

   Input: text
*/

void print_dword(DWORD text)
{
  // We print 64-bit numbers in hexadecimal representation.

  printf ("%04X%04X%04X%04X", (UINT)(text>>48), (UINT)((text>>32)&0x0000FFFF),
	  (UINT)((text>>16)&0x00000000FFFF), (UINT)(text&0x000000000000FFFF));
}

/* The following function prints a pair

   Input: pair
*/

void print_pair(PAIR pair)
{
  printf("(");
  print_word(pair.x);
  printf(", ");
  print_word(pair.y);
  printf(")");
}

/* The following function randomly initializes COCONUT98 constants, S boxes and keys */

void coconut98_random_init()
{
  // Initializing COCONUT98's constant, S box and keys with random values. 
  
  coconut98_c = 0xB7E15162;

  for(WORD i=0; i<256; i++) coconut98_S[i] = lrand48()%16777216; // 16777216 = 2^24

  // COCONUT98's key schedule algorithm.

  random_word(coconut98_K1);
  random_word(coconut98_K2);
  random_word(coconut98_K3);
  random_word(coconut98_K4);
  random_word(coconut98_K5);
  random_word(coconut98_K6);
  random_word(coconut98_K7);
  random_word(coconut98_K8);
  
  coconut98_k1 = coconut98_K1;
  coconut98_k2 = coconut98_K1^coconut98_K3;
  coconut98_k3 = coconut98_K1^coconut98_K3^coconut98_K4;
  coconut98_k4 = coconut98_K1^coconut98_K4;
  coconut98_k5 = coconut98_K2;
  coconut98_k6 = coconut98_K2^coconut98_K3;
  coconut98_k7 = coconut98_K2^coconut98_K3^coconut98_K4;
  coconut98_k8 = coconut98_K2^coconut98_K4;

  // Fast decorrelation module support.
  
#ifdef DECORRELATION_MODULE
  coconut98_K5_K6_pair = PAIR(coconut98_K5, coconut98_K6);
  coconut98_K7_K8_pair = PAIR(coconut98_K7, coconut98_K8);
  coconut98_K5_K6_pair_inv = pair_fast_mul(coconut98_K5_K6_pair,coconut98_K7_K8_pair);
  coconut98_K7_K8_pair_inv = pair_inv(coconut98_K7_K8_pair); 
#endif
}

/* The following function prints COCONUT98's constants, S boxes and keys */

void coconut98_print_config()
{
  // Printing COCONUT98's configuration

  printf("COCONUT98's configuration is as follows:\n");
  printf("S box:");
  for(WORD i=0; i<256; i++)
    {
      if(!(i%16)) NEWLINE;
      else SPACE;
      printf("%02X%04X", (UINT)(coconut98_S[i]>>16), (UINT)(coconut98_S[i]&0x0000FFFF));
    }
  NEWLINE;

  printf("K1: ");
  print_word(coconut98_K1);
  printf(" K2: ");
  print_word(coconut98_K2);
  printf(" K3: ");
  print_word(coconut98_K3);
  printf(" K4: ");
  print_word(coconut98_K4);
  printf(" K5: ");
  print_word(coconut98_K5);
  printf(" K6: ");
  print_word(coconut98_K6);
  printf(" K7: ");
  print_word(coconut98_K7);
  printf(" K8: ");
  print_word(coconut98_K8);
  
  NEWLINE;
  
  printf("k1: ");
  print_word(coconut98_k1);
  printf(" k2: ");
  print_word(coconut98_k2);
  printf(" k3: ");
  print_word(coconut98_k3);
  printf(" k4: ");
  print_word(coconut98_k4);
  printf(" k5: ");
  print_word(coconut98_k5);
  printf(" k6: ");
  print_word(coconut98_k6);
  printf(" k7: ");
  print_word(coconut98_k7);
  printf(" k8: ");
  print_word(coconut98_k8);

  NEWLINE;
}

/* The following function implements COCONUT98's F function

   Input: text (argument to the F function)

   Output: F(text)
*/

WORD coconut98_F(WORD text)
{
  // Applying first psi.

  DWORD tmp = text+256*coconut98_S[text%256];
  tmp &= 0x00000000FFFFFFFF;

  // Left rotation by 11 bits. 

  tmp = (tmp << 11)^(tmp >> 21);
  tmp &= 0x00000000FFFFFFFF;

  // Adding coconut98_c.

  tmp += coconut98_c;
  tmp &= 0x00000000FFFFFFFF;

  // Applying second psi.

  tmp += 256*coconut98_S[tmp%256];
  
  return tmp&0x00000000FFFFFFFF;
}

/* The following function implements one encryption round of COCONUT98

   Inputs: plaintext - plaintext before the encryption
           key - key of this round

   Output: ciphertext
*/

PAIR coconut98_encryption_round(PAIR plaintext, WORD key)
{
  return PAIR(plaintext.y, plaintext.x^coconut98_F(plaintext.y^key));
}

/* The following function implements one decryption round of COCONUT98

   Inputs: ciphertext - ciphertext before the decryption
           key - key of this round

   Output: plaintext
*/

PAIR coconut98_decryption_round(PAIR ciphertext, WORD key)
{
  return PAIR(ciphertext.y^coconut98_F(ciphertext.x^key), ciphertext.x);
}

/* The following function calculates probabilities of characteristics for the F function of COCONUT98

   Inputs: delta - input difference
           nabla - output difference
	   trials - number of trials in the calculation

   Ouput: probability of the characteristic 	   
*/

double coconut98_calc_approx_prob_of_characteristic_for_F(WORD delta, WORD nabla, UINT trials)
{
  // We calculate here approximately probability of charactersitic delta-->nabla for function F.
  WORD counter = 0;
  for(WORD i=0; i<trials; i++)
    {
      WORD P;
      random_word(P);
      WORD P_tag = P^delta;
      WORD output_difference = coconut98_F(P)^coconut98_F(P_tag);
      if(output_difference == nabla) counter++;
    }
  return (double)(counter) / (double)(trials);
}

/* The following function implements the Psi function of COCONUT98

   Inputs: plaintext - plaintext before the encryption
           k1,k2,k3,k4 - keys used in first, second, third and fourth rounds respectively

   Output: ciphertext
*/   

PAIR coconut98_Psi(PAIR& plaintext, WORD& k1, WORD& k2, WORD& k3, WORD& k4)
{
  // We implement here four rounds of COCONUT98. Instead of doing

  WORD out2 = plaintext.x^coconut98_F(plaintext.y^k1); 
  WORD out3 = plaintext.y^coconut98_F(out2^k2);
  WORD out4 = out2^coconut98_F(out3^k3);
  WORD out5 = out3^coconut98_F(out4^k4);
  return PAIR(out5, out4);
}

/* The following function generates null (zero) polynomial

   Input/Ouput: a - polynomial (represented as an array)
*/

void poly_null(UINT a[])
{
  for(WORD i = 0; i < MAX_POLY_SIZE; i++) a[i] = 0; 
}

/* The following function prints polynomial

   Input: a - polynomial
*/

void poly_print(UINT a[])
{
  for(long i = MAX_POLY_SIZE-1; i >= 0; i--)
    {
      printf("%u", a[i]);
      if(i%100 == 0) NEWLINE;
    }
}

/* The following function generates one (1) polynomial

   Input/Output: a - polynomial
*/

void poly_one(UINT a[])
{
  poly_null(a);
  a[0] = 1;
}

/* The following function performs copying operation between polynomials

   Input: src - source polynomial

   Input/Output: dst - destination polynomial
*/

void poly_copy(UINT src[], UINT dst[])
{
  for(WORD i = 0; i < MAX_POLY_SIZE; i++) dst[i] = src[i]; 
}

/* The following function calculates degree of a polynomial

   Input: a - polynomial

   Output: degree of a
*/

long poly_deg(UINT a[])
{
  for(long i = MAX_POLY_SIZE-1; i >= 0; i--)
    if(a[i]) return i;
  return -1;
}

/* The following function performs addition operation between polynomials

   Inputs: a, b - polynomials

   Input/Output: c = a + b - result of the addition
*/

void poly_add(UINT a[], UINT b[], UINT c[])
{
  for(WORD i = 0; i < MAX_POLY_SIZE; i++) c[i] = a[i]^b[i];
}

/* The following function performs multiplication operation between polynomials

   Inputs: a, b - polynomials

   Input/Output: c = a * b - result of the multiplication
*/

void poly_mul(UINT a[], UINT b[], UINT c[])
{
  poly_null(c);
  for(WORD i = 0; i < MAX_POLY_SIZE; i++)
      for(WORD j = 0; j <= i; j++) c[i] ^= (a[j]&b[i-j]);
}

/* The following function performs shift of a polynomial

   Input/Output: a - polynomial

   Input: shift - some integer shift
*/ 	    

void poly_shift(UINT a[], WORD shift)
{
  //printf("degree of a before %d ", poly_deg(a));
  UINT aa[MAX_POLY_SIZE], elementary[MAX_POLY_SIZE];
  poly_null(elementary);
  elementary[shift] = 1;
  poly_copy(a, aa);
  poly_mul(aa, elementary, a);
  //printf("degree of a after %d\n", poly_deg(a));
}

/* The following function performs division with remainder operation between polynomials

   Inputs: a - polynomial we want to divide
           b - polynomial by which we want to divide

   Inputs/Outputs: q, r - s.t. a = q * b + r and deg(r) < deg(b)
*/

long poly_div(UINT a[], UINT b[], UINT q[], UINT r[])
{
  //printf("Beg\n");
  UINT aa[MAX_POLY_SIZE], aaa[MAX_POLY_SIZE], bb[MAX_POLY_SIZE], bq[MAX_POLY_SIZE];
  if(poly_deg(b) == -1) return -1;
  poly_null(q);
  poly_null(r);
  poly_copy(a, aa);
  //printf("At the beggining: deg(aa):%d deg(b):%d\n", poly_deg(aa), poly_deg(b));
  while(poly_deg(aa) >= poly_deg(b))
    {
      poly_copy(b, bb);
      poly_shift(bb, poly_deg(aa)-poly_deg(b));
      q[poly_deg(aa)-poly_deg(b)] = 1;
      //printf("First: deg(aa):%d deg(b):%d deg(bb):%d\n", poly_deg(aa), poly_deg(b), poly_deg(bb));
      poly_add(aa, bb, aaa);
      //printf("aa[199]: %d\n", aa[199]);
      //printf("bb[199]: %d\n", bb[199]);
      //printf("aaa[199]: %d\n", aaa[199]);
      poly_copy(aaa, aa);
      //  printf("Second: deg(aa):%d deg(b):%d deg(bb):%d deg(aaa): %d\n",
      //	     poly_deg(aa), poly_deg(b), poly_deg(bb), poly_deg(aaa));
    }
  poly_mul(b, q, bq);
  poly_add(a, bq, r);
  //printf("End\n");
  return 0;
}

/* The following function implements Extended Euclidean algorithm

   Inputs: a, b - polynomials

   Inputs/Outputs: d, x, y - polynomials satsifying a * x + b * y = d = gcd(a,b)
*/

void poly_ext_euclidean(UINT a[], UINT b[], UINT d[], UINT x[], UINT y[])
{
  //printf("1\n");
  
  UINT aa[MAX_POLY_SIZE], bb[MAX_POLY_SIZE], u[MAX_POLY_SIZE], v[MAX_POLY_SIZE], xa[MAX_POLY_SIZE],
    yb[MAX_POLY_SIZE], q[MAX_POLY_SIZE], r[MAX_POLY_SIZE], qu[MAX_POLY_SIZE], qv[MAX_POLY_SIZE],
    xqu[MAX_POLY_SIZE], yqv[MAX_POLY_SIZE];
  if(poly_deg(a) < poly_deg(b)) // In this case just applying this procedure will cause infinite loop.
    {
      poly_ext_euclidean(b,a,d,y,x);
      return;
    }
  if(poly_div(a,b,q,r) == -1) // b is zero polynomial.
    {
      poly_copy(a,d); // MUST CHECK
      return;
    }
  // printf("2\n");
  poly_one(x);
  poly_one(v);
  poly_null(y);
  poly_null(u);
  poly_copy(a, aa);
  poly_copy(b, bb);
  //printf("3\n");
  do{
    //   printf("4\n");
    poly_div(aa,bb,q,r);
    //printf("4a\n");
    poly_mul(q,u,qu);
    poly_add(x,qu,xqu);
    //printf("4b\n");
    poly_mul(q,v,qv);
    poly_add(y,qv,yqv);
    poly_copy(bb,aa);
    poly_copy(r,bb);
    //printf("4c\n");
    poly_copy(u,x);
    poly_copy(v,y);
    poly_copy(xqu,u);
    //printf("4d\n");
    poly_copy(yqv,v);
  }while(poly_deg(bb) != -1);
  // printf("5\n");
  poly_mul(x,a,xa);
  poly_mul(y,b,yb);
  poly_add(xa,yb,d);
  //printf("6\n");
}

/* The following function performs transformation from DWORD representation to polynomial representation

   Input: dword

   Input/Output: poly
*/

void dword_to_poly(DWORD& dword, UINT poly[])
{
  poly_null(poly);
  for(WORD i = 0; i < 64; i++) poly[i] = (dword>>i)&1; 
}

/* The following function performs transformation from polynomial representation to DWORD representation

   Input: poly

   Input/Output: dword
*/

void poly_to_dword(UINT poly[], DWORD& dword)
{
  DWORD aux = 1;
  dword = 0;
  for(WORD i = 0; i < 64; i++) if(poly[i]) dword ^= aux<<i; 
}

/* The following function performs transformation from DWORD representation to PAIR representation

   Input: dword

   Input/Output: pair
*/

void dword_to_pair(DWORD& dword, PAIR& pair)
{
  DWORD copy = dword;
  pair.y = copy&0xFFFFFFFF;
  copy >>= 32;
  pair.x = copy&0xFFFFFFFF;
}

/* The following function performs transformation from PAIR representation to DWORD representation

   Input: pair

   Input/Output: dword
*/

void pair_to_dword(PAIR& pair, DWORD& dword)
{
  dword = pair.x;
  dword <<= 32;
  dword ^= pair.y;
}

/* The following function performs transformation from polynomial representation to PAIR representation

   Input: poly

   Input/Output: pair
*/

void poly_to_pair(UINT poly[], PAIR& pair)
{
  DWORD aux;
  poly_to_dword(poly, aux);
  dword_to_pair(aux, pair);
}

/* The following function performs transformation from PAIR representation to polynomial representation

   Input: pair

   Input/Output: poly
*/

void pair_to_poly(PAIR& pair, UINT poly[])
{
  DWORD aux;
  pair_to_dword(pair, aux);
  dword_to_poly(aux, poly);
}

/* The following function performs pair multiplication operation

   Input: pair_a - first pair, pair_b - second pair

   Ouput: pair_a * pair_b
*/

PAIR pair_mul(PAIR pair_a, PAIR pair_b)
{
  PAIR pair_c;
  UINT a[MAX_POLY_SIZE], b[MAX_POLY_SIZE], m[MAX_POLY_SIZE],
    q[MAX_POLY_SIZE], r[MAX_POLY_SIZE], res[MAX_POLY_SIZE];
  poly_null(m);
  m[0] = m[1] = m[2] = m[11] = m[64] = 1;
  pair_to_poly(pair_a, a);
  pair_to_poly(pair_b, b);
  poly_mul(a,b,res);
  poly_div(res,m,q,r);
  poly_to_pair(r,pair_c);
  return pair_c;
}

/* The following function performs fast pair multiplication operation

   Input: pair_a - first pair, pair_b - second pair

   Ouput: pair_a * pair_b
*/

PAIR pair_fast_mul(PAIR pair_a, PAIR pair_b)
{
  PAIR pair_c;
  DWORD a, aa, aaa = 0, b, c, m = 2055; // m[0] = m[1] = m[2] = m[11] = 1; 1+2+4+2048 = 2055
  pair_to_dword(pair_a, a);
  pair_to_dword(pair_b, b);

  for(WORD i = 0; i < 64; i++)
    if((b>>i)&1)
      {
	// c +=  xtime(...(xtime(a))...), i times

	aa = a;
	for(WORD j = 0; j < i; j++)
	  {
	    if((aa>>63)&1)
	      {
		aa <<= 1;
		aa ^= m;
	      }
	    else aa <<= 1;
	  }
	aaa ^= aa;
      }
  
  dword_to_pair(aaa,pair_c);
  return pair_c;
}

/* The following function performs pair invertion operation

   Input: pair_a

   Output: (pair_a)^{-1}
*/

PAIR pair_inv(PAIR pair_a)
{
  PAIR pair_inv_a; 
  UINT a[MAX_POLY_SIZE], inv_a[MAX_POLY_SIZE], m[MAX_POLY_SIZE], d[MAX_POLY_SIZE],
    x[MAX_POLY_SIZE], y[MAX_POLY_SIZE], q[MAX_POLY_SIZE];
  poly_null(m);
  m[0] = m[1] = m[2] = m[11] = m[64] = 1;
  pair_to_poly(pair_a, a);
  poly_ext_euclidean(a,m,d,x,y);
  poly_div(x,m,q,inv_a);
  poly_to_pair(inv_a, pair_inv_a);
  return pair_inv_a;
}

/* The following function checks performance of polynomial arithmetic operations */

void poly_check()
{
  /*PAIR a, aa, aaa;
  random_word(a.x);
  random_word(a.y);
  aa = pair_inv(a);
  aaa = pair_inv(aa);
  print_pair(a); NEWLINE;
  print_pair(aa); NEWLINE;
  print_pair(aaa); NEWLINE;
  NEWLINE;
  NEWLINE;
  NEWLINE;*/
  /*WORD K5, K6, K7, K8;
  random_word(K5);
  random_word(K6);
  random_word(K7);
  random_word(K8);
  PAIR K5_K6_pair(K5, K6);
  PAIR K7_K8_pair(K7, K8);
  PAIR K5_K6_pair_inv = pair_mul(K5_K6_pair,K7_K8_pair);
  PAIR K7_K8_pair_inv = pair_inv(K7_K8_pair); 
  PAIR middle;
  random_pair(middle);
  print_pair(middle); NEWLINE;
  middle = middle + K5_K6_pair;
  PAIR aux = pair_mul(middle, K7_K8_pair);
  middle = aux;
  middle = middle + K5_K6_pair_inv;
  aux = pair_mul(middle, K7_K8_pair_inv);
  middle = aux;
  print_pair(middle); NEWLINE;*/
  /*coconut98_random_init();
  PAIR middle;
  random_pair(middle);
  print_pair(middle); NEWLINE;
  middle = middle + coconut98_K5_K6_pair;
  PAIR aux = pair_mul(middle, coconut98_K7_K8_pair);
  middle = aux;
  middle = middle + coconut98_K5_K6_pair_inv;
  aux = pair_mul(middle, coconut98_K7_K8_pair_inv);
  middle = aux;
  print_pair(middle); NEWLINE;*/
  /*PAIR a,b,c,ab,ac,bc,abc1, abc2;
  random_pair(a);
  random_pair(b);
  random_pair(c);
  ab = a+b;
  ac = pair_mul(a,c);
  bc = pair_mul(b,c);
  abc1 = pair_mul(ab,c);
  abc2 = ac+bc;
  print_pair(abc1); NEWLINE;
  print_pair(abc2); NEWLINE;*/
  PAIR a, b, c1, c2;
  random_pair(a);
  random_pair(b);
  c1 = pair_mul(a,b);
  c2 = pair_fast_mul(a,b);
  print_pair(c1); NEWLINE;
  print_pair(c2); NEWLINE;
}

/* The following function performs encryption operation of COCONUT98

   Inputs: plaintext - plaintext before the encryption
           k1, k2, k3, k4, k5, k6, k7, k8, K5, K6, K7, K8 - keys used for the encryption

   Output: ciphertext	    
*/

PAIR coconut98_encrypt(PAIR& plaintext, WORD& k1, WORD& k2, WORD& k3, WORD& k4, WORD& k5,
		       WORD& k6, WORD& k7, WORD& k8, WORD& K5, WORD& K6, WORD& K7, WORD& K8)
{
  PAIR middle = coconut98_Psi(plaintext, k1, k2, k3, k4);

#ifdef DECORRELATION_MODULE
  PAIR aux1(K5, K6), aux2(K7, K8), aux3;
  middle = middle + aux1;
  aux3 = pair_fast_mul(middle,aux2);
  middle = aux3;
#endif  
  
  return coconut98_Psi(middle, k5, k6, k7, k8);
}

/* The following function performs decryption operation of COCONUT98

   Inputs: ciphertext - ciphertext before the decryption
           k1, k2, k3, k4, k5, k6, k7, k8, K5, K6, K7, K8 - keys used for the decryption

   Output: plaintext	    
*/

PAIR coconut98_decrypt(PAIR& ciphertext, WORD& k1, WORD& k2, WORD& k3, WORD& k4, WORD& k5,
		       WORD& k6, WORD& k7, WORD& k8, WORD& K5, WORD& K6, WORD& K7, WORD& K8)
{
  PAIR middle = coconut98_Psi(ciphertext, k8, k7, k6, k5);
  
#ifdef DECORRELATION_MODULE
  PAIR aux1(K5, K6), aux2(K7, K8), aux3, aux4, aux5;
  aux3 = pair_fast_mul(aux1,aux2);
  aux4 = pair_inv(aux2);
  middle = middle + aux3;
  aux5 = pair_fast_mul(middle,aux4);
  middle = aux5;
#endif
  
  return coconut98_Psi(middle, k4, k3, k2, k1);
}

/* The following function performs encryption trial of the attacker

   Input: plaintext 

   Ouput: ciphertext
*/

PAIR coconut98_encryption_trial(PAIR& plaintext)
{
  PAIR middle = coconut98_Psi(plaintext, coconut98_k1, coconut98_k2, coconut98_k3, coconut98_k4);

#ifdef DECORRELATION_MODULE
  middle = middle + coconut98_K5_K6_pair;
  PAIR aux = pair_fast_mul(middle, coconut98_K7_K8_pair);
  middle = aux;
#endif  
  
  return coconut98_Psi(middle, coconut98_k5, coconut98_k6, coconut98_k7, coconut98_k8);
}

/* The following function performs decryption trial of the attacker

   Input: ciphertext 

   Ouput: plaintext
*/

PAIR coconut98_decryption_trial(PAIR& ciphertext)
{
  PAIR middle = coconut98_Psi(ciphertext, coconut98_k8, coconut98_k7, coconut98_k6, coconut98_k5);
  
#ifdef DECORRELATION_MODULE
  middle = middle + coconut98_K5_K6_pair_inv;
  PAIR aux = pair_fast_mul(middle, coconut98_K7_K8_pair_inv);
  middle = aux;
#endif
  
  return coconut98_Psi(middle, coconut98_k4, coconut98_k3, coconut98_k2, coconut98_k1);
}

/* The following function performs encryption trial of the attacker of one round reduced cipher

   Inputs: plaintext - plaintext before the encryption
           key - subkey used as k1

   Ouput: ciphertext
*/

PAIR coconut98_encryption_trial_2345678(PAIR& plaintext, WORD key)
{
  PAIR tmp = coconut98_decryption_round(plaintext, key);
  return coconut98_encryption_trial(tmp);
}

/* The following function performs decryption trial of the attacker of one round reduced cipher

   Inputs: ciphertext - ciphertext before the decryption
           key - subkey used as k1

   Ouput: plaintext
*/

PAIR coconut98_decryption_trial_2345678(PAIR& ciphertext, WORD key)
{
  PAIR tmp = coconut98_decryption_trial(ciphertext);
  return coconut98_encryption_round(tmp, key);
}

/* The following function performs encryption trial of the attacker of two round reduced cipher

   Inputs: plaintext - plaintext before the encryption
           key1, key2 - subkeys used as k1, k2 respectively

   Ouput: ciphertext
*/

PAIR coconut98_encryption_trial_345678(PAIR& plaintext, WORD key1, WORD key2)
{
  PAIR tmp = coconut98_decryption_round(plaintext, key2);
  PAIR tmp1 = coconut98_decryption_round(tmp, key1);
  return coconut98_encryption_trial(tmp1);
}

/* The following function performs decryption trial of the attacker of two round reduced cipher

   Inputs: ciphertext - ciphertext before the decryption
           key1, key2 - subkeys used as k1, k2 respectively

   Ouput: plaintext
*/

PAIR coconut98_decryption_trial_345678(PAIR& ciphertext, WORD key1, WORD key2)
{
  PAIR tmp = coconut98_decryption_trial(ciphertext);
  PAIR tmp1 = coconut98_encryption_round(tmp, key1);
  return coconut98_encryption_round(tmp1, key2);
}

/* The following function checks that COCONUT98's encryption and decryption operations
are inverses of one another */

void coconut98_check_inverse()
{
  // Checking that COCONUT98's decryption is inverse of its encryption.
  
  PAIR plaintext, ciphertext;
  random_pair(plaintext);
  print_pair(plaintext);
  NEWLINE;
  ciphertext = coconut98_encryption_trial(plaintext);
  print_pair(ciphertext);
  NEWLINE;
  PAIR back_plaintext = coconut98_decryption_trial(ciphertext);
  print_pair(back_plaintext);
  NEWLINE;
}

/* The following function implements boomerang attack on COCONUT98 */

void coconut98_boomerang_attack()
{
  // Here we implement Boomerang attack on COCONUT98. We use functions coconut98_enryption_trial and
  // coconut98_decryption_trial as encryption and decryption quiries of the cipher respectively.  
  
  PAIR P, P_tag, Q, Q_tag, Q_diff, C, C_tag, D, D_tag, delta[COCONUT98_CHARS], nabla[COCONUT98_CHARS],  
    right_k1_P[COCONUT98_CHARS][COCONUT98_RIGHT_PAIRS], right_k1_P_tag[COCONUT98_CHARS][COCONUT98_RIGHT_PAIRS],
    right_k8_P[COCONUT98_CHARS][COCONUT98_RIGHT_PAIRS], right_k8_P_tag[COCONUT98_CHARS][COCONUT98_RIGHT_PAIRS],
    right_k2_P[COCONUT98_CHARS_SECOND][COCONUT98_RIGHT_PAIRS_SECOND],
    right_k2_P_tag[COCONUT98_CHARS_SECOND][COCONUT98_RIGHT_PAIRS_SECOND],
    right_k3_P[COCONUT98_CHARS_THIRD][COCONUT98_RIGHT_PAIRS_THIRD],
    right_k3_P_tag[COCONUT98_CHARS_THIRD][COCONUT98_RIGHT_PAIRS_THIRD],
    p1, pp1, c1, cc1, p2, pp2, c2, cc2, pp12, cc12, pp12_inv, K5_K6, K7_K8, inv_K7_K8, cc1_inv_K7_K8;  
  
  WORD right_pairs, after_F_diff, x_diff, i, key, iter, j, k1, k2, k3, k4,
    k5, k6, k7, k8, K1, K2, K3, K4, K5, K6, K7, K8, queries_counter = 0;
  bool is_correct_key, found_key;

#ifndef SEARCH_RANGE
  printf("Attack should take around 2 hours...\n\n");
#endif
  
  // Characteristics.
  
#ifdef DECORRELATION_MODULE
  delta[0].x = 0x00080000; 
  delta[0].y = 0x00000100; // 0000 0000 0000 0000 0000 0001 0000 0000 bit 8 
  nabla[0] = delta[0];

  delta[1].x = 0x00100000; 
  delta[1].y = 0x00000200; // 0000 0000 0000 0000 0000 0010 0000 0001 bit 9 
  nabla[1] = delta[1];

  delta[2].x = 0x20000000; 
  delta[2].y = 0x00040000; // 0000 0000 0000 0100 0000 0000 0000 0000 bit 18  
  nabla[2] = delta[2];

  delta[3].x = 0x40000000; 
  delta[3].y = 0x00080000; // 0000 0000 0000 1000 0000 0000 0000 0000 bit 19
  nabla[3] = delta[3]; 

  delta[4].x = 0x80000000; 
  delta[4].y = 0x00100000; // 0000 0000 0001 0000 0000 0000 0000 0000 bit 20 
  nabla[4] = delta[4];

  delta[5].x = 0x00000100; 
  delta[5].y = 0x20000000; // 0010 0000 0000 0000 0000 0000 0000 0000 bit 29
  nabla[5] = delta[5];

  delta[6].x = 0x00000200; 
  delta[6].y = 0x40000000; // 0100 0000 0000 0000 0000 0000 0000 0000 bit 30 
  nabla[6] = delta[6];
  
  delta[7].x = 0x00000400; 
  delta[7].y = 0x80000000; // 1000 0000 0000 0000 0000 0000 0000 0000 bit 31 
  nabla[7] = delta[7];

  // These two take two much trials
  
  //delta[6].x = 0x00000200; 
  //delta[6].y = 0x40000000; // 0100 0000 0000 0000 0000 0000 0000 0000 bit 30 
  //nabla[6] = delta[6];

  //delta[7].x = 0x40000000; 
  //delta[7].y = 0x00080000; // 0000 0000 0000 1000 0000 0000 0000 0000 bit 19
  //nabla[7] = delta[7]; 
#else
  delta[2].x = 0x10000000;
  delta[2].y = 0x00020000; // 0000 0000 0000 0010 0000 0000 0000 0000 bit 17
  nabla[2] = delta[2]; 
  
  delta[1].x = 0x01000000; 
  delta[1].y = 0x00002000; // 0000 0000 0000 0000 0010 0000 0000 0000 bit 13 
  nabla[1] = delta[1]; 
  
  delta[0].x = 0x02000000; 
  delta[0].y = 0x00004000; // 0000 0000 0000 0000 0100 0000 0000 0000 bit 14 
  nabla[0] = delta[0]; 

  delta[3].x = 0x04000000; 
  delta[3].y = 0x00008000; // 0000 0000 0000 0000 1000 0000 0000 0000 bit 15
  nabla[3] = delta[3]; 

  delta[4].x = 0x08000000; 
  delta[4].y = 0x00010000; // 0000 0000 0000 0001 0000 0000 0000 0000 bit 16
  nabla[4] = delta[4]; 

  delta[5].x = 0x00800000; 
  delta[5].y = 0x00001000; // 0000 0000 0000 0000 0001 0000 0000 0000 bit 12
  nabla[5] = delta[5]; 
 
  delta[6].x = 0x00400000; 
  delta[6].y = 0x00000800; // 0000 0000 0000 0000 0000 1000 0000 0000 bit 11 
  nabla[6] = delta[6];
#endif

  // Generating right quartets for k1 and k8 recovery.

  int ccc;
  
  for(iter = 0; iter < COCONUT98_CHARS; iter++)
    {
      //printf("Iter: %d\n", iter);
      right_pairs = 0;
      ccc = 0;

      while(right_pairs < COCONUT98_RIGHT_PAIRS) 
	{
	  ccc++;
	  
	  // Generating quartet.
	
	  random_pair(P);
	  P_tag =  P+delta[iter];
	  C = coconut98_encryption_trial(P);
	  C_tag = coconut98_encryption_trial(P_tag);
	  D = C+nabla[iter];
	  D_tag = C_tag+nabla[iter];
	  Q = coconut98_decryption_trial(D);
	  Q_tag = coconut98_decryption_trial(D_tag);
	  Q_diff = Q + Q_tag;

	  // Updating queries counter.
	  
	  queries_counter += 4;
	  
	  // Checking whether this quartet is right.

	  if(Q_diff.x != delta[iter].x || Q_diff.y != delta[iter].y) continue;
 
	  // Filling array of right pairs for k1 and k8.

	  right_k1_P[iter][right_pairs] = P;
	  right_k8_P[iter][right_pairs] = C;
	  right_k1_P_tag[iter][right_pairs] = P_tag;
	  right_k8_P_tag[iter][right_pairs++] = D;
	  right_k1_P[iter][right_pairs] = Q;
	  right_k8_P[iter][right_pairs] = C_tag;
	  right_k1_P_tag[iter][right_pairs] = Q_tag;
	  right_k8_P_tag[iter][right_pairs++] = D_tag;

#ifdef PRINT
	  PAIR middle,aux;
	  
	  WORD out2 = P.x^coconut98_F(P.y^coconut98_k1);    
      	  WORD out3 = P.y^coconut98_F(out2^coconut98_k2);   
	  WORD out4 = out2^coconut98_F(out3^coconut98_k3);  
	  WORD out5 = out3^coconut98_F(out4^coconut98_k4);
	  WORD out2_tag = P_tag.x^coconut98_F(P_tag.y^coconut98_k1); 
	  WORD out3_tag = P_tag.y^coconut98_F(out2_tag^coconut98_k2);
	  WORD out4_tag = out2_tag^coconut98_F(out3_tag^coconut98_k3);
	  WORD out5_tag = out3_tag^coconut98_F(out4_tag^coconut98_k4);

	  print_word(out2); SPACE; print_word(out2_tag); SPACE; print_word(out2^out2_tag); NEWLINE;
	  print_word(out3); SPACE; print_word(out3_tag); SPACE; print_word(out3^out3_tag); NEWLINE;
	  print_word(out4); SPACE; print_word(out4_tag); SPACE; print_word(out4^out4_tag); NEWLINE;
	  print_word(out5); SPACE; print_word(out5_tag); SPACE; print_word(out5^out5_tag); NEWLINE; NEWLINE;  

#ifdef DECORRELATION_MODULE	  
	  middle = PAIR(out5,out4);
	  middle = middle + coconut98_K5_K6_pair;
	  aux = pair_fast_mul(middle, coconut98_K7_K8_pair);
	  out5 = aux.x;
	  out4 = aux.y;
	  
	  middle = PAIR(out5_tag,out4_tag);
	  middle = middle + coconut98_K5_K6_pair;
	  aux = pair_fast_mul(middle, coconut98_K7_K8_pair);
	  out5_tag = aux.x;
	  out4_tag = aux.y;
#endif
	  
       	  WORD out6 = out5^coconut98_F(out4^coconut98_k5);    
	  WORD out7 = out4^coconut98_F(out6^coconut98_k6);   
	  WORD out8 = out6^coconut98_F(out7^coconut98_k7);  
	  WORD out9 = out7^coconut98_F(out8^coconut98_k8);
	  WORD out6_tag = out5_tag^coconut98_F(out4_tag^coconut98_k5);    
	  WORD out7_tag = out4_tag^coconut98_F(out6_tag^coconut98_k6);   
	  WORD out8_tag = out6_tag^coconut98_F(out7_tag^coconut98_k7);  
	  WORD out9_tag = out7_tag^coconut98_F(out8_tag^coconut98_k8);
	  PAIR W, W_tag;
	  W.x = out9 ^ nabla[iter].x;
	  W.y = out8 ^ nabla[iter].y;
	  W_tag.x = out9_tag ^ nabla[iter].x;
	  W_tag.y = out8_tag ^ nabla[iter].y;
	  
	  WORD out10 = W.x^coconut98_F(W.y^coconut98_k8); 
	  WORD out11 = W.y^coconut98_F(out10^coconut98_k7);
	  WORD out12 = out10^coconut98_F(out11^coconut98_k6);
	  WORD out13 = out11^coconut98_F(out12^coconut98_k5);
	  WORD out10_tag = W_tag.x^coconut98_F(W_tag.y^coconut98_k8); 
	  WORD out11_tag = W_tag.y^coconut98_F(out10_tag^coconut98_k7);
	  WORD out12_tag = out10_tag^coconut98_F(out11_tag^coconut98_k6);
	  WORD out13_tag = out11_tag^coconut98_F(out12_tag^coconut98_k5);

	  print_word(out6); SPACE; print_word(out6_tag); SPACE; print_word(out6^out6_tag); NEWLINE; 
	  print_word(out7); SPACE; print_word(out7_tag); SPACE; print_word(out7^out7_tag); NEWLINE;
	  print_word(out8); SPACE; print_word(out8_tag); SPACE; print_word(out8^out8_tag); NEWLINE;
	  print_word(out9); SPACE; print_word(out9_tag); SPACE; print_word(out9^out9_tag); NEWLINE; NEWLINE;
	  print_word(out10); SPACE; print_word(out10_tag); SPACE; print_word(out10^out10_tag); NEWLINE;
	  print_word(out11); SPACE; print_word(out11_tag); SPACE; print_word(out11^out11_tag); NEWLINE;
	  print_word(out12); SPACE; print_word(out12_tag); SPACE; print_word(out12^out12_tag); NEWLINE;
	  print_word(out13); SPACE; print_word(out13_tag); SPACE; print_word(out13^out13_tag); NEWLINE; NEWLINE;

#ifdef DECORRELATION_MOFULE
	  middle = PAIR(out13,out12);
	  middle = middle + coconut98_K5_K6_pair_inv;
	  aux = pair_fast_mul(middle, coconut98_K7_K8_pair_inv);
	  out13 = aux.x;
	  out12 = aux.y;
	  
	  middle = PAIR(out13_tag,out12_tag);
	  middle = middle + coconut98_K5_K6_pair_inv;
	  aux = pair_fast_mul(middle, coconut98_K7_K8_pair_inv);
	  out13_tag = aux.x;
	  out12_tag = aux.y;
#endif
	  
	  WORD out14 = out13^coconut98_F(out12^coconut98_k4); 
	  WORD out15 = out12^coconut98_F(out14^coconut98_k3);
	  WORD out16 = out14^coconut98_F(out15^coconut98_k2);
	  WORD out17 = out15^coconut98_F(out16^coconut98_k1);
	  WORD out14_tag = out13_tag^coconut98_F(out12_tag^coconut98_k4); 
	  WORD out15_tag = out12_tag^coconut98_F(out14_tag^coconut98_k3);
	  WORD out16_tag = out14_tag^coconut98_F(out15_tag^coconut98_k2);
	  WORD out17_tag = out15_tag^coconut98_F(out16_tag^coconut98_k1);

	  print_word(out14); SPACE; print_word(out14_tag); SPACE; print_word(out14^out14_tag); NEWLINE; 
	  print_word(out15); SPACE; print_word(out15_tag); SPACE; print_word(out15^out15_tag); NEWLINE;
	  print_word(out16); SPACE; print_word(out16_tag); SPACE; print_word(out16^out16_tag); NEWLINE;
	  print_word(out17); SPACE; print_word(out17_tag); SPACE; print_word(out17^out17_tag); NEWLINE;
#endif
	}

      //print_word(delta[iter].y);
      //printf(" used %d quartets\n", (int)(ccc/(COCONUT98_RIGHT_PAIRS/2.0)));
      //printf("%d & ", (int)(ccc/(COCONUT98_RIGHT_PAIRS/2.0)));
    }

  //printf("in average used %d quartets\n", (int)(queries_counter/4.0/COCONUT98_CHARS/(COCONUT98_RIGHT_PAIRS/2.0)));
  //printf("%d\n\n", (int)(queries_counter/4.0/COCONUT98_CHARS/(COCONUT98_RIGHT_PAIRS/2.0)));
  //printf("queries counter: %d\n", queries_counter);

  // k1 recovery.

  printf("Recovering k1...\n\n");

  found_key = false;
  
#ifdef SEARCH_RANGE
  //for(key = coconut98_k1; key <= coconut98_k1; key++)
  for(key = max(coconut98_k1 - RANGE,0); key < min(coconut98_k1 + RANGE,0xFFFFFFFF); key++)
#else
    for(key = 0; key <= 0xFFFFFFFF; key++)
    //for(key = 0xFFFFFFFF; key >= 0; key--)
#endif
    { 
      is_correct_key = true;
      
      for(iter = 0; iter < COCONUT98_CHARS; iter++)
	{
	  for(j = 0; j < COCONUT98_RIGHT_PAIRS; j++)
	    {
	      after_F_diff = coconut98_F(key^right_k1_P[iter][j].y)^
		coconut98_F(key^right_k1_P_tag[iter][j].y);
	      x_diff = right_k1_P[iter][j].x^right_k1_P_tag[iter][j].x;
	      
	      if(after_F_diff != x_diff) 
		{
		  // Key is wrong.
		  
		  is_correct_key = false;
		  break;
		}
	    }

	  if(!is_correct_key) break;
	}

      if(is_correct_key)
	{
	  printf("Attack found following value for k1: ");
	  print_word(key);
	  NEWLINE;
	  NEWLINE;
	  k1 = key;
	  found_key = true;
	  break;
	}
    }

  if(!found_key)
    {
      printf("Attack failed to find k1\n");
      return;
    }

  // k8 recovery.

  printf("Recovering k8...\n\n"); 

  found_key = false;

#ifdef SEARCH_RANGE  
  for(key = max(coconut98_k8 - RANGE,0); key < min(coconut98_k8 + RANGE,0xFFFFFFFF); key++)
#else
    //for(key = 0; key <= 0xFFFFFFFF; key++)
    for(key = 0xFFFFFFFF; key >= 0; key--) 
#endif    
    {
      is_correct_key = true;
      
      for(iter = 0; iter < COCONUT98_CHARS; iter++)
	{
	  for(j = 0; j < COCONUT98_RIGHT_PAIRS; j++)
	    {
	      after_F_diff = coconut98_F(key^right_k8_P[iter][j].y)^
		coconut98_F(key^right_k8_P_tag[iter][j].y);
	      x_diff = right_k8_P[iter][j].x^right_k8_P_tag[iter][j].x;
	      if(after_F_diff != x_diff)
		{
		  // Key is wrong.
		  
		  is_correct_key = false;
		  break;
		}
	    }

	  if(!is_correct_key) break;
	}

      if(is_correct_key)
	{
	  printf("Attack found following value for k8: ");
	  print_word(key);
	  NEWLINE;
	  NEWLINE;
	  k8 = key;
	  found_key = true;
	  break;
	}
    }

  if(!found_key)
    {
      printf("Attack failed to find k8\n");
      return;
    }
  
  // Generating right quartets for k2 recovery.

  for(iter=0; iter < COCONUT98_CHARS_SECOND; iter++)
    {
      right_pairs = 0;
      while(right_pairs < COCONUT98_RIGHT_PAIRS_SECOND)
	{
	  // Generating quartet.
	  
	  random_pair(P);
	  P_tag =  P+delta[iter];
	  D = coconut98_encryption_trial_2345678(P,k1)+nabla[iter];
	  Q = coconut98_decryption_trial_2345678(D,k1);
	  D_tag = coconut98_encryption_trial_2345678(P_tag,k1)+nabla[iter];
	  Q_tag = coconut98_decryption_trial_2345678(D_tag,k1);
	  Q_diff = Q + Q_tag;

	  // Updating queries counters.

	  queries_counter += 4;
	  
	  // Checking whether this quartet is right.

	  if(Q_diff.x != delta[iter].x || Q_diff.y != delta[iter].y) continue;
	  
	  // Filling array of right pairs.
	  
	  right_k2_P[iter][right_pairs] = P; 
	  right_k2_P_tag[iter][right_pairs++] = P_tag;
	  right_k2_P[iter][right_pairs] = Q;
	  right_k2_P_tag[iter][right_pairs++] = Q_tag; 
	}
    }

  // k2 recovery.

  printf("Recovering k2...\n\n");

  found_key = false;

#ifdef SEARCH_RANGE
  for(key = max(coconut98_k2 - RANGE,0); key < min(coconut98_k2 + RANGE,0xFFFFFFFF); key++)
#else
  for(key = 0; key <= 0xFFFFFFFF; key++)
#endif
    {
      is_correct_key = true;
      
      for(iter = 0; iter < COCONUT98_CHARS_SECOND; iter++)
	{
	  for(j = 0; j < COCONUT98_RIGHT_PAIRS_SECOND; j++)
	    {
	      after_F_diff = coconut98_F(key^right_k2_P[iter][j].y)^
		coconut98_F(key^right_k2_P_tag[iter][j].y);
	      x_diff = right_k2_P[iter][j].x^right_k2_P_tag[iter][j].x;
	      
	      if(after_F_diff != x_diff) 
		{
		  // Key is wrong.
		  
		  is_correct_key = false;
		  break;
		}
	    }

	  if(!is_correct_key) break;
	}

      if(is_correct_key)
	{
	  printf("Attack found following value for k2: ");
	  print_word(key);
	  NEWLINE;
	  NEWLINE;
	  k2 = key;
	  found_key = true;
	  break;
	}
    }

  if(!found_key)
    {
      printf("Attack failed to find k2\n");
      return;
    }
  
  // Generating right quartets for k3 recovery.

  for(iter=0; iter < COCONUT98_CHARS_THIRD; iter++)
    {
      right_pairs = 0;
      while(right_pairs < COCONUT98_RIGHT_PAIRS_THIRD)
	{
	  // Generating quartet.
	  
	  random_pair(P);
	  P_tag =  P+delta[iter];
	  D = coconut98_encryption_trial_345678(P,k1,k2)+nabla[iter];
	  Q = coconut98_decryption_trial_345678(D,k1,k2);
	  D_tag = coconut98_encryption_trial_345678(P_tag,k1,k2)+nabla[iter];
	  Q_tag = coconut98_decryption_trial_345678(D_tag,k1,k2);
	  Q_diff = Q + Q_tag;

	  // Updating queries counters.

	  queries_counter += 4;
	  
	  // Checking whether this quartet is right.

	  if(Q_diff.x != delta[iter].x || Q_diff.y != delta[iter].y) continue;
	  
	  // Filling array of right pairs.
	  
	  right_k3_P[iter][right_pairs] = P; 
	  right_k3_P_tag[iter][right_pairs++] = P_tag;
	  right_k3_P[iter][right_pairs] = Q;
	  right_k3_P_tag[iter][right_pairs++] = Q_tag;
	}
    }
  
  // k3 recovery.

  printf("Recovering k3...\n\n"); 

  found_key = false;

#ifdef SEARCH_RANGE
  for(key = max(coconut98_k3 - RANGE,0); key < min(coconut98_k3 + RANGE,0xFFFFFFFF); key++)
#else
  for(key = 0; key <= 0xFFFFFFFF; key++)
#endif
    {
      is_correct_key = true;
      
      for(iter = 0; iter < COCONUT98_CHARS_THIRD; iter++)
	{
	  for(j = 0; j < COCONUT98_RIGHT_PAIRS_THIRD; j++)
	    {
	      after_F_diff = coconut98_F(key^right_k3_P[iter][j].y)^
		coconut98_F(key^right_k3_P_tag[iter][j].y);
	      x_diff = right_k3_P[iter][j].x^right_k3_P_tag[iter][j].x;
	      
	      if(after_F_diff != x_diff) 
		{
		  // Key is wrong.
		  
		  is_correct_key = false;
		  break;
		}
	    }
	  
	  if(!is_correct_key) break;
	}
      
      if(is_correct_key)
	{
	  printf("Attack found following value for k3: ");
	  print_word(key);
	  NEWLINE;
	  NEWLINE;
	  k3 = key;
	  found_key = true;
	  break;
	}
    }

  if(!found_key)
    {
      printf("Attack failed to find k3\n");
      return;
    }

  // After recovering k1 k2 k3 k8 we can seacrh for right keys k3 and k8 in little Hamming distance from keys k3 and
  // k8 that we found. For each choice k1 k2 _k3 _k8 we can recover rest keys via weaknesses in key schedule algorithm
  // and check the result by encrypting random plaintext and comparing result to true ciphertext.

  printf("Trying to variate the keys k3 and k8 a little in order to get the right keys...\n\n");

  WORD _k3, _k8;
  PAIR plaintext, ciphertext, ciphertext2;
  random_pair(plaintext);
  ciphertext = coconut98_encryption_trial(plaintext);
  random_pair(p1);
  random_pair(p2);
  c1 = coconut98_encryption_trial(p1);
  c2 = coconut98_encryption_trial(p2);
  
  for(WORD bit1 = 0; bit1 <= 32; bit1++)
      for(WORD bit2 = 0; bit2 <= 32; bit2++)
	for(WORD bit3 = 0; bit3 <= 32; bit3++)
	  for(WORD bit4 = 0; bit4 <= 32; bit4++)
	    {
	      _k3 = k3;
	      _k8 = k8;
	      if(bit1 != 32) _k3 ^= 1<<bit1;
	      if(bit2 != 32) _k3 ^= 1<<bit2;
       	      if(bit3 != 32) _k8 ^= 1<<bit3;
	      if(bit4 != 32) _k8 ^= 1<<bit4;
		  
	      K1 = k1;
	      K3 = k1^k2;
	      K4 = k2^_k3;
	      K2 = _k8^K4;
	      
	      k4 = K1^K4;
	      k5 = K2;
	      k6 = K2^K3;
	      k7 = k6^K4;

#ifdef DECORRELATION_MODULE
	      // Breaking decorrelation module.
	     
	      pp1 = coconut98_Psi(p1, k1, k2, _k3, k4);
	      pp2 = coconut98_Psi(p2, k1, k2, _k3, k4);
	      cc1 = coconut98_Psi(c1, _k8, k7, k6, k5); 
	      cc2 = coconut98_Psi(c2, _k8, k7, k6, k5);
	      cc12 = cc1+cc2;
	      pp12 = pp1+pp2;
	      pp12_inv = pair_inv(pp12);
	      K7_K8 = pair_fast_mul(cc12,pp12_inv);
	      inv_K7_K8 = pair_inv(K7_K8);
	      cc1_inv_K7_K8 = pair_fast_mul(cc1,inv_K7_K8);
	      K5_K6 = cc1_inv_K7_K8+pp1;
	      K5 = K5_K6.x;
	      K6 = K5_K6.y;
	      K7 = K7_K8.x;
	      K8 = K7_K8.y;
#endif
	      
	      // Checking correctness of the keys. 
	      
	      ciphertext2 = coconut98_encrypt(plaintext, k1, k2, _k3, k4, k5, k6, k7, _k8, K5, K6, K7, K8);
	      
	      if(ciphertext.x == ciphertext2.x && ciphertext.y == ciphertext2.y) goto Out;	     
	    }
  
 Out:
  
  k3 = _k3;
  k8 = _k8;
  
  K1 = k1;
  K3 = k1^k2;
  K4 = k2^k3;
  K2 = k8^K4;
  
  k4 = K1^K4;
  k5 = K2;
  k6 = K2^K3;
  k7 = k6^K4;

  printf("Attack used %d encryption/decryption trials and found:\n", queries_counter);
  printf("K1: ");
  print_word(K1);
  printf(" K2: ");
  print_word(K2);
  printf(" K3: ");
  print_word(K3);
  printf(" K4: ");
  print_word(K4);
  printf("\nk1: ");
  print_word(k1);
  printf(" k2: ");
  print_word(k2);
  printf(" k3: ");
  print_word(k3);
  printf(" k4: ");
  print_word(k4);
  printf(" k5: ");
  print_word(k5);
  printf(" k6: ");
  print_word(k6);
  printf(" k7: ");
  print_word(k7);
  printf(" k8: ");
  print_word(k8);

  ciphertext2 = coconut98_encrypt(plaintext, k1, k2, k3, k4, k5, k6, k7, k8, K5, K6, K7, K8);
  if(ciphertext.x == ciphertext2.x && ciphertext.y == ciphertext2.y)
    printf("\n\nResults of the attack are correct with very high probability.");
  else
  printf("\n\nResults of the attack are wrong.");     
}

/* The following function initializes randomly S boxes of Khufu */ 

void khufu_random_init()
{
  WORD i, j, tmp, khufu_S1_byte0[256], khufu_S1_byte1[256], khufu_S1_byte2[256], khufu_S1_byte3[256],
    khufu_S2_byte0[256], khufu_S2_byte1[256], khufu_S2_byte2[256], khufu_S2_byte3[256];  

  // Generating random S(256) permutations khufu_S[1,2]_byte[1,2,3,4].

  khufu_S1_byte0[0] = khufu_S1_byte1[0] = khufu_S1_byte2[0] = khufu_S1_byte3[0] =
  khufu_S2_byte0[0] = khufu_S2_byte1[0] = khufu_S2_byte2[0] = khufu_S2_byte3[0] = 0;
  
  for(i=1; i<256; i++)
    {
      j = random()%(i+1);
      khufu_S1_byte0[i] = khufu_S1_byte0[j];
      khufu_S1_byte0[j] = i;
      j = random()%(i+1);
      khufu_S1_byte1[i] = khufu_S1_byte1[j];
      khufu_S1_byte1[j] = i;
      j = random()%(i+1);
      khufu_S1_byte2[i] = khufu_S1_byte2[j];
      khufu_S1_byte2[j] = i;
      j = random()%(i+1);
      khufu_S1_byte3[i] = khufu_S1_byte3[j];
      khufu_S1_byte3[j] = i;
      j = random()%(i+1);
      khufu_S2_byte0[i] = khufu_S2_byte0[j];
      khufu_S2_byte0[j] = i;
      j = random()%(i+1);
      khufu_S2_byte1[i] = khufu_S2_byte1[j];
      khufu_S2_byte1[j] = i;
      j = random()%(i+1);
      khufu_S2_byte2[i] = khufu_S2_byte2[j];
      khufu_S2_byte2[j] = i;
      j = random()%(i+1);
      khufu_S2_byte3[i] = khufu_S2_byte3[j];
      khufu_S2_byte3[j] = i;
    }

  // Using khufu_S[1,2]_byte[1,2,3,4] generating khufu_S1 and khufu_S2.

  for(i=0; i<256; i++)
    {
      khufu_S1[i] = khufu_S1_byte0[i] + 256*khufu_S1_byte1[i] + 65536*khufu_S1_byte2[i] + 16777216*khufu_S1_byte3[i];
      khufu_S2[i] = khufu_S2_byte0[i] + 256*khufu_S2_byte1[i] + 65536*khufu_S2_byte2[i] + 16777216*khufu_S2_byte3[i];
    }
}  

/* The following function prints Khufu S boxes */

void khufu_print_config()
{
  printf("Khufu's configuration is as follows:\n");
  printf("S box 1:");
  for(WORD i=0; i<256; i++)
    {
      if(!(i%12)) NEWLINE;
      else SPACE;
      print_word(khufu_S1[i]);
    }
  printf("\nS box 2:");
  for(WORD i=0; i<256; i++)
    {
      if(!(i%12)) NEWLINE;
      else SPACE;
      print_word(khufu_S2[i]);
    }
}

/* The following function performs one round of Khufu's encryption operation

   Inputs: plaintext - plaintext before the encryption
           S - S box used in this round

   Output: ciphertext
*/

PAIR khufu_encryption_round(PAIR plaintext, WORD S[])
{
  PAIR ciphertext;
  ciphertext.x = S[plaintext.x&0x000000FF]^plaintext.y;
  ciphertext.y = ((plaintext.x&0x000000FF)<<8) ^ ((plaintext.x&0x0000FF00)<<8) ^
    ((plaintext.x&0x00FF0000)<<8) ^ ((plaintext.x&0xFF000000)>>24);
  return ciphertext;
}

/* The following function performs one round of Khufu's decryption operation

   Inputs: ciphertext - ciphertext before the decryption
           S - S box used in this round

   Output: plaintext
*/

PAIR khufu_decryption_round(PAIR ciphertext, WORD S[])
{
  PAIR plaintext;
  plaintext.x = ((ciphertext.y&0x000000FF)<<24) ^ ((ciphertext.y&0x0000FF00)>>8) ^
    ((ciphertext.y&0x00FF0000)>>8) ^ ((ciphertext.y&0xFF000000)>>8);
  plaintext.y = S[plaintext.x&0x000000FF]^ciphertext.x;
  return plaintext;
}

/* The following function performs Khufu's encryption operation

   Inputs: plaintext - plaintext before the encryption
           S1, S2 - S boxes used

   Output: ciphertext
*/

PAIR khufu_encrypt(PAIR& plaintext, WORD S1[], WORD S2[])
{
  int i;
  PAIR tmp(plaintext.x, plaintext.y);
  
  for(i=0; i<8; i++) tmp = khufu_encryption_round(tmp, S1);
  for(i=0; i<8; i++) tmp = khufu_encryption_round(tmp, S2);

  return PAIR(tmp.x, tmp.y);
}

/* The following function performs Khufu's decryption operation

   Inputs: ciphertext - ciphertext before the decryption
           S1, S2 - S boxes used 

   Output: plaintext
*/

PAIR khufu_decrypt(PAIR& ciphertext, WORD S1[], WORD S2[])                                       
{
  int i;
  PAIR tmp(ciphertext.x, ciphertext.y);
  
  for(i=0; i<8; i++) tmp = khufu_decryption_round(tmp, S2);
  for(i=0; i<8; i++) tmp = khufu_decryption_round(tmp, S1);

  return PAIR(tmp.x, tmp.y);
}

/* The following function performs attacker's encryption trial operation

   Input: plaintext

   Output: ciphertext
*/

PAIR khufu_encryption_trial(PAIR& plaintext)
{
  return khufu_encrypt(plaintext, khufu_S1, khufu_S2);
}

/* The following function performs attacker's decryption trial operation

   Input: ciphertext

   Output: plaintext
*/

PAIR khufu_decryption_trial(PAIR& ciphertext)
{
  return khufu_decrypt(ciphertext, khufu_S1, khufu_S2);
}

/* The following function checks that Khufu's encryption and decryption operations
   are inverses of one another */

void khufu_check_inverse()
{
  // Checking that Khufu's decryption is inverse of its encryption.
  
  PAIR plaintext, ciphertext;
  random_pair(plaintext);
  print_pair(plaintext);
  NEWLINE;
  ciphertext = khufu_encryption_trial(plaintext);
  print_pair(ciphertext);
  NEWLINE;
  PAIR back_plaintext = khufu_decryption_trial(ciphertext);
  print_pair(back_plaintext);
  NEWLINE;
}

/* The following function implements first variant of the partition algorithm */ 

WORD Partition(PAIR arr[], PAIR snd_arr[], WORD p, WORD r)
{
  WORD x = arr[p].x;
  WORD i = p-1, j = r+1;
  for(;;)
    {
      do{ j--; } while(arr[j].x > x);
      do{ i++; } while(arr[i].x < x);
      if(i < j)
	{
	  PAIR tmp = arr[i];
	  arr[i] = arr[j];
	  arr[j] = tmp;
	  tmp = snd_arr[i];
	  snd_arr[i] = snd_arr[j];
	  snd_arr[j] = tmp;
	}
      else return j;
    }
}

/* The following function implements first variant of the quicksort algorithm */

void Quicksort(PAIR arr[], PAIR snd_arr[], WORD p, WORD r)
{
  if(p < r)
    {
      WORD q = Partition(arr, snd_arr, p, r);
      Quicksort(arr, snd_arr, p, q);
      Quicksort(arr, snd_arr, q+1, r);
    }
}

/* The following function implements first variant of the binary search algorithm */

long BinarySearch(PAIR arr[], WORD val, WORD p, WORD r)
{
  if(r <= p+1)
    {
      if(arr[p].x == val) return p;
      if(arr[r].x == val) return r;
      return -1;	
    }
  
  WORD mid = (long)((double)(p+r)/2.0);

  if(arr[mid].x == val) return (long)(mid);
  if(val < arr[mid].x) return BinarySearch(arr, val, p, mid);
  if(val > arr[mid].x) return BinarySearch(arr, val, mid, r);
}

/* The following function implements second variant of the partition algorithm */

WORD PartitionSnd(PAIR arr[], PAIR snd_arr[], WORD p, WORD r)
{
  WORD x = arr[p].y;
  WORD i = p-1, j = r+1;
  for(;;)
    {
      do{ j--; } while(arr[j].y > x);
      do{ i++; } while(arr[i].y < x);
      if(i < j)
	{
	  PAIR tmp = arr[i];
	  arr[i] = arr[j];
	  arr[j] = tmp;
	  tmp = snd_arr[i];
	  snd_arr[i] = snd_arr[j];
	  snd_arr[j] = tmp;
	}
      else return j;
    }
}

/* The following function implements second variant of the quicksort algorithm */

void QuicksortSnd(PAIR arr[], PAIR snd_arr[], WORD p, WORD r)
{
  if(p < r)
    {
      WORD q = PartitionSnd(arr, snd_arr, p, r);
      QuicksortSnd(arr, snd_arr, p, q);
      QuicksortSnd(arr, snd_arr, q+1, r);
    }
}

/* The following function implements second variant of the binary search algorithm */

long BinarySearchSnd(PAIR arr[], WORD val, WORD p, WORD r)
{
  if(r <= p+1)
    {
      if(arr[p].y == val) return p;
      if(arr[r].y == val) return r;
      return -1;	
    }
  
  WORD mid = (long)((double)(p+r)/2.0);

  if(arr[mid].y == val) return (long)(mid);
  if(val < arr[mid].y) return BinarySearchSnd(arr, val, p, mid);
  if(val > arr[mid].y) return BinarySearchSnd(arr, val, mid, r);
}

/* The following function implements simple variant of the partition algorithm */ 
   
WORD PartitionSimple(WORD arr[], WORD p, WORD r)
{
  WORD x = arr[p];
  WORD i = p-1, j = r+1;
  for(;;)
    {
      do{ j--; } while(arr[j] > x);
      do{ i++; } while(arr[i] < x);
      if(i < j)
	{
	  WORD tmp = arr[i];
	  arr[i] = arr[j];
	  arr[j] = tmp;
	}
      else return j;
    }
}

/* The following function implements simple variant of the quicksort algorithm */ 

void QuicksortSimple(WORD arr[], WORD p, WORD r)
{
  if(p < r)
    {
      WORD q = PartitionSimple(arr, p, r);
      QuicksortSimple(arr, p, q);
      QuicksortSimple(arr, q+1, r);
    }
}

/* The following function implements simple variant of the binary search algorithm */ 

long BinarySearchSimple(WORD arr[], WORD val, WORD p, WORD r)
{
  if(r <= p+1)
    {
      if(arr[p] == val) return p;
      if(arr[r] == val) return r;
      return -1;	
    }
  
  WORD mid = (long)((double)(p+r)/2.0);

  if(arr[mid] == val) return (long)(mid);
  if(val < arr[mid]) return BinarySearchSimple(arr, val, p, mid);
  if(val > arr[mid]) return BinarySearchSimple(arr, val, mid, r);
}

/* Following function checks birthday paradox

   Inputs: pool - pool size
           trials - number of trials in the calculation

   Output: estimation of probability of the birthday paradox	   
*/

double birthday_paradox_prob_calc(WORD pool, WORD trials)
{
  WORD iter, i, num_success = 0, P[MAX_POOL_SIZE], Q[MAX_POOL_SIZE];
  
  for(iter = 0; iter < trials; iter++)
    {
      for(i = 0; i < pool; i++) 
	{
	  random_word(P[i]);
	  random_word(Q[i]);
	}
      
      QuicksortSimple(Q, 0, pool-1);
      
      for(i = 0; i < pool; i++)
	if(BinarySearchSimple(Q, P[i], 0, pool-1) != -1)
	  {
	    num_success++;
	    break;
	  }
    }

  return (double)(trials - num_success) / (double)(trials);
}

/* The following function implements boomerang attack on Khufu */

void khufu_boomerang_attack()
{
  // Here we implement Boomerang attack on Khufu. We use functions khufu_enryption_trial and
  // khufu_decryption_trial as encryption and decryption quiries of the cipher respectively.  
  
  WORD a, L, L_tag, i, j, x_diff, y_diff, pair_number, k, alpha, tagged_Q, iter = 0;
  PAIR P[KHUFU_POOL_SIZE], P_tag[KHUFU_POOL_SIZE], C, C_tag, D, D_tag, _Q, _Q_tag,
    Q[KHUFU_POOL_SIZE], Q_tag[KHUFU_POOL_SIZE], nabla, alpha_pair, CC[KHUFU_POOL_SIZE],   
    right_P[KHUFU_ITERS][KHUFU_RIGHT_PAIRS], right_P_tag[KHUFU_ITERS][KHUFU_RIGHT_PAIRS],
    right, right_tag, new_P, new_P_tag, candidate_P[KHUFU_RIGHT_PAIRS],
    candidate_P_tag[KHUFU_RIGHT_PAIRS];
  WORD S1[256], S2[256], num_found, old_num_found, candidate_num, queries_counter = 0;  
  bool recovered[256], right_candidate;

  // Preparing to recover S1.

  printf("\nPreparing to recover S1...\n");

  random_word(L);
  
  for(i = 0; i < KHUFU_POOL_SIZE; i++) 
    {
      // First pool of 2^16 plaintexts/ciphertexts.  
      
      P[i].x = L; 
      random_word(P[i].y); 
      CC[i] = khufu_encryption_trial(P[i]);
       
      queries_counter++;
    }
    
  for(;;)
    { 
    Beg:      
      if(iter == KHUFU_ITERS) break;

      // Generating random a, L, L_tag, alpha and alpha_pair.

      a = 1<<iter;
      nabla = PAIR(a,0);
      L_tag = L^a;

      // Generating pool of plaintexts and ciphertexts.
      
      for(i = 0; i < KHUFU_POOL_SIZE; i++) 
	{
	  // First pool of 2^16 plaintexts/ciphertexts.  
	   
          CC[i] = khufu_encryption_trial(P[i]);	  
	  D = CC[i]+nabla; 
	  Q[i] = khufu_decryption_trial(D);
		  
	  // Second pool of 2^16 plaintexts/ciphertexts.
	  
	  P_tag[i].x = L_tag;
	  random_word(P_tag[i].y);
	  C_tag = khufu_encryption_trial(P_tag[i]); 
	  D_tag = C_tag+nabla; 
	  Q_tag[i] = khufu_decryption_trial(D_tag);

	  queries_counter += 3;
	}
      
      // Quicksort of Q_tag[j]'s and P_tag[j]'s for j in [0, KHUFU_POOL_SIZE) w.r.t. Q_tag[j].x's.
      
      Quicksort(Q_tag, P_tag, 0, KHUFU_POOL_SIZE-1);
      
      // Seeking for candidadte right pair via binary search algorithm.

      pair_number = 0;
      
      for(i = 0; i < KHUFU_POOL_SIZE; i++)
	{
	  tagged_Q = Q[i].x^a;
	  if(BinarySearch(Q_tag, tagged_Q, 0, KHUFU_POOL_SIZE-1) == -1) continue;
	  
	  j = BinarySearch(Q_tag, tagged_Q, 0, KHUFU_POOL_SIZE-1);
	  
	  candidate_P[pair_number] = P[i]; 
	  candidate_P_tag[pair_number++] = P_tag[j];	  
	}
      
      if(pair_number == 0)
	{
	  printf("Failing to find any candidate.\n");
	  goto Beg;
	}
       
      // KHUFU_RIGHT_PAIRS more right pairs generation.

      for(candidate_num = 0; candidate_num < pair_number; candidate_num++)
	{
	  right_candidate = true;
	  for(k = 0; k < KHUFU_RIGHT_PAIRS; k++)
	    {
	      random_word(alpha);
	      alpha &= 0xFFFFFF00;
	      alpha_pair = PAIR(alpha, 0);
	      new_P = candidate_P[candidate_num]+alpha_pair;
	      C = khufu_encryption_trial(new_P);
	      D = C+nabla;
	      _Q = khufu_decryption_trial(D);
	      new_P_tag = candidate_P_tag[candidate_num]+alpha_pair;
	      C_tag = khufu_encryption_trial(new_P_tag);
	      D_tag = C_tag+nabla;
	      _Q_tag = khufu_decryption_trial(D_tag);
	      right_P[iter][k] = _Q;
	      right_P_tag[iter][k] = _Q_tag;
	      
	      // Checking this candidate.
	      
	      new_P = _Q+alpha_pair;
	      C = khufu_encryption_trial(new_P);
	      D = C+nabla;
	      _Q = khufu_decryption_trial(D);
	      new_P_tag = _Q_tag+alpha_pair;
	      C_tag = khufu_encryption_trial(new_P_tag);
	      D_tag = C_tag+nabla;
	      _Q_tag = khufu_decryption_trial(D_tag);
	      x_diff = _Q.x^_Q_tag.x;

	      queries_counter += 8;
	      
	      if(x_diff != a)
		{
		  right_candidate = false;
		  break;
		}		  	      	      
	  }

	  if(right_candidate)
	    {
	      printf("Found right candidate.\n");	      
	      iter++;
	      goto Beg;
	    }
	}

      printf("All candidates are wrong.\n");
    }
  
  // Recovering S1.
  // We use the fact that
  // S1[right_P[iter][i].x%256] ^ S1[right_P_tag[iter][i].x%256] = right_P[iter][i].y^right_P_tag[iter][i].y for every
  // iter in [0, KHUFU_ITERS) and i in [0, KHUFU_RIGHT_PAIRS).
  
  printf("\nRecovering S1...\n");
  
  for(i = 0; i < 256; i++) recovered[i] = false;
  recovered[0] = true;  
  S1[0] = 0; // here we can put arbitrary value

  old_num_found = 0;
  
  for(;;) 
    {
      num_found = 0; 
      for(i = 0; i < 256; i++) if(recovered[i]) num_found++;

      if(num_found == 256) break; 
      
      if(num_found == old_num_found)
	{
	  printf("\nAttack failed to find S1.");
	  return;
	}

      for(j = 0; j < 256; j++) 
	if(recovered[j]) 
	  for(iter = 0; iter < KHUFU_ITERS; iter++)
	    for(i = 0; i < KHUFU_RIGHT_PAIRS; i++)		  
	      {
		if(right_P[iter][i].x%256 == j) 
		  {
		    recovered[right_P_tag[iter][i].x%256] = true; 
		    S1[right_P_tag[iter][i].x%256] = S1[j] ^ right_P[iter][i].y ^ right_P_tag[iter][i].y; 
		  }
		if(right_P_tag[iter][i].x%256 == j) 
		  {
		    recovered[right_P[iter][i].x%256] = true; 
		    S1[right_P[iter][i].x%256] = S1[j] ^ right_P[iter][i].y ^ right_P_tag[iter][i].y; 
		  }
	      }
      
      old_num_found = num_found;
    }

  // Preparing to recover S2.

  iter = 0;
  
  printf("\nPreparing to recover S2...\n");

  random_word(L);

  for(i = 0; i < KHUFU_POOL_SIZE; i++) 
    { 
      // First pool of 2^16 plaintexts/ciphertexts.  
      
      P[i].y = L; 
      random_word(P[i].x); 
      CC[i] = khufu_decryption_trial(P[i]);

      queries_counter++;
    }
  
  for(;;)
    {
    BegSnd:
      if(iter == KHUFU_ITERS) break;

      // Generating random a, L, L_tag, alpha and alpha_pair.

      a = 0x100<<iter;
      a &= 0x0000FF00;
      nabla = PAIR(0,a);
      L_tag = L^a;

      // Generating pool of plaintexts and ciphertexts.
      
      for(i = 0; i < KHUFU_POOL_SIZE; i++) 
	{ 
	  // First pool of 2^16 plaintexts/ciphertexts.  
	  
	  D = CC[i]+nabla; 
	  Q[i] = khufu_encryption_trial(D);
	  
	  // Second pool of 2^16 plaintexts/ciphertexts.
	  
	  P_tag[i].y = L_tag;
	  random_word(P_tag[i].x);
	  C_tag = khufu_decryption_trial(P_tag[i]); 
	  D_tag = C_tag+nabla; 
	  Q_tag[i] = khufu_encryption_trial(D_tag);

	  queries_counter += 3;
	}
      
      // Quicksort of Q_tag[j]'s and P_tag[j]'s for j in [0, KHUFU_POOL_SIZE) w.r.t. Q_tag[j].y's.
      
      QuicksortSnd(Q_tag, P_tag, 0, KHUFU_POOL_SIZE-1);
      
      // Seeking for candidadte right pair via binary search algorithm.
      
      pair_number = 0;
      
      for(i = 0; i < KHUFU_POOL_SIZE; i++)
	{
	  tagged_Q = Q[i].y^a;
	  if(BinarySearchSnd(Q_tag, tagged_Q, 0, KHUFU_POOL_SIZE-1) == -1) continue;
	  
	  j = BinarySearchSnd(Q_tag, tagged_Q, 0, KHUFU_POOL_SIZE-1);
	  candidate_P[pair_number] = P[i]; 
	  candidate_P_tag[pair_number++] = P_tag[j];	  
	}
      
      if(pair_number == 0)
	{
	  printf("Failing to find any candidate.\n");
	  goto BegSnd;
	}
       
      // KHUFU_RIGHT_PAIRS more right pairs generation.

      for(candidate_num = 0; candidate_num < pair_number; candidate_num++)
	{
	  right_candidate = true;
	  for(k = 0; k < KHUFU_RIGHT_PAIRS; k++)
	    {
	      random_word(alpha);
	      alpha &= 0xFFFF00FF;
	      alpha_pair = PAIR(0, alpha);
	      new_P = candidate_P[candidate_num]+alpha_pair;
	      C = khufu_decryption_trial(new_P);
	      D = C+nabla;
	      _Q = khufu_encryption_trial(D);
	      new_P_tag = candidate_P_tag[candidate_num]+alpha_pair;
	      C_tag = khufu_decryption_trial(new_P_tag);
	      D_tag = C_tag+nabla;
	      _Q_tag = khufu_encryption_trial(D_tag);
	      right_P[iter][k] = _Q;
	      right_P_tag[iter][k] = _Q_tag;
	      
	      // Checking this candidate.
	      
	      new_P = _Q+alpha_pair;
	      C = khufu_decryption_trial(new_P);
	      D = C+nabla;
	      _Q = khufu_encryption_trial(D);
	      new_P_tag = _Q_tag+alpha_pair;
	      C_tag = khufu_decryption_trial(new_P_tag);
	      D_tag = C_tag+nabla;
	      _Q_tag = khufu_encryption_trial(D_tag);
	      y_diff = _Q.y^_Q_tag.y;

	      queries_counter += 8;
	      
	      if(y_diff != a)
		{
		  right_candidate = false;
		  break;
		}		  	      	      
	  }

	  if(right_candidate)
	    {
	      printf("Found right candidate.\n");	      
	      iter++;
	      goto BegSnd;
	    }
	}

      printf("All candidates are wrong.\n");
    }
  
  // Recovering S2.
  // We use the fact that
  // S2[(right_P[iter][i].y>>8)%256] ^ S2[(right_P_tag[iter][i].x>>8)%256] =
  // right_P[iter][i].x^right_P_tag[iter][i].x for every
  // iter in [0, KHUFU_ITERS) and i in [0, KHUFU_RIGHT_PAIRS).
  
  printf("\nRecovering S2...\n");
  
  for(i = 0; i < 256; i++) recovered[i] = false;
  recovered[0] = true;  
  S2[0] = 0; // here we can put arbitrary value

  old_num_found = 0;
  
  for(;;) 
    { 
      num_found = 0; 
      for(i = 0; i < 256; i++) if(recovered[i]) num_found++;
      
      if(num_found == 256) break;
      
      if(num_found == old_num_found)
	{
	  printf("\nAttack failed to find S2.");
	  return;
	}

      for(j = 0; j < 256; j++) 
	if(recovered[j]) 
	  for(iter = 0; iter < KHUFU_ITERS; iter++)
	    for(i = 0; i < KHUFU_RIGHT_PAIRS; i++)		  
	      {
		if((right_P[iter][i].y>>8)%256 == j) 
		  {
		    recovered[(right_P_tag[iter][i].y>>8)%256] = true; 
		    S2[(right_P_tag[iter][i].y>>8)%256] = S2[j] ^ right_P[iter][i].x ^ right_P_tag[iter][i].x; 
		  }
		if((right_P_tag[iter][i].y>>8)%256 == j) 
		  {
		    recovered[(right_P[iter][i].y>>8)%256] = true; 
		    S2[(right_P[iter][i].y>>8)%256] = S2[j] ^ right_P[iter][i].x ^ right_P_tag[iter][i].x; 
		  }
	      }

       old_num_found = num_found;
    }

  // Printing results.

  printf("\nAttack used %d encryption/decryption trials and found S boxes of the cipher up to 32-bit shift.",
	 queries_counter); 
    
  // Checking results of the attack.

  bool right_result = true;

  for(i = 0; i<256; i++)
    {
      WORD our_S1_diff = S1[i]^S1[0], true_S1_diff = khufu_S1[i]^khufu_S1[0],
	our_S2_diff = S2[i]^S2[0], true_S2_diff = khufu_S2[i]^khufu_S2[0];
      if(our_S1_diff != true_S1_diff || our_S2_diff != true_S2_diff) right_result = false;
    }
  
  if(right_result) 
    printf("\n\nResults of the attack are correct."); 
  else 
  printf("\n\nResults of the attack are wrong.");  
}

/* The following function calculates probability of right quartet of Khufu

   Input: trials - number of trials used in the calculation

   Output: probability of right quartet
*/

double khufu_right_quartet_prob_calc(WORD trials)
{
  // Here we implement Boomerang attack on Khufu. We use functions khufu_enryption_trial and
  // khufu_decryption_trial as encryption and decryption quiries of the cipher respectively.  
  
  WORD a, L, L_tag, i, j, x_diff, y_diff, pair_number, k, alpha, tagged_Q, iter = 0;
  PAIR P[KHUFU_POOL_SIZE], P_tag[KHUFU_POOL_SIZE], C, C_tag, D, D_tag, _Q, _Q_tag,
    Q[KHUFU_POOL_SIZE], Q_tag[KHUFU_POOL_SIZE], nabla, alpha_pair, CC[KHUFU_POOL_SIZE],   
    right_P[KHUFU_ITERS][KHUFU_RIGHT_PAIRS], right_P_tag[KHUFU_ITERS][KHUFU_RIGHT_PAIRS],
    right, right_tag, new_P, new_P_tag, candidate_P[KHUFU_RIGHT_PAIRS],
    candidate_P_tag[KHUFU_RIGHT_PAIRS];
  WORD S1[256], S2[256], num_found, old_num_found, candidate_num, queries_counter = 0, jjj, ccc = 0;  
  bool recovered[256], right_candidate;

  // Preparing to recover S1.

  random_word(L);
  
  for(i = 0; i < KHUFU_POOL_SIZE; i++) 
    {
      // First pool of 2^16 plaintexts/ciphertexts.  
      
      P[i].x = L; 
      random_word(P[i].y); 
      CC[i] = khufu_encryption_trial(P[i]);
       
      queries_counter++;
    }
    
  for(jjj = 0; jjj < trials; jjj++)
    {      
      // Generating random a, L, L_tag, alpha and alpha_pair.

      a = 1<<(jjj%8);
      nabla = PAIR(a,0);
      L_tag = L^a;

      // Generating pool of plaintexts and ciphertexts.
      
      for(i = 0; i < KHUFU_POOL_SIZE; i++) 
	{
	  // First pool of 2^16 plaintexts/ciphertexts.  
	   
          CC[i] = khufu_encryption_trial(P[i]);	  
	  D = CC[i]+nabla; 
	  Q[i] = khufu_decryption_trial(D);
		  
	  // Second pool of 2^16 plaintexts/ciphertexts.
	  
	  P_tag[i].x = L_tag;
	  random_word(P_tag[i].y);
	  C_tag = khufu_encryption_trial(P_tag[i]); 
	  D_tag = C_tag+nabla; 
	  Q_tag[i] = khufu_decryption_trial(D_tag);

	  queries_counter += 3;
	}
      
      // Quicksort of Q_tag[j]'s and P_tag[j]'s for j in [0, KHUFU_POOL_SIZE) w.r.t. Q_tag[j].x's.
      
      Quicksort(Q_tag, P_tag, 0, KHUFU_POOL_SIZE-1);
      
      // Seeking for candidadte right pair via binary search algorithm.

      pair_number = 0;
      
      for(i = 0; i < KHUFU_POOL_SIZE; i++)
	{
	  tagged_Q = Q[i].x^a;
	  if(BinarySearch(Q_tag, tagged_Q, 0, KHUFU_POOL_SIZE-1) == -1) continue;
	  
	  j = BinarySearch(Q_tag, tagged_Q, 0, KHUFU_POOL_SIZE-1);
	  
	  candidate_P[pair_number] = P[i]; 
	  candidate_P_tag[pair_number++] = P_tag[j];	  
	}
      
      if(pair_number == 0) continue;
       
      // KHUFU_RIGHT_PAIRS more right pairs generation.

      for(candidate_num = 0; candidate_num < pair_number; candidate_num++)
	{
	  right_candidate = true;
	  for(k = 0; k < KHUFU_RIGHT_PAIRS; k++)
	    {
	      random_word(alpha);
	      alpha &= 0xFFFFFF00;
	      alpha_pair = PAIR(alpha, 0);
	      new_P = candidate_P[candidate_num]+alpha_pair;
	      C = khufu_encryption_trial(new_P);
	      D = C+nabla;
	      _Q = khufu_decryption_trial(D);
	      new_P_tag = candidate_P_tag[candidate_num]+alpha_pair;
	      C_tag = khufu_encryption_trial(new_P_tag);
	      D_tag = C_tag+nabla;
	      _Q_tag = khufu_decryption_trial(D_tag);
	      right_P[iter][k] = _Q;
	      right_P_tag[iter][k] = _Q_tag;
	      
	      // Checking this candidate.
	      
	      new_P = _Q+alpha_pair;
	      C = khufu_encryption_trial(new_P);
	      D = C+nabla;
	      _Q = khufu_decryption_trial(D);
	      new_P_tag = _Q_tag+alpha_pair;
	      C_tag = khufu_encryption_trial(new_P_tag);
	      D_tag = C_tag+nabla;
	      _Q_tag = khufu_decryption_trial(D_tag);
	      x_diff = _Q.x^_Q_tag.x;

	      queries_counter += 8;
	      
	      if(x_diff != a)
		{
		  right_candidate = false;
		  break;
		}		  	      	      
	  }

	  if(right_candidate)
	    {
      	      ccc++;
	      break;
	    }
	}
    }

  return (double)(ccc) / (double)(trials);
}