// QuadratSieb20.java
// Implementierung, wie Morrsion und Brillhard.
// Beginn: August 2008
// Versucht, einen quadratischen Sieb zu implementieren. Testzahl knnte z.B. 4999574933, 56312533736935564420935284701870003 sein !
// Nimm nur die Primzahlen p in die Faktorbasis auf mit Jac(n%p, p) == +1; Zahlen -1 und -2 immer drin.

// Bei groen Zahlen ergeben sich Probleme:
// a) Die b(Gre Wurzel(n)) knnen nicht komplett zerlegt werden, oder nur mit groem Rest
// b) Die Gre der Faktorbasis steigt exponentiell an, wenn auch langsam.
//    dann mssen natrlich zwangslufig viele Relationen gesammelt werden; hier ist gute Auswahl wichtig
// c) Der Versuch, die Reste (nicht in der Faktorbasis) mit ECM zu zerlegen, brachte fast nie Erfolg,
//    nur kleine Teiler konnten gefunden werden, plus jeweils immer groem Primrest
// d) Die Darstellung in Summen vom Typ x^2 = y^2 + z^2 mod n hat nicht weitergefhrt, weil man auch hier
// Wurzel zu ziehen hat von der Grenordnung O(n)
// e) ein ggT von BNeu und BUralt lieferte nur selten einen Primfaktor > max der Primbasis. Lohnt sich nicht.
//    ebenso wie (bNeu / BUralt) = (aNeu / aAlt)^2 leider selten Erfolg brachte
// f) Versuche, kleine Quadrate zu erzeugen mit [k * n]. Diese sind aber nicht immer ber der Faktorbasis
//    zerlegbar, z.B. bei 40633 7 oder 17 ! Die Ursache, warum die Quadrate aus Kettenbruch anders sein sollen,
//    ist mir nicht klar geworden.
// 10.01.2010 wegen Version 1.1 von QuadratSieb5.java nach QuadratSieb11.java umbenannt
// g) QuadratSieb20: Asymptotisch besserern Teilertest bei Zahl 3 mod 4 eingebaut

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

public final class QuadratSieb20 extends Thread {

  private BigInteger n = null;
  private BigInteger zwein; // ganzzahliges Vielfachen
  private BigInteger viern;
  private BigInteger achtn;
  public BigInteger a,b, w1, w2;
  private BigInteger wurzel;

  private static final BigInteger SIEBEN = BigInteger.valueOf(7);
  private static final BigInteger DREI = BigInteger.valueOf(3);
  private static final BigInteger ZWEI = BigInteger.valueOf(2);
  private static final BigInteger FUENF = BigInteger.valueOf(5);
  // Rein fr Ausfiltern von Quadraten
  private static final BigInteger ELF = BigInteger.valueOf(11);
  private static final BigInteger DREIZEHN = BigInteger.valueOf(13);

  // fr Faktor-Basis
  private BigInteger wurzeln [];  // Feld mit zu zerlegenden Quadraten
  private BitSet rel [];
  private int anzahlRel;
  private boolean DEBUG = false;  // zeigt die Matrix und Verrechnungen an
  private boolean DEBUG2 = false;  // zeigt nicht geglckte Zerlegungen an
  private boolean DEBUG3 = true;  // zeigt geglckte Zerlegungen an !

  private BigInteger primBasis []; // ein paar Zahlen, um die Periode besser bestimmen zu knnen
  private BigInteger invers [] = new BigInteger[50];  // modulare Inverse der untersten (hufigsten) Primteiler
  boolean istQuadrat = false;
  public boolean fertig;
  public int anzahlGeschafft;
  public int anzahlNichtGeschafft;
  private boolean ist1Mod4;
  private BigInteger nqTestzahl = null; // erste Primzahl, die wegen Jacobi nicht in Basis ist

  // Hier noch die feste Primgrenze fr RSA-704
  // private static BigInteger primGrenze = new BigInteger("860189464792");
  // private static BigInteger primQuadrat = primGrenze.multiply(primGrenze);


  // Konstruktor
  public QuadratSieb20 (String args[]) {

    boolean fehlerEingabe = false;
    a = null; b = null; w1 = null; w2 = null;
    if (args.length >= 1) // Test: Eingabe ber Kommandozeile
    {
      try
      {
        if (args.length >= 2)
        {
          int n1 = Integer.parseInt(args[0]);
          int n2 = Integer.parseInt(args[1]);
          if (n1 < 4 || n2 < 4) fehlerEingabe = true;  // mind. diese Bitzahl pro Faktor
          else n = PrimGenerator.hole_rsa_zahl(n1, n2);  // Aus 2 Primzahlen Eingangszahl generieren
          if (n == null) { System.out.println ("Es konnte keine RSA-Zahl generiert werden"); System.exit(1); }
        }
        else
        {
          // n = BigInteger.ONE.shiftLeft(512).add(BigInteger.ONE);  // 7/9-te Fermat Zahl
          // System.out.println (n);
          n = new BigInteger(args[0]);  // 1 Eingangszahl
        }
      } 
      catch (Exception ex) { fehlerEingabe = true; }
    }
    else
    {
      System.out.println ("Aufruf: QuadratSieb20 [[Zahl] | [Bitgroesse1 Bitgroesse2]]");
      System.out.println ("Faktorisiert die Zahl, oder eine aus den Bitlaengen ermittelte schwere Testzahl "
        + "(RSA-Typ)");
      System.exit(1);
    }


    // Wertebereiche prfen
    if (fehlerEingabe || n == null || (n.compareTo(BigInteger.valueOf(15)) < 0) || 
        (a != null && a.signum() < 0) || (b != null && b.signum() < 0) || 
        (w1 != null && w1.signum() < 0) || (w2 != null && w2.signum() < 0))
    {
      fehlerEingabe = true;
    }
    if (fehlerEingabe)
    {
      System.out.println ("Aufruf: QuadratSieb20 [[Zahl] | [Bitgroesse1 Bitgroesse2]]");
      System.out.println ("Alle Eingaben sind natrliche Zahlen, Zahl >= 15, Bitgroesse >= 4");
      System.exit(1);
    }
    System.out.println ("Zu faktorisierende Zahl ist " + n.toString() );

    // Optionale Wert u.U. noch setzen !
    if (a == null) a = BigInteger.ZERO;
    if (b == null) b = BigInteger.ONE;
    if (w1 == null) w1 = BigInteger.ZERO;
    if (w2 == null) w2 = BigInteger.ONE;

    // Bei Primzahl macht das hier wenig Sinn ... Anzahl Tests steigt mit Gre
    if (n.isProbablePrime(2))
    {
      System.out.println ("Es ist wahrscheinlich eine Primzahl !"); 
      System.exit(0);
    }

    zwein = n.shiftLeft(1);
    viern = zwein.shiftLeft(1);
    achtn = viern.shiftLeft(1);
    wurzel = holeWurzel(n, false);  // unten angenäherter Wert
    ist1Mod4 = (n.and(DREI).intValue() == 1);

    anzahlGeschafft = 0; anzahlNichtGeschafft = 0;

    // folgende Anweisungen fr das Betriebssystem
    fertig = false;
    this.setPriority(Thread.MIN_PRIORITY);  // damit man nebenher noch was machen kann

    // Haken, dass beim Herunterfahren (oder CTRL + C) die aktuellen Werte nicht verloren gehen
    Runtime.getRuntime().addShutdownHook(new QuadratSieb20Haken(this));
  } // Ende Konstruktor



  // 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 nur die aufgenommen mit Jac(n%p, p) == +1. Diejenigen mit == -1 können im Quadrat
  // vorkommen, dies wird hier aber anders erkannt und beseitigt
  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;
    }
    // Beim Kettenbruchverfahren haben die zu faktorisierenden Zahlen maximal halbe n-Bitlänge
    int anzahl = PrimGenerator.holeAnzahlPrimzahlen(n.bitLength()/* >> 1*/);  // wie viele Primzahlen
    if (anzahl < 2) anzahl = 2;  // mindestanzahl
    // Optimierung mit Jacobi-Symbol
    boolean dreiMod4 = n.testBit(1) && n.testBit(0);
    if (!dreiMod4) anzahl >>= 1;  // n == 1 mod 4: die mit Jacobi -1 raus
    else anzahl = (int) Math.round(anzahl * 0.75);  // Faktoren 3 mod 4 stets, die anderen mit Jacobi -1 raus
    System.out.println ("Anzahl Primzahlen: " + anzahl);

    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); fertig = true; return false;
      }
      if (!prim_mr(aktuell)) continue;  // wenn zusammengesetzt, dann sicher nicht verwenden
      // Achtung: Manchmal ist die 9 durchgeschlpft, es knnen also zus. Zahlen in der Basis stehen !!!
      if (dreiMod4 && ((zahl & 3) == 3)) {}  // Primfaktor stets übernehmen
      else
      {
        // Jacobi prüfen: wenn -1, dann nicht im Primbasis. Achtung: Quadrat x * x kann immer vorkommen !
        int erg = PrimGenerator.jacobiSymbol (n.remainder(aktuell), aktuell);
        if (erg == 0) { System.out.println ("Teiler ist " + zahl); fertig = true; return false; }
        else if (erg < 0) { System.out.println ("Ign.: " + zahl); if (nqTestzahl == null) nqTestzahl = aktuell; continue; }
      }

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

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

    return true;  // Primfaktor-Basis wurde aufgebaut
  }


  // sucht solange, bis die 1 da ist, oder der ggT einen nicht-trivialen Wert gefunden hat
  public final void run ()
  {
    int weg, weg2, i;

    // die folgenden beiden Aufrufe knnte man auch im Konstruktor machen - aber das sie das ganze
    // beenden knnen sollen, also hier in der Methode - Reihenfolge ist wichtig.

    // Test auf Quadratwurzel hier sinnvoll zu testen, sonst evt. Probleme mit Kettenbruch
    {
      BigInteger wur = this.istWurzel(n, false /* kein Zusatztest */);
      if (wur != null) { System.out.println ("Quadratzahl zu " + wur); fertig = true; return; }
    }

    // Noch Basis von Primfaktoren gebraucht
    boolean erzeugt = primBasisBerechnen ();
    if (! erzeugt) { fertig = true; System.exit(0); }

    long zeitVorher = System.currentTimeMillis(); int restsekunden = -1;
    System.out.println ("Suchen ...");
    BigInteger dar [];   // Division und Rest in einem Schritt
    BigInteger bNeu, aNeu, wNeu, erg;
    while (true )
    {
      dar = (wurzel.add(a)).divideAndRemainder(b);
      erg = dar[0];
      aNeu = wurzel.subtract(dar[1]); dar = null;
      bNeu = (n.subtract(aNeu.multiply(aNeu))).divide(b);  // b(neu)
      if (DEBUG) System.out.println ("a=" + aNeu.toString() + ", b=" + bNeu.toString());
      
      // Wurzel noch aktualisieren
      wNeu = erg.equals(BigInteger.ONE) ? w2 : w2.multiply(erg).remainder(n);
      wNeu = wNeu.add(w1); if (wNeu.compareTo(n) >= 0) wNeu = wNeu.subtract(n);

      // Suche noch nach Periodenmitte
      if (ist1Mod4)
      {
        if (bNeu.equals(b))
        {
          System.out.println ("\n@@@\nPeriodenmitte ohne Mittelelement\n@@@\n"); 
          // Nur Wurzel von -1 findbar
          // System.out.println (w2 + " und " + wNeu);
          BigInteger inv = wNeu.modInverse(n);
          BigInteger tmp = w2.multiply(inv).remainder(n);
          testeWurzelMinusEins(tmp);
          zeige_statistik();
          fertig = true; return;
        }
      }
      // das nchste kann wahrs. bei 1&3 mod 4 auftreten
      {
        if (aNeu.equals(a))
        {
          System.out.print ("\n@@@\nPeriodenmitte mit Mittelelement " + b.toString()); 
          // hier knnte man noch Wurzel vergleichen ...
          if (w1.equals(wNeu) || w1.equals(n.subtract(wNeu)) )
          {
            System.out.println (" -->Keine weitere Info");
          }
          else
          {
            BigInteger test = (w1.subtract(wNeu)).gcd(n);
            if (! test.equals(BigInteger.ONE))
            {
              System.out.println (" -->Teiler ist " + test);
            }
          }
          System.out.println("\n@@@\n");
          zeige_statistik();
          fertig = true; return;
        }
      }


      //
      // Hier konsistent zuweisen, damit es keine Inkonsistenzen gibt !
      //
      if (! fertig)
      {
        w1 = w2; w2 = wNeu; b = bNeu; a = aNeu;
        if (DEBUG) System.out.println ("Paar: " + w2.toString() + " --> " + b.toString() );
      }
      else
      {
        try { Thread.sleep(1000); } catch (InterruptedException ex) {}
      }


      // Eintragen in die Faktorbasis
      boolean geklappt = zerlegenEintragen(b, w2, istQuadrat);  // Aktualisiert die Faktor-Basis
      if (geklappt && anzahlRel >= 5)
      {
        // Extrapoliere Die Zeit nach oben
        long zeitNachher = System.currentTimeMillis();
        int anzahlRelevant = (int) ((primBasis.length - 1) * 0.90);  // -1 nicht mitzhlen; etwa 10% nicht gebraucht
        int anzahlErreicht = anzahlRel;
        if (anzahlErreicht > anzahlRelevant) anzahlErreicht = anzahlRelevant;
        // Heuristik: ungefhr Anteil bei oberen Primfaktoren, nicht genutzt
        // Methode length()-1 von BitSet gibt das hchste gesetzte Bit zurck. Wenn man jetzt
        // eine gute erste Relation hat, dann nimm diese zur Berechnung der Zeit, vor allem am
        // "Ende" kann hier die Restzeit besser geschtzt werden.
        double erledigt = (double) anzahlErreicht / anzahlRelevant;  
        long diff = zeitNachher - zeitVorher;
        // System.out.println("Bisher: " + diff + " und Anteil " + erledigt);
        long ganz = (long) (diff / erledigt);
        int dauerneu = Math.round((ganz - diff) / 1000);
        if (dauerneu != restsekunden)  // Anzeige nur, wenn Zeit verndert
        {
          System.out.println ("Restsekunden: " + dauerneu );  // gerundet, damit nicht so schwankt
          restsekunden = dauerneu;
        }
      }

      istQuadrat = !istQuadrat;
    }

  }


  // diese Methode ist fr die Effizenz sehr entscheidend - Versucht, in den ganzen Zahlen Z eine Wurzel zu ziehen
  // 30.09.2008: Berechnet jetzt nur eine Wurzel, andere am Schluss
  // 28.09.2010: da es fr RSA sinnlos ist, n zu addieren, hier modifiziert: Zuerst Test, ob 
  // dies in Z berhaupt ein Quadrat sein kann; wenn ja, dann mod. Heron-Verfahren mitnutzen
  // Bei Der Zahl 1512732637186330689679100908197028741080553 war die neue Lsung mit 17 min 11 sek
  // etwas schneller als die alte Lsung mit 18 min 7 sek.
  // 10.01.2010 Noch mod 3, 5, 7 betrachten; Zusatztest Vorbedingung: nicht durch p^2 teilbar
  private static final BigInteger istWurzel(BigInteger para, boolean zusatztest)
  {
    if ((para.getLowestSetBit() & 1) == 1) return null; // 2, 6 etc.
    // ist ungerade oder durch 4 teilbar
    int weg = para.getLowestSetBit();
    BigInteger q = para.shiftRight(weg);
    weg = weg >> 1;
    if (1 != q.and(SIEBEN).intValue()) return null;  //  kein Quadrat

    // Heuristik mit 3 ,5 und 7: Rest muss Quadat mod p sein. Annahme: nicht durch p^2 teilbar !!!
    // Anm.: der frher hier stehende Primzahltest hat sich als zu teuer erwiesen !
    if (zusatztest)
    {
      for (short i = 0; i < 5; ++i)
      {
        BigInteger teile = null;
        switch (i)
        {
          case 0: teile = DREI; break;
          case 1: teile = FUENF; break;
          case 2: teile = SIEBEN; break;
          case 3: teile = ELF; break;
          case 4: teile = DREIZEHN;
        }
        int rr = q.remainder(teile).intValue();
        if (rr == 0) continue; // Achtung: sollte eigentlich nicht auftreten
        if ((rr == 1) || (rr == 4) || (rr == 9)) continue;  // sicher Quadrat 
        // Weitere nicht-triviale Tests
        if ((i == 2) && (rr != 2)) return null;  // 7
        else if ((i == 3) && (rr != 3 && rr != 5)) return null;  // 11
        else if ((i == 4) && (rr != 3 && rr != 10 && rr != 12)) return null;  // 13
      }
    } // Ende Zusatztest

    // Jetzt heron mitbenutzen
    return holeWurzel(para, true);  // != null, wenn Wurzel ganzzahlig in Z
  } 


  // Ruft die Zerlegung von Quadraten auf, und trgt das Ergebnis bei Erfolg in eine Matrix ein
  // Eine Relation wird von rechts her gesehen, d.h. vom hchsten und damit seltensten Primteiler
  // Diese hchsten dienen zur Verrechnung von relationen, da ja mglichst kleine Primteiler
  // aufteilen sollen, bzw. die Relation eine Wurzel der -+1 erzeugen sollen
  // Liefert wahr zurck, wenn eine neue Relation gefunden wurde, sonst falsch
  public final boolean zerlegenEintragen (BigInteger q, BigInteger w, boolean istQuadrat)
  {
    int i, j, k, l;

    BitSet bel = new BitSet(primBasis.length);   // die aktuelle quadratfreie Faktorisierung
    if (! istQuadrat) bel.set(0);

    ParaBig p1 = new ParaBig(); p1.zahl = q;
    ParaBig p2 = new ParaBig(); p2.zahl = w;
    boolean gefunden = zerlege_quadrat (p1, p2, bel);  // alle Parameter mit Rckgabe

    if (! gefunden)
    {
      // konnte nicht komplett zerlegt werden, aber vielleicht ist Quadrat jetzt kleiner
      if (DEBUG2) System.out.println ("\t" + q.toString() + " konnte nicht komplett zerlegt werden");
      anzahlNichtGeschafft++;
      return false;
    }
    else  // fertig, komplett zerlegt
    {
      // Nachzustand: p1.zahl, p2.zahl und bel sind verndert; die Belegung hat mindestens einen ungerade Koeff.,
      // da am Ende schon Abprfung in zerlege_quadrat stattfand !
      w = p2.zahl;

      if (DEBUG3) 
      {
        System.out.print ("\t" + w + " = ");
        for (k = 0; k < primBasis.length; ++k)
        {
          if (bel.get(k)) System.out.print("," + primBasis[k]);
        }
        System.out.println("");
      }
      anzahlGeschafft++;

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

      int letzteVariable1 = letzteVariable(bel);

      // Schaue in die Matrix der Belegungen von den schlechten zu den guten rein
      // die Matrix ist dann geordnet, d.h. es ist eine Rechtsdreiecksmatrix
      for (j = anzahlRel - 1; j >= 0; --j)
      {
        if (DEBUG)
        {
          // Testen Quadrat (aus Belegung) und gespeicherte Wurzel !
          System.out.print("Vorh:\t");
          BigInteger vergleich = BigInteger.ONE;
          for (k = 0; k < primBasis.length; ++k)
          {
            // Hier Ausgabe der Belegungstafel !
            if (DEBUG) System.out.print(rel[j].get(k) ? "1" : "0");
            if (rel[j].get(k)) 
            {
              if (k == 0) vergleich = n.subtract(vergleich);
              else vergleich = vergleich.multiply(primBasis[k]).remainder(n);
            }
          }
          if (DEBUG) System.out.println("");
          BigInteger tmp = wurzeln[j].modPow(ZWEI, n);  // Quadrat aus Wurzel berechnet
          if (! vergleich.equals(tmp))
          {
            System.out.println ("Vorhandene Relation an " + j + " von " + anzahlRel + " ist falsch");
            System.out.println ("W^2 ist " + tmp.toString() + ", Vergleich ist " + vergleich.toString() ); 
            warten();
          }
        }  // Ende DEBUG

        // welches ist die letzte benutzte Primzahl dieser Relation ?
        int letzteVariable2 = letzteVariable(rel[j]);
        if (letzteVariable2 > letzteVariable1) 
        {
          if ((j == 0) || (letzteVariable1 > letzteVariable(rel[j-1])))
          {
            break;  // an Stelle j einfgen, ohne zu verrechnen
          }
          else
          {
            continue;  // weitersuchen
          }
        }
        else if (letzteVariable2 == letzteVariable1)
        {
          // Neue verbessern
          w = w.multiply(wurzeln[j]).remainder(n);
          if (DEBUG) System.out.println ("Neu verrechnen mit " + j);
          for (k = 1; k < primBasis.length; ++k)
          {
            if (bel.get(k) && rel[j].get(k))
            {
              // nochmals Quadrate rauswerfen, ausser bei -1; Performanz-Behandlung bei 2 an [1]
              if (k == 1)
              {
                if (w.testBit(0)) w = w.add(n);
                w = w.shiftRight(1);
              }
              else
              {
                BigInteger durch;
                if ((k < invers.length) && (invers[k] != null))
                {
                  durch = invers[k];
                }
                else
                {
                  durch = primBasis[k].modInverse(n);  // ggT != 1 s.o. bei Berechnung Primbasis
                  if (k < invers.length) invers[k] = durch;  // sichern fr weitere Verwendung
                }
                w = w.multiply(durch).remainder(n);
              }
              bel.clear(k);  // Potenz wieder 0
            }
            else if (rel[j].get(k))
            {
              bel.set(k);  // es kommt eine Potenz dazu ! 
            }
          }
          // Vorzeichen getrennt (an Index 0) !!!
          if (bel.get(0) != rel[j].get(0)) bel.set(0); else bel.clear(0);
          letzteVariable1 = letzteVariable(bel);  // neu berechnen !
        }
        else if (letzteVariable2 < letzteVariable1)
        {
          ++j;break;  // schlechter als das aktuelle: danach einfgen
        }
      } // 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 (1^2 == (leer) )
      if (letzteVariable1 <= 0)
      {
        // if (anzahlRel > 0) System.out.println ("Neue Relation konnte komplett aufgelst werden: " + w);
        if (! bel.get(0)) teste_erfolg(w); // Wurzel der +1
        else testeWurzelMinusEins (w);
      }
      else if (anzahlRel == primBasis.length)
      {
        // hier muss ein Fehler vorliegen:
        System.out.println("Fehler in der Verrechnung: Gleichung nicht mehr speicherbar");
        for (i = 0; i < anzahlRel; ++i)
        {
          System.out.println (rel[i]);
        }
        System.out.println ("An " + j + " kommt " + bel);
        System.exit(1); 
      }
      else
      {
        // neue Relation bel in Zeile j abspeichern
        if (j < 0) j = 0;  // es gab bisher nichts
        if (j == 0) System.out.println ("--> Neuer Spitzenreiter: " + bel); // zur Info
        if (DEBUG) belegung_testen (w, bel);  // nur zum Test !

        if (j < anzahlRel)
        {
          if (DEBUG) System.out.println("Einfgen an " + j);
          for (k = anzahlRel; k > j; --k)
          {
            rel[k] = rel[k-1];
            wurzeln[k] = wurzeln[k-1]; 
          }
        }
        if (DEBUG)
        {
          System.out.println ("Belegen bei #" + j + " und Wurzel " + w.toString());
          for (k = 0; k < primBasis.length; ++k) 
          {
            System.out.print(bel.get(k) ? "1" : "0");
          }
        }
        rel[j] = (BitSet) bel.clone();  // ist clone() hier notwendig ?
        if (DEBUG) System.out.println("");
        wurzeln[j] = w;
        ++anzahlRel;
      }

      return true;
    }


  }


  // sucht nach der letzten Variable
  public int letzteVariable (BitSet bel)
  {
    for (int k = primBasis.length-1; k >= 1; --k)  // bis zur 2 an
    {
      if (bel.get(k)) return k;
    }
    return -1;
  }

  // prft, ob Teiler findbar (w ist eine Zahl, deren Quadrat +-1 ergibt), wenn ja, dann Ende
  private final void teste_erfolg(BigInteger w)
  {
    if (DEBUG2) System.out.println ("Nach " + anzahlRel + " Relationen:");
    if (w.equals(BigInteger.ONE) || w.equals(n.subtract(BigInteger.ONE)) )
    {
      if (DEBUG2) System.out.println ("Triviale Faktorisierung der 1");
    }
    else
    {
      System.out.println (w.toString() + " +- 1 liefert Teiler"); 
      BigInteger erg = (w.subtract(BigInteger.ONE)).gcd(n);
      if (! erg.equals(BigInteger.ONE))
      {
        System.out.println ("Faktor1 ist " + erg.toString() ); fertig = true;
      }
      erg = (w.add(BigInteger.ONE)).gcd(n);
      if (! erg.equals(BigInteger.ONE))
      {
        System.out.println ("Faktor2 ist " + erg.toString() ); fertig = true;
      }
      if (! fertig)
      {
        System.out.println ("Wurzel der 1 brachte keinen Erfolg");
        System.exit(1);
      }
      zeige_statistik();
      System.exit(0);
    }
  }


  // zeigt an, welche vorher gespeicherte Primfaktoren fr den Faktorisierungsprozess nicht gebraucht
  // wurden
  // Es scheint folgende Regel zu gelten: wenn Jacobi(2, n) = -1, dann braucht es nicht vorkommen
  // die anderen nicht verwendeten scheinen auch Jac(a, n) == 1 zu haben.
  private void zeige_statistik()
  {
    // noch ermitteln, welche Faktoren wirklich gebraucht wurden (ist ja jetzt tatschlich faktorisiert)
    // aber maximal 10, damit von andere Ausgabe noch was bleibt
    int gezeigt = 0;
    System.out.println ("Zahl ist " + n + ", die 10 kleinsten nicht gebrauchten Primzahlen:");
    for (int i = 0; i < primBasis.length; ++i)
    {
      boolean gef = false;
      for (short j = 0; j < anzahlRel; ++j) if (rel[j].get(i)) { gef = true; break; }
      if (!gef && gezeigt < 10) { 
        System.out.print (primBasis[i]+","); if ((gezeigt % 5) == 4) System.out.println(""); }
      if (! gef) ++gezeigt;
    }
    if (gezeigt > 0) System.out.println ("\n" + gezeigt + " von " + primBasis.length + " Faktoren nicht fr Matrix gebraucht");
  }


  // Testet die Belegung, und gibt einen Fehler aus, wenn nicht stimmt
  private final void belegung_testen (BigInteger w, BitSet bel)
  {
    // Testen, ob q zu vergleich passt  ...
    BigInteger vergleich = w.modPow(ZWEI, n); 
    BigInteger lauf = BigInteger.ONE;
    for (int j = 1; j < primBasis.length; ++j) // -1 nicht dabei !
    {
      if (bel.get(j)) lauf = lauf.multiply(primBasis[j]).remainder(n);
    }
    if (bel.get(0)) lauf = n.subtract(lauf);

    if (! vergleich.equals(lauf)) 
    { 
      System.out.println(w + " --> " + vergleich + " zu " + bel ); 
      System.exit(1);  // Inkonsistenter Zustand !
    }
  }


  // Zerlegt ein Quadrat über eine Faktorbasis. Es wird hierbei nicht der Kontext betrachtet !
  // Es wird ein Quadrat und eine Vergleichswurzel übergeben, und über der Primfaktor-Basis dargestellt
  // Im Gegensatz zu früheren Quadrat-Sieben erfolgt hier aber keine rekursive Bearbeitung.
  // arbeitet modifiziert mit m3mod4, das brachte 5% Ersparnis zur seq. Teilersuche
  private final boolean zerlege_quadrat (ParaBig para1, ParaBig para2 /* wegen Rckgabe */, BitSet bel)
  {
    if ((para1.zahl == null) || (para2.zahl == null) || (bel == null)) return false;
    int i, j;

    if (DEBUG) 
      System.out.println ("Aufruf zerlege_quadrat fuer " + para2.zahl.toString() + " --> " + para1.zahl.toString() );

    BigInteger q = para1.zahl; BigInteger w = para2.zahl;
    BigInteger rueck = q;  // rueck muss zu w passen, wobei q auch durch Faktoren der Potenz^1 div. wird !
    BigInteger dAR [];

    if (q.equals(BigInteger.ZERO)) // um bei Eingabe von Quadraten Endlosschleifen zu vermeiden
    {
      para2.zahl = BigInteger.ZERO;
      return true;
    }

    // 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)
    {
      if ((weg & 1) == 1) bel.set(1);  // eine 2 brig
      q = q.shiftRight(weg);  // ist nun ungerade
      if ((weg & 1) == 1) --weg;
      rueck = rueck.shiftRight(weg); // ungerade oder durch 2 teilbar (nun gerade Potenzen weg)
      weg = weg >> 1;  // fr die Wurzel nur die Hlfte weg
      while (weg > 0)
      {
        int moegl = w.getLowestSetBit();
        if (weg <= moegl) { w = w.shiftRight(weg); weg = 0;}
        else { w = w.shiftRight(moegl); weg -= moegl; w = w.add(n); }
      }
    }

    // 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 (DEBUG) System.out.println ("Vorzeitiger Abbruch wegen Primzahl " + q);
      return false;
    }

    // alle Potenzen von Primfaktoren raushauen ...
    boolean m3mod4 = q.testBit(1); int start3 = 2; int start1 = 2;
    for (j = 2; j < primBasis.length; ++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)
      {
        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
      }
      if (wieoft >= 2) // fr alle Quadrate die Wurzel anpassen, mit dem Inversen
      {
        // schaue nach, ob das Inverse schon da ist
        BigInteger inv;
        if ((j < invers.length) && (invers[j] != null))
        {
          inv = invers[j];
        }
        else
        {
          inv = f.modInverse(n);  // ggT != 1 s.o. bei Berechnung Primbasis
          if (j < invers.length) invers[j] = inv;  // sichern fr weitere Verwendung
        }
        while (wieoft >= 2)
        {
          if (DEBUG) System.out.println ("Weg Quadrat " + f.toString() );
          w = w.multiply(inv).remainder(n); wieoft -= 2;
          rueck = rueck.divide(f.multiply(f));
        }
      }
      if (wieoft > 0)
      {
        // kann eigentlich nur 1 sein, da alle 2er-Potenzen abgezogen wurden
        bel.set(j);
      }
      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 (DEBUG) System.out.println ("Abbruch " + j + " von " + primBasis.length); 
          return 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
    // Rckgabewerte (um Quadrate bereinigt): Achtung: q weicht von rueck ab, q <= rueck !!!
    
    if (! q.equals(BigInteger.ONE)) 
    {
      // eigentlich gehrt bei folgendem auch DEBUG dazu, damit nicht im Leistungsbereich aktiviert ist
      if ((DEBUG || DEBUG2 || DEBUG3) &&
          q.compareTo(primBasis[primBasis.length-1]) <= 0 /* && q.isProbablePrime(1) */) 
      { 
        // Fehlersituation, auch wenn q nicht tatschlich prim !
        System.out.println ("Quadrat " + q + " gefunden, das nicht in Basis ist"); System.exit(1);
      }
      // teste Quadrat von q
      BigInteger wur = this.istWurzel(q, true /* Zusatztest */);
      if (wur != null)
      {
        // Glücksfall: Aus dem Rest konnte noch Quadrat gezogen werden
        System.out.println ("Glcksfall: Rest ist Quadrat " + wur);
        BigInteger inv = wur.modInverse(n);
        w = w.multiply(inv).remainder(n); rueck = rueck.divide(wur.multiply(wur)); q = BigInteger.ONE;
        // Bitbelegung nicht verndern
      }
      else
      {
        if (DEBUG2) System.out.println ("Nicht zerlegbarer Rest: " + q + " ");
        // Test mit der kleinsten nicht aufgenommen Primzahl, ob Primbasis vollstndig:
        if (nqTestzahl != null && q.remainder(nqTestzahl).equals(BigInteger.ZERO))
        {
          System.out.println ("Fehler: Quadrat " + q + ", nicht aufgenommener Teiler " + nqTestzahl); System.exit(1);
        }
        return false;   // konnte nicht komplett zerlegt werden
      }
    }

    // Wenn nur gerade Potenzen und Wurzel zur Verfgung, dann hier evt. Faktorisierung mglich
    // falls nur -1, ist auch eine Faktorisierung mglich.
    if (q.equals(BigInteger.ONE))
    {
      for (j = 1; j < primBasis.length; ++j) if (bel.get(j)) break;  // -1 nicht dabei
      if (j == primBasis.length)
      {
        // kann man rueck faktorisieren ?
        if (!bel.get(0)) teste_erfolg(w); // weil es keine Erkenntnis bringt
        else testeWurzelMinusEins(w);
      }
    }

    // Rckgabewert ist die evt. modifizierte Wurzel, und die Belegung der Primfaktor-Basis
    para1.zahl = rueck; para2.zahl = w;
    return true;
  }

  // Berechnet mithilfe der Gauss'schen Zahlen einen ggT mit der Wurzel von -1
  // Bildet ggt fr (n, w + imag). Mathematisch sinnvolle Vorbedingung: n, p und q jeweils 1 mod 4
  private final void testeWurzelMinusEins(BigInteger w)
  {
    System.out.println ("Wurzel der -1");
    {
      BigInteger vergl = w.modPow(ZWEI, n);
      if (! vergl.equals(n.subtract(BigInteger.ONE))) { System.out.println ("stimmt nicht"); System.exit(1); }
    }

    BigInteger a = n; BigInteger b = BigInteger.ZERO;
    BigInteger c = w; BigInteger d = BigInteger.ONE;
    BigInteger e = null, f = null;
    BigInteger betr1 = n.multiply(n); BigInteger betr2 = (w.multiply(w)).add(BigInteger.ONE);
    while (betr1.compareTo(betr2) > 0)
    {
      // Wegen Wurzelberechnung (Betrge liegen nur im Quadrat vor) mit int rechnen, sofern mglich
      // System.out.println (a + " und " + b + " durch " + c + " und " + d);
      int erg = (betr1.divide(betr2)).intValue();
      if (erg < 0) { System.out.println ("berlauf bei int in testeWurzelMinusEins"); System.exit(1); } // Mist, klappt nicht
      // else System.out.println("Erg ist " + erg);
      BigInteger faktor = BigInteger.valueOf((int) Math.sqrt(erg));  // wird autom. abgerundet, wird auch so gebraucht
      e = a.subtract(faktor.multiply(c)); f = b.subtract(faktor.multiply(d));
      // neue Werte zuweisen
      a = c; b = d; c = e; d = f;
      betr1 = betr2; betr2 = (e.multiply(e)).add(f.multiply(f));
    }
    if (betr1.equals(betr2))
    {
      // n konnte in Summe zweier Quadrate zerlegt werden. Bei Primzahlen 1 mod 4 gegeben
      // Achtung: a oder b knnen hier negativ sein, ist aber fr weiteres hoffentlich kein
      // Problem
      System.out.println ("a = " + a + ", b = " + b); 
      // Bei Zusammengesetzt knnte man mit ECM ggf. eine Gruppenordnung erraten, und damit
      // ein Element mit einem nicht-trivialen ggT erzeugen

      // Sind die Primteiler p und q beide 3 mod 4, dann lassen sich die Primteiler hier nach der
      // einfachen Formel gewinnen:  a = (p + q)/2, b = (p - q) / 2
      BigInteger versuch = a.add(b).gcd(n);
      if (! versuch.equals(BigInteger.ONE))
      {
        System.out.println ("Faktor1 ist " + versuch); fertig = true;
      }
      versuch = a.subtract(b).gcd(n);
      if (! versuch.equals(BigInteger.ONE))
      {
        System.out.println ("Faktor2 ist " + versuch); fertig = true;
      }
      if (fertig)
      {
        zeige_statistik();
        System.exit(0);
      }
      else
      {
        System.out.println ("Entweder Primzahl oder beide Faktoren nicht 3 mod 4");
      }
    }
    else
    {
      System.out.println ("Verfahren gescheitert"); System.exit(1);
    }
  }


  // berechnet Wurzel mit Heron-Verfahren.
  // wenn angegeben nurExakt, wird Wurzel nur bei exakter Gleichheit zurckgeliefert
  // x[n+1] = (x[n] + q/x[n]) / 2
  private final static BigInteger holeWurzel (BigInteger eingabe, boolean nurExakt)
  {
    // System.out.println ("Wurzel der " + eingabe + " exakt " + nurExakt);
    if (eingabe.signum() < 0) return null;
    if (eingabe.compareTo(ZWEI) < 0) return eingabe;  // 0, 1
    // Fr Startwert bilde die obersten 2 Bit ab (auf die Hlfte des Exponenten)
    BigInteger wert = BigInteger.ONE;
    for (int i = eingabe.bitLength() -1, tmp = 0; (i >= 2) && (tmp < 2); --i)
    {
      if (eingabe.testBit(i) && !wert.testBit(i >> 1))
      {
        wert = wert.setBit(i >> 1); ++tmp;
      }
    }
    BigInteger erg = null, vorgaenger = null;
    BigInteger rueck [];

    // Schleife durchlaufen. Divsision arbeitet ganzzahlig
    // bei 35: Hier Oszillation zwischen 6 --> 5 --> 6 ) !
    while (true) {
      rueck = eingabe.divideAndRemainder(wert);
      if ((rueck[0].equals(wert)) &&
          (rueck[1].equals(BigInteger.ZERO)) ) return wert;  // exakt gefunden
      // nchsten Folgewert berechnen
      erg = rueck[0].add (wert);
      erg = erg.shiftRight(1);
      // System.out.println (erg.toString() );  // fr Testzwecke, zum Sehen der Approximation
      if (erg.equals(wert)) 
        return (nurExakt ? null : wert);  // keine Vernderung mehr, aber Rest

      if ((vorgaenger != null) && (erg.equals(vorgaenger)))
      {
        if (nurExakt) return null;
        if (wert.compareTo(erg) < 0) return wert; else return erg;  // immer das kleinere nehmen
      }
      vorgaenger = wert; wert = erg;
    }
  } // Ende Methode holeWurzel

  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) {}
  }

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


}

// einfache, venderbare Hllklasse
class ParaBig
{
 public BigInteger zahl;
}

// wegen Shutdown
class QuadratSieb20Haken extends Thread
{
  private QuadratSieb20 w;
  QuadratSieb20Haken(QuadratSieb20 p) { w = p; }
  public void run()
  {
    if (! w.fertig)
    {
      w.fertig = true;  // damit oben nicht mehr weitergemacht wird !

      System.out.println ("Abbruchstand:\na= " + w.a.toString() + "\nb= " + w.b.toString() );
      System.out.println ("w1= " + w.w1.toString() + "\nw2= " + w.w2.toString() );
    }
    if ((w.anzahlGeschafft > 0) || (w.anzahlNichtGeschafft > 0))
    {
      System.out.println ("Zelegt: " + w.anzahlGeschafft + ", nicht " + w.anzahlNichtGeschafft);
      System.out.println ("Erfolgsfaktor " + ((float) w.anzahlGeschafft / (w.anzahlGeschafft + w.anzahlNichtGeschafft)));
    }
  }
}

/*
Statistik (ber nicht gebrauchte Primfaktoren, was den Faktor 90% erklren soll)

101 Bit
54 / 302 = 17.88 %
39 / 316 = 12.34 %
49 / 302 = 16,22 %

111 Bit
60 / 499 = 12,02 %
86 / 499 = 17,12 %

121 Bit
97 / 745 = 13,02 %
79 / 779 = 10,14 %

131 Bit
158 / 1999 = 7,90 %
165 / 1149 = 14,36 %

141 Bit
238 / 1750 = 13,6 %
252 / 1824 = 13,81 %
*/