// LagrangeFit.java
//
//   L_n,j(x) = product i=0,n i!=j (x-x_i)/(x_j-x_i)  
//   let 0 < x < 1 for points x_0, x_1, ... , x_n
//   let h = x_i - x_i-1,  x_i = i * h,  h = 1/n
//   L_n(x) = sum j=0,n  f(x_j) L_n,j(x)
//   This L_n(x) fits f(x) with a n th order polynomial
 
import java.awt.*;
import java.awt.event.*;
import java.text.*;

class LagrangeFit extends Frame
{
  
  LagrangeFit()
  {
    setTitle("LagrangeFit(x) 0 to 2 of e^x");
    setSize(450,500);
    setBackground(Color.white);
    setForeground(Color.black);
    addWindowListener(new WindowAdapter()
    {
      public void windowClosing(WindowEvent e)
      {
        System.exit(0);
      }
    });
    setVisible(true);
  }

  double f(double x) // "exact" function e^x 
  {
    return Math.exp(x);
  }

  void L(int n, int j, double p[]) // Lagrange polynomial generation
  {
    double termx, termc, h;
    double q[] = new double[20]; // working polynomial coefficients
    
    System.out.println("computing terms of L_"+n+","+j);
    h = 1.0/(double)n;
    p[0] = 1.0; // identity polynomial to start product
    for(int i=1; i<=n; i++) p[i] = 0.0;
    for(int i=0; i<=n; i++)
    {
      if(i!=j)
      {
        termx = 1.0/((j-i)*h); // (x    )/(x_j-x_i)
        termc = -h*i*termx;    // ( -x_i)/(x_j-x_i)
        // do product p[] * (termx * x + termc)
        for(int k=0; k<=n; k++) q[k] = p[k] * termc;
        for(int k=1; k<=n; k++) q[k] = q[k] + p[k-1]*termx;
        for(int k=0; k<=n; k++) p[k] = q[k];
      }
    }
    // print L_n,j(x) = p[...]
    System.out.println("L_"+n+","+j+"(x) = "+p[0]+" x^0 + ");
    for(int i=1; i<=n; i++)
      System.out.println("           "+p[i]+" x^"+i+" + ");
    System.out.println(" ");
  }

  public void paint(Graphics g)
  {
    // set up function to be approximated
    final int n=100;            // number of samples for plotting
    final double nf=(double)n;
    double fx[] = new double[n]; // f(x)
    double ax[] = new double[n]; // approximation(x)
    final double xmin=0.0;
    final double xmax=2.0;
    double ymin=1.0;            // for automatic scaling some day
    double ymax=7.0;
    double dx;                  // x step increment
    double xstart;              // initial value
    double x;                   // stepped value xmin to xmax
    
    xstart=xmin;
    dx = (xmax-xmin)/(double)n;
    // do NOT use x=x+dx, very inaccurate!
    for(int i=0; i<n; i++) fx[i]=f(xstart+(double)i*dx); // get "exact" function
    ymin=fx[0];
    ymax=fx[0];
    for(int i=0; i<n; i++)
    {
      ymin=Math.min(ymin,fx[i]);
      ymax=Math.max(ymax,fx[i]);
    }
    System.out.println("computed ymin="+ymin+", ymax="+ymax);
    // may fudge ymin and ymax for extra room for approximation
    ymin=1.0;
    ymax=7.0;
 
    // set up plot area
    final int hw=400;   // graph area height and width
    final double xsca=hw/(xmax-xmin);  // x data scale 0 to 2
    final double ysca=hw/(ymax-ymin);  // y data scale 1 to 7
    final int f1xOff=25;  // offsets
    final int f1yOff=70;
    final int f1xC=f1xOff+hw/2; // centers
    final int f1yC=f1yOff+hw/2;
    
    g.drawRect(f1xOff, f1yOff, hw, hw);
    g.drawLine(f1xOff, f1yOff+hw/3, f1xOff+hw, f1yOff+hw/3);
    g.drawLine(f1xOff, f1yOff+2*hw/3, f1xOff+hw, f1yOff+2*hw/3);
    g.drawLine(f1xC,   f1yOff, f1xC, f1yOff+hw);
    g.drawString("7",  f1xOff-14,   f1yOff+4);
    g.drawString("5",  f1xOff-14,   f1yOff+hw/3+4);
    g.drawString("3",  f1xOff-14,   f1yOff+2*hw/3+4);
    g.drawString("e",  f1xOff-14,   f1yOff+2*hw/3+24);
    g.drawString("1",  f1xOff-14,   f1yOff+hw+4);
    g.drawString("0",  f1xOff-4,    f1yOff+hw+15);
    g.drawString("1",  f1xC-4,      f1yOff+hw+15);
    g.drawString("2",  f1xOff+hw-4, f1yOff+hw+15);
    g.drawString("y",  f1xOff,      f1yOff-12);
    g.drawString("x",  f1xOff+hw+6, f1yOff+hw+4);

    // plot function
    double x1new, y1new;
    double x1old=xstart;
    double y1old=fx[0];
    g.setColor(Color.black);
    for(int i=1; i<n; i++) // n-1 lines
    {
      x1new=xstart+(double)i*dx;
      y1new=fx[i];
      g.drawLine((int)(f1xOff+x1old*xsca), (int)(f1yOff+hw-(y1old-ymin)*ysca),
                 (int)(f1xOff+x1new*xsca), (int)(f1yOff+hw-(y1new-ymin)*ysca));
      x1old=x1new;
      y1old=y1new;
    }
    
    // compute approximation
 
    Color col[]={Color.yellow, Color.green, Color.red, Color.blue};
    int t=0;                     // power actually used
    double avgerr, maxerr, rmserr, err, y=0.0;
    double La[] = new double[6]; // the coefficients of the fit
    double p[]  = new double[6]; // coefficients of Lagrange polynomial
    DecimalFormat f1 = new DecimalFormat("0.000E00");
 
    for(int k=0; k<4; k++)
    {
      // generate the t_th approximation
      t=k+2;
      for(int i=0; i<=t; i++) La[i] = 0.0;
      for(int j=0; j<=t; j++)
      {
        L(t,j,p); // generate Lagrange polynomial coefficients  
        // perform summation (integration not needed)
        for(int i=0; i<=t; i++) La[i] = La[i] +f(j*(1.0/t))*p[i];
      }
      System.out.println("Lagrange fit coefficients for pwr="+t);
      for(int i=0; i<=t; i++) System.out.println("x^"+i+"* La["+i+"]="+La[i]);

      // print La[i]
      avgerr=0.0;
      maxerr=0.0;
      rmserr=0.0;
      for(int i=0; i<n; i++) // n samples to plot
      {
        x1new=xstart+(double)i*dx;
        // y1new is the fit
        y1new=La[t]*x1new;
        for(int j=t-1; j>0; j--)y1new=(La[j]+y1new)*x1new;
        ax[i]=y1new+La[0];
        err=Math.abs(fx[i]-ax[i]);
        avgerr=avgerr+err;
        maxerr=Math.max(maxerr, err);
        rmserr=rmserr+err*err;
        
      }
      avgerr=avgerr/nf;
      rmserr=Math.sqrt(rmserr/nf);
      System.out.println("for "+t+" power: avgerr="+avgerr);
      System.out.println("            maxerr="+maxerr);
      System.out.println("            rmserr="+rmserr);

      // plot approximation
      g.setColor(col[k]);
      x1old=xstart;
      y1old=ax[0];
      for(int i=1; i<n; i++) // n-1 lines
      {
        x1new=xstart+(double)i*dx;
        y1new=ax[i];
        g.drawLine((int)(f1xOff+x1old*xsca), (int)(f1yOff+hw-(y1old-ymin)*ysca),
                   (int)(f1xOff+x1new*xsca), (int)(f1yOff+hw-(y1new-ymin)*ysca));
        x1old=x1new;
        y1old=y1new;
      }
      g.drawString(t+" power,  avg error ="+f1.format(avgerr), f1xC+20, f1yOff+hw-10-15*k); 
    }
    g.setColor(Color.black);
    g.drawString(n+" samples", f1xOff+20, f1yOff+hw-20);
    g.setColor(Color.red);
    Font cur18 = new Font("courier", Font.BOLD, 18); 
    g.setFont(cur18);
    g.drawString("LagrangeFit(x) 0 to 2 of e^x", 10, 50);
  }
  
  public static void main(String args[])
  {
    System.out.println("LagrangeFit.java running");
    new LagrangeFit();
  }
}