import java.applet.*;
import java.awt.*;

public class Example2_12 extends Applet
{
public static void Mrcheng(double[][] a,double[][] b,double[][] c,int m,int l,int n)
{double[][] d=new double[m][n];
 int i,j,k;
 for(i=0;i<m;i++)
   for(j=0;j<n;j++)
    {d[i][j]=0;
     for(k=0;k<l;k++)d[i][j]+=a[i][k]*b[k][j];}
 for(i=0;i<m;i++)
   for(j=0;j<n;j++)c[i][j]=d[i][j];
}

public static void  svd(double[][] a,double[][] U,double[][] V, int m,int n )
{int i,j,k,l,loop=0;
 double c,s,t,u,v,w,x,y;
 double[] b=new double[m]; double[] d=new double[m];
 boolean flag= true;
 for(i=0;i<m;i++)
    {for(j=0;j<m;j++)U[i][j]=0;
     U[i][i]=1;}
 for(i=0;i<n;i++)
	{for(j=0;j<n;j++)V[i][j]=0;
     V[i][i]=1;}
  //化為雙對角部分
  for(k=0;k<m-1;k++)
  { //下三角部分
	c=0;
    for(l=k+1;l<m;l++)c+=a[l][k]*a[l][k];
    if(c!=0)
    {c=Math.sqrt(c+a[k][k]*a[k][k]);
    if(a[k][k]<0)c=-c;
    s=(c+a[k][k])*c;
    for(i=0;i<m;i++)
     {b[i]=U[i][k]*c;
      for(l=k;l<m;l++)b[i]+=U[i][l]*a[l][k];}
    for(i=0;i<m;i++)
    {U[i][k]-=b[i]/c;
	 for(j=k+1;j<m;j++)U[i][j]-=b[i]*a[j][k]/s;}
     for(j=0;j<n;j++)
     {b[j]=c*a[k][j];
      for(l=k;l<m;l++)b[j]+=a[l][k]*a[l][j];}
    for(j=k+1;j<n;j++)
     {a[k][j]-=b[j]/c;
      for(i=k+1;i<m;i++)a[i][j]-=a[i][k]*b[j]/s;}
    for(l=k+1;l<m;l++)a[l][k]=0;
   a[k][k]=-c;}
   //上三角次對角線外部分
   if(k>=n-2)continue;
   c=0;
   for(l=k+2;l<n;l++)c+=a[k][l]*a[k][l];
   if(c!=0)
   {c=Math.sqrt(c+a[k][k+1]*a[k][k+1]);
    if(a[k][k+1]<0)c=-c;
    s=(c+a[k][k+1])*c;
    for(j=0;j<n;j++)
    {b[j]=V[k+1][j]*c;
     for(l=k+1;l<n;l++)b[j]+=V[l][j]*a[k][l];}
    for(j=0;j<n;j++)
    {V[k+1][j]-=b[j]/c;
	 for(i=k+2;i<n;i++)V[i][j]-=a[k][i]*b[j]/s;}
    for(i=0;i<m;i++)
    {b[i]=a[i][k+1]*c;
     for(l=k+1;l<n;l++)b[i]+=a[i][l]*a[k][l];}
     for(i=k+1;i<m;i++)
     {a[i][k+1]-=b[i]/c;
     for(j=k+2;j<n;j++)a[i][j]-=b[i]*a[k][j]/s;}
     for(j=k+2;j<n;j++)a[k][j]=0;
    a[k][k+1]=-c;}
   }
  //對角化部分
   while(loop++<100)
   {for(k=0;k<n-1;k++)
    {u=Math.abs(a[k][k+1]);v=Math.abs(a[k][k])+Math.abs(a[k+1][k+1]);
     if(u+v==v)a[k][k+1]=0;
     else break;}
     if(k==n-1)break;
      t=a[k][k+1]/a[k][k];c=1/Math.sqrt(1+t*t);s=c*t;
      u=a[k][k]*c+a[k][k+1]*s;v=s*a[k+1][k+1];w=-c*a[k+1][k+1];
      a[k][k]=u;a[k][k+1]=0;a[k+1][k]=v;a[k+1][k+1]=w;
      for(j=0;j<n;j++){b[j]=c*V[k][j]+s*V[k+1][j];d[j]=s*V[k][j]-c*V[k+1][j];}
      for(j=0;j<n;j++){V[k][j]=b[j];V[k+1][j]=d[j];}
     for(i=k;i<n-2;i++)
     { if(a[i+1][i]!=0)
      {t=a[i+1][i]/a[i][i];c=1/Math.sqrt(1+t*t);s=c*t;
       u=c*a[i][i]+s*a[i+1][i];v=c*a[i][i+1]+s*a[i+1][i+1];w=s*a[i+1][i+2];
       x=s*a[i][i+1]-c*a[i+1][i+1];y=-c*a[i+1][i+2];
       a[i][i]=u;a[i][i+1]=v;a[i][i+2]=w;a[i+1][i]=0;a[i+1][i+1]=x;a[i+1][i+2]=y;
       for(j=0;j<m;j++){b[j]=U[j][i]*c+U[j][i+1]*s;d[j]=U[j][i]*s-U[j][i+1]*c;}
       for(j=0;j<m;j++){U[j][i]=b[j];U[j][i+1]=d[j];}}
      if(a[i][i+2]!=0)
     {t=a[i][i+2]/a[i][i+1];c=1/Math.sqrt(1+t*t);s=c*t;
      u=a[i][i+1]*c+a[i][i+2]*s;v=a[i+1][i+1]*c+a[i+1][i+2]*s;w=s*a[i+2][i+2];
      x=a[i+1][i+1]*s-c*a[i+1][i+2];y=-c*a[i+2][i+2];
      a[i][i+1]=u;a[i][i+2]=0;a[i+1][i+1]=v;a[i+1][i+2]=x;a[i+2][i+1]=w;a[i+2][i+2]=y;
      for(j=0;j<n;j++){b[j]=V[i+1][j]*c+s*V[i+2][j];d[j]=V[i+1][j]*s-c*V[i+2][j];}
      for(j=0;j<n;j++){V[i+1][j]=b[j];V[i+2][j]=d[j];}}
    }
   if(a[n-1][n-2]!=0)
   {t=a[n-1][n-2]/a[n-2][n-2];c=1/Math.sqrt(1+t*t);s=c*t;
    u=c*a[n-2][n-2]+s*a[n-1][n-2];v=c*a[n-2][n-1]+s*a[n-1][n-1];w=s*a[n-2][n-1]-c*a[n-1][n-1];
    a[n-2][n-2]=u;a[n-2][n-1]=v;a[n-1][n-2]=0;a[n-1][n-1]=w;
    for(j=0;j<m;j++){b[j]=U[j][n-2]*c+U[j][n-1]*s;d[j]=U[j][n-2]*s-c*U[j][n-1];}
    for(j=0;j<m;j++){U[j][n-2]=b[j];U[j][n-1]=d[j];}}
  }

   for(i=0;i<n;i++)
   {if(a[i][i]<0)
    {a[i][i]=-a[i][i];
     for(j=0;j<n;j++)V[i][j]=-V[i][j];}
   }
   for(i=0;i<n-1;i++)
   {k=i;
    for(j=i;j<n;j++){if(a[k][k]<a[j][j])k=j;}
    if(k!=i)
    {t=a[k][k];a[k][k]=a[i][i];a[i][i]=t;
     for(j=0;j<n;j++){t=V[k][j];V[k][j]=V[i][j];V[i][j]=t;}
     for(j=0;j<m;j++){t=U[j][k];U[j][k]=U[j][i];U[j][i]=t;}}
   }
}

 public void MInv(double[][] a,double[][] aa,int m,int n)
 {double[][] U=new double[m][m];
  double[][] V=new double[n][n];
  double[][] bb=new double[n][n];
  int i,j,k,p;double t;
  svd(a,U,V,m,n);
  for(k=n-1;k>=0;k--){if(a[k][k]!=0)break;}
  p=k+1;
  for(k=0;k<p;k++){bb[k][k]=1/a[k][k];}
  for(i=0;i<m;i++)
   for(j=i+1;j<m;j++)
   {t=U[i][j];U[i][j]=U[j][i];U[j][i]=t;}
  for(i=0;i<n;i++)
      for(j=i+1;j<n;j++)
   {t=V[i][j];V[i][j]=V[j][i];V[j][i]=t;}
  Mrcheng(bb,U,aa,p,p,m);Mrcheng(V,aa,aa,n,p,m);
 }

 public  void paint(Graphics g)
  {double[][] a={{1,2,3,4},{6,7,8,9},{1,2,13,0},{16,17,8,9},{2,4,3,4}};
   double[][] aa=new double[4][5];
   int i,j;
   g.drawString("原矩陣A=",10,10);
   for(i=0;i<5;i++)
   g.drawString(""+a[i][0]+"  "+a[i][1]+"  "+a[i][2]+"  "+a[i][3],10,20+10*i);
   MInv(a,aa,5,4);
   g.drawString("其廣義逆A+=",10,70);
   for(i=0;i<4;i++)
   g.drawString(""+aa[i][0]+"  "+aa[i][1]+"  "+aa[i][2]+"  "+aa[i][3]+"  "+aa[i][4],10,80+10*i);

  }
}