/*
  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 Gauss per la determinazione di un generatore del gruppo
  Z_p^*, dove "p" è un numero primo.
  Si veda il Paragrafo 6.7 del libro citato.
*/

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

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

/*
  Determina l'ordine di "a" modulo "p", memorizzando le potenze
  successive di "a" nel vettore "potenze_di_a".
  Alla fine si ordina il vettore "potenze_di_a" per semplificare la
  ricerca di un elemento di Z_p^* che non sia una potenza di "a" (vedi
  sotto).
  La funzione ritorna l'ordine di "a" in Z_p^*.
*/
long long
calcola_potenze_mod_p(long long p, long long a,
		      std::vector <long long> &potenze_di_a) {

  long long potenza = 1;
  long long r_a = 0;
  do {
    // Appendi la nuova potenza di "a" in fondo al vettore
    // "potenze_di_a"
    potenze_di_a.push_back(potenza);
    potenza = moltiplicazione_mod_n(potenza, a, p);
    ++r_a;
  } while (potenza != 1);
  // Ordina le potenze di "a", ridotte modulo "p", contenute nel
  // vettore "potenze_di_a". NB. È meglio usare la funzione "sort"
  // del linguaggio perché è sicuramente piú efficiente
  sort(potenze_di_a.begin(), potenze_di_a.end());

  return(r_a);
}

/*
  Decomponi "v = MCM(r_a, r_b) = r_a * r_b / mcd(r_a, r_b)"
  come "v = n * m", dove
  1. n divide r_a;
  2. m divide r_b;
  3. mcd(n, m) = 1.
  Strategia: sia d = mcd(r_a, r_b), e sia p un fattore primo di d.
  Se p^a || r_a, e p^b || r_b:
    se a >= b assegna il fattore p^a ad n;
    se a <  b assegna il fattore p^b ad m.
  Gli altri fattori primi di "v" dividono uno solo fra "r_a" ed "r_b" e
  vengono assegnati ad "n" o "m" rispettivamente.
*/

long long
decomponi(long long r_a, long long r_b, long long &n, long long &m) {

  long long d = mcd(r_a, r_b);
  // "v = mcm(r_a, r_b)" è il valore ritornato dalla funzione
  long long v = r_a * r_b / d;
  n = 1;
  m = 1;

  // Molto inefficiente: sostituire con qualcosa di meglio
  for(long long p = 2; p <= d; ++p) {
    if ((d % p) == 0) {
      // "p" è un fattore primo di "d":
      // determina la potenza "m_a" con cui "p" divide "r_a"
      long long m_a = 1;
      while ((r_a % p) == 0) {
	r_a /= p;
	m_a *= p;
      }
      // Determina la potenza "m_b" con cui "p" divide "r_b"
      long long m_b = 1;
      while ((r_b % p) == 0) {
	r_b /= p;
	m_b *= p;
      }
      // Elimina il primo "p" dal numero "d"
      while ((d % p) == 0)
	d /= p;
      // Aggiorna il valore di "n" o di "m" con la potenza corretta
      // del primo "p"
      if (m_a >= m_b)
	n *= m_a;
      else
	m *= m_b;
    }
  }
  // A questo punto "mcd(r_a, r_b) = 1":
  // assegna i fattori primi residui ad "n" o "m"
  // NB: non è necessario scomporre "r_a", "r_b" in fattori
  n *= r_a;
  m *= r_b;
  return(v);
}

/*
  Algoritmo di Gauss per la determinazione di un generatore mod "p"
  N. B.: Questo algoritmo NON determina necessariamente il minimo
  intero nell'intervallo [2, p-1] che genera Z_p^*
*/
long long
gauss(long long p) {

  // Se "p = 2", allora 1 è un generatore
  if (p == 2)
    return(1);
  // Da qui in poi, "p" è un numero primo dispari

  // Cominciamo sempre dal valore "a = 2" perché conviene lavorare
  // con interi piccoli: qui sfruttiamo il fatto che "p" è dispari
  long long a = 2;
  long long v = 1;

  // Ripetiamo il ciclo fino a quando non si trova un intero di ordine
  // "p - 1".
  // NB: se "p" non è un numero primo, questo ciclo viene ripetuto
  // infinite volte perché nessun elemento di Z_p^* ha ordine "p - 1".
  while (v < p-1) {
    std::vector <long long> potenze_di_a;
    // "r_a" è l'ordine di "a" in Z_p^*
    long long r_a = calcola_potenze_mod_p(p, a, potenze_di_a);
    // Se "a" è un generatore abbiamo finito
    if (r_a == p-1)
      return(a);

    long long b = 0;
    // Cerca il primo intero "b > 1" che non sia una potenza di "a".
    // Il vettore "potenze_di_a[]" è stato riempito con le potenze
    // successive di "a", ridotte modulo "p", e poi ordinato, nella
    // funzione "calcola_potenze_mod_p()".
    while (potenze_di_a[b] == b + 1) 
      ++b;
    ++b;
    // A questo punto "b" è un elemento di Z_p^* che non è potenza
    // di "a": in effetti, è il piú piccolo intero nell'intervallo
    // [1, p-1] con questa proprietà.
    // Stiamo sfruttando il fatto che il vettore "potenze_di_a" è
    // ordinato.

    std::vector <long long> potenze_di_b;
    // "r_b" è l'ordine di "b" in Z_p^*
    long long r_b = calcola_potenze_mod_p(p, b, potenze_di_b);
    // Se "b" è un generatore abbiamo finito
    if (r_b == p-1)
      return(b);

    // Se si raggiunge questo punto, né "a" né "b" sono generatori.
    // Da qui in poi si determina un altro elemento di Z_p^* con
    // ordine piú grande sia di "a" che di "b".
    long long n, m;
    v = decomponi(r_a, r_b, n, m);
    long long e_a = r_a / n;
    long long e_b = r_b / m;
    long long nuovo_a = potenze_mod_n(a, e_a, p);
    long long temp = potenze_mod_n(b, e_b, p);
    nuovo_a = moltiplicazione_mod_n(nuovo_a, temp, p);

    // "nuovo_a" ha ordine "v > max(r_a, r_b)"
    a = nuovo_a;
    // Se il nuovo valore di "a" è un generatore abbiamo finito
    if (v == p-1)
      return(a);
    // In caso contrario, continuiamo con il ciclo while, con un nuovo
    // valore di "a", avente un ordine piú grande dell'ordine del
    // valore precedente di "a".
  }
  // Questa istruzione dovrebbe essere irraggiungibile, ma è
  // richiesta dal compilatore.
  return(-1);
}

int main() {

  long long unsigned
    primes[] = {3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, \
		59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, \
		113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, \
		179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, \
		239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, \
		307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, \
		373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, \
		439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, \
		503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, \
		587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, \
		647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, \
		727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, \
		809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, \
		877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, \
		953, 967, 971, 977, 983, 991, 997
    };

  // Come prova, calcoliamo i generatori modulo i primi fino a 1000
  long long p;
  for (unsigned i = 0; i < sizeof(primes)/sizeof(p); ++i) {
    p = primes[i];
    std::cout << "Gauss (" << p << "):   --->    " << gauss(p) << std::endl;
  }

  // Calcoliamo generatori modulo altri primi 
  p = 65537;
  std::cout << "Gauss (" << p << "):   --->    " << gauss(p) << std::endl;

  p = 65543;
  std::cout << "Gauss (" << p << "):   --->    " << gauss(p) << std::endl;

  p = 999983;
  std::cout << "Gauss (" << p << "):   --->    " << gauss(p) << std::endl;
}