// ChebyshevFit.java
//
// approximate f(x) ~ a0 + a1*T1(x) + a2*T2(x) + ...
// ai = 2/pi * int -1 to 1 f(x)*Ti(x)/sqrt(1-x^2) dx
// a0 = a0/2


import java.awt.*;
import java.awt.event.*;
import java.text.*;

class ChebyshevFit extends Frame
{
  ChebyshevFit()
  {
    setTitle("ChebyshevFit(x) -1 to 1 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);
  }

  double T(int i, double x) // Chebyshev polynomial evaluation
  {
    if(i<=0) return 1.0;
    if(i==1) return x;
    return 2.0*x*T(i-1,x)-T(i-2,x);
  }

  public void paint(Graphics g)
  {
    final double pi=Math.PI;

    // set up function to be approximated
    final int n=100;            // number of samples
    final double nf=(double)n;
    double fx[] = new double[n]; // f(x)
    double ax[] = new double[n]; // approximation(x)
    double bx[] = new double[n]; // alternate approximation(x)
    final double xmin=-1.0;
    final double xmax=1.0;
    double ymin=0.0;            // for automatic scaling some day
    double ymax=3.0;
    double dx;                  // x step increment
    double xstart;              // initial value
    if(n%2==0) // even, center samples, not end points or center
    {
      dx=(xmax-xmin)/nf;
      xstart=xmin+dx/2.0;
    } 
    else // odd, end points are samples, center point in center
    {
      dx=(xmax-xmin)/(double)(n-1);
      xstart=xmin;
    }
    // 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]);
    }
    // may fudge ymin and ymax for extra room for approximation
    ymin=0.0;
    ymax=3.0;
 
    // set up plot area
    final int hw=400;   // graph area height and width
    final double xsca=hw/(xmax-xmin);  // x data scale -1 to 1
    final double ysca=hw/(ymax-ymin);  // y data scale 0 to 3
    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("3",  f1xOff-14,   f1yOff+4);
    g.drawString("e",  f1xOff-14,   f1yOff+400-362+4);
    g.drawString("2",  f1xOff-14,   f1yOff+hw/3+4);
    g.drawString("1",  f1xOff-14,   f1yOff+2*hw/3+4);
    g.drawString("1/e", f1xOff-24,  f1yOff+400-49+4);
    g.drawString("0",  f1xOff-14,   f1yOff+hw+4);
    g.drawString("-1", f1xOff-4,    f1yOff+hw+15);
    g.drawString("0",  f1xC-4,      f1yOff+hw+15);
    g.drawString("1",  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, x;
    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)(f1xC+x1old*xsca), (int)(f1yOff+hw-y1old*ysca),
                 (int)(f1xC+x1new*xsca), (int)(f1yOff+hw-y1new*ysca));
      x1old=x1new;
      y1old=y1new;
    }
    
    // compute approximation
    final int nt=10;
    double a[] = new double[nt]; // coefficients
    double b[] = new double[nt];
    
    // compute Chebyshev coefficients a's and b's 
    for(int i=0; i<nt; i++)
    {
      a[i]=0.0;
      b[i]=0.0;
      for(int j=0; j<n; j++) // n samples
      {
        x=xstart+(double)j*dx;
        a[i]=a[i]+fx[j]*T(i,x)/Math.sqrt(1.0-x*x);
        b[i]=b[i]+fx[j]*T(i,x);
      }
      a[i]=(2.0/pi)*(xmax-xmin)*a[i]/nf; // scale numerical integration
      b[i]=(2.0/pi)*(xmax-xmin)*b[i]/nf;
      if(i==0)a[0]=a[0]/2.0; // constant term special
      if(i==0)b[0]=b[0]/2.0;
      System.out.println("a["+i+"]="+a[i]+"   b="+b[i]);
    }
    // override with Maple values
    a[0]=1.2660658777;
    a[1]=1.1303182079;
    a[2]=0.2714953395;
    a[3]=0.0443368498;
    a[4]=0.0054742404;
    a[5]=0.0005429263;
    a[6]=0.0000449773;
    a[7]=0.0000031984;
    System.out.println("values from Maple");
    for(int i=0; i<nt; i++)System.out.println("a["+i+"]="+a[i]);
 
    Color col[]={Color.yellow, Color.green, Color.red, Color.blue};
    int t=0;                     // power actually used
    double avgerr, maxerr, rmserr, err;
    double avgerrb, maxerrb, rmserrb, errb, y=0.0;
    DecimalFormat f1 = new DecimalFormat("0.000E00");
 
    for(int k=0; k<4; k++)
    {
      // generate the k_th approximation
      t=k+3;
      avgerr=0.0;
      maxerr=0.0;
      rmserr=0.0;
      avgerrb=0.0;
      maxerrb=0.0;
      rmserrb=0.0;
      for(int i=0; i<n; i++) // n samples
      {
        x1new=xstart+(double)i*dx;
        y1new=a[0];
        for(int j=1; j<t; j++)y1new=y1new+a[j]*T(j,x1new);
        ax[i]=y1new;
        err=Math.abs(fx[i]-ax[i]);
        avgerr=avgerr+err;
        maxerr=Math.max(maxerr, err);
        rmserr=rmserr+err*err;
        
        y=b[0];
        for(int j=1; j<t; j++)y=y+b[j]*T(j,x1new);
        bx[i]=y/Math.sqrt(1.0-x1new*x1new);
        errb=Math.abs(fx[i]-bx[i]);
        avgerrb=avgerrb+errb;
        maxerrb=Math.max(maxerrb, errb);
        rmserrb=rmserrb+errb*errb;
      }
      avgerr=avgerr/nf;
      rmserr=Math.sqrt(rmserr/nf);
      System.out.println("for "+(t-1)+" power: avgerr="+avgerr);
      System.out.println("            maxerr="+maxerr);
      System.out.println("            rmserr="+rmserr);
      avgerrb=avgerrb/nf;
      rmserrb=Math.sqrt(rmserrb/nf);
      System.out.println("for "+(t-1)+" power: avgerrb="+avgerrb);
      System.out.println("            maxerrb="+maxerrb);
      System.out.println("            rmserrb="+rmserrb);

      // 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)(f1xC+x1old*xsca), (int)(f1yOff+hw-y1old*ysca),
                   (int)(f1xC+x1new*xsca), (int)(f1yOff+hw-y1new*ysca));
        x1old=x1new;
        y1old=y1new;
      }
      g.drawString((t-1)+" 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("ChebyshevFit(x) -1 to 1 of e^x", 10, 50);
  }
  
  public static void main(String args[])
  {
    new ChebyshevFit();
  }
}