// Repräsentiert ein Polynom im Bereich von Z, also nicht modulo gerechnet. Bei modulo
// siehe PolynomModZ
// Siehe das gute Buch:
// Jörg Bewersdorf, "Algebra für Einsteiger",
// Vieweg Verlag, 3. Auflage, Fussnote 32 Kapitel 6

// Änd. Nov. 14: kleinere Fehlerkorrekturen, Sortierung der NS,
// Kombinatorik schneller; Test bis P'Grad 60

  // Hilfsklasse für PolynomModPZ
  class PTypZ implements Cloneable
  {
    public int grad;
    public Komplex [] poly;

    public PTypZ () { grad = -1; poly = null; }
    public PTypZ (int p_grad, long [] p_poly) throws Exception
    {
      String fehler = null;
      if (p_grad < -1) fehler = "Grad < 0";
      else if (p_grad > 15876)  // 4 * 63 * 63
    	  throw new Exception("Zu großer Grad " + p_grad);
      else if (p_grad >= 0 &&
         (p_poly == null || p_poly.length <= p_grad || p_poly[p_grad] == 0))
	    fehler = "kein Polynom, Grad passt nicht";
      if (fehler != null) throw new Exception("Fehler(" + fehler + ") bei Grad " + p_grad + " und Polynom " + p_poly);
      grad = p_grad;
      poly = new Komplex[grad+1];
      for (int i = 0; i < p_poly.length; ++i)
      {
    	poly[i] = new Komplex(p_poly[i]);
      }
    }

    public int grad_bestimmen()
    {
      int i = -1;
      for (i = grad; i >= 0; --i) if (! poly[i].istNull()) break;
      return i;
    }
    // tiefe Kopie (des Feldes)
    @Override
    protected Object clone ()
    {
      PTypZ erg = new PTypZ();
      erg.grad = this.grad;
      erg.poly = new Komplex[this.grad+1];
      for (int i = 0; i <= this.grad; ++i) erg.poly[i] = this.poly[i].clone();

      return erg;
    }

    // Gibt Polynom absteigend aus (hohe Koeffiz. zuerst), z.B. 1x^2 3x 0
    @Override
    public String toString()
    {
      StringBuilder rueckg = new StringBuilder(); // "Grad " + new Integer(grad).toString() + "; ";
      if (grad == -1) return "0";
      for (int i = grad; i >= 0; --i)
      {
        if (!poly[i].istNull())
        {
        	rueckg.append (" +(" + poly[i].toString() + ")");
        	if (i > 0) rueckg.append("x"); 
        	if (i > 1) rueckg.append( "^" + i);
        }
      }
      return rueckg.toString();
    }
  }

  // kapselt Logik für komplexe Zahlen
  class Komplex implements Cloneable
  {
    private final static double EPS = 1e-7;  // große Gen.
    public final static double EPS2 = 1e-4;  // kleine Gen.

	public double a;
	public double b;
	public Komplex (long real) { a = real; }
	public Komplex (double real) { a = real; }
	public Komplex (double real, double imag) { a = real; b = imag; }
    @Override
	public Komplex clone()
    {
    	Komplex rueck = new Komplex(a);
    	rueck.b = b;
    	return rueck;
    }
    @Override
    public boolean equals(Object vergl)
    {
    	if (vergl == null || ! (vergl instanceof Komplex)) return false;
    	Komplex was = (Komplex) vergl;
        if (Math.abs(a - was.a) > EPS) return false;
        if (Math.abs(b - was.b) > EPS) return false;

        return true;
    }
    @Override
    public String toString()
    {
      String ausg = new Double(a).toString();
      if (Math.abs(b) > EPS2) ausg += " + i * " + new Double(b);
      return ausg; 
    }
    public boolean istNull()
    {
      return Math.abs(a) < EPS2 && Math.abs(b) < EPS2;
    }
	public void add (Komplex vergl) {
		a += vergl.a; b += vergl.b;
	}
	public void div (Komplex vergl) {
		double betrag = vergl.a * vergl.a +
		  vergl.b * vergl.b;
		double real = (a * vergl.a + b * vergl.b);
		double imag = (b * vergl.a - a * vergl.b);
		a = real / betrag; b = imag / betrag;
	}
	public void mult (Komplex vergl) {
		double neu = (a * vergl.a - b * vergl.b);
		double neu2 = (a * vergl.b +  b * vergl.a);
		a = neu; b = neu2;
	}
	public void sub (Komplex vergl) {
		a -= vergl.a; b -= vergl.b;
	}
	
  // Ermittelt (mit red. Genauigkeit), ob ein anderes Element das komplex konjugierte ist
	public boolean istKonjugiert(Komplex wert2)
	{
		return (Math.abs(a - wert2.a) < Komplex.EPS2 &&
				Math.abs(b + wert2.b) < Komplex.EPS2);
	}
  }

// Liest Polynom mit Koeff. ein
  public class PolynomModZ {

		final static long poly [][] = {{4,11,+2,-4,-2,1},
				              {-2,8,-4,-4,1}, {-1, +1, +2,0,0,1},
				              {4,11,2,-4,-2,1}};      //
		final static boolean DEBUGKONFIG = false; // zeigt Komb. der Wurzeln an
		

	public static void main (String [] args) throws Exception
	{
		PolynomModZ pz = new PolynomModZ();

		// Entweder Eingabe ueber Kommandozeile oder eins von oben
    long [] ak = null;
		if (args.length >= 1)
    {
    	ak = new long[args.length];  // führende Ziffern zuerst
    	for (int i = args.length-1; i >= 0 ; --i)
    		ak[args.length-i-1] = Long.parseLong(args[i]);
    }
		else
		{
			ak = poly[0];  // zuf. eins ziehen
		}
		PTypZ neu = new PTypZ(ak.length-1, ak);
		System.out.println("Polynom: " + neu);

		long start = System.currentTimeMillis();
    Komplex [] rueck = pz.alleNullstellen(neu);
    for (int i = 1; i <= ak.length>>1; ++i)
    {
    	System.out.println("Grad" + i);
      PTypZ erg = pz.faktorisieren(rueck, i, neu);
      if (erg != null) { System.out.println("Teiler: " + erg); break; }
    }
		long ende = System.currentTimeMillis();
		System.out.println("Dauer sek: " + ((ende-start) / 1000) );
	}

    // Rechenfunktionen

	// Berechnet die Ableitung des Polynoms. Wichtig fuer Newton-Verfahren
	public static PTypZ ableitung (PTypZ eingabe) throws Exception
	{
	  if (eingabe.grad < 0) return eingabe;
	  PTypZ rueck = (PTypZ) eingabe.clone();
	  for (int i = 1; i <= rueck.grad; ++i)
	  {
	    rueck.poly[i-1].a = i * rueck.poly[i].a;
	    rueck.poly[i-1].b = i * rueck.poly[i].b;
	  }
	  rueck.poly[rueck.grad].a = 0;
	  rueck.poly[rueck.grad].b = 0;
	  rueck.grad = eingabe.grad - 1;

	  return rueck;
	}

	// wertet ein Polynom an der Stelle x aus
	public static Komplex auswerten (Komplex x,
			PTypZ poly)
	{
	  if (poly.grad < 0) return new Komplex(0);
	  Komplex ak = poly.poly[poly.grad].clone();
	  for (int i = poly.grad-1; i >= 0; --i)
	  {
		  ak.mult(x); ak.add(poly.poly[i]);
	  }
	  return ak;
	}

	// wertet ein Polynom an der Stelle x aus
	// Vorbedingung: Funktionswert an x ist 0.
	public static PTypZ abdividieren (Komplex x,
			PTypZ poly) throws Exception
	{
	  if (poly.grad < 0) return poly;
    if (poly.grad == 0) throw new Exception("Konstantes Polynom");

    long feld [] = new long [poly.grad]; feld[poly.grad-1] = 1;
    PTypZ rueck = new PTypZ(poly.grad-1, feld);
    // der höchste Koeffizient bleibt
    rueck.poly[poly.grad-1] = poly.poly[poly.grad];
    // System.out.println(rueck.poly[rueck.grad]);
	  for (int i = poly.grad-2; i >= -1; --i)
	  {
		  Komplex ak = rueck.poly[i+1].clone();
		  ak.mult(x); ak.add(poly.poly[i+1]);
		  // System.out.println(i + ": " + ak);
		  if (i >= 0) rueck.poly[i] = ak;
		  else if (!ak.istNull()) 
		  	System.out.println("Warnung: " + ak + " nach Abdividieren"); // sollte 0 sein
	  }

	  return rueck;
	}


	// Sucht eine komplexe Nullstelle mit Newton
	public static Komplex nullstelleNewton(
			Komplex start,
			PTypZ p1, PTypZ p2) throws Exception
	{
		if (p1.grad < 0 || p1.poly[0].istNull()) return new Komplex(0);
	  if (p1.grad == 0) throw new Exception("Konstante hat keine NS"); 

	  Komplex wo = start.clone();
	  Komplex alt = null; int i = 0;

	  // Es muss keine Konvergenz geben, z.B. bei doppelten
	  // Nullstellen, deshalb hier Zahl an Näherungen,
	  // Quadrat des Polynoms, da quadr. Konvergenz von Newton
	  final int maxIter = (p1.grad <= 5 ? 100 : 4 * p1.grad*p1.grad);
	  for (i = 0; i < maxIter; ++i)
    {
    	alt = wo.clone();
		  Komplex oben = auswerten(wo, p1);
      //System.out.println("Oben: " + oben);
      if (oben.istNull()) break;
      Komplex unten = auswerten(wo, p2);
      if (unten.istNull()) return null;  // Extremwert !
      oben.div(unten); wo.sub(oben); // Diff. bestimmen
    }

	  if (i == maxIter) 
	  {
	  	System.out.println("Warnung: Keine Konvergenz!");
	  	return null; // bei ungenauen Werten kann es zu Folgefehler kommen
	  }
    if (alt == null) return wo;
    // Gewichte besseren Wert doppelt so wie alten
	  wo.a = (2 * wo.a + alt.a) / 3;
	  wo.b = (2 * wo.b + alt.b) / 3;
	  return wo;  // Appr.
	}

    // Löst das Polynom über den komplexen Zahlen auf, d.h.
	// liefert alle Nullstellen, ggf. Probleme bei mehrfachen
	Komplex [] alleNullstellen(PTypZ poly) throws Exception
	{
	  if (poly.grad < 0) return new Komplex[0];
	  if (poly.grad == 0) throw new Exception("konstantes Polynom");

	  Komplex [] rueck = new Komplex[poly.grad];
    Komplex start = new Komplex(1.0, 1.0);  // wichtig: Imaginaerteil
    
    // das ursprüngliche Polynom und dessen Ableitung behalten
    PTypZ ursprung = poly;
    PTypZ ursprung_abl = ableitung(poly);

	  // Nur abdiv., solange Grad > 1
    for (int i = 0; i < rueck.length-1; ++i)
	  {
		  PTypZ abl = ableitung(poly);

		  Komplex wo = nullstelleNewton(start, poly, abl);
		  if (wo == null)
		  {
			  start.a += 1.0; --i; continue;
		  }
		  else
		  {
		  	// bei stimmten Polynomen wie x^24-9=0 kam es zu Rundungsproblemen, bei -1,058 wurde EPS2 verletzt, so dass das Polynom
		  	// nicht faktorisiert werden konnte. Folgendes wg. Genauigkeit, also nachiterieren:
		  	Komplex weiter = nullstelleNewton(wo, ursprung, ursprung_abl);
		  	if (weiter != null)
		  	{
		  	  if (Math.abs(wo.a - weiter.a) < Komplex.EPS2 &&
		  		  	Math.abs(wo.b - weiter.b) < Komplex.EPS2)
		  	  {
		  		  wo = weiter;  // sollte genauer im Orginal sein
		  	  }
		  	  else System.out.println("Hinweis: Nachiterieren erfolglos " + wo + " und " + weiter);
		  	}  
		  	
		  }
		  System.out.println(i + "-te Ns ist " + wo);
		  rueck[i] = wo;
		  poly = abdividieren(wo, poly);
		  start.a = wo.a; start.b = -wo.b; // dann komplex konj.
	  }
    // Rest Grad==1
    if (! poly.poly[1].istNull())
      poly.poly[0].div(poly.poly[1]);
    else
    	throw new Exception("Grad passt nicht !");
    poly.poly[0].a = -poly.poly[0].a;
    poly.poly[0].b = -poly.poly[0].b;
    rueck[rueck.length-1] = poly.poly[0];
	  System.out.println("letzte Ns ist " + rueck[rueck.length-1]);

	  sort(rueck); // sortieren wg. komplex konjugiert

    // for (int i = 0; i < rueck.length; ++i)
    // 	System.out.println(rueck[i]);
	  return rueck;
	}

	// Sortieren: erst komplex konjugierte, dann komplexe, dann rin
	// reelle. Gibt zurück, ab welcher Stelle es nur
	// reelle Werte gibt
	private static void sort(Komplex [] feld) throws Exception
	{
		int einfEinz = 0; boolean anders = false; // soll: i >= einfEinz !
		for (int i = 0; i < feld.length; ++i)
		{
			if (Math.abs(feld[i].b) < Komplex.EPS2 && i >= einfEinz) {} // reell
			else
			{
				// Suche, ob es komplex konj. dazu gibt
				int j;
				for (j = i+1; j < feld.length; ++j)
				{
					if (feld[i].istKonjugiert(feld[j])) break;
				}
				if (j < feld.length)
				{
					// wenn die NS aufeinander folgend sind, muss ggf. nichts gemacht werden
					if ((j == i + 1) && (i == einfEinz)) {}
					else
					{
					  // Elemente i,j zu einfEinz
					  Komplex hilfe = feld[i];
					  for(int k = i; k > einfEinz; --k) feld[k] = feld[k-1];
					  feld[einfEinz] = hilfe; 
					  hilfe = feld[j];
					  for(int k = j; k > einfEinz+1; --k) feld[k] = feld[k-1];
					  feld[einfEinz+1] = hilfe; anders = true;
					}  
					einfEinz += 2; ++i; 
				}
				else 
				{
          System.out.println("Warnung: kein konj. zu " + feld[i]);
          if (i < einfEinz) // ohne komplex konj. !
          {
          	throw new Exception("Unklarer Fall beim Einsortieren");
          }  
				}
			}
		}
		if (anders)
		{
      //System.out.println("Sortiertes Feld: ");
			//for (int i = 0; i < feld.length; ++i)
      //	System.out.println(feld[i]);
		}
	}

	// Ermittelt, ob es mit den uebergebenen komplexen
	// Nullstellen ein Teilerpolynom vom Grad grad gibt, sonst null. Es werden Nullstellen durchkombiniert,
	// die Nullstellen sind sortiert, ein paar komplex konjugierter steht hintereinander, das wird benutzt.
	public PTypZ faktorisieren (Komplex [] ns, int grad,
			PTypZ poly) throws Exception
	{
      if (grad < 1 || grad == poly.grad) return poly;
      if (grad > poly.grad) throw new Exception("Falscher Grad");

      // Feld mit Indizes, welche Lösungen benützt werden; komplex konjugierte stehen paarweise am Anfang
      int feld [] = new int[grad];
      for (int i = 0; i < feld.length; ++i) feld[i] = i;

      // ermitteln, bis wohin die Nullstellen komplex konjugiert sind (dann nur paarweise verwenden)
      int abEinzeln = 0;
      for (int j = 0; j < ns.length; j += 2)
      {
      	if ((j+1 == ns.length) || !ns[j].istKonjugiert(ns[j+1])) break;
      	else abEinzeln += 2;
      }
      if (DEBUGKONFIG) System.out.println("Abeinzeln: " + abEinzeln);
      
      // ein Teilerpolynom muss wieder ganzzahlige Koeff.
      // haben, ohne Imaginaerteil
      boolean fertig;
      do
      {
        // bei ungeradem Grad, falls letztes auf einem einzelnen eines komp. konj. Paares steht
        if ((grad & 1) == 1 && feld[grad-1] < abEinzeln)
        {
      		if (abEinzeln == ns.length) { System.out.println("Bei Grad " + grad + ": keine singuläre Ns"); return null; }
      		feld[grad-1] = abEinzeln;
        }
      	
      	// multiplizieren(Runde 1) und addieren(Runde 0) muss auf Ganzzahl führen; Addition ist schneller.
        if (DEBUGKONFIG)
        {
          System.out.print("Konfig: ");
          for (int i = 0; i < grad; ++i) System.out.print(feld[i]+"|");
          System.out.println("");
        }  

        boolean verletzt = false;
        for (short runde = 0; runde < 2 && !verletzt; ++runde)
        {
          Komplex ak = ns[feld[0]].clone();
          for (int i = 1; i < feld.length && !verletzt; ++i)
    	    {
    	      if (feld[i] >= ns.length) throw new Exception(feld[i] + ": falsch gestellt");
        	  if (runde == 1) ak.mult(ns[feld[i]]);
    	      else ak.add(ns[feld[i]]);
    	    }
          double realGer = Math.round(ak.a);
  	      // double imagGer = Math.round(ak.b);
          if (Math.abs(ak.a - realGer) < Komplex.EPS2 &&
              Math.abs(ak.b/* - imagGer*/) < Komplex.EPS2)
    	    {
   		      // System.out.println(runde == 0 ? ("Kandidat Mult. " + ak.a + " und " + ak.b) : ("Kandidat Add. " + ak.a + " und " + ak.b)); // passt !
    	    }
          else verletzt = true;
        }

        // wenn Multiplikation und Addition zielführend, dann Polynom zusammensetzen, bei Erfolg fertig
        if (!verletzt)
        {
          PTypZ erg = zusammensetzen(ns, feld);
          if (erg != null) return erg;
        }

    	  // weiterstellen, verschiedene Positionen; komplexe Werte werden möglichst zusammen genutzt wg. Laufzeit
        // 5 Nullstellen, Grad 3:
        // (0,1,4), (2,3,4)
    	  int wo = grad - 1; fertig = false;
    	  while (true /* bis gültige Konf. oder kein Weiterstellen mehr */)
    	  {
      	  if (wo < 0) { fertig = true; break; }
    	  	int index = feld[wo]; boolean mod = false;
    	    if (index >= abEinzeln)
    	    {
    	    	++feld[wo]; mod = true;
    	    }
    	    else if (index < abEinzeln && ((index & 1) == 0))
    	    {
    	    	feld[wo] += 2; mod = true;
    	    }
    	    else if (index < abEinzeln && ((index & 1) == 1))
    	    {
    	    	--wo; feld[wo] += 2; mod = true; 
    	    }
    	    // Wenn modifziert und eines der Indizes überläuft, dann weiter zurück
    	    if (mod)
    	    {
    	    	int j;
    	    	for (j = wo; j < grad; ++j)
    	    	{
    	    		if (j > wo) feld[j] = feld[j-1]+1; // nächstes der vorherg. 
    	    		if (feld[j] >= ns.length) { --wo; break; } // zuviel
    	    	}
    	    	if (j == grad) break;  // gültige Konf.
    	    }
    	    else { --wo; }
    	  } // ende while weiterstellen
    	}
      while (!fertig); // ende teilerpolynome komb.

	  return null;
	}

	// Setzt aus Nullstellen ein Polynom zusammen,
	// und gibt es zurück, wenn es aus Z[x] ist
	static PTypZ zusammensetzen (Komplex [] ns, int feld []) throws Exception
	{
      if (feld.length == 0 || feld.length > ns.length)
    	  throw new Exception("keine passenden Eingabedaten");

      Komplex zw1 [] = new Komplex[feld.length+1];
      Komplex zw2 [] = new Komplex[feld.length+1];
      // Setze die erste Nullstelle initial, baue dann dazu
      zw1[0] = ns[feld[0]].clone(); zw1[0].a = -zw1[0].a; zw1[0].b = -zw1[0].b;
      zw1[1] = new Komplex(0); zw1[1].a = 1.0; zw1[1].b = 0.0;

      // jeweils eine NS dazu
      for (int i = 1; i < feld.length; ++i)
      {
    	  Komplex c = ns[feld[i]].clone();
    	  c.a = -c.a; c.b = -c.b;    // (z - c)
    	  for (int j = 0; j <= i; ++j)
    	  {
    	    zw2[j+1] = zw1[j].clone();
    	  }
    	  zw2[0] = new Komplex(0);
    	  for (int j = 0; j <= i; ++j)
    	  {
    	    Komplex neu = zw1[j].clone();
    	    neu.mult(c); zw2[j].add(neu);
    	  }
        // Tausche wieder zw1 und zw2, danach zw1 aktuell
    	  {
    	    Komplex t [] = zw1; zw1 = zw2; zw2 = t;
    	  }
      }

      // Test, ob alle Koeff ganzzahlig sind !
      for (int i = 0; i < feld.length; ++i)
      {
    	  double realGer = Math.round(zw1[i].a);
    	  if (Math.abs(zw1[i].a - realGer) > Komplex.EPS2) return null;
    	  //double imagGer = Math.round(zw1[i].b);
    	  if (Math.abs(zw1[i].b/* - imagGer*/) > Komplex.EPS2) return null;
      }

      // bastle ein Polynom daraus (ist nach obigem reell) !
      long neuPoly [] = new long[zw1.length];
      for (int i = 0; i < zw1.length; ++i)
      {
        neuPoly[i] = (long) Math.round(zw1[i].a); // TODO klappt dies ?
      }
      return new PTypZ(feld.length, neuPoly);
	}
}