/**
 * Diese Klasse enthaelt eine Reihe nuetzlicher mathematischer Methoden, die
 * von Sun nicht mitgeliefert werden.
 *
 * @author Olaf Müller
 * @version 1.0 25.08.2003
 */
public class MiscMath {
  public static final float EPS = 0.001f;

  /**
   * Liefert den cosinus hyperbolicus des uebergebenen Wertes value.
   */
  public static double cosh(double value) {
    return 0.5*(Math.exp(value) + Math.exp(-value));
  }

  /**
   * Liefert den sinus hyperbolicus des uebergebenen Wertes value.
   */
  public static double sinh(double value) {
    return 0.5*(Math.exp(value) - Math.exp(-value));
  }

  /**
   * Liefert den area cosinus hyperbolicus des uebergebenen Wertes value.
   * Value muss >= 1 sein!
   * Der area cosinus hyperbolicus ist die Umkehrfunktion des cosinus 
   * hyperbolicus. Dieser hat zwei Aeste. Geliefert wird der Wert, der sich im
   * Monotonieintervall [0, infinity) der Funktion value = cosh x ergibt.
   * Die Methode benutzt: acosh(value) = ln(value+sqrt(value*value-1)).
   *
   * Den Wert des anderen Astes erhaelt man durch Umkehrung des Vorzeichens. 
   * Dieser Wert entspricht dann demjenigen, der sich im Monotonieintervall
   * (-infinity, 0] von value = cosh x ergibt. 
   * Er wird mit acosh(value) = ln(value - sqrt(value*value -1)) errechnet.
   */
  public static double acosh(double value) {
    if(value < 1) {
      throw new RuntimeException("MiscMath.acosh: value ="+value+". Function is not defined for values < 1!");
    }
    return Math.log(value + Math.sqrt(value*value - 1));
  }

  /**
   * Liefert den area sinus hyperbolicus des uebergebenen Wertes value.
   */
  public static double asinh(double value) {
    return Math.log(value + Math.sqrt(value*value + 1));
  }


  /**
   * Loest die durch das Array coeff gegebene Gleichung 4. Grades
   *
   * 0 = coeff[0] + coeff[1]*x + coeff[2]*x^2 + coeff[3]*x^3 + coeff[4]*x^4
   * 
   * und liefert
   * im den Arrays re und im die Loesungen getrennt nach Imaginaer- und 
   * Realteil zurueck. Die Arrays re und im muessen je 4 Elementen Platz 
   * bieten.
   */
  public static void solveQuartic(double[] coeff, double[] re, double[] im) {
    double a,b,c,d,e,rl,alpha,beta,probe;
    double h,h1,h2;

    if(coeff.length < 5) {
      System.out.println("MiscMath.solveQuartic: array coeff must have at least 5 elements!");
      return;
    }
    if(re.length < 4) {
      System.out.println("MiscMath.solveQuartic: array re must have at least 4 elements!");
      return;
    }
    if(im.length < 4) {
      System.out.println("MiscMath.solveQuartic: array im must have at least 4 elements!");
      return;
    }
    if(coeff[4] == 0.0) {
      System.out.println("MiscMath.solveQuartic: quartic coefficient == 0 -> calling solveCubic()!");
      solveCubic(coeff, re, im);
      re[3] = im[3] = 0;
      return;
    }
    // normieren
    a = coeff[3]/coeff[4]; 
    b = coeff[2]/coeff[4]; 
    c = coeff[1]/coeff[4]; 
    d = coeff[0]/coeff[4];
    //System.out.println("MiscMath.solveQuartic: a="+a+"; b="+b+"; c="+c+"; d="+d);

    // reduzierte Gleichung der kubischen Resolvente bestimmen
    h1 = 0.25*a*c - b*b/12.0 - d;
    h2 = a*b*c/24.0 - 0.125*c*c - b*b*b/108.0 - 0.125*a*a*d + b*d/3.0;
    //System.out.println("MiscMath.solveQuartic: h1="+h1+"; h2="+h2);
  
    solveCubic(new double[] {h2, h1, 0, 1}, re, im);
    //System.out.println("MiscMath.solveQuartic: Loesungen der kubischen Resolvente:");
    //System.out.println("MiscMath.solveQuartic: "+re[0]+"+i*"+im[0]);
    //System.out.println("MiscMath.solveQuartic: "+re[1]+"+i*"+im[1]);
    //System.out.println("MiscMath.solveQuartic: "+re[2]+"+i*"+im[2]);

    // alert("Lösgen der Resolv= "+lr_1+" "+lk_1+" / "+lr_2+" "+lk_2+" / "+lr_3+" "+lk_3)
  
    rl = re[0];
    h1 = 0.25*a*a - 2.0*b/3.0 + 2.0*rl + 1e-15;
    h2 = b*b/36.0 + rl*rl + b*rl/3.0 - d + 1e-15;
    if (Math.abs(h1) < 1e-10) { 
      h1 = 0.0;
    }
    if (Math.abs(h2) < 1e-10) { 
      h2 = 0.0;
    }
    alpha = Math.sqrt(h1);
    beta = Math.sqrt(h2);
    //System.out.println("MiscMath.solveQuartic: alpha="+alpha+"; beta="+beta);

    probe = 2.0*alpha*beta - (a*b/6.0 + a*rl - c);
    if (Math.abs(probe) > 0.000001) {
      alpha = -alpha;
      probe = 2.0*alpha*beta - (a*b/6.0 + a*rl - c);
    }
    if (probe > 0.000001) { 
      System.out.println("MiscMath.solveQuartic: Probe nicht gelungen !");
    } 
    h1 = 0.5*a + alpha;
    h2 = b/6.0 + rl + beta;
    //System.out.println("MiscMath.solveQuartic: h1="+h1+"; h2="+h2);
    solveQuadric(new double[] {h2, h1, 1}, re, im);
    //System.out.println("MiscMath.solveQuartic: Loesungen der quadratischen NR:");
    //System.out.println("MiscMath.solveQuartic: "+re[0]+"+i*"+im[0]);
    //System.out.println("MiscMath.solveQuartic: "+re[1]+"+i*"+im[1]);
    // Loesungen an einen sicheren Ort kopieren
    re[2] = re[0];
    im[2] = im[0];
    re[3] = re[1];
    im[3] = im[1];

    h1 = 0.5*a - alpha;
    h2 = b/6.0 + rl - beta;
    //System.out.println("MiscMath.solveQuartic: h1="+h1+"; h2="+h2);
    solveQuadric(new double[] {h2, h1, 1}, re, im);
  }



  /**
   * Loest die durch das Array coeff gegebene kubische Gleichung
   *
   * 0 = coeff[0] + coeff[1]*x + coeff[2]*x^2 + coeff[3]*x^3
   * 
   * und liefert
   * im den Arrays re und im die Loesungen getrennt nach Imaginaer- und 
   * Realteil zurueck. Die Arrays re und im muessen je 3 Elementen Platz 
   * bieten.
   */
  public static void solveCubic(double[] coeff, double[] re, double[] im) {
    double r,s,t,p,q,u,v;
    double det,help,phi;
    double y1,y2,y3;

    if(coeff.length < 4) {
      System.out.println("MiscMath.solveCubic: array coeff must have at least 4 elements!");
      return;
    }
    if(re.length < 3) {
      System.out.println("MiscMath.solveCubic: array re must have at least 3 elements!");
      return;
    }
    if(im.length < 3) {
      System.out.println("MiscMath.solveCubic: array im must have at least 3 elements!");
      return;
    }
    if(coeff[3] == 0.0) {
      System.out.println("MiscMath.solveCubic: cubic coefficient == 0 -> calling solveQuadric()!");
      solveQuadric(coeff, re, im);
      re[2] = im[2] = 0;
      return;
    }

    // normieren
    r=coeff[2]/coeff[3];
    //System.out.println("MiscMath.solveCubic: r="+r);
    s=coeff[1]/coeff[3]; 
    //System.out.println("MiscMath.solveCubic: s="+s);
    t=coeff[0]/coeff[3];
    //System.out.println("MiscMath.solveCubic: t="+t);

    r /= 3.0;
    // Hilfsgroessen fuer reduzierte Gleichung
    p=s-3.0*r*r;
    //System.out.println("MiscMath.solveCubic: p="+p);
    q=2.0*r*r*r-r*s+t;
    //System.out.println("MiscMath.solveCubic: q="+q);

    // Entscheidungsvariable fuer Loesungsverhalten
    det=q*q/4+p*p*p/27;
    //System.out.println("MiscMath.solveCubic: det="+det);

    // Falls 3 relle Loesungen existieren
    // Der Fall det==0 koennte auch hier abgehandelt werden, aber die Methode
    // Math.pow liefert dann NaN, was falsch ist. Also wird dieser Fall im
    // else-Zweig bearbeitet
    if (det<0) {
      //System.out.println("MiscMath.solveCubic: det<=0 -> 3 real solutions!");
      double R = 2.0*Math.sqrt(Math.abs(p/3.0));
      //help=-0.5*q/Math.pow(R,3.0);
      help=-4.0*q/Math.pow(R,3.0);
      //System.out.println("MiscMath.solveCubic: help="+help);
      phi=Math.acos(help);
      //System.out.println("MiscMath.solveCubic: phi="+phi);
      y1= R*Math.cos(phi/3.0);
      //System.out.println("MiscMath.solveCubic: y1="+y1);
      y2=-R*Math.cos((phi-Math.PI)/3.0);
      //System.out.println("MiscMath.solveCubic: y2="+y2);
      y3=-R*Math.cos((phi+Math.PI)/3.0);
      //System.out.println("MiscMath.solveCubic: y3="+y3);
    
      // Uebergang von reduzierter zu Originalgleichung
      re[0]=y1-r;
      re[1]=y2-r;
      re[2]=y3-r;
      im[0] = 0;
      im[1] = 0;
      im[2] = 0;
    }
    else {
    // Ggf. existieren komplexe Loesungen
      help = -0.5*q+Math.sqrt(det);
      if (help<0) { 
        u = -Math.pow(-help, 1.0/3.0);
      }
      else {
        u = Math.pow(help, 1.0/3.0);
      }

      help = -0.5*q-Math.sqrt(det);
      if (help<0) {
        v = -Math.pow(-help, 1.0/3.0);
      }
      else {
        v = Math.pow(help, 1.0/3.0);
      }

      y1 = u+v;
      y2 = -y1/2.0;
      y3 = (-u+v)*Math.sqrt(0.75);

      // Uebergang von reduzierter zu Originalgleichung
      re[0] = y1-r;
      re[1] = y2-r;
      re[2] = re[1];
      im[0] = 0;
      im[1] = y3;
      im[2] = -im[1];
    }
  }


  public static void solveQuadric(double[] coeff, double[] re, double[] im) {
    if(coeff.length < 3) {
      System.out.println("MiscMath.solveQuadric: array coeff must have at least 3 elements!");
      return;
    }
    if(re.length < 2) {
      System.out.println("MiscMath.solveQuadric: array re must have at least 2 elements!");
      return;
    }
    if(im.length < 2) {
      System.out.println("MiscMath.solveQuadric: array im must have at least 2 elements!");
      return;
    }
    if(coeff[2] == 0.0) {
      System.out.println("MiscMath.solveQuadric: quadric coefficient == 0 -> calling solveLinear()");
      solveLinear(coeff, re, im);
      re[1] = im[1] = 0;
      return;
    }

    double p,q,det;
    p=coeff[1]/coeff[2];
    q=coeff[0]/coeff[2];
    det=p*p/4-q;
    if (det >= 0) {
      // 2 reelle Loesungen
      re[0] = -p/2+Math.sqrt(det);
      re[1] = -p/2-Math.sqrt(det);
      im[0] = 0;
      im[1] = 0;
    }
    else {
      // 2 konjugiert komplexe Loesungen
      re[0] = -p/2;
      re[1] = re[0];
      im[0] = Math.sqrt(-det);
      im[1] = -im[0];
    }
  }

  public static void solveLinear(double[] coeff, double[] re, double[] im) {
    if(coeff.length < 2) {
      System.out.println("MiscMath.solveLinear: array coeff must have at least 2 elements!");
      return;
    }
    if(re.length < 1) {
      System.out.println("MiscMath.solveLinear: array re must have at least 1 element!");
      return;
    }
    if(im.length < 1) {
      System.out.println("MiscMath.solveLinear: array im must have at least 1 element!");
      return;
    }

    if(coeff[1] == 0) {
      // Konstante Gleichung -> Loesung steht fest
      re[0] = im[0] = 0;
    }
    else {
      // Lineare Gleichung -> Loesung berechnen
      re[0] = -coeff[0]/coeff[1];
      im[0] = 0;
    }
  }

  public static void main(String[] argv) {
    System.out.println("MiscMath.main: Testing Methods!");
    System.out.println("\nsolveCubic: Using coefficients:");
    double[] c = {-36, 60, -25, 0, -5.3};
    //double[] c = {-6, 11, -6, 1};
    System.out.println("0="+c[4]+"x^4+"+c[3]+"x^3+"+c[2]+"x^2+"+c[1]+"x+"+c[0]);
    double[] r = new double[4];
    double[] i = new double[4];

    long start = System.currentTimeMillis();
    int j=0;
    //for(; j<1000000; j++) {
    for(; j<1; j++) {
      solveQuartic(c, r, i);
    }
    System.out.println("Solving "+j+" times costs: "+(System.currentTimeMillis()-start)+"ms");

    System.out.println("Solutions:");
    System.out.println(r[0]+"+i*"+i[0]);
    System.out.println(r[1]+"+i*"+i[1]);
    System.out.println(r[2]+"+i*"+i[2]);
    System.out.println(r[3]+"+i*"+i[3]);

  }
}