/* DiskreterLog_bsgs - Soll versuchen, diskrete Logarithmen zu "knacken"
   a ^ x = b modulo p   ; p ist primzahl, a ist Primitivwurzel
   Autor: J.Gamenik
   Beginn: Anfang 2008
   15.02.2012 vereinfacht, aus experimentellen Status geholt
   Aufruf: java DiskreterLog_bsgs <Basis> <Wert> <Primzahl>
   ACHTUNG: kann mit OutOfMemoryError abbrechen, weil Ergebnisses des Einzelschritts gespeichert werden
*/

import java.io.IOException;   // wegen "warten"
import java.math.BigInteger;  // für Primtest, modPow
import java.util.HashMap;  // speichern im Einzelschritt

/*
  relativ schwere Eingabewerte sind z.B.
2 7 8589942539
5 7 8589942863
5 7 8589942983
5 7 8589943127
13 5 8589944303

5 7 137438953427
*/

public final class DiskreterLog_bsgs extends Thread {

  private static int primBasis [] = {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}; // eine Hand voll Primzahlen
  private static int basisExp [];
  private static BigInteger phiPrimRest;  // falls Wert über Basis bis auf einen Primrest bestimmbar
  private BigInteger n;  // "p"
  private BigInteger phi;
  private long phi_long;
  private BigInteger basis;  // "a"
  private BigInteger wert;  // "b"

  // zur Anzeige schwerer Primzahlen ist folgendes benutzbar:
  static
  {
    // suche_primzahlen(25, 60);
  }

  // Konstruktor. Übernimmt Eingabewerte
  public DiskreterLog_bsgs (String args[])
  {
    if (args.length < 3) { 
      System.out.println ("Syntax: java DiskreterLog_bsgs Basis Wert Primzahl"); System.exit(1); }

    // Gleich Bedingungen prüfen:
    n = new BigInteger(args[2]); phi = n.subtract(BigInteger.ONE);
    phi_long = (phi.bitLength() < 64) ? phi.longValue() : -1;  // zur Ganzzahl-Überlaufbehandlung
    if (! n.isProbablePrime(2)) { System.out.println ("keine Primzahl"); System.exit(1); }
    wert = new BigInteger(args[1]);
    if (wert.compareTo(BigInteger.ONE) < 0 ||
        wert.compareTo(n) >= 0) { System.out.println ("Wert ausserhalb Bereich"); System.exit(1); }
    basis = new BigInteger(args[0]);
    if ((basis.compareTo(BigInteger.ONE) < 0) ||
        (basis.compareTo(n) >= 0)) { System.out.println ("Basis ausserhalb Bereich"); System.exit(1); }
    // Folgendes Kritierium ist notwendig, aber nicht hinreichend für Primitivwurzel
    if (PrimGenerator.jacobiSymbol(basis, n) != -1)
      { System.out.println ("Basis ist keine Primitivwurzel"); System.exit(1); }

    if (!ordnung_zerlegen()) { System.out.println ("Phi konnte nicht komplett zerlegt werden"); System.exit(1); }
    System.out.println ("Primzahl hat " + n.bitLength() + " Bits");

    this.setPriority(Thread.MIN_PRIORITY);  // damit man nebenher noch was machen kann
  } // Ende Konstruktor

  // wird über start() beim Thread aufgerufen
  public void run()
  {
    long rueck = -1;
    try
    {
      rueck = logmodulo();
    } catch (Exception ex) { System.out.println ("Fehler erkannt: " + ex); }
    if (rueck != -1) System.out.println("Logarithmus ist " + rueck);
  }
  
  // Hauptmethode zum Berechnen der Logarithmen. zerlegt phi, teilt also die Größe auf.
  // löst die (einfacheren) Teilprobleme und setzt das Ergebnis mit dem chinesichen Restsatz zusammen
  // Vorbedingung: Ordnung phi war komplett zerlegbar !
  public long logmodulo() throws Exception
  {
    // Annahme: basis a hat die Ordnung phi = p - 1, d.h. jeder Wert b ist erzeugbar.

    // Einfache Werte abtesten
    if (wert.equals(BigInteger.ONE)) return 0;
    else if (wert.equals(BigInteger.ZERO)) throw new Exception ("Kein Logarithmus verfügbar");
    else if (basis.equals(wert)) return 1;

    // Teilprobleme anhand der primBasis erzeugen
    long teilLoesung [] = new long[primBasis.length + 1]; int nrProb = 1;
    for (int i = 0; i < basisExp.length; ++i)
    {
      if (basisExp[i] == 0) continue;
      BigInteger maxOrdnung = BigInteger.valueOf(primBasis[i]).pow(basisExp[i]);
      BigInteger hoch = phi.divide(maxOrdnung);
      BigInteger basisNeu = basis.modPow(hoch, n);
      BigInteger wertNeu = wert.modPow(hoch, n);
      System.out.println (nrProb + ": " + basisNeu + " ^ x = " + wertNeu + " mod " + maxOrdnung);
      long rueck = logmodulo_bsgs(basisNeu.longValue(), wertNeu.longValue(), maxOrdnung.longValue());
      teilLoesung[i] = rueck;  ++nrProb;
    }
    if (phiPrimRest != null)
    {
      BigInteger hoch = phi.divide(phiPrimRest);
      BigInteger basisNeu = basis.modPow(hoch, n);
      BigInteger wertNeu = wert.modPow(hoch, n);
      System.out.println (nrProb + ": " + basisNeu + " ^ x = " + wertNeu + " mod " + phiPrimRest);
      long rueck = logmodulo_bsgs(basisNeu.longValue(), wertNeu.longValue(), phiPrimRest.longValue());
      teilLoesung[primBasis.length] = rueck;
    }

    // mit chinesischen Restesatz zusammensetzen aus Teillösungen
    long gesamtLoesung = 0; long moderg = 0; long anzLoesung = 0;
    for (int i = 0; i < teilLoesung.length; ++i)
    {
      boolean neu = false; long dazu = 0; long mod = 0;
      if (i >= basisExp.length)
      {
        if (phiPrimRest == null) break;
        neu = true; System.out.println("mod " + phiPrimRest + " --> " + teilLoesung[i]);
        dazu = teilLoesung[i]; mod = phiPrimRest.longValue();
      }
      else if (basisExp[i] != 0)
      {
        neu = true; System.out.println("mod " + primBasis[i] + " --> " + teilLoesung[i]);
        dazu = teilLoesung[i]; mod = BigInteger.valueOf(primBasis[i]).pow(basisExp[i]).longValue();
      }
      if (neu)
      {
        if (anzLoesung == 0) { gesamtLoesung = dazu; moderg = mod; }
        else 
        {
          long erg [] = chinRestesatz(gesamtLoesung, moderg, dazu, mod);
          if (erg != null) { gesamtLoesung = erg[0]; moderg = erg[1]; }
        }
        ++anzLoesung;
      }
    }
    System.out.print ("Errechnete Lösung: " + gesamtLoesung + " modulo " + moderg);

    // prüfen, ob tatsächlich passt
    if (basis.modPow(BigInteger.valueOf(gesamtLoesung), n).equals(wert)) 
    { System.out.println (" ... stimmt"); return gesamtLoesung; }
    else { System.out.println (" ... falsch"); return -1; }
  } 


  // Versucht, phi über der Faktorbasis zu zerlegen. Wenn dies nicht gelingt, kann im Programm
  // nicht richtig gerechnet werden. Ist einmalig pro Basis machbar, ist ja unabhängig von Werten b
  // liefert wahr, wenn Primzerlegung von Phi bekannt, sonst falsch.
  private boolean ordnung_zerlegen()
  {
    basisExp = new int[primBasis.length];
    BigInteger lauf = phi; phiPrimRest = null;
    int weg = lauf.getLowestSetBit();  // Sonderfall 2
    if (weg > 0)
    {
      basisExp[0] = weg; lauf = lauf.shiftRight(weg);
    }
    for (int i = 1; i < basisExp.length; ++i)
    {
      while (true) // solange, bis 1 oder zum Faktor nicht teilbar.
      {
        BigInteger teil = BigInteger.valueOf(primBasis[i]);
        BigInteger [] rueck = lauf.divideAndRemainder(teil);
        if (rueck[1].equals(BigInteger.ZERO))
        {
          basisExp[i]++; lauf = rueck[0]; if (lauf.equals(BigInteger.ONE)) break;
        }
        else break; // nicht teilbar
      } 
    }
    if (lauf.equals(BigInteger.ONE)) return true;  // fertig zerlegt
    // jetzt noch Hoffnung, wenn Quadrat (hier nicht geprüft) oder Primzahl
    if (lauf.isProbablePrime(2)) { phiPrimRest = lauf; return true; }

    System.out.println ("Ordnung hat hier nicht faktorisierbaren Rest " + lauf);
    return false; // Ordnung nicht zerlegt.
  }


  // Versucht, einen Logaritmus mit giant-step, baby-step zu finden: einfaches Suchverfahren
  // gedacht ist es, dass die Methode nur gerufen wird, wenn maxOrdnung prim ist (oder Primpotenz)
  private final long logmodulo_bsgs (long a, long b, long maxOrdnung) throws Exception
  {
    int i, j;
    BigInteger basisNeu = BigInteger.valueOf(a);
    BigInteger wertNeu = BigInteger.valueOf(b);
    // Einfache Werte abtesten
    if (wertNeu.equals(BigInteger.ONE)) return 0;
    else if (wertNeu.equals(BigInteger.ZERO)) throw new Exception ("Kein Logarithmus verfügbar");
    else if (basisNeu.equals(wertNeu)) return 1;

    // Sei Q die kleinste natürliche Zahl mit q*q >= maxOrdnung
    long q = (long) Math.sqrt(maxOrdnung);
    long summe = (q-1) * q + (q-1); 
    if (summe < 0) { System.out.println ("Interner Fehler: Überlauf"); System.exit(1); }
    while (summe < maxOrdnung) { ++q; summe = (q-1) * q + (q-1); }
    // Der gesuchte Exponent läßt sich damit schreiben k * q + l, also
    System.out.println ("q ist " + q);

    // Kleiner Schritt. Speichere (wert, index) in einer Map, am besten Hashmap. 
    // Anzahl Einträge ist genau bekannt. Der 0.75 ist der stand. Ladefaktor.
    long startZeit = System.currentTimeMillis();
    //BigInteger speicher [] = new BigInteger[(int)q];
    HashMap speicher = new HashMap((int)(q / 0.75 + 1));  // q ~ SQRT(maxPeriode), dies ist hart !
    BigInteger invers = basisNeu.modInverse(n);
    //speicher[0] = wertNeu;
    speicher.put (wertNeu,new Integer(0)); 
    BigInteger lauf = wertNeu;
    for (i = 1; i < q; ++i)
    {
      lauf = lauf.multiply(invers).remainder(n); // lauf hprof diese Stelle sehr teuer.
      //speicher[i] = lauf;
      speicher.put (lauf, new Integer(i)); // in map speichern
    }   

    // Großer Schritt
    lauf = basisNeu.modPow(BigInteger.valueOf(q), n); BigInteger aktuell = BigInteger.ONE;
    boolean gef = false;
    for (i = 0, j=0; i < q && !gef; ++i)
    {
      // suchen: schon da ? durch HashMap geht es hier schnell
      if (speicher.containsKey(aktuell)) { j = ((Integer) speicher.get(aktuell)).intValue(); gef = true; break; }
      //for (j = 0; j < q && !gef; ++j) if (speicher[j].equals(aktuell)) { gef = true; break; }
      if (gef) break;
      aktuell = aktuell.multiply(lauf).remainder(n);
    }
    System.out.println ("Zeitdauer " + (System.currentTimeMillis() - startZeit) / 1000.00);
    if (! gef) { System.out.println ("Es gibt keinen Logarithmus"); return -1; }
    else return (i * q + j);
  }


  // Realisiert den Chinesischen Restesatz, aber nur speziell für den Fall ggT(n,m) = 1
  // und speziell a, b beide >= 0. Macht Tests auf Zahlenüberlauf (z.B. pos. Zahlen gehen ins neg.)
  private final static long [] chinRestesatz (long a, long n, long b, long m) throws Exception
  {
    System.out.println ("Chin. Restesatz: " + a + " m " + n + ", " + b + " m " + m);
    // x == a mod n, x == b mod m

    if ((n == 0) || (m == 0) || (a < 0) || (b < 0)) throw new Exception ("Ungültiger Aufruf von chinRestesatz"); // Div durch 0 vermeiden

    long rueck [] = new long [2];  // um Rückgabewerte zu verpacken
    if (a >= n) a = a % n;       // Zahlen klein halten
    if (b >= m) b = b % m;
    // drei einfache Fälle, wo sich nichts verändert
    if (n == 1) { rueck[0] = b; rueck[1] = m; return rueck; }
    else if (m == 1) { rueck[0] = a; rueck[1] = n; return rueck; }
    else if (a == b) { rueck[0] = a; rueck[1] = n * m; return rueck; }  // Erkl. siehe Formel unten

    // wegen modinvers: Annahme n < m sicherstellen, evt. tauschen
    if (n > m)
    {
      long tmp = a; a = b; b = tmp;
      tmp = n; n = m; m = tmp;
    }
    long y = -1;
    try
    {
      y = BigInteger.valueOf(n).modInverse(BigInteger.valueOf(m)).longValue();   // n^-1 mod m
      if (y < 0) { System.out.println ("Überlauf bei Inversem"); System.exit(1); }
    }
    catch (Exception ex)  // ggT(n,m) > 1 hier als Ausnahme erkennen
    { System.out.println ("CRT nur speziell bei relativ prim"); System.exit(1); }
    // Rechenformel laut Lehrbuch:
    // rueck[0] = a - (n * y) * (a - b);  // 2) Gesuchte Zahl. Wird ggf. < 0, wenn (a > b)
    long teilMult = modmul(y ,a - b, m);
    rueck[0] = a - n * teilMult;  // teilMult < m
    rueck[1] = m * n;                  // 1) Eindeutigkeitsbereich
    if (rueck[1] < 0) { System.out.println ("Überlauf bei neuem Modul"); System.exit(1); }
    rueck[0] = rueck[0] % rueck[1];  
    if (rueck[0] < 0) rueck[0] = rueck[0] + rueck[1];  // soll immer positiv sein

    // Test:
    if (((rueck[0] - a) % n != 0) || ((rueck[0] - b) % m != 0))
    {
      System.out.println ("Fehler bei chinRestesatz: ermittelt " + rueck[0] + " in " + rueck[1]); 
      System.exit(1);
    }
    return rueck;
  }


  // finde schwere Primzahlen, d.h. solche, welche auch p- 1 prim sind. Dazu werden zwei Bitzahlen
  // angegeben, dann die nächst größere Bestimmt
  public static void suche_primzahlen (int b1, int b2)
  {
    if (b1 <= 0 || b1 >= 62 || b2 <= 0 || b2 >= 62)
    { System.out.println ("Bitzahlen für long nicht geeignet"); return; }

    final int stufe = 3;
    for (int laenge = b1; laenge <= b2; ++laenge)
    {
      long start = 1l << (laenge-1);
      // Suche nun nächste passende Primzahl
      if ((start & 1) == 0) ++start ; // immer ungerade
      for (long lauf = start; lauf <= Long.MAX_VALUE; lauf += 2)
      {
        if ((lauf % 3) == 0) continue;
        if ((lauf % 7) == 0) continue;
        BigInteger neu = BigInteger.valueOf(lauf);
        if (neu.isProbablePrime(stufe))
        {
          BigInteger zusatz = neu.shiftRight(1);
          if (zusatz.isProbablePrime(stufe)) { System.out.println(lauf); break; }
        }
      }
    }

  }

  // modulare Multipl. a und b sollen zwischen 0 und n-1 liegen
  private final static long modmul (long a, long b, long n)
  {
    long summe = 0;
    boolean neg = false;
    if (a < 0) { neg = !neg; a = -a; }
    if (b < 0) { neg = !neg; b = -b; }
    if (b < a) { long temp = a; a = b; b = temp; }
    while (a != 0)
    {
      if ((a <= Integer.MAX_VALUE) && (b <= Integer.MAX_VALUE))
      {
        summe += (a * b); summe = summe % n; break;  // kein Überlauf
      }
      if ((a & 1) == 1)
      {
        summe += b; if (summe >= n) summe -= n;
      }
      a >>= 1;
      b <<= 1; if (b >= n) b -= n;
    }
    return (!neg ? summe : -summe);
  }


  // Haupt-Methode
  public static void main (String[] args) throws Exception {
    DiskreterLog_bsgs bsgs = new DiskreterLog_bsgs (args);
    bsgs.start();   
  }


}