/*
  Questo programma fa parte di una libreria messa a disposizione dei
  lettori del libro
  Alessandro Languasco & Alessandro Zaccagnini,
  Introduzione alla Crittografia, Ulrico Hoepli Editore, Milano, 2004.

  Algoritmo di Shanks per la determinazione del logaritmo discreto
  nel gruppo Z_p^*, p primo.
  Si veda il Paragrafo 6.8.1 del libro citato.
*/

#include <iostream>
#include <cmath>
#include <vector>
#include <algorithm> 
#include "aritmetica.cc"

/*
  Questo dato contiene una coppia del tipo "(g^a mod p, a)", dove "g"
  è un generatore di Z_p^*, "a" è un intero nell'intervallo
  [0, p-1].
*/
typedef struct coppia_potenza_esponente {
  long long potenza;
  long long esponente;
};

/*
  Modificare il valore di questa costante in relazione alla quantità
  di memoria disponibile
*/
static const long long unsigned
MAX_PRIME = 1000000;

/*
  Questa funzione stabilisce quale sia "minore" fra due oggetti di tipo
  "coppia_potenza_esponente": è rilevante solo il valore nel campo
  "potenza".
  Questa funzione, passata come argomento all'algoritmo di "sort",
  permette un ordinamento efficiente delle coppie generate nei "Baby
  steps".  In sostanza, le coppie "(a^n mod p, n)", sono ordinate
  lessicograficamente, tenendo conto solo del primo elemento della
  coppia.
*/
bool
minore(coppia_potenza_esponente a, coppia_potenza_esponente b) {

  if (a.potenza < b.potenza)
    return(true);
  else
    return(false);
}

/*
  Ricerca binaria del valore "x" nel vettore ordinato "coppie"
*/
long long
ricerca_binaria(long long x,
		std::vector <coppia_potenza_esponente> &coppie) {

  long long unsigned dim = coppie.size();
  long long unsigned inizio = 0;
  long long unsigned fine = dim - 1;
  // NB: se "inizio = fine - 1" si entra in un ciclo infinito a causa
  // degli arrotondamenti
  while (inizio < fine - 1) {
    long long unsigned meta = (inizio + fine) / 2;
    // Se si trova l'elemento desiderato, si ritorna la sua posizione
    if (coppie[meta].potenza == x)
      return(meta);
    // Se non si trova l'elemento cercato, si elimina la metà del
    // vettore che certamente non lo contiene
    if (coppie[meta].potenza < x)
      inizio = meta;
    else
      fine = meta;
  }
  // Se si esce dal ciclo while, allora "inizio = fine", oppure
  // "inizio = fine - 1"

  if (coppie[inizio].potenza == x)
    return(inizio);
  if (coppie[fine].potenza == x)
    return(fine);
  // Se l'elemento non compare, restituisci -1
  return(-1);
}

/*
  Calcolo del logaritmo discreto modulo "p" per mezzo dell'Algoritmo
  di Shanks. Qui
  1. "p" è un numero primo;
  2. "g" è un generatore di Z_p^*;
  3. "x" è l'elemento di Z_p^* di cui si cerca il logaritmo discreto.
  La funzione ritorna "y" in Z_p^* tale che g^y = x (in Z_p^*).
  NB: La funzione NON controlla che "p" sia primo né che "g" sia un
  generatore.
*/
long long
shanks(long long p, long long g, long long x) {

  if (p < 0)
    p *= -1;
  // "m" è il numero di Baby steps necessari
  long long unsigned m = 1 + radice_quadrata(p);
  if (m >= MAX_PRIME)
    return(-1);

  // Memorizziamo le coppie "(g^i mod p, i)", per i = 0, ..., m-1
  // nel vettore valori
  std::vector <coppia_potenza_esponente> valori;
  long long potenze_di_g = 1;

  // Baby steps: calcolo delle potenze successive del generatore "g"
  for(unsigned i = 0; i < m; ++i) {
    coppia_potenza_esponente valore;
    valore.potenza = potenze_di_g;
    valore.esponente = i;
    valori.push_back(valore);
    if(potenze_di_g == x)
      return(i);
    potenze_di_g = moltiplicazione_mod_n(potenze_di_g, g, p);
  }
  // A questo punto "potenze_di_g = g^m mod p"

  if (x == potenze_di_g)
    return(m);

  std::cout << "Baby steps, prima dell'ordinamento ";
  for(unsigned i = 0; i < m; ++i)
    std::cout << valori[i].potenza << " (" << valori[i].esponente << ") ";
  std::cout << std::endl;

  // Invochiamo l'algoritmo sort sul vettore "valori", usando la
  // funzione di confronto "minore" definita sopra
  sort(valori.begin(), valori.end(), minore);

  std::cout << "Baby steps,  dopo    l'ordinamento ";
  for(unsigned i = 0; i < m; ++i)
    std::cout << valori[i].potenza << " (" << valori[i].esponente << ") ";
  std::cout << std::endl;

  // Preparazione dei Giant steps
  long long lambda;
  long long mi;
  long long MCD = mcd_est(potenze_di_g, p, lambda, mi);
  std::cout << "mcd come combinazione lineare ";
  std::cout << MCD << " = " << potenze_di_g << " * " << lambda << " + ";
  std::cout << p << " * " << mi << std::endl;
  // Correggi "lambda" se "lambda < 0", in modo che appartenga
  // all'intervallo [0, p-1]
  if (lambda < 0)
    lambda += p;
  // Ora "lambda" è l'inverso di "potenze_di_g mod p"

  // Giant steps
  long long potenze_alte = moltiplicazione_mod_n(lambda, x, p);
  for(unsigned i = 0; i < m; ++i) {
    std::cout << "x * g^" << -(i+1) * m << " = " << potenze_alte << std::endl;
    // Ricerca binaria
    long long j = ricerca_binaria(potenze_alte, valori);
    if (j != -1)
      return((i+1) * m + valori[j].esponente);
    potenze_alte = moltiplicazione_mod_n(potenze_alte, lambda, p);
  }

  // Ritorna -1 in caso di errore; per esempio, se "g" non fosse un
  // generatore, il ciclo qui sopra potrebbe non dare una risposta
  return(-1);
}

void
SHANKS(long long p, long long g, long long x) {

  long long log_discreto = shanks(p, g, x);
  std::cout << "Il logaritmo discreto di " << x;
  std::cout << " rispetto al generatore " << g;
  std::cout << " in Z/" << p << "^* e' " << log_discreto << std::endl;
  std::cout << "---------------------------------------";
  std::cout << "---------------------------------------" << std::endl;
}

int
main() {

  long long p;
  long long g;
  long long x;

  // Alcuni esempi di uso della funzione "SHANKS"
  p = 101;
  g = 3;
  x = 30;
  SHANKS(p, g, x);

  p = 9973;
  g = 2;
  x = 5;
  SHANKS(p, g, x);

  p = 65537;
  g = 3;
  x = 5;
  SHANKS(p, g, x);

  p = 61;
  g = 2;
  x = 15;
  SHANKS(p, g, x);

  p = 61;
  g = 2;
  x = 24;
  SHANKS(p, g, x);

  p = 43;
  g = 3;
  x = 5;
  SHANKS(p, g, x);

  p = 31;
  g = 3;
  x = 7;
  SHANKS(p, g, x);

  x = 1;
  for(int i = 1; i < 31; ++i) {
    SHANKS(p, g, x);
    x *= g;
    x = x % p;
  }

  // 2 non è un generatore
  p = 31;
  g = 2;
  x = 6;
  SHANKS(p, g, x);
}