/*
  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 dell'intero "n": Algoritmo di
  Fermat-Lehman.
  Si veda il Paragrafo 6.5.2 del libro citato.
*/

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

/*
  Variabile per il conteggio del numero di iterazioni necessarie
  (serve per confrontare fra loro i vari algoritmi proposti)
*/
long long int
iterazioni;

/*
  Se l'intero "n" è divisibile per l'intero "m", restituisci l'esponente
  "k" per cui esiste un intero "a" tale che "n = m^k * a", dove "a" non è
  divisibile per "n". All'uscita della funzione, "n" vale "a".
*/

long long int
prova_il_fattore(long long int& n, const long long int m) {

  ++iterazioni;
  int k = 0;
  while (n % m == 0) {
    ++k;
    n /= m;
  }
  return k;
}

/*
  Prova a dividere "n" per tentativi fino ad "m".
  Stampa gli eventuali fattori primi trovati con i loro esponenti
  (memorizzandoli nei vettori "basi" ed "esponenti" rispettivamente).
  Restituisci la parte di "n" eventualmente non fattorizzata.
*/

long long int
dividi_per_tentativi(long long int nn, long long int m,
		     std::vector<long long int>& basi,
		     std::vector<int>& esponenti) {

  long long int n = nn;
  if (n < 0)
    n *= -1;
  // Se "n" è pari, rimuovi la giusta potenza di 2
  int k = prova_il_fattore(n, 2);
  if (k > 0) {
    basi.push_back(2);
    esponenti.push_back(k);
    std::cout << " 2 ";
    if (k > 1)
      std::cout << "^ " << k;
    std::cout << "* ";
  }
  // Verifica l'eventuale divisibilità di "n" per i primi dispari
  // "piccoli", cioè minori di "m"
  for (long long int i = 3; i <= m; i += 2) {
    k = prova_il_fattore(n, i);
    if (k > 0) {
      basi.push_back(i);
      esponenti.push_back(k);
      std::cout << i;
      if (k > 1)
	std::cout << "^ " << k;
      std::cout << " * ";
    }
    // Se "i" è grande, esci dal ciclo
    if (i * i > n)
      i = m;
  }
  return(n);
}

/*
  Algoritmo di Lehman
*/

void
lehman(long long n) {

  // Prima fase: divisione per tentativi fino ad "n^(1/3)"
  std::vector<long long int> basi;
  std::vector<int> esponenti;

  std::cout << std::endl << "Numero da fattorizzare: " << n << std::endl;
  iterazioni = 0;
  long long int s = int(exp(log((double) n) / 3.0));
  long long m = dividi_per_tentativi(n, s, basi, esponenti);

  std::cout << std::endl;
  if (m > 1)
    std::cout << "Parte non fattorizzata: " << m << std::endl;
  else
    std::cout << "Il numero dato è stato completamente fattorizzato" << std::endl;
  std::cout << "Divisione per tentativi: " << iterazioni << " iterazioni" << std::endl;

  // Seconda fase: scrivi "n" nella forma "x^2- y^2"
  n = m;
  long long R = int(exp(log((double) n)/3.0));

  for(long long k = 1, r = 1; k < R + 1; ++k) {
    long long a = 1 + 2 * radice_quadrata(k * n);
    long long d = radice_quadrata(R / k) / 4;
    for(long long x = a; x < a + d + 1 && r > 0; ++x) {
      ++iterazioni;
      long long y = radice_quadrata(x * x - 4 * k * n);
      long long z = x * x - y * y - 4 * k * n;
      if (z == 0) {
	r = 0;
	a = mcd(x + y, n);
	long long b = mcd(x - y, n);
	std::cout << "k = " << k << "  x = " << x << "  y = " << y;
	std::cout << std::endl;
	std::cout << "a = " << a << "  b = " << b << std::endl;
      }
    }
  }
  std::cout << "Numero di iterazioni (totale): " << iterazioni << std::endl;
}

int main() {

  // Come numero di prova usiamo il "numero di Jevons" definito a p. 175
  long long int n = 8616460799LL;
  lehman(n);

  // Facciamo un'altra prova con un numero con molti fattori
  n = 510510000LL;
  lehman(n);

  // Un altro numero interessante
  n = 408704709982LL;
  lehman(n);
}