//================================================C
//  Newton法による複素多項式の根の計算             C
//   多項式f(z)=(z^n-1)*(z^m-c)                   C
//     c = ar^3 + ai^3・i                         C
//   収束回数に色を付けて描画(未収束は黒)          C
//   表示ボタンをクリックすると描画                C
//   多項式の次数nは3,4,5,6,7,8,9,10から選択   　  C
//              mは2,3,4,5,6,7,8,9,10から選択     C
//   ar,aiは-1.0,-0.9,...,1.0から選択             C
//   拡大縮小比(sl)は0.2,0.4,...,2.0から選択       C
//   2009年5月16日　後 保範（東京工芸大学)         C
//================================================C
import java.applet.*;
import java.awt.*;
import java.awt.event.*;

public class newton1 extends Applet
{ Color col[] = { new Color(255,0,0), new Color(130,0,0),
            new Color(0,255,0), new Color(0,130,0),
            new Color(0,0,255), new Color(0,0,130),
            new Color(200,200,0), new Color(120,0,120),
            new Color(0,160,160) };
  int n, m, leng = col.length;
  double ar, ai, sl;
  Button btn;                     //ボタンの設定.
  Choice n_v;                     //次数nのアイテム設定
  Choice m_v;                     //次数mのアイテム設定
  Choice ar_v;                    //複素係数aのアイテム設定
  Choice ai_v;
  Choice sl_v;                    //倍率slのアイテム設定
// 初期値の設定
 public void init() 
  { Label La1=new Label("次数n,m=");
    n_v = new Choice();             //次数nの選択
    m_v = new Choice();             //次数mの選択
    ar_v = new Choice();            //係数aの選択
    ai_v = new Choice();
    sl_v = new Choice();            //倍率slの選択
    btn = new Button("表示");       //表示ボタン
    for (int i=0; i<=7; i++)        //nの選定リスト
      { n_v.addItem(""+(i+3)); }
    for (int i=0; i<=8; i++)        //mの選定リスト
      { m_v.addItem(""+(i+2)); }
    Label La2=new Label("複素a=");
    for (int i=-10; i<=10; i++)     //arの選定リスト
      { ar_v.addItem(""+(i/10.0)); }
    for (int i=-10; i<=10; i++)     //aiの選定リスト
      { ai_v.addItem(""+(i/10.0)); }
    Label La3=new Label("倍率=");
    for (int i=0; i<=9; i++)        //slの選定リスト
      { sl_v.addItem(""+((i+1)/5.0)); }
    add(La1);  add(n_v);            //パネルへnを追加 
    add(m_v);                       //パネルへmを追加
    add(La2);  add(ar_v);           //パネルへaを追加
    add(ai_v);
    add(La3);  add(sl_v);           //パネルへslを追加
    add(btn);                       //パネルへボタン追加
    btn.addActionListener(new ActionAdp()); //ボタン処理
    n_v.addItemListener(new ItemAdp());     //アイテム追加
    n_v.select(2); n = 5;           //n=5(2番目)で表示
    m_v.addItemListener(new ItemAdp());     //アイテム追加
    m_v.select(3); m = 5;           //n=5(3番目)で表示
    ar_v.addItemListener(new ItemAdp());    //アイテム追加
    ar_v.select(20); ar = 1.0;      //ar=1.0(20番目)で表示
    ai_v.addItemListener(new ItemAdp());    //アイテム追加
    ai_v.select(11); ai = 0.1;      //ai=0.1(11番目)で表示
    sl_v.addItemListener(new ItemAdp());    //アイテム追加
    sl_v.select(3); sl = 0.8;       //sl=0.8(3番目)で表示
  }
// ボタンクリックの実装クラス
 class ActionAdp implements ActionListener
 { public void actionPerformed(ActionEvent e)
    { if (e.getSource() == btn)        //「表示」クリック
       { repaint(); }                  //再描画
    }
 }
// アイテム(nの選定)の実装クラス
 class ItemAdp implements ItemListener
 { public void itemStateChanged(ItemEvent e)
    { Object source = e.getSource();
      if (source == n_v)                   //nの選定
        { n = n_v.getSelectedIndex() + 3; }
      if (source == m_v)                   //mの選定
        { m = m_v.getSelectedIndex() + 2; }
      if (source == ar_v)                  //arの選定
        { ar = (ar_v.getSelectedIndex()-10)/10.0; }
      if (source == ai_v)                  //aiの選定
        { ai = (ai_v.getSelectedIndex()-10)/10.0; }
      if (source == sl_v)                  //slの選定
        { sl = (sl_v.getSelectedIndex()+1)/5.0; }
    }
 }
 public void paint(Graphics g) 
  { int xs = 500, ys = 500;               //描画領域
    double X[] = new double[2];
    double Z[] = new double[2];
    g.setColor(Color.black);              //黒で埋める
    g.fillRect(0, 0, xs, ys);
    for (int y=0; y<ys; y++)
     { for( int x=0; x<xs; x++)
        { double XR = (2*x-xs)/(xs*sl);   //Xの実部(x対応)
          double XI = -(2*y-ys)/(ys*sl);  //Xの虚部(y対応)
          for (int k=0; k<=80; k++)       //収束テスト回数
           { X[0] = XR;  X[1] = XI;
             double X0 = XR;  double X1 = XI;
             fnewt(X, Z);                 //X=X-f(X)/f'(X)
             XR = Z[0];  XI = Z[1];
             if((XR-X0)*(XR-X0)+(XI-X1)*(XI-X1) < 1.0e-8)
              { g.setColor(col[k%leng]);   //色の設定
                g.drawLine(x, y, x, y);    //点を描く
                break; }
           }
        }
     }
  }
//複素数のn乗計算
 public void CXN(int l, double x[], double z[])
  { double W0, W1;
    W0 = x[0];  W1 = x[1];
    for (int k=1; k<l; k++)
     { z[0]=x[0]*W0-x[1]*W1;
       z[1]=x[0]*W1+x[1]*W0; 
       W0 = z[0];  W1 = z[1]; }
  }
//複素数z=x-f(x)/f'(x)の計算
 public void fnewt(double x[], double z[])
  { double W[] = new double[2];
    double V[] = new double[2];
    double U[] = new double[2];
    double S[] = new double[2];
    double T[] = new double[2];
    CXN(n, x, U);
    U[0] -= 1.0;                //U=x^n - 1
    CXN(m, x, V);
    V[0] -= ar*ar*ar;           //V=x^m - (ar^3+ai^3・i)
    V[1] -= ai*ai*ai;
    cmult(U, V, W);             //W=f(x)
    CXN(m-1, x, S);  
    S[0] *= m;  S[1] *= m;      //S=mx^(m-1)
    cmult(S, U, T);             //T=mx^(m-1)*(x^n-1)
    CXN(n-1, x, S);
    S[0] *= n;  S[1] *= n;      //S=nx^(n-1)
    cmult(S, V, U);             //U=nx^(n-1)*(x^m-i)
    U[0] += T[0];               //U=f'(x)
    U[1] += T[1];
    Cdiv(W, U, z);              //z=f(x)/f'(x)
    z[0] = x[0] - z[0];         //z=x-z
    z[1] = x[1] - z[1]; 
  }
//複素数z=x*yの計算
 public void cmult(double x[], double y[], double z[])
  { z[0] = x[0]*y[0] - x[1]*y[1];
    z[1] = x[0]*y[1] + x[1]*y[0]; }
//複素数z=x/yの計算
 public void Cdiv(double x[], double y[], double z[])
  { double d;
    d = y[0]*y[0] + y[1]*y[1];
    if(d < 1.0e-30) { d = 1.0-30; }
    z[0] = (x[0]*y[0] + x[1]*y[1])/d;
    z[1] = (x[1]*y[0] - x[0]*y[1])/d;
  }
}