/* DiskreterLog_Faktor.java
   Soll versuchen, diskrete Logarithmen modulo einer Primzahl zu ermitteln. 
 
   Dazu werden zur (Grund)basis eine Reihe von Quadrat-Potenzen modulo berechnet, dann
   diese faktorisiert und in eine Tabelle eingetragen. Damit wird versucht, den
   diskreten Logarithmus zu zum Wert zu finden. Dies setzt natürlich voraus, dass auch
   der Wert oder eine bekannte kleine Potenz davon über der Primbasis faktorisierbar ist.
   Autor: J.Gamenik
   Beginn: 30.12.2011
   Aufruf: java DiskreterLog_Faktor <Basis> <Wert> <Primzahl> 
   z.B. 7 2 41 ergibt 14
     Basis ^ x == wert modulo Primzahl. Was ist x ?
*/

/*
   Testsequenzen z.B. 7 5 1100303
   2 5 11003059
   5 1 11003077
   3 2 1100301211
   (3,17,31) 5 1100308133
   (2,59,67, 83, 101, 113, 327) 2 11003012203
   2 6 1100301220307
   (5,7,11,19, 43,  ) 2 1100301220307
   53 2 110030122030729
*/

import java.math.BigInteger;
import java.util.BitSet;

public final class DiskreterLog_Faktor extends Thread {
  private BigInteger n;  // Primzahl
  private BigInteger phi;  // Phi(n), also n-1
  private BigInteger basis;  // Basiszahl
  private BigInteger wert;
  private boolean fertig;
  private BigInteger exponenten [];  // Potenz der Basiszahl
  private long rel [][];  // Exponenten, leider nicht binär
  private int anzahlRel;
  private BigInteger primBasis []; // hier alle verwertbaren Primzahlen 
  private long phi_long;
  long schritteWert = 0; // Quadrierungen des Zielwertes
  long belWert [] = null; // Relation des ggf. transformierten Zielwertes

  private static final BigInteger ZWEI = BigInteger.valueOf(2);
  private static final BigInteger DREI = BigInteger.valueOf(3);
  private boolean DEBUG = false;  // zeigt die Matrix und Verrechnungen an
  private boolean DEBUG2 = false;  // zeigt nicht geglückte Zerlegungen an
  private boolean DEBUG3 = false;  // zeigt geglückte Zerlegungen an !
  private boolean DEBUG4 = false;  // zeigt Matrix zwischendrin an.

  // Konstruktor. Übernimmt Eingabewerte
  public DiskreterLog_Faktor (String args[])
  {
    if (args.length < 3) { 
      System.out.println ("Syntax: java DiskreterLog_Faktor 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); }

    BigInteger teil = phi.shiftRight(1);
    if (!teil.testBit(0) || !teil.isProbablePrime(2))
    {
      // Hinweis ausgeben, dass es dafür auf jeden Fall bessere Verfahren gibt
      System.out.println ("\n\tHINWEIS:\n\tPhi der Primzahl ist nicht das Doppelte einer " +
        "anderen Primzahl.\n\tFuer diesen Fall gibt es bessere Verfahren.\n");
    }

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


  // Hauptmethode für Thread. Werte den Status aus 0=Vorarbeit, 1=Zahl faktorisieren,
  // 2=Basis füllen, 3 = Restarbeiten zum Ermitteln des Exponenten
  public final void run ()
  {
    int status = 0;
    boolean erfolg;
    do
    {
      if (status == 0)
      {
        // Noch Basis von Primfaktoren gebraucht
        erfolg = primBasisBerechnen ();
        if (! erfolg) { fertig = true; System.exit(0); }

        // Nochmals einfacher Test, ob überhaupt Primitivwurzel
        zusatztest_primitivwurzel();
        status = 1;
      }
      else if (status == 1)
      {
        // Zahl 'wert' hoffentlich faktorisierbar, auch mit Transformationen wie Quadrat
        erfolg = run_zahlFaktorisieren();
        if (erfolg) status = 2; else break;
      }
      else if (status == 2)
      {
        // Basis hoffentlich füllbar
        erfolg = run_basisFuellen();
        if (erfolg) status = 3; else break;
      }
      else if (status == 3)
      {
        // Matrix in Diagonalmatrix umwandeln, dann Exponent zu evt. transformierten 'wert'
        // berechnen, dann ggf. rückrechnen/-transformieren.
        erfolg = restarbeiten(); break;
      }
    } 
    while (true);
    // System.out.println ("Ende mit Status " + status); 
  }


  /*
    Hier die Methoden für einzelne Teilschritte
  */

  // PrimBasis noch berechnen (Wenn Zahl zu gross, dann evt. zu viele Primzahlen). Wenn einer
  // der gespeicherten Primfaktoren die Zahl bereits teilt, dann hier gar nicht mehr weitermachen
  // Es werden alle Primzahlen aufgenommen. Grund: Quadrate sind nicht einfach zu beseitigen, wegen
  // der festen Basis.
  private final boolean primBasisBerechnen()
  {
    final int obergrenze = 100000;
    if (n.getLowestSetBit() > 0)
    {
      System.out.println ("Zahl ist gerade, hat also den Faktor 2"); return false;
    }
    int anzahl = PrimGenerator.holeAnzahlPrimzahlen(n.bitLength());  // wie viele Primzahlen
    System.out.println ("Anzahl Primzahlen: " + anzahl);
    if (anzahl < 2) anzahl = 2;  // mindestanzahl
    /* Mit dem Jacobi-Symbol kann hier im Gegensatz zur Faktorisierung keine Optimierung gemacht werden
       Grund ist, dass ja eine Gruppe hier zyklisch erzeugt wird, also alles zwischen [1; p-1]
       Zwar werden im folgenden Potenzen von g^2 untersucht, also Jac(g^2,p) = +1, scheinbar können
       aber doch alle Primfaktoren vorkommen, nicht nur Jac(p % g^2, p) == +1, heuristisch ermittelt
       z.B. in <7, 41> und <5, 43>.
       Es sind also (leider) alle Primzahlen unter einer Schranke zu nehmen
    */

    // wegen der festen Basis nicht rausgeschmissen werden kann.
    if (anzahl > obergrenze)
    {
      System.out.println ("Reduziere Primzahlbasis auf " + obergrenze);
      anzahl = obergrenze;
    }

    anzahl ++;  // wegen -1
    primBasis = new BigInteger [anzahl]; primBasis[0] = BigInteger.valueOf(-1); primBasis[1] = ZWEI;
    int zahl = 3;
    for (int i = 2; i < anzahl; zahl += 2 /* i s.u. */)  // Start bei i=2, weil -1 und 2 drin
    {
      // Jac(n%p, p) == 1 und Primzahl, dann aufnehmen !
      BigInteger aktuell = BigInteger.valueOf(zahl);
      // Hier zur Sicherheit Teilertest - weil sonst evt. die Primzahl nicht invertiert werden kann !
      if (n.remainder(aktuell).equals(BigInteger.ZERO))
      {
         System.out.println ("ggT ist " + aktuell); return false;
      }
      if (!prim_mr(aktuell)) continue;  // wenn zusammengesetzt, dann sicher nicht verwenden
      // Achtung: Manchmal ist die 9 durchgeschlpft, es können also zus. Zahlen in der Basis stehen !!!
      // Achtung: Anders als beim Quadratischen Sieb kann hier nicht mit Jacobi-Symbol die Primbasis
      // verringert werden, dies zeigt sich schon bei 41 und 42.

      // Diese Primzahl für Basis eintragen
      primBasis[i] = aktuell; ++i; // System.out.println ("Primzahl ist " + zahl);
    }

    // Faktor-Basis
    exponenten = new BigInteger[anzahl];
    rel = new long [anzahl][anzahl];
    // die Bitsets in Tab. rel werden erst bei Bedarf angelegt, also wenn neue Relation gefunden wurde !
    System.out.println ("Größte Primzahl ist " + primBasis[primBasis.length-1]);
    anzahlRel = 0;

    return true;  // Primfaktor-Basis wurde aufgebaut
  }


  // versucht, die Basis zu füllen. Gibt zurück, ob fertig
  private final boolean run_basisFuellen ()
  {
    int k;

    System.out.println ("Suchen ...");
    BigInteger teil = basis.modPow(ZWEI, n);
    BigInteger lauf = ZWEI;
    BigInteger aktuell = teil;
    boolean erfolg = false;
    long zeitVorher = System.currentTimeMillis();
    while (lauf.compareTo(phi) < 0 )  // Exponent phi entspricht 0
    {
      // -1 kann erzeugt werden, +1 aber nicht. Also hier Test, sonst sind FolgeFehler möglich
      if (aktuell.equals(BigInteger.ONE)) { 
        System.out.println ("Exponent " + lauf + ": Basis ist keine Primitivwurzel"); System.exit(1); }
      boolean geklappt = zerlegenEintragen(aktuell, lauf);
      // schon alles zerlegt ?
      if (geklappt)
      {
        if (anzahlRel == primBasis.length-1)
        {
          // sind die höchsten Exponenten alle 1, oder bei faktorisiertem Wert <= alle dort
          // vorkommenden Exponenten ?
          for (k = 0; k < primBasis.length-1; ++k) if (rel[k][k+1] != 1) break;
          if (k == primBasis.length-1) { System.out.println("... Faktorisierung fertig"); erfolg=true; break; }
        }
      }
      aktuell = aktuell.multiply(teil).remainder(n);
      lauf = lauf.add(ZWEI);
    }
    System.out.println ("Basis erstellen in Sekunden: " + (System.currentTimeMillis() - zeitVorher) / 1000.0 + "\n");

    return erfolg;
  }


  // versucht, die Zahl zu faktorisieren. Gibt zurück, ob fertig
  // es wird ausgehend von wert eine Zahl gesucht, welche über der Primbasis komplett zerlegt werden kann.
  // Dann wird deren Exponent ermittelt, um dann zurück zu rechnen.
  private final boolean run_zahlFaktorisieren ()
  {
    BigInteger einf = null;  // einfache Fälle hier abtesten ...
    if (wert.equals(BigInteger.ONE)) einf = BigInteger.ZERO;
    else if (wert.equals(basis)) einf = BigInteger.ONE;
    else if (wert.equals(n.subtract(basis))) einf = phi.shiftRight(1).add(BigInteger.ONE);
    if (einf != null) { System.out.println("Log = " + einf); System.exit(0); }

    // irgendwie eine ableitbare Zahl (Quadrat, inverses) finden, das faktorisiert werden kann.
    BigInteger start = wert; boolean gefunden = false;
    int schritteMax = n.bitLength();  // so oft dürfte maximal quadriert werden, wichtig, falls n nicht prim !
    long zeitVorher = System.currentTimeMillis();
    for (int i = 0; i < schritteMax; ++i, ++schritteWert)
    {
      belWert = new long[primBasis.length];
      gefunden = zerlege_quadrat(start, belWert, true);
      if (gefunden) break;
      start = start.modPow(ZWEI, n);  // ggf. wird -1|+1 erzeugt, was einfach zu zerlegen ist
    }
    if (!gefunden) { System.out.println ("Zu " + wert + " kann kein Logarithmus hiermit berechnet werden"); System.exit(1); }
    if (schritteWert > 0) System.out.println ("Achtung: Es sind noch " + schritteWert + " Quadrierungen notwendig"); 
    System.out.println ("Wert zerlegen in Sekunden: " + (System.currentTimeMillis() - zeitVorher) / 1000.0 + "\n");

    return true;
  }

  // da die Verwendung von Nicht-Primitivwurzeln zu Problemen bei Berechnung führen kann
  // bzl. erst spät erkannt wird, hier nochmals Test (einfach). basis^(Ordnung/2) ergibt -1, ungerade
  // Teiler davon dürfen dann nicht auch die -1 erzeugen.
  private void zusatztest_primitivwurzel ()
  {
    // boolean raus = true; if (raus) return;

    // basis ^ Phi/2 == -1 ist hier Voraussetzung. Suche nach ungeraden Teilern in der Primbasis,
    // welche eventuell die -1 auch schon erzeugen
    BigInteger vergl = basis.modPow(phi.shiftRight(1), n);
    BigInteger minusEins = n.subtract(BigInteger.ONE);
    if (!vergl.equals(minusEins)) { System.out.println ("Jacobi passt nicht"); System.exit(1); }

    BigInteger phiHalbe = phi.shiftRight(1);

    for (int i = 2; i < primBasis.length; ++i)
    {
      if (phiHalbe.remainder(primBasis[i]).equals(BigInteger.ZERO))
      {
        vergl = basis.modPow(primBasis[i], n);
        if (vergl.equals(minusEins)) { System.out.println (basis + " ist keine Primitivwurzel"); System.exit(1); }
        vergl = basis.modPow(phiHalbe.divide(primBasis[i]), n);
        if (vergl.equals(minusEins)) { System.out.println (basis + " ist keine Primitivwurzel"); System.exit(1); }
      }
    }
    System.out.println ("Primitivwurzel-Test zu Ende");
  }


  // zerlegt das Quadrat und arbeitet es in Relationenmatrix ein
  private boolean zerlegenEintragen (BigInteger quadrat, BigInteger hoch)
  {
    int i, j, k, l;
    long bel [] = new long[primBasis.length];
    boolean gefunden = zerlege_quadrat(quadrat, bel, true);

    if (! gefunden)
    {
      // konnte nicht komplett zerlegt werden, aber vielleicht ist Quadrat jetzt kleiner
      if (DEBUG2) System.out.println ("\t" + quadrat.toString() + " konnte nicht komplett zerlegt werden");
      return false;
    }
    else  // fertig, komplett zerlegt
    {
      // Sonderfall: wenn alle Exponenten gerade (hoch ist es ohnehin), dann hier Kürzung möglich
      if (DEBUG4) 
      {
        System.out.println ("Beginn aktuelle Relationenmatrix");
        for (k = 0; k < anzahlRel; ++k) relation_ausgeben(exponenten[k], rel[k]);
        System.out.println ("Ende aktuelle Relationenmatrix");
      }
      int letzteVariable1 = letzteVariable(bel);  // Grad der akt. Relation von oben
      hoch = normieren_kuerzen(hoch, bel, letzteVariable1);

      if (DEBUG3) 
      {
        relation_ausgeben (hoch, bel);
      }

      // Eine Belegung gilt als erfolgsversprechend, wenn a) wenige Varablen belegt sind
      // oder b) die größte Variable mglichst klein ist.
      // d.h. es wird in den Belegungen von den "seltenen" zu den "häufigen" Primfaktoren gesucht.

      // Schaue in die Matrix der Belegungen von den schlechten zu den guten rein
      // Lediglich aus historischen Gründen stehen die guten mit niedrigem Grad bei niedrigen Indizes,
      // die schlechten mit hohem Grad bei hohen Indizes drin. Es ist also bildlich eine
      // Linksdreiecksmatrix, die in folgender for-Schleife von unten nach oben durchlaufen wird,
      // in der Hoffnung, eine Relation einzufügen (zum Grad noch nichts da), eine
      // vorhandene zu verbessern (Exponenten teilen nicht ganzzahlig) oder die Relation
      // um einen Grad zu verringern (Exponenten teilen ganzzahlig).
      for (j = anzahlRel - 1; j >= 0; --j)
      {
        // welches ist die letzte benutzte Primzahl dieser Relation ?
        int letzteVariable2 = letzteVariable(rel[j]);   // Grad der Relation aus Tabelle
        //System.out.println ("Aktuelle Zeile " + j + " mit " + letzteVariable2 + 
        //   ", laufend " + letzteVariable1);
        // wenn besser als aktuelle TabellenZeile, dann einfügen
        if (letzteVariable2 > letzteVariable1) 
        {
          // erst später entscheiden ...
        }
        // verrechnen. Wenn Grad des einen modulo grad des anderen 0 ist, kann eine Variable reduz. werden
        else if (letzteVariable2 == letzteVariable1)
        {
          // System.out.println ("Gleicher Grad " + letzteVariable2);
          // relation_ausgeben(exponenten[j], rel[j]);  // gibt Zeile der Matrix aus

          long belNeu [] = new long[primBasis.length];
          BigInteger hochNeu = null;
          long grad2 = rel[j][letzteVariable2];
          long grad1 = bel[letzteVariable1];  // Annahme (s. normiert_gekuerzt(): Grad1/2 nicht < 0)
          if (grad1 <= 0 || grad2 <= 0)
          {
            System.out.println ("Grade nicht positiv: " + grad1 + " und " + grad2); System.exit(1);
          }
          long mal = (grad2 >= grad1) ? (grad2 / grad1) : (grad1 / grad2);
          if (mal == 1)
          {
            // Sonderfall. einfach abziehen, 2 - 1
            hochNeu = exponenten[j].subtract(hoch);
            if (hochNeu.signum() < 0) hochNeu = hochNeu.add(phi);
            for (k = 0; k < primBasis.length; ++k)
            {
              belNeu[k] = rel[j][k] - bel[k];
              if (phi_long >= 0) belNeu[k] = (belNeu[k] % phi_long);  // Überlauf behandeln
            }
          }
          else
          {
            BigInteger malB = BigInteger.valueOf(mal);
            hochNeu = ( (grad2 > grad1) ?
                         exponenten[j].subtract(hoch.multiply(malB)) :
                         hoch.subtract(exponenten[j].multiply(malB)) ).remainder(phi);
            if (hochNeu.signum() < 0) hochNeu = hochNeu.add(phi);
            for (k = 0; k < primBasis.length; ++k)
            {
              belNeu[k] = (grad2 > grad1) ?
                          (rel[j][k] - modmul(mal, bel[k], phi_long)) :
                          (bel[k] - modmul(mal, rel[j][k], phi_long));  // Überlauf behandeln
              //belNeu[k] = (grad2 > grad1) ?
              //            (rel[j][k] - mal * bel[k]) :
              //            (bel[k] - mal * rel[j][k]);
              if (phi_long >= 0) belNeu[k] = (belNeu[k] % phi_long);  // Überlauf behandeln
            }
          }
          belNeu[0] = ((belNeu[0] & 1) == 1) ? 1 : 0;  // Sonderfall Vorzeichen
          letzteVariable1 = letzteVariable(belNeu);
          //for (k = 0; k < belNeu.length; ++k) System.out.print(belNeu[k]+",");
          //System.out.println ("Vor Normieren: " + hochNeu);
          hochNeu = normieren_kuerzen (hochNeu, belNeu, letzteVariable1);
          //System.out.println("Ergebnis: ");
          // relation_ausgeben (hochNeu, belNeu);

          // Was ist mit der neuen Relation zu machen ?
          if (((grad1 % grad2) == 0) || ((grad2 % grad1) == 0))
          {
            // es kommt eine neue Relation ohne diesen Teiler, die Tabelle wird ggf. verbessert,
            // die neue Relation wandert weiter
            if (grad1 < grad2)
            {
              exponenten[j] = hoch; if (DEBUG) System.out.println("Bessere übernehmen"); // in Matrix, da besser als dortige
              // if (j == 0 && DEBUG3) System.out.println ("--> Neuer Spitzenreiter: "); // zur Info
              for (k = 0; k < primBasis.length; ++k)
              {
                rel[j][k] = bel[k];
              }
            }
            if (letzteVariable1 <= 0)
            {
              return true;  // hat keine Information mehr
            } 
            hoch = hochNeu;
            for (k = 0; k < primBasis.length; ++k)
            {
              bel[k] = belNeu[k];
            }
            if (letzteVariable1 == letzteVariable2)
            {
              System.out.println ("Fehler beim Verrechnen: Grad nicht verringert");  // kann bei negat. Exponenten auftreten
              relation_ausgeben(hoch, bel); System.exit(1);
            }
            if (DEBUG) System.out.println(" -Mit besserer Relation weiter");
          }
          else
          {
            if (letzteVariable1 <= 0)
            {
              return true;  // hat keine Information mehr
            } 
            // es bleibt dieser Teiler mit einem niedrigeren Exponenten übrig. Dieser ersetzt die
            // Relation in der Tabelle, dann stop.
            exponenten[j] = hochNeu;
            for (k = 0; k < primBasis.length; ++k)
            {
              rel[j][k] = belNeu[k];
            }       
            // if (j == 0 && DEBUG3) System.out.println ("--> Neuer Spitzenreiter: "); // zur Info
            if (DEBUG) System.out.println(" -Relation ersetzt bisherige. Ende");
            return true;
          }
        }
        // einfügen
        else if (letzteVariable2 < letzteVariable1)
        {
          if (DEBUG) System.out.println(" -einfügen");
          ++j;break;  // schlechter als das aktuelle: danach einfügen
        }
      } // alle Relationen von unten nach oben durchlaufen


      // Was ergibt die Belegung ? Gibt es keine Quadratfaktoren mehr, dann auswerten,
      // wobei die Triviale Belegung ja auch vorkommen kann (b^0 == (leer), b^((p-1)/2) == -1 )
      if (anzahlRel == primBasis.length-1)  // wegen -1 eine Relation weniger als Primteiler
      {
        // hier muss ein Fehler vorliegen:
        System.out.println("Fehler in der Verrechnung: Gleichung nicht mehr speicherbar");
        System.exit(1); 
      }
      else
      {
        // neue Relation bel in Zeile j abspeichern
        if (j < 0) j = 0;  // es gab bisher nichts
        // if (j == 0 && DEBUG3) System.out.println ("--> Neuer Spitzenreiter: "); // zur Info

        if (j < anzahlRel)
        {
          if (DEBUG) System.out.println("Einfügen an " + j);
          for (k = anzahlRel; k > j; --k)
          {
            rel[k] = rel[k-1];
            exponenten[k] = exponenten[k-1]; 
          }
        }
        if (DEBUG)
        {
          System.out.println ("Belegen bei #" + j + " und Hoch " + hoch.toString());
          for (k = 0; k < primBasis.length; ++k) 
          {
            System.out.print(bel[k]);
          }
        }
        rel[j] = (long []) bel.clone();  // ist clone() hier notwendig ?
        if (DEBUG) System.out.println("");
        exponenten[j] = hoch;
        ++anzahlRel;
      }

      return true;
    } // komplett zerlegt

  }


  // Gibt Exponent und nicht-triviale Faktoren in einer Zeile aus. Bei Schalter noch Testrechnung
  // Im Debug-Modus wird noch muehsam nachgerechnet
  public void relation_ausgeben(BigInteger hoch, long [] bel)
  {
    if (hoch.signum() < 0 || hoch.compareTo(phi) >= 0) {
       System.out.println (hoch + " entspricht nicht der Norm"); System.exit(1); }
    System.out.print ("\t" + hoch + " = ");
    BigInteger tats = DEBUG ? basis.modPow(hoch,n) : null;
    BigInteger vergl = DEBUG ? BigInteger.ONE : null;
    if (DEBUG) System.out.print (tats + " = ");
    for (int k = 0; k < primBasis.length; ++k)
    {
      if (bel[k] != 0) {
        System.out.print("," + primBasis[k]); if (bel[k] != 1) System.out.print ("^" + bel[k]);
      }
      if (DEBUG && bel[k] != 0) { // nachrechnen
        BigInteger faktor = primBasis[k].modPow(BigInteger.valueOf(Math.abs(bel[k])), n);
        if (bel[k] > 0) vergl = vergl.multiply(faktor).remainder(n); 
        else tats = tats.multiply(faktor).remainder(n);
        // System.out.println ("Zusatz " + faktor);
      }
    }
    System.out.println("");
    if (DEBUG && !vergl.equals(tats)) { 
      System.out.println ("Fehler Belegung: Aus pos. Koeff " + vergl + ", aus Exp. und neg. Koeff " + tats); System.exit(1); }
  }


  // versucht, 1) den höchsten Koeffizienten auf positiv zu normen
  // 2) alle Exponenten um 2 zu kürzen. Sinnvoll, damit Vielfaches von Einträgen nicht gesp. werden
  // Der Wert des höchsten Grades hoechsteV ist wichtig für Vorzeichenbestimmung
  private BigInteger normieren_kuerzen(BigInteger hoch, long [] bel, int hoechsteV)
  {
    int i;
    if (hoch.signum() < 0 || hoch.compareTo(phi) >= 0) { System.out.println ("Hoch " + hoch + " ist nicht normiert"); System.exit(1); }
    if (hoechsteV < 0) 
    {
      if (hoch.equals(BigInteger.ZERO)) return hoch;  // trivial
      // stimmt das überhaupt ?
      BigInteger tmp = basis.modPow(hoch, n);
      if (! tmp.equals(BigInteger.ONE)) { System.out.println ("Rechenfehler für Exp. " + hoch); System.exit(1); }
      // dann ist Basis aber keine Primitivwurzel !
      System.out.println (basis + " ist keine Primitivwurzel, Ordnung " + hoch); System.exit(1);
    }
    if (bel[hoechsteV] == 0) { System.out.println (hoechsteV + " stimmt nicht"); System.exit(1); }

    // 1) Norm: höchste Koeffizienten positiv machen, dann auch Hoch-Faktor
    if (bel[hoechsteV] < 0)
    {
      for (int k = 1; k < primBasis.length; ++k) bel[k] = - bel[k];  // -1 muss nicht
      hoch = phi.subtract(hoch); // System.out.println ("Hoch genormt");
    }

    // 2) Gerade Koeff. kürzen
    do  // Vorsicht: Endlosschleife, wenn alle Exponenten 0
    {
      // das vorzeichen an [0] behandeln
      if (bel[0] != 0)
      {
        hoch = hoch.add(phi.shiftRight(1));  // phi muss gerade sein, Potenz der Hälfte gibt -1
        bel[0] = 0; if (hoch.compareTo(phi) >= 0) hoch = hoch.subtract(phi);
      }
      for (i = 1; i < bel.length; ++i) if ((bel[i] & 1) == 1) break;
      if ((i < bel.length) || hoch.testBit(0)) 
      {
        // zweite Chance alle durch 3 teilbar testen; andere wie z.B. 5 wird zu unwahrsch. sein.
        for (i = 1; i < bel.length; ++i) if ((bel[i] % 3) != 0) break;
        if (i < bel.length) break;  // auch das nicht
        // jetzt muesste noch hoch irgendwie durch 3 teilbar sein.
        if ((hoch.remainder(DREI).equals(BigInteger.ZERO)) ||
            ((phi.subtract(hoch)).remainder(DREI).equals(BigInteger.ZERO)) )
        {
          boolean DEBUGSOND3 = false;  // hier lokal Ausgabe anstellen; Sachverhalt unten nicht ganz klar.
          if (DEBUGSOND3) System.out.println ("Sonderfall der 3"); 
          // Hier gibt es Chance, dass durch 3 geteilt werden kann. Test mit 2 5 11003059
          // rückgängig gemacht werden muss das ganze im Fall x^3 == y^3, aber x != +- y,
          // hier weiss ich keine Abhilfe.
          
          // Exponenten ganzzahlig teilen
          if (DEBUGSOND3)
          {
            System.out.print (hoch.toString() + '|');
            for (i = 0; i < bel.length; ++i) System.out.print(bel[i] + " ");
            System.out.println (""); 
          }
          BigInteger hochSicher = hoch;
          for (i = 1; i < bel.length; ++i) if (bel[i] != 0) bel[i] /= 3;
          if (hoch.remainder(DREI).equals(BigInteger.ZERO))
          {
            hoch = hoch.divide(DREI);
          }
          else { 
             hoch = phi.subtract( (phi.subtract(hoch)).divide(DREI) ); // phi=11002, h=7561, h' = 9855
          }
          if (DEBUGSOND3)
          {
            System.out.print (hoch.toString() + '|');
            for (i = 0; i < bel.length; ++i) System.out.print(bel[i] + " ");
            System.out.println ("");
          }

          boolean test_ok = true;
          { // Testsequenz
            BigInteger tmp = basis.modPow(hoch, n);
            BigInteger pos = BigInteger.ONE;            BigInteger neg = BigInteger.ONE;
            for (i = 1; i < primBasis.length; ++i)
            {
              BigInteger zw = primBasis[i].modPow(BigInteger.valueOf(Math.abs(bel[i])), n); 
              if (bel[i] < 0)
                neg = neg.multiply(zw).remainder(n);
              else if (bel[i] > 0)
                pos = pos.multiply(zw).remainder(n);
            }
            if (DEBUGSOND3) System.out.println ("Nachher: " + tmp + ", positiv " + pos + ", negativ " + neg);
            if (! tmp.multiply(neg).remainder(n).equals(pos) )
            {
              // x^3 == y^3, aber x != y; weiss hier nicht, wie weiter, also Rolle Rückwärts
              if (DEBUGSOND3) System.out.println ("Teilung durch 3 war nicht erfolgreich");
              hoch = hochSicher; test_ok = false;
              for (i = 1; i < bel.length; ++i) if (bel[i] != 0) bel[i] = bel[i] * 3;
            }
          }
          if (DEBUGSOND3) warten();
          if (test_ok) continue; else break;  // abhängig, ob Vergleich geklappt hat
        }
        else break;  // auch keine Teilung durch 3 möglich
      }
      // System.out.println ("Kürzung möglich !");
      hoch = hoch.shiftRight(1);
      for (i = 1; i < bel.length; ++i) bel[i] >>= 1;
      // Das Vorzeichen kann nun anderes sein, da es faktisch einer Wurzel-Ziehung entspricht
      // z.B. mod 179: 11^12 == 4 mod 179, 11^6 == -2 mod 179, +2 wäre hier falsch !
      // folgende Lösung ist aufwendig - gibt es nicht besseres; Test auf höchsten Primteiler Rest 0 nicht zuverl. !
      BigInteger tats = basis.modPow(hoch, n);  // leider TEUER !!!
      BigInteger vergl = BigInteger.ONE;
      for (i = 1; i < primBasis.length; ++i)
      {
        if (bel[i] != 0) { // nachrechnen
          BigInteger faktor = primBasis[i].modPow(BigInteger.valueOf(Math.abs(bel[i])), n); 
          if (bel[i] > 0) vergl = vergl.multiply(faktor).remainder(n); 
          else tats = tats.multiply(faktor).remainder(n);
        }
      }
      if (vergl.equals(tats)) {} // OK
      else if ((n.subtract(vergl)).equals(tats))
      {
        bel[0] = 1; // System.out.println ("Vorzeichen minus"); 
      }
      else { System.out.println ("Rechenfehler bei Vorzeichenwechsel erkannt"); 
         relation_ausgeben (hoch, bel); System.exit(1); }
    }
    while (true);

    return hoch;  // zurückgeben, da ggf. verändert
  }


  // Zerlegt quadrat über der Faktorbasis
  // einmalig mit Flagge anfang wird ein nicht fakt. Rest Invertiert und versucht zu zerlegen
  // arbeitet modifiziert, erst alle 3er suchen
  // Anm.: Die Erweiterung des herkömmlichen Siebes besteht aus drei Teilen (Dekl., Teiler, weiter)
  private boolean zerlege_quadrat (BigInteger q, long [] bel, boolean anfang)
  {
    int i, j;
    if (q.equals(BigInteger.ZERO)) { System.out.println ("Null hier nicht erlaubt"); System.exit(1); }

    BigInteger dAR [];

    // da die 2er eine Sonderstellung einnehmen, aus Performanzgrnden hier rausziehen
    // Anm: Bei MBH kommt scheinbar die 2 nur in gerade Potenz vor, wenn Jac(q, n) = -1
    int weg = q.getLowestSetBit();
    if (weg > 0)
    {
      q = q.shiftRight(weg);  // ist nun ungerade
      bel[1] += weg;
    }

    BigInteger qsicher = q;  // nur zum Test
    boolean pruefen = true;
    // Es es eine einzige Primzahl ist, die grer als das max. Element der Basis ist, dann nichts tun
    if (q.compareTo(primBasis[primBasis.length-1]) > 0 && q.isProbablePrime(1) ) 
    {
      // Achtung: Hier werden auch zusammengesetzte Zahlen ggf. als Prim erkannt, z.B. auch 289
      if (DEBUG2) System.out.println ("Vorzeitiger Abbruch wegen Primzahl " + q);
      pruefen = false;
    }

    // alle Potenzen von Primfaktoren raushauen ...
    boolean m3mod4 = q.testBit(1); int start3 = 2; int start1 = 2;
    for (j = 2; j < primBasis.length && pruefen; ++j)   // -1, 2 nicht beachten
    {
      BigInteger f = primBasis[j];
      if (m3mod4 && !f.testBit(1)) continue;  // 1 mod 4 überspringen
      if (!m3mod4 && j < start3 && f.testBit(1)) continue;  // schon geteilte 3 überspringen
      int wieoft = 0; boolean geteilt = false;
      while (true)
      {
        // nur diejenigen, die auch mit einfacher Potenz vorkommen. Rest unten mit Quadrattest ermitteln
        dAR = q.divideAndRemainder(f);  // diese Stelle ist laut hprof besonders teuer
        if (! dAR[1].equals(BigInteger.ZERO))
        {
          break;  // Rest nicht 0
        }
        dAR[1] = null;
        q = dAR[0] /* Divisionsergebnis  */  ; ++wieoft; geteilt = true; // ganzzahlig geteilt
      }
      bel[j] += wieoft;
      if (geteilt && q.equals(BigInteger.ONE)) break; // schn, konnte komplett zerlegt werden

      // Da Primzahltest teuer ist, wird er hier erst gemacht, wenn ein gewisser Teil der 
      // Basis durch ist (hier einfach mal 40 ) !
      if ((j == 40) || (j > 40 && geteilt))
      {
        if (q.compareTo(primBasis[primBasis.length-1]) > 0 && q.isProbablePrime(1) ) 
        { 
          // Achtung: Hier kann auch bei zusammenges. Zahlen abgebrochen werden
          if (DEBUG2) System.out.println ("Abbruch " + j + " von " + primBasis.length); 
          pruefen = false; // die verbleibende Zahl ist prim und kann nicht weiter zerlegt werden
        }
      }

      if (wieoft > 0)
      {
        // wie geht es weiter ?
        boolean n3 = q.testBit(1);
        if (m3mod4 && !n3) 
        { 
          start3 = j+1; j = start1-1; // Wechsel 3 --> 1
        }
        if (!m3mod4 && n3)   // Wechsel 1 --> 3
        {
          if (j >= start3) start3 = j+1; // 3 er können weiter liegen, dann nicht überschr.
          start1 = j+1;
          j = start3-1;
        }
        m3mod4 = n3; 
      }

    } // Ende alle Primfaktoren durchmachen

    if (anfang && !q.equals(BigInteger.ONE))
    {
      // Abbruch obiger Teilertests: einmalig mit dem inversen
      // einmalig, um Rekursion zu Vermeiden. Einfachstes Beispiel: q ist prim, das inverse auch,
      // dann würde es dazwischen hin- und herspringen
      BigInteger invers = q.modInverse(n);
      if (DEBUG2) System.out.println ("Weiter mit dem inversen " + invers);
      long belNeu [] = new long[primBasis.length];
      boolean jetzt = zerlege_quadrat(invers, belNeu, false /* zweiter Versuch */);
      if (jetzt)
      {
        // System.out.println ("... Inverses bringt Erfolg");
        for (j = 0; j < primBasis.length; ++j) bel[j] -= belNeu[j];
        q = BigInteger.ONE;  // anzeigen, dass doch komplett zerlegt; höchster Koeff ggf. < 0 !
      }
    }
  
    return q.equals(BigInteger.ONE);
  }


  // Matrix in Dialogonalmatrix umwandeln, Exponent zu tranf. 'wert' berechnen, dann
  // rückrechnen.
  private boolean restarbeiten ()
  {
    int j, k;

    if (DEBUG4)
    {
      System.out.println ("Aktuelle Linksmatrix der Relationen");
      for (k = 0; k < anzahlRel; ++k) relation_ausgeben(exponenten[k], rel[k]);
      System.out.println ("Ende Linksmatrix der Relationen");
    }

    // Es liegt eine Linksdreiecksmatrix vor, d.h. die guten Relationen mit wenig
    // Primfaktoren stehen oben. 0-te Spalte ist Vorzeichen, sollte immer + sein.
    //Jetzt in Diagonalmatrix umwandeln
    for (j = 0; j < anzahlRel; j++)
    {
      // welches ist die letzte benutzte Primzahl dieser Relation ?
      int letzteVariable1 = letzteVariable(rel[j]);   // Grad der Relation aus Tabelle
      if (rel[j][0] != 0) { System.out.println ("Vorzeichen noch übrig"); System.exit(1); }
      if (letzteVariable1 != (j+1)) { System.out.println ("Fehler: keine Linksdreiecksmatrix"); System.exit(1); }
      if (rel[j][letzteVariable1] != 1) { System.out.println ("Mittelelement nicht 1"); System.exit(1); }
      // alle Spalten in nachfolgenden Zeilen auf 0 setzen, da glatt teilbar.
      for (k = j+1; k < anzahlRel; ++k)
      { 
        long faktor = rel[k][letzteVariable1]; if (faktor == 0) continue;
        BigInteger zusatz = exponenten[j].multiply(BigInteger.valueOf(Math.abs(faktor))).remainder(phi);
        if (faktor > 0)
        {
          exponenten[k] = exponenten[k].subtract(zusatz);
          if (exponenten[k].signum() < 0) exponenten[k] = exponenten[k].add(phi);
          rel[k][letzteVariable1] -= modmul(rel[j][letzteVariable1], faktor, phi_long);
        }
        else
        {
          exponenten[k] = exponenten[k].add(zusatz);
          if (exponenten[k].compareTo(phi) >= 0) exponenten[k] = exponenten[k].subtract(phi);
          rel[k][letzteVariable1] += modmul(rel[j][letzteVariable1], -faktor, phi_long);
        }
        if (rel[k][letzteVariable1] != 0) { System.out.println ("Verrechnungsfehler, nicht 0"); System.exit(1); }
      }
    }

    {
      System.out.println ("Exponenten zur Basis (links) fuer Primfaktor (rechts)");
      for (k = 0; k < anzahlRel; ++k) relation_ausgeben(exponenten[k], rel[k]);
      System.out.println ("Ende Exponenten fuer Primfaktor");
    }

    // Test, ob die einzelnen Koeffizienten passen
    for (k = 0; k < anzahlRel; ++k)
    {
      BigInteger w1 = primBasis[k+1].modPow(BigInteger.valueOf(rel[k][k+1]), n);
      BigInteger w2 = basis.modPow(exponenten[k], n);
      if (! w1.equals(w2)) { System.out.println (k + ": Differenz " + w1 + " und " + w2); System.exit(1); } 
    }

    // Exponent von dem ggf. transformierten Wert bestimmen
    BigInteger erg = BigInteger.ZERO;
    for (k = 0; k < primBasis.length; ++k)
    {
      if (belWert[k] != 0)
      {
        BigInteger mal = BigInteger.valueOf(Math.abs(belWert[k]));
        mal = mal.multiply(exponenten[k-1]).remainder(phi);
        if (belWert[k] > 0) erg = erg.add(mal);
        else erg = erg.subtract(mal);
      }
    }
    erg = erg.remainder(phi); if (erg.signum() < 0) erg = erg.add(phi);
    System.out.println ("... Exponent " + erg);
    if (schritteWert == 0)
    {
      BigInteger vergleich = basis.modPow(erg, n);
      if (vergleich.equals(vergleich)) System.out.println ("--> Treffer");
    }
    else
    {
      System.out.println ("--> Ausgehend davon noch " + schritteWert + " Wurzeln ziehen");
    }

    System.out.println ("Restarbeiten erledigt");

    return true;
  }


  /*
    Methoden ohne viel Änderungsfrequenz
  */

  // Haupt-Methode
  public static void main (String[] args) {
    DiskreterLog_Faktor w = new DiskreterLog_Faktor (args);
    w.start ();   
  }


  // sucht nach der letzten Variable. Gibt die Position im Feld zurück oder -1
  private int letzteVariable (long [] bel)
  {
    for (int k = primBasis.length-1; k >= 1; --k)  // bis zur 2 an
    {
      if (bel[k] != 0) return k;
    }
    return -1;
  }


  // Multpliziert zwei Exponenten unter Berücksichtigung von modulo,
  // um Zahlenüberlauf zu vermeiden, aus Klasse PTyp kopiert, aber modulo muss nicht bek. sein, 
  // die Eingabewerte können auch negativ sein, Vorzeichen wird also gesichert
  private long modmul (long e1, long e2, long modulo)
  {
    boolean minus = ((e1 < 0) != (e2 < 0)); 
    if (e1 < 0) e1 = -e1; if (e2 < 0) e2 = -e2;  // danach beide positiv
    long summe = 0;  // Überlaufbehandlung
    while (true)
    {
      if (e2 < e1) { long hilfe = e2; e2 = e1; e1 = hilfe; }  // am besten das kleinere, Performanz
      if (e1 == 0) break;
      if ((e1 < Integer.MAX_VALUE) && (e2 < Integer.MAX_VALUE))
      {
        long erg = (e1 * e2); // passt jetzt in long-Wert
        summe += erg;
        if ((modulo >= 0) && (summe >= modulo)) summe = summe % modulo; 
        break;
      }
      if ((e1 & 1) == 1) { summe = (summe + e2); if ((modulo >= 0) && (summe >= modulo)) summe -= modulo; }
      e1 = e1 >> 1; e2 = e2 << 1; if ((modulo >= 0) && (e2 >= modulo)) e2 -= modulo;
    }
    return (minus ? -summe : summe);
  }


  private final static int prim100 [] = { 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
    59, 61, 67, 71, 73, 79, 83, 89, 97};

  // Eigenimplementierung, da isProbablePrime oft zusammengesetzte auch durchgehen lies
  // wird momentan aber nur bei Erzeugung der Primzahlbasis verwendet.
  private final static boolean prim_mr(BigInteger tmp)
  {
    if (!tmp.testBit(0) && !tmp.equals(ZWEI)) return false;  // gerade, aber nicht 2

    // 1) fr int schnellen Test, sofern mglich
    if (tmp.bitLength() <= 31)
    {
      int aktuell = tmp.intValue();
      for (int lauf = 0; lauf < prim100.length; ++lauf)
      {
        if (aktuell == prim100[lauf]) return true;
        if ((aktuell % prim100[lauf]) == 0) return false;
      }
    }

    // 2) Miller-Rabin mit BigInt

    // 3) normale Primtest
    return tmp.isProbablePrime(1);  // irgendwie doch viel durchgeschlpft
  }

  // Wartet auf Eingabe des Benutzers
  public static void warten()
  {
    int rueck;
    try
    {
      while ((rueck = System.in.read()) != '\n');
    }
    catch (Exception ex) {}
  }

}