/*
  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.

  Scomposizione in fattori primi di un intero "n"
  Crivello quadratico di Pomerance
  Si veda il Paragrafo 6.6.1 del libro citato.
*/

#include <stdio.h>
#include <iostream>
#include <cmath>
#include "aritmetica.cc"

/*
  Il numero dei numeri primi precalcolati a disposizione
*/
static const
unsigned primi_prec = 40;

/*
  Il numero massimo dei numeri primi nella base di fattori
*/
static const
unsigned dim_bf = 15;

/*
  Il numero massimo di congruenze attese
*/
static const
unsigned max_num_congr = 5000;

/*
  Il numero di valori del polinomio ausiliario esplorati
*/
static const
unsigned max_num_val_poli = 15000;

/*
  Se "verbose = true", durante l'esecuzione vengono stampate
  informazioni dettagliate, e vengono creati alcuni files di dati.
*/
static const
bool verbose = false;

/*
  "primi_prec" numeri primi piccoli precalcolati per comodit�.
  In un'applicazione reale si dovrebbe usare il Crivello di Eratostene.
*/
static const
long long unsigned
primi[]={2,  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};

/*
  Alcuni file in cui vengono scritti i risultati delle varie fasi del
  calcolo.
*/
FILE *f_fb, *f_vec, *f_factors, *f_sol;

/*
  Prima del calcolo si verifica la condizione necessaria di Fermat
  "a^{n-1} = 1 mod n" per vari valori di "a".
  Se "n" soddisfa tutte queste condizioni allora � un "probabile primo".
  Pi� precisamente, "n" � uno pseudoprimo in tutte le basi prime
  memorizzate in "primi[]"; quindi � probabilmente opportuno che sia
  sottoposto ad un criterio di primalit� prima di tentare una sua
  scomposizione in fattori.
  NB: Le verifiche nelle basi composte sono superflue (se "n" �
  pseudoprimo in base 2 allora � pseudoprimo in base 4; se "n" �
  pseudoprimo in base 2 e 3, allora � pseudoprimo anche in base 6,
  ecc.) Viceversa, se "n" _non_ � pseudoprimo in base 2 o in base 3,
  il calcolo si interrompe senza raggiungere il valore j = 4, ecc.

  La funzione restituisce
  - "false" se trova un valore "a" nell'elenco dato sopra per cui
     a^(n-1) != 1 mod n,
  - "true" se a^(n-1) = 1 _per tutti_ gli "a" di cui sopra, cio� se
    "n" � uno pseudoprimo in tutte le basi esaminate.
*/
bool
criterio_di_fermat(long long int n) {

  for (unsigned j = 0; j < primi_prec; ++j) {
    long long unsigned p = primi[j];
    long long i = potenze_mod_n(p, n - 1, n);
    // Se i != 1 allora "n" non � primo
    if (i != 1) {
      std::cout << "N non supera il criterio di Fermat in base " << p;
      std::cout << ": " << p << "^" << n - 1 << " = " << i;
      std::cout << " mod "<< n << std::endl;
      return(false);
    }
  }
  // "n" � un probabile primo
  return(true);
}

/*
  Determinazione della base di fattori, cio� dei primi "p" "piccoli"
  (quelli memorizzati nell'array "primi[]") per cui il simbolo di
  Legendre (n | p) = 1, cio� per cui l'equazione "x^2 = n mod p" ha
  soluzione.
  La funzione restituisce il numero di elementi della base di fattori.
  Le soluzioni dell'equazione "Q(A) = n mod p" (quelle che determinano
  le classi di congruenza da esplorare) sono memorizzate in
  "soluzioni[][]". I primi nella base di fattori sono copiati
  nell'array "primi_bf[]".
*/
long long unsigned
base_di_fattori(long long unsigned n, long long unsigned s,
		long long unsigned primi_bf[primi_prec],
		long long unsigned soluzioni[primi_prec][2]) {

  // Il primo 2 appartiene sempre alla base di fattori
  primi_bf[0] = 2;
  // "k" contiene il numero di elementi della base di fattori trovati
  // finora
  long long unsigned k = 1;
  for (unsigned i = 1; i < primi_prec; ++i) {
    long long unsigned p = primi[i];
    long long unsigned m = n % p;
    // Data la simmetria delle soluzioni di "x^2 = n mod p" quando "p"
    // � dispari, � sufficiente esaminare i valori di "x"
    // nell'intervallo [0, (p - 1) / 2]
    for (unsigned j = 0; 2 * j < p ; ++j)
      if ((j*j) % p == m) {
	// Se si entra qui vuol dire che "n" � un residuo quadratico
	// modulo "p", e quindi "p" fa parte della "base di fattori".
	// Le due soluzioni di "Q(A) = n mod p" si possono calcolare a
	// partire dalle soluzioni di "x^2 = n mod p".  Bisogna
	// accertarsi che la soluzione calcolata cada nell'intervallo
	// [0, p - 1], e per questo sono necessarie le operazioni di
	// modulo.
	soluzioni[k][0] = (p - ((s - j) % p)) % p;
	soluzioni[k][1] = (p - ((j + s) % p)) % p;
	// Memorizza il k-esimo elemento della base di fattori e
	// incrementa "k"
	primi_bf[k] = p;
	++k;
      }
  }
  if (verbose) {
    for (unsigned i = 0; i < k; ++i)
      printf("%5lld", primi_bf[i]);
    std::cout << std::endl;
  }

  // C'� una sola classe da esplorare modulo 2 durante la fase di
  // crivello, perch�, se "N" � dispari, l'equazione "Q(A) = 0 mod 2"
  // diventa "(A + s) = 1 mod 2", la cui soluzione � "A = (s + 1) mod 2".
  soluzioni[0][0] = (s + 1) % 2;
  soluzioni[0][1] = (s + 1) % 2;
  if (verbose) {
    fprintf(f_sol, "&       &");
    for (unsigned i = 0; i < dim_bf; ++i)
      fprintf(f_sol, "&  %3lld  &", primi_bf[i]);
    fprintf(f_sol, "\\cr\n");
    for (int j = 0; j < 2; ++j) {
      fprintf(f_sol, "&       &");
      for (unsigned i = 0; i < dim_bf; ++i)
	fprintf(f_sol, "&  %3lld  &", soluzioni[i][j]);
      fprintf(f_sol, "\\cr\n");
    }
  }
  return(k);
}

/*
  Stampa le fattorizzazioni in varie forme in vari files ausiliari
*/
void
stampa_fattori(long long int a, long long int b, unsigned k,
	       long long int c, long long unsigned primi_bf[primi_prec],
	       long long unsigned vettore[max_num_congr][dim_bf]) {

  fprintf(f_fb,      "A=%6lld; Q(A)=%11lld =", a, b);
  fprintf(f_factors, "& %6lld && %11lld &&", a, b);
/*    fprintf(f_vec,     "%5lld ", c + 1); */
  fprintf(f_vec,     "&%6lld &&%11lld && (", a, b);
  bool flag = false;
  for (unsigned i = 0; i < k; ++i) {
    if (vettore[c][i] > 0) {
      if (flag) {
	fprintf(f_fb, " *");
	fprintf(f_factors, "\\cdot");
      }
      else
	flag = true;
      fprintf(f_fb, "%3lld", primi_bf[i]);
      fprintf(f_factors, "%3lld", primi_bf[i]);
      if(vettore[c][i] > 1) {
	fprintf(f_fb,      "^%1lld", vettore[c][i]);
	fprintf(f_factors, "^%1lld", vettore[c][i]);
      }
      else {
	fprintf(f_fb, "  ");
	fprintf(f_factors, "  ");
      }
    }
    fprintf(f_vec,"%1lld", vettore[c][i]%2);
    if (i < k - 1)
      fprintf(f_vec, ",");
  }
  fprintf(f_fb, "\n");
  fprintf(f_factors, " & \\cr\n");
  // fprintf(f_vec," -- %5lld %3lld %3lld", c, wt[c][0], wt[c][1]);
  fprintf(f_vec,") & \\cr\n");
}

void
somma_mod_2(long long int a, long long int b, long long int c,
	    long long unsigned vettore_mod_2[max_num_congr][dim_bf]) {

  for (long i = 0; i < c; ++i) {
    long j = vettore_mod_2[a][i] + vettore_mod_2[b][i];
    vettore_mod_2[b][i] = j % 2;
  }
}

void
somma(long long unsigned a, long long unsigned b, long long unsigned c,
      long long unsigned vettore[max_num_congr][dim_bf]) {

  for (unsigned i = 0; i < c; ++i)
    vettore[b][i] = vettore[a][i] + vettore[b][i];
}

void
somma_bis(long long unsigned a, long long unsigned b, long long unsigned c,
	  unsigned dep_bis[max_num_congr][200]) {

  for (unsigned i = 0; i < c; ++i)
    dep_bis[b][i] = (dep_bis[a][i] + dep_bis[b][i]) % 2;
}

/*
  Funzione ausiliaria
*/
void
get_wt(long long unsigned a, long long unsigned b,
       long long unsigned vettore_mod_2[max_num_congr][dim_bf],
       long long unsigned wt[max_num_congr][2]) {

  wt[a][0] = 0;
  wt[a][1] = b;
  for (long i = b - 1; i > -1; --i) {
    wt[a][0] += vettore_mod_2[a][i];
    if (vettore_mod_2[a][i])
      wt[a][1] = i;
  }
}

/*
  Questa funzione realizza il crivello.
  Sia "Q(A) = (A + s)^2 - N", dove "N" � il numero da scomporre in
  fattori, ed "s" � la sua radice quadrata intera (approssimata per
  difetto).
  "num_primi" � il numero di elementi nella base di fattori.
  1. Si memorizzano i valori di "Q(A)" per "A = 0", ...,
     "max_num_val_poli - 1" nel vettore "valori_pol_aus[]".
  2. Per ogni primo "p" nella "base di fattori" determinata
     dall'apposita funzione, si esplorano le 2 classi di congruenza
     mod p (una sola se "p = 2"), memorizzate nel vettore
     "soluzioni[][]" per cui "Q(A) = 0 mod p" se e solo se "A" � in
     una delle due classi, e si divide ripetutamente "Q(A)" per "p",
     finch� ci� � possibile, memorizzando l'esponente di "p"
     corrispondente nel vettore "fattor_pol_aus[][]".
  3. I valori di "A" per cui "Q(A)" si fattorizza completamente sulla
     base di fattori, e le corrispondenti fattorizzazioni di "Q(A)",
     sono memorizzati per essere usati nella fase della ricerca delle
     dipendenze lineari.
  La funzione ritorna il numero di fattorizzazioni complete determinate.
*/
long long unsigned
crivello(long long unsigned n, long long unsigned s,
	 long long unsigned num_primi,
	 long long unsigned primi_bf[primi_prec],
	 long long unsigned soluzioni[primi_prec][2],
	 long long unsigned vettore_mod_2[max_num_congr][dim_bf],
	 long long unsigned wt[max_num_congr][2],
	 long long unsigned vettore[max_num_congr][dim_bf],
	 long long unsigned valori[max_num_congr],
	 long long unsigned dep[max_num_congr][2],
	 unsigned dep_bis[max_num_congr][200]) {

  // Vettore in cui memorizziamo i valori del polinomio ausiliario
  long long valori_pol_aus[max_num_val_poli] = {0};

  // Calcolo dei valori di "Q(i)"
  long long unsigned i;
  for (i = 0; i < max_num_val_poli; ++i) {
    valori_pol_aus[i] = s * s - n + i * (i + 2 * s);
  }

  // Vettore in cui memorizziamo le fattorizzazioni degli interi nel
  // vettore "valori_pol_aus[]" rispetto alla base di fattori.
  unsigned fattor_pol_aus[max_num_val_poli][dim_bf] = {0};

  for (unsigned j = 0; j < num_primi; ++j) {
    unsigned p = primi_bf[j];
    // Vi sono 2 classi di congruenza da esplorare per tutti i primi
    // dispari, ma una sola per "p = 2"
    unsigned num_classi = 2;
    if (p == 2)
      num_classi = 1;
    for (unsigned k = 0; k < num_classi; ++k) {
      for (i = soluzioni[j][k]; i < max_num_val_poli; i += p) {
	long long unsigned n = valori_pol_aus[i];
	while (n % p == 0) {
	  ++fattor_pol_aus[i][j];
	  n /= p;
	}
	valori_pol_aus[i] = n;
      }
    }
  }
  // Qui termina il crivello

  if (verbose) {
    std::cout << "Per questi valori il polinomio si fattorizza";
    std::cout << " completamente sulla base di fattori";
  }
  unsigned num_fatt_compl = 0;
  for (i = 0; i < max_num_val_poli; ++i)
    if (valori_pol_aus[i] == 1) {
      if (verbose)
	std::cout << i << "  ";
      // L'obiettivo di questa parte del codice � costruire la
      // matrice che contiene gli esponenti (rispetto ai primi nella
      // base di fattori) delle fattorizzazioni dei valori di "Q(A)"
      // che si fattorizzano completamente sulla base di fattori, e
      // un'altra matrice che contiene gli stessi esponenti ridotti
      // modulo 2. � su quest'ultima matrice che verr� eseguita
      // l'eliminazione di Gauss.

      // Per l'eliminazione di Gauss sono utili due ulteriori
      // quantit� per ciascuna riga della matrice costruita come
      // detto sopra:
      // - il suo "peso" in bit, cio� il numero di bit = 1,
      //   memorizzato in "wt[][0]";
      // - la posizione del bit pi� a sinistra in ciascuna riga,
      //   memorizzato in "wt[][1]".
      // Per questo motivo, ciascuna riga � scorsa da destra verso
      // sinistra e ad ogni bit = 1 viene incrementato il valore di
      // "wt[][0]" e aggiornato il valore di "wt[][1]".
      wt[num_fatt_compl][1] = num_primi;
      for (int j = num_primi - 1; j > -1; --j) {
	long long unsigned l             = fattor_pol_aus[i][j] % 2;
	vettore      [num_fatt_compl][j] = fattor_pol_aus[i][j];
	vettore_mod_2[num_fatt_compl][j] = fattor_pol_aus[i][j] % 2;
	wt[num_fatt_compl][0] += l;
	if (l > 0)
	  wt[num_fatt_compl][1] = j;
	//    printf("%6lld %6lld\n", i, num_fatt_compl);
      }

      // Queste quantit� serviranno per costruire la congruenza
      // "X^2 = Y^2 mod n" alla fine del crivello quadratico
      valori[num_fatt_compl] = i;
      dep[num_fatt_compl][0] = i + s;
      long long unsigned q = s * s - n + i * (i + 2 * s);
      dep[num_fatt_compl][1] = q;
      dep_bis[num_fatt_compl][num_fatt_compl] = 1;
      //   printf("%8lld  %8lld\n", i, num_fatt_compl);
      if (verbose)
	stampa_fattori(i, q, num_primi, num_fatt_compl, primi_bf, vettore);
      ++num_fatt_compl;
    }
  if (verbose)
    std::cout << std::endl;
  std::cout << "Ho trovato " << num_fatt_compl << " fattorizzazioni complete";
  std::cout << std::endl << std::endl;
  return(num_fatt_compl);
}

/*
  Questa funzione realizza il crivello quadratico. Si noti che, nel
  caso sia in grado di produrre un fattore di "N", non c'� garanzia
  che questo fattore sia un numero primo, e quindi � opportuno
  chiamare ricorsivamente la stessa funzione sui fattori primi, fino a
  determinarne la fattorizzazione completa.
*/
bool
crivello_quadratico(long long unsigned N) {

  std::cout << "Numero da fattorizzare: N = " << N << std::endl;
  // 1. Ricerca di eventuali fattori primi "piccoli" di "N"
  std::cout << "Passo 1: ricerca di eventuali fattori primi piccoli";
  std::cout << std::endl;
  for (unsigned i = 0; i < primi_prec; ++i) {
    long long unsigned p = primi[i];
    if (N % p == 0) {
      unsigned alpha = 0;
      while (N % p == 0) {
	N /= p;
	++alpha;
      }
      std::cout << "Ho trovato il fattore " << p;
      if (alpha > 1)
	std::cout << "^" << alpha;
      std::cout << std::endl;
    }
  }
  long long unsigned p = primi[primi_prec - 1];
  // Se abbiamo gi� rimosso tanti fattori primi, pu� darsi che la
  // parte non fattorizzata di "N" sia un numero primo, ed allora non
  // � necessario eseguire il resto del programma.
  if (p * p > N) {
    if (N > 1)
      std::cout << "Ho trovato il fattore " << N << std::endl;
    std::cout << "Fattorizzazione completata" << std::endl;
    return(true);
  }

  std::cout << "Numero da fattorizzare: N = " << N << std::endl;
  long long s = radice_quadrata(N);
  std::cout << "Radice quadrata:        s = " << s << std::endl << std::endl;

  // 2. Verifica della condizione necessaria per la primalit�
  std::cout << "Passo 2: verifica della condizione di Fermat" << std::endl;
  bool fermat = criterio_di_fermat(N);
  if (fermat)
    std::cout << N << " � un probabile primo" << std::endl;
  std::cout << std::endl;

  // 3. Determinazione dei primi piccoli "p" per cui "( N | p ) = 1",
  // cio� dei primi "p" per i quali l'equazione "x^2 = N mod p" ha
  // soluzione
  std::cout << "Passo 3: determinazione della base di fattori" << std::endl;

  // "primi_prec" coppie di eventuali soluzioni di "x^2 = N mod p"
  long long unsigned soluzioni[primi_prec][2] = {0};

  // I numeri primi nella "base di fattori".
  long long unsigned primi_bf[primi_prec] = {0};

  long long unsigned num_primi = base_di_fattori(N, s, primi_bf, soluzioni);
  std::cout << "La base di fattori contiene " << num_primi << " numeri primi";

  // Permettiamo al massimo "dim_bf" numeri primi nella base di
  // fattori per due motivi: occupazione di memoria, e numero di
  // dipendenze lineari necessarie
  if (num_primi > dim_bf) {
    num_primi = dim_bf;
    std::cout << ". Consideriamo solo i primi " << dim_bf;
  }
  std::cout << std::endl << std::endl;

  // Array in cui si memorizza la scomposizione di "Q(A)" in fattori
  // primi appartenenti alla base di fattori.
  long long unsigned vettore[max_num_congr][dim_bf] = {0};

  // Valori di "A" per cui "Q(A)" si fattorizza completamente sulla
  // base di fattori.
  long long unsigned valori[max_num_congr] = {0};

  // Array in cui si memorizza la scomposizione di "Q(A)" in fattori
  // primi appartenenti alla base di fattori, con esponenti ridotti
  // mod 2
  long long unsigned vettore_mod_2[max_num_congr][dim_bf] = {0};

  // Qui sono contenute le dipendenze lineari trovate alla fine del
  // crivello quadratico.
  long long unsigned dep[max_num_congr][2] = {0};

  unsigned dep_bis[max_num_congr][200] = {0};

  // Matrici ausiliarie in uso nella ricerca delle dipendenze lineari
  // - wt[i][0] contiene il peso in bit (il numero di bit di valore 1)
  //   dell'i-esima riga;
  // - w[i][1] la posizione del bit = 1 pi� a sinistra nell'i-esima
  //   riga.
  long long unsigned wt[max_num_congr][2] = {0};

  // 4. Crivello
  std::cout << "Passo 4: il crivello" << std::endl;
  long long unsigned fb = crivello(N, s, num_primi, primi_bf, soluzioni,
				   vettore_mod_2, wt, vettore, valori, dep,
				   dep_bis);

  // A questo punto
  // - "dep[j][0]" contiene il (j+1) esimo valore di "A"
  //   per cui "Q(A)" si fattorizza sulla base di fattori;
  // - "dep[j][1]" contiene "Q(A)";
  // - "wt[j][0]" contiene il "peso" della fattorizzazione di "Q(A)",
  //   cio� il numero di esponenti dispari nella fattorizzazione;
  // - "wt[j][1]" contiene la posizione del bit pi� a sinistra = 1.

  // 5. algebra lineare: eliminazione di Gauss
  std::cout << "Passo 5: eliminazione di Gauss" << std::endl;
  for (unsigned i = 0; i < num_primi; ++i) {
    unsigned j;
    // Cerchiamo il primo vettore che abbia il bit pi� a sinistra
    // esattamente nella posizione i-esima
    for (j = 0; (j < fb) & (wt[j][1] != i); ++j);
    // Usiamo questo vettore per eliminare un eventuale bit = 1 nella
    // posizione i-esima da tutti i vettori successivi. Teniamo
    // traccia di queste eliminazioni (e quindi delle dipendenze
    // lineari che troviamo) aggiornando i valori degli array
    // "dep[][]"
    for (unsigned k = j + 1; k < fb; ++k) {
      if (vettore_mod_2[k][i]) {
	somma_mod_2(j, k, num_primi, vettore_mod_2);
	somma(j, k, num_primi, vettore);
	somma_bis(j, k, fb, dep_bis);
	dep[k][0] = moltiplicazione_mod_n(dep[j][0], dep[k][0], N);
	dep[k][1] = moltiplicazione_mod_n(dep[j][1], dep[k][1], N);
	get_wt(k, num_primi, vettore_mod_2, wt);
      }
    }
  }

  long unsigned successi = 0;
  // I vettori che sono stati "eliminati" hanno peso = 0
  for (unsigned i = 0; i < fb; i++)
    if (wt[i][0] == 0) {
      long long j = 1;
      long long l;
      long long unsigned m;
      for (unsigned k = 0; k < num_primi; ++k) {
	l = vettore[i][k] / 2;
	m = potenze_mod_n(primi_bf[k], l, N);
	j = moltiplicazione_mod_n(j, m, N);
      }
      long long X = dep[i][0];
      long long Y = j;
      printf("** %12lld^2 = %12lld^2 mod N;", X, Y);
      // 5. calcolo del MCD
      m = mcd(X + Y, N);
      l = N / m;
      printf("  N = %12lld * %12lld\n", m, l);
      if ((m > 1) & (m < N))
	++successi;
    }

  std::cout << "Fattorizzazione di N = " << N;
  std::cout << ". Numero di successi: " << successi << std::endl;
  if (verbose) {
    fprintf(f_vec, "\n\n");
    for (unsigned i = 0; i < fb; ++i)
      if (wt[i][0] == 0) {
	printf("%4u --- ", i);
	fprintf(f_vec, "& ");
	for (unsigned k = 0; k < fb; ++k)
	  if (dep_bis[i][k] > 0) {
	    printf("%3u,", k);
	    fprintf(f_vec, "\\vec v(%2lld) + ", valori[k]);
	  }
	printf("\n");
	fprintf(f_vec, "\\equiv 0 \\bmod 2 \\\\\n");
      }
  }
  std::cout << "Fine dell'esecuzione" << std::endl;
  return (successi > 0);
}

int
main() {

  if (verbose) {
    f_fb      = fopen("FactorBase.dat", "w");
    f_factors = fopen("Fattorizzazioni.dat", "w");
    f_sol     = fopen("Soluzioni.dat", "w");
    f_vec     = fopen("Vettori.dat", "w");
  }

/*
  Alcuni interi interessanti con cui provare il codice: l'asterisco
  indica i valori del parametro "dim_bf" per cui il crivello
  quadratico ha successo (e cio� scompone "N" in due fattori non
  banali, ma non necessariamente primi), fermi restando i valori di
  tutti gli altri parametri.

  N            = fattorizzazione completa                     dim_bf
                                                           10   15   20
  65537        = primo
  654079       = 593 * 1103                                 *    *    *
  67108865     = 5 * 53 * 157 * 1613                        *    *    *
  220424623    = 337 * 593 * 1103                           *    *    *
  536813567    = 8191 * 65537                                         *
  1073741825   = 5^2 * 13 * 41 * 61 * 1321 = 2^{30} + 1     *    *    *
  4294967297LL = 641 * 6700417 = 2^{32} + 1                 *    *    *
  8616460799LL = 89681 * 96079                                   *    *
  9049465682LL = 2 * 28879 * 156679                         *    *    *

  408704709982LL = 2 * 399601 * 511391
  Per quest'ultimo intero � necessario porre dim_bf = 40,
  max_num_val_poli = 40000.
*/

  long long unsigned test[] = {65537, 654079, 67108865, 220424623,
			       536813567, 1073741825, 4294967297LL,
			       8616460799LL, 9049465682LL,
			       408704709982LL};

  for (unsigned i = 0; i < sizeof(test) / sizeof(long long); ++i) {
    long long unsigned N = test[i];
    bool status = crivello_quadratico(N);
    if (!status) {
      std::cout << "Il crivello quadratico non ha trovato fattorizzazioni";
      std::cout << std::endl;
    }
    std::cout << "Premi un tasto per continuare...";
    char c;
    std::cin >> c;
    std::cout << std::endl;
  }

  if (verbose) {
    fclose(f_fb);
    fclose(f_factors);
    fclose(f_sol);
    fclose(f_vec);
  }
}