// Erstelldatum: März 2011
// Fachliche Anmerkungen:
// Arbeitet mit Polynomen über endlichen Körper (vorerst nur mod 2). Das Verfahren heisst Berlekamp
// und baut intern eine Matrix auf. Beginnend vom grad n/2+1 bis n-1 werden Polynome konstruiert.
// Liegt eine lineare Abhängigkeit vor, kann ein Teilerpolynom gefunden werden.

// Technische Anmerkungen:
// In einem long-Wert passen 64 Bit rein. Polynome werden durch die Potenzen eingeben, getrennt von leer oder Komma.
// Da die Lösung mit dem linearen Gleichungssystem bestimmte Gleichungen (z.B. x^2 + 1 = 0, x^6 + x^2 + 1 = 0)
// nur mit b(x) == 1 lösen kann, wird hier noch der Sonderfall zur Erkennung von Quadratpolynomen gemacht
// (nur gerade Potenzen). Ob dann damit immer eine Lösung gefunden werden kann, ist mir unbekannt.
// Januar 2014: Durchsicht, kleinere technische Verbesserungen
// Literatur: Wolfram Koepf, "Computeralgebra", Kapitel 8.2, Springer Verlag 2006

import java.awt.*;
import java.util.StringTokenizer;

  // Hilfsklasse für PolynomMod2
  class PTyp2
  {
    public int grad;
    public long poly;
    public PTyp2 () { grad = -1; poly = 0; }
    public PTyp2 (int p_grad, long p_poly)
    { 
      String fehler = null;
      if (p_grad < -1) fehler = "Grad < 0";
      else if (p_grad > 62) fehler = "Grad > 62";
      else if (p_grad >= 0 && (p_poly & (1L << p_grad)) == 0) fehler = "Grad passt nicht zu Polynom";
      if (fehler != null) { System.out.println("Fehler(" + fehler + ") bei Grad " + p_grad + " und Polynom " + p_poly); 
        System.exit(1); }
      grad = p_grad; poly = p_poly; 
    }
    public int grad_bestimmen()
    {
      int wert = -1; long ak = poly;
      if (poly != 0) while (ak != 0) { ++wert; ak = ak >>> 1; }
      return wert;
    }
    protected Object clone ()
    {
      PTyp2 erg = new PTyp2();
      erg.grad = this.grad; erg.poly = this.poly;
      return erg;
    }
  }

public class PolynomMod2 extends Frame
{
  private final boolean mitQuadrattest = true;  // siehe Beschreibung oben
  private boolean FAKDEBUG = true;   // Ausführliche Infomationen z.B. Polynomstufe, Verrechnungsmatrix

  private Label leing;
  private TextField teing;
  private Checkbox cMitMatrix;
  private Checkbox cbox;
  private Checkbox cMitDebug;
  private Checkbox cAusgabe;
  private Label lausg;
  private TextField tausg;
  private Button bstart;

  // Konstruktor, legt das Fenster an
  public PolynomMod2(boolean keineGUI)
  {
    if (keineGUI) return;
  	setTitle("Polynom mod 2 - 64 Bit - Faktorisierer");
    setLayout(new GridLayout(3,3));

    // 1. Zeile
    leing = new Label("Eingabe (Liste Potenz)"); this.add(leing);
    teing = new TextField("0 4 5"); this.add(teing);
    cMitMatrix = new Checkbox("Matrix-Ausgabe"); this.add(cMitMatrix);

    // 2. Zeile    
    cbox = new Checkbox("Kreisteilungspolynom"); this.add(cbox);
    cMitDebug = new Checkbox("mit Konsole-Meldungen"); this.add(cMitDebug);
    cAusgabe = new Checkbox("Ausgabe nur Potenzen"); this.add(cAusgabe);

    // 3. Zeile
    lausg = new Label("Ausgabe"); this.add(lausg);
    tausg = new TextField(); tausg.setEditable (false); this.add(tausg);
    bstart = new Button("starten"); this.add(bstart);

    /*this.setSize (new Dimension(450, 100));*/ this.pack();
    this.setVisible(true);
  }


  // AWT-Ereignisbehandlung 1.0 - Achtung: läuft ohne Threads
  public boolean handleEvent (Event evt)
  {
    if (evt.id == Event.WINDOW_DESTROY) { System.exit(0); return true; }
    return super.handleEvent(evt);   
  }


  // Hier werden nur Ereignisse vom Typ "Action" angesteuert
  public boolean action (Event evt, Object obj)
  {
    if (evt.target == bstart && evt.id == 1001)
    {
      // Checkboxen abfragen
      FAKDEBUG = cMitDebug.getState();
      boolean nurPotenzen = cAusgabe.getState();

      String text = teing.getText().trim(); tausg.setText("");
      long poly = 0; int grad = -1; boolean fehler = false;
      // Kreisteilungspolynom, ein Koeffizient
      if (cbox.getState())
      {
        try
        {
          grad = Integer.parseInt(text);
          // kurzer Test
          if (grad < 0 || grad > 62) fehler = true;  // Potenz zu groß/klein
        } catch (NumberFormatException ex) { fehler = true; }
        if (fehler) { tausg.setText("ungültiger Koeffizient"); return true; }
        poly = ((1L << grad) | 1); teing.setText(text + " 0"); cbox.setState(false);
      }
      else
      {
        StringTokenizer st = new StringTokenizer(text, ", ");  // Eingabe durch Leerzeichen oder Komma getrennt
        while (st.hasMoreTokens())
        {
          String teil = st.nextToken();
          try
          {        
             int pot = Integer.parseInt(teil);
             // kurzer Test
             if (pot < 0 || pot > 62) fehler = true;  // Potenz zu groß/klein
             else if ((poly & (1L << pot)) != 0) fehler = true;  // schon gesetzt
             else { poly = poly | (1L << pot); if (pot > grad) grad = pot; }
          } catch (NumberFormatException ex) { fehler = true; }
        }
        if (fehler) { tausg.setText("ungültiges Polynom"); return true; }
      }


      // Polynom erzeugen, und faktorisieren
      System.out.println ("Zu tun: Grad " + grad + " und Polynom " + poly);
      PTyp2 pol = new PTyp2(grad, poly);
      PTyp2 ziel = faktorisieren(pol);
      // Noch Koeffizienten der Ausgabe auflisten
      if (ziel != null) 
      {
        String neu = ""; long lauf = 1; int i = 0;
        while ((lauf >= 0) && (lauf <= ziel.poly))
        {
          if ((ziel.poly & lauf) != 0) 
          {
            if (!nurPotenzen)
            { if (neu.length() > 0) neu += "+ "; neu += "x^"; }
            neu = (neu + (new Integer(i)).toString() + " ");
          }
          lauf <<= 1; ++i;
        }
        tausg.setText (neu);
      }
      else tausg.setText ("irreduzibel");

      return true;
    }

    return false; // super.handleEvent(evt);
  }

  // multiplizieren xor 2 ohne modulo
  public static PTyp2 multiplizieren (PTyp2 p1, PTyp2 p2)
  {
    if (p1.grad < 0) return p2;
    else if (p2.grad < 0) return p1;

    // Grad ist Summe der Ausgangswerte
    long summe = 0;
    long a = p1.poly; long b = p2.poly;
    if (a > b) { long hilfe = a; a = b; b = hilfe; }  // ggf. tauschen, um schneller zu sein
    while (a > 0)
    {
      if ((a & 1) == 1) summe = summe ^ b;
      a >>>= 1; b <<= 1; if (b < 0) { System.out.println("Zahlenüberlauf"); System.exit(1); }
    }
    PTyp2 erg = new PTyp2(p1.grad + p2.grad, summe);

    return erg;
  }


  // teilen, und ermittelt den Rest
  public static PTyp2[] teilenUndRest (PTyp2 p1, PTyp2 p2)
  {
    if (p2.grad == -1) return null;  // keine Div. durch 0

    PTyp2 rueck [] = new PTyp2[2];  // erste Komp.: Teilung, zweite der Rest
    // Ist ist was zu dividieren
    int diff = p1.grad - p2.grad;
    if (diff < 0) { rueck[0] = new PTyp2(); rueck[1] = p1; return rueck; }

    long setzen = 1;  // zum Setzen der Bits für das Ergebnispolynom
    PTyp2 erg = new PTyp2();
    PTyp2 b = (PTyp2) p2.clone();
    if (diff > 0)
    {
      b.grad += diff; b.poly = b.poly << diff; setzen = setzen << diff;
    }
    // haben jetzt gleichen Grad
    PTyp2 a = (PTyp2) p1.clone();
    // while (a.grad >= p2.grad)
    for (int i = 0; i <= diff; ++i)
    {
      if (a.grad == b.grad) // vergleiche a und dem geschobenen Polynom b
      {
        a.poly = a.poly ^ b.poly; a.grad = a.grad_bestimmen(); erg.poly += setzen;
        if (a.grad > p1.grad) { System.out.println ("Logikfehler"); System.exit(1); }
      }
      b.grad--; b.poly = b.poly >>> 1; setzen = setzen >>> 1;
    }

    erg.grad = erg.grad_bestimmen();
    rueck[0] = erg; rueck[1] = a;

    return rueck;
  }

  // ggT
  public static PTyp2 ggT (PTyp2 p1, PTyp2 p2)
  {
    if (p1.grad < p2.grad)
    {
      PTyp2 hilfe = p1; p1 = p2; p2 = hilfe;  // vertauschen
    }
    while (p2.grad > 0)
    {
      PTyp2 neu [] = teilenUndRest (p1, p2); 
      if (neu == null) return null;
      if (neu[1].grad == -1) return p2;
      p1 = p2; p2 = neu[1];
    }
    return p2;
  }

  public static PTyp2 kreisteilungspolynom(int p_grad)
  {
    if (p_grad > 62) { System.out.println ("Zu großes Polynom"); return null; }
    PTyp2 neu = new PTyp2(); neu.grad = p_grad; 
    neu.poly = (1L << p_grad) | 1;  // Fehlerquelle: 1 << p_grad bleibt (neg.) int !
    if (neu.poly < 0) { System.out.println ("Überlauf"); System.exit(1); }
    return neu;
  }

  // Polynom vom Grad n, Teilerpolynom vom Grad m, Faktorisierung mit Matrix
  // (b(x)^m - b(x)) kongruent 0 mod a(x), a(x) gegeben, b(x) gesucht. Rückgabe null, wenn irreduzibel.
  // Ist nicht rekursiv, d.h. wenn 2 Teile erzeugt wurden, wird dieser zurückgeliefert.
  public PTyp2 faktorisieren(PTyp2 eingabe)
  {
    int i, j, k; long sp;

    // Feld: zum Speichern der m nicht-trivialen Polynome. Matrix: n Zeilen pro x-Poly, zur Best. der Koeff bx
    // System.out.println ("Eingabe: " + eingabe.grad + " und Polynom " + eingabe.poly);
    if (eingabe.grad <= 0) return eingabe;

    if (mitQuadrattest)
    {
      // wenn alle Potenzen gerade sind, ist Faktorisierung mod 2 einfach
      long versuch = 0; long muster = 1; long setzen = 1; int stufe = 0; boolean moeglich = true;
      while (moeglich && (muster >= 0) && (muster <= eingabe.poly))  // muster >= 0, damit bei Umschalten 62 auf 63 Ende
      {
        if ((eingabe.poly & muster) != 0)
        {
          if ((stufe & 1) == 1) moeglich = false;  // ungerade Potenz
          else versuch = versuch | setzen;  // gerade Potenz, für Ergebnis übernehmen
        }
        muster = muster << 1; ++stufe; if ((stufe & 1) == 1) setzen = setzen << 1;
      }
      if (moeglich) { PTyp2 rueck = new PTyp2(stufe >> 1, versuch); return rueck; }
    }

    int n = eingabe.grad;
    long mat [] = new long[n];
    int m = (n >> 1); if ((n & 1) == 1) ++m; // n = 6 --> 3, 5 --> 3, mind. einmal mod rechnen
    PTyp2 feld [] = new PTyp2[n-1 /* m */];  // hier bereits max. dimension.

    // Es scheint nicht immer ab m eine Faktorisierung möglich, wohl aber bei höhreren, also m bis n-1 gehen
    for (; m < n; ++m)
    {
      if (FAKDEBUG) System.out.println("\nZum Grad " + m);

      // a) Teilpolynome bestimmen: Das Polyom für b0(=1) nicht abspeichern, also x^2p, x^(2p-2), ... x^(2)
      int berAb = 0; PTyp2 ak = null;
      if ((m > 1) && (feld[m-2] != null))
      {
        // inkrementell, [m-2] bereits vorhanden
        ak = (PTyp2) feld[m-2].clone(); berAb = m - 1;
      }
      else
      {
        // komplett von vorne
        ak = new PTyp2(0, 1);
      }
      for (i = berAb; i < m; ++i)
      {
        ak.grad += 2; ak.poly <<= 2;
        PTyp2 [] erg = teilenUndRest(ak, eingabe);
        if (FAKDEBUG) System.out.println("x^" + ((i+1) << 1) + ": Grad " + erg[1].grad + ", Polynom " + erg[1].poly);
        feld[i] = (PTyp2) erg[1].clone();  // wichtig, da teilenUndRest ggf. this zurückliefert
      }
      // for (i = 0; i < m; ++i) System.out.println("Grad " + feld[i].grad + ", Polynom " + feld[i].poly);    

      // Matrix für die Koeffizienten bestimmen
      // mat[0] = x^0, mat[1] = x, mat[2] = x^2, mat[n-1] = x^(n-1)
      // b) Die Koeffizienten aus den Gleichungen rausholen. Erste zu b1, zweite zu b2 u.s.w.
      // Durchlaufe alle m erzeugten Polynome, schaue nach, welche x-Potenz, und dann in der Zeile eintragen
      for (i = 0, sp = 2 /* b1 */; i < m; ++i, sp <<= 1)
      {
        long lauf = 1;
        // System.out.println(feld[i].poly + " und Spalte " + sp);
        for (j = 0; j <= feld[i].grad; ++j, lauf <<= 1)
        {
          if ((feld[i].poly & lauf) != 0) mat[j] = mat[j] ^ sp;
        }
      }

      // c) noch minus b(x) = b^m*x^m + b^m-1*x^m-1 + b1 * x + b0, b0 nicht gespeichert
      long lauf = 2;  // in Matrix diagonal laufen, ist aber nicht quadr.
      for (i = 1; i <= m; ++i, lauf <<= 1) mat[i] = mat[i] ^ lauf;

      // if (FAKDEBUG) for (i = 0; i < n; ++i) System.out.println (mat[i]);

      // d) Matrix verrechnen, und Untersuchung abschliessen. m Spalten (b0 egal), n zeilen, also Redundanz
      for (i = 0, lauf=2 /* b1 */; i < m; ++i, lauf <<= 1)  // lauf = Spalte, i,j = Zeile
      {
        // alle Zeilen durchsuchen, bis das erste Bit gesetzt ist
        for (j = i; j < n; ++j) if ((mat[j] & lauf) != 0) break;
        if (j < n)
        {
          // gib diese Zeile nach oben, ggf. Zeile j mit Zeile i tauschen
          if (j != i) { long hilfe = mat[i]; mat[i] = mat[j]; mat[j] = hilfe; }
          // xore alle Zeilen k drunter, falls nötig
          for (k = j+1; k < n; ++k)
          {
            if ((mat[k] & lauf) != 0) mat[k] = mat[k] ^ mat[i];
          }
        }
      }

      if (cMitMatrix != null && cMitMatrix.getState())  // gibt die Matrix aus
      {
        System.out.println ("Matrix verrechnet: x^k/b1 ... bm");
        for (i = 0; i < n; ++i) // zwar gibt es n-x-Potenzen, aber nur m relevante Zeilen !
        {
          lauf = 2; String zeile = "x^" + i + " "; // bei b1 anfangen, b0 immer 1
          for (j = 1; j <= m; ++j, lauf <<= 1) zeile += ( ((mat[i] & lauf) != 0) ? "1": "0");
          System.out.println (zeile);
        }
      }

      // Da auf der rechten Seite 0er stehen, muss es mind. eine lin. Abhängigkeit geben.
      // Bei m Zeilen (b0 ist immer frei wählbar) gibt es nur die 0-Lösung, bei weniger eine nicht-triviale
      boolean gefunden = (mat[m-1] == 0);
      // System.out.print ("Faktorisierung ");
      if (gefunden)
      {
        // Koeffizienten b bestimmen, von b(m) bis b(1), b(0) = 1
        lauf = (1L << m) /* L ist wichtig ! */; long hoechst = lauf; long polyneu = 1 | hoechst;
        // Letzte vorhandene Gleichung ermitteln, mat[m-1] ist ja leer !
        for (i = m-2; i >= 0; --i)
        {
          if (mat[i] == 0) continue;  // Leerpolynome unten in der Matrix überspringen
          boolean bisher = false;  // 0, nun die bisher ermittelten benützen
          long bis = 2; while ((mat[i] & bis) == 0) bis <<= 1;  // minimum b in dieser Gleichung bestimmen
          for (long lauf2 = hoechst; lauf2 > lauf; lauf2 >>>= 1) 
            if (((polyneu & lauf2) != 0) && ((mat[i] & lauf2) != 0)) bisher = !bisher;
          // zu diesem Polynom: Gibt es nächsten tieferen Koeffizienten auch darin, dann lege den aktuellen mit 1 fest
          // if (FAKDEBUG) System.out.println ("Gl " + i + ": " + lauf + " bis " + bis);
          while (lauf > bis) { polyneu = polyneu | lauf; if ((mat[i] & lauf) != 0) bisher = !bisher; lauf >>>= 1; }
          // lauf liegt nun auf dem minimalsten Koeffizienten bis
          if (bisher) polyneu = polyneu | lauf; 
          lauf >>>= 1;
        }
        PTyp2 pneu = new PTyp2(m, polyneu);
        System.out.print ("Ermittelt bei Grad " + m + ": ");
        System.out.println (polyneu);
        // Polynom erzeugt, jetzt muss noch ggT geprüft werden, da b(x) ggf. zusätzlichen Faktor haben kann; b0 bestimmmen
        PTyp2 erg = ggT(eingabe, pneu);
        if (erg.grad > 0) { System.out.println ("Teiler: Grad " + erg.grad + " und Polynom " + erg.poly + " gefunden"); return erg; }
        // 2te Möglichkeit: b0 = 0
        pneu.poly &= 0xfffffffe;
        erg = ggT(eingabe, pneu);
        if (erg.grad > 0) { System.out.println ("Teiler: Grad " + erg.grad + " und Polynom " + erg.poly + " gefunden"); return erg; }        
        System.out.println ("versucht, aber ohne Erfolg"); return null;
      }
      // System.out.println("nicht gefunden");
      for (i = 0; i < n; ++i) mat[i] = 0;  // für nächsten Durchlauf wieder löschen
    }

    return null; // nichts gefunden
  }


  public static void main (String[] args)
  {
    boolean keineGUI = false; // true;
  	PolynomMod2 p = new PolynomMod2(keineGUI); 
    if (keineGUI) p.test();
  }

  public void test()
  {
    boolean einzeltest = true;
  	if (einzeltest)
    {	
  	  PTyp2 p1 = new PTyp2(9, 769); // s.u.
      PTyp2 p2 = new PTyp2(8, 425); // s.u.
      PTyp2 erg  = multiplizieren(p1, p2);
      PTyp2 ziel = faktorisieren(erg);
      if (ziel != null) System.out.println ("Ergebnis: " + ziel.grad + " und " + ziel.poly);
      else System.out.println ("Fehler, nicht faktorisiert !");      
    }
  	else
  	{	
      // Lauf von relativ weit oben
  	  for (long lauf = Long.MAX_VALUE; lauf > Long.MAX_VALUE - 20000; lauf -= 2)
      {
        PTyp2 erg = new PTyp2(62, lauf);
        // System.out.println (erg.grad + " " + erg.poly);
        System.out.println("");
        PTyp2 ziel = faktorisieren(erg);
        if (ziel != null) System.out.println ("Ergebnis: " + ziel.grad + " und " + ziel.poly);
        else System.out.println ("ist irreduzibel");
      }  
    }
  }

}

/*
Hier sind die Codierungen für einige Polynome (Grad, Dualwert)

Irreduzible zu Grad 8 ist 369
Irreduzible zu Grad 8 ist 375
Irreduzible zu Grad 8 ist 379
Irreduzible zu Grad 8 ist 391
Irreduzible zu Grad 8 ist 395
Irreduzible zu Grad 8 ist 397
Irreduzible zu Grad 8 ist 415
Irreduzible zu Grad 8 ist 419
Irreduzible zu Grad 8 ist 425

Irreduzible zu Grad 9 ist 661
Irreduzible zu Grad 9 ist 665
Irreduzible zu Grad 9 ist 675
Irreduzible zu Grad 9 ist 677
Irreduzible zu Grad 9 ist 687
Irreduzible zu Grad 9 ist 695
Irreduzible zu Grad 9 ist 701
Irreduzible zu Grad 9 ist 719
Irreduzible zu Grad 9 ist 721
Irreduzible zu Grad 9 ist 731
Irreduzible zu Grad 9 ist 757
Irreduzible zu Grad 9 ist 761
Irreduzible zu Grad 9 ist 769
Irreduzible zu Grad 9 ist 787

*/