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

public class newton0 extends Applet
{ Color col[] = { Color.pink, Color.blue, Color.orange,
                  Color.green, Color.lightGray, Color.red,
                  Color.yellow, Color.cyan, Color.magenta };               
  int n, leng = col.length;
  Button btn;                       //ボタンの設定.
  Choice n_v;                       //係数nのアイテム設定
// 初期値の設定
 public void init() 
  { Label La1=new Label("n=");
    n_v = new Choice();             //実数nの選択
    btn = new Button("表示");       //表示ボタン
    for (int i=0; i<=7; i++)        //nの選定リスト
      { n_v.addItem(""+(i+3)); }
    add(La1);  add(n_v);            //パネルへnを追加 
    add(btn);                       //パネルへボタン追加
    btn.addActionListener(new ActionAdp()); //ボタン処理
    n_v.addItemListener(new ItemAdp());     //アイテム追加
    n_v.select(0); n = 3;           //n=3(0番目)で表示
  }
// ボタンクリックの実装クラス
 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; }
    }
 }
 public void paint(Graphics g) 
  { int xs = 500, ys = 500;                //描画領域
    double sl = 0.5;                       //倍率
    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<=50; k++)       //収束テスト回数
           { X[0] = XR;  X[1] = XI;
             double X0 = XR;  double X1 = XI;
             fnewt(n, 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 n, double x[], double z[])
  { double W0, W1;
    W0 = x[0];  W1 = x[1];
    for (int k=1; k<n; 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-(x^n-1)/(nx^(n-1))の計算
 public void fnewt(int n, double x[], double z[])
  { double W[] = new double[2];
    double U[] = new double[2];
    CXN(n, x, W);
    W[0] -= 1.0;                //W=x^n - 1
    CXN(n-1, x, U);  
    U[0] *= n;  U[1] *= n;      //U=nx^(n-1)
    Cdiv(W, U, z);              //z=W/U
    z[0] = x[0] - z[0];         //z=x-z
    z[1] = x[1] - z[1]; 
  }
//複素数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;
  }
}