コンピュータの最近のブログ記事

「コンピュータによる計算」というサブジェクトでブログ記事を書こうとしているのだが、これほど、何がいいたいの?と曖昧な疑問を呈させるサブジェクトもないので、説明する。

すなわち、最近、ポナンザ囲碁AI新時代など、将棋AIや囲碁AIが強くなって人間が勝てないという議論がよく行われているが、私はつい、そうはいうけど、もともとコンピュータは人間に比べて化け物のような計算量を誇るので、将棋AIや囲碁AIの人間の負かし方は相当手ぬるいと思ってしまう。

例えば、Lucas-Lehmer判定法であるが、Lucasはほぼ一生を費やして手計算で、2127-1がメルセンヌ素数であることを発見したが、それが手計算の限界で、次の2521-1は、それから70年以上かかってコンピュータを使うことで発見された。もし2521-1を手計算でLucas-Lehmer判定により素数であることを発見しようとしたら、もしかして1000年くらいかかるかもしれない。しかし、私が使っているパソコンでも、その計算にかかる時間はせいぜい10分くらいで、その差は2000万倍程度である。

もしかして、21000-1のコラッツ数列の計算なら、まだしもその差は4万倍程度で、それほど大差ではないかもしれない。

それではこれならどうだろう、すなわち、3体問題のルンゲクッタの数値計算を人間が手計算すると考えると、Lucas-Lehmer判定やコラッツ数列と違って、平方根の計算も含み、1ステップだけでも下手したらまる一日くらいかかりそうな面倒な浮動小数点計算を少なくとも100万ステップは計算しなくては意味がなく、それは発狂的とさえいえる。

これがナビエストークスのような流体力学の偏微分方程式の数値計算となると、発狂的の累乗で、尻尾を巻いて逃げるよりない。

私がいいたいことは、将棋AIや囲碁AIは相当に人間に近い土俵で、コンピュータが戦っているのだというただけそれだけである。

Visual Studio

user-pic
0

私は実は昔からマイクロソフトのCコンパイラを使うのが憧れで、しかしとても高価なので仕方なく、ほぼ1万円程度で買えるボーランドのTurbo C++を購入して、プログラミングを勉強した。最初に買ったTurbo C++は安いのはいいが、日本語処理に難があって、2バイトコードがきちんと処理できないので、そこを前処理するプログラムを書いたりした。

ボーランドへの慣れで、GUIベースのC++コンパイラも、ボーランドのC++ Builderを買った。当時は私は予算権限があり、仕事に使うツールなら私の一存で買えたので、最早マイクロソフトのVisual C++を買うことができたのだが、C++ Builderを使い続けた。

それで、外資系の会社を辞めたときは、自宅のWindows 7パソコンに、手持ちのC++ Builder 3.0は導入できなかったので、仕方なくボーランドの後継であるエンバカデロC++ Builder XE4を購入して使った。5月にソニーVAIOのWindows 7パソコンから、今のレノボのWindows10パソコンに移行したときも、エンバカデロC++ Builder XE4を移行した。

ここまでは順調だったのだが、Windows10パソコンでエンバカデロC++ Builder XE4が、LME288のリンクエラーを頻発して、信頼が下がってきた。それでアマゾンでいろいろ探していたら、「アプリを作ろう! Visual C++入門」WINGSプロジェクト著という本が見つかった。これによると、無料で使えるVisual Studio Community 2015でのプログラミングが学べるという。そこで「アプリを作ろう! Visual C++入門」を買ってみると、ウェブからVisual Studio Community 2015をダウンロードするのだという。

ところが、「アプリを作ろう! Visual C++入門」はたった1年前の本なのに、ウェブでVisual Studio Community 2015はみつからず、代わりにVisual Studio Community 2017が出てきた。うーん、と考えて仕方なくVisual Studio Community 2017を導入したら、「アプリを作ろう! Visual C++入門」の画面と見比べると、そんなにインターフェースは違わないようで一応安心した。

しかし、Visual Studio Community 2017は操作といい、C++の文法といい、あまりにも違うのでカルチャーショックを受けるほどだった。そもそもが、C++ Builderはエントリポイントを自動生成してくれるが、Visual Studio Community 2017のVisual C++は、自分でエントリポイントを作成しなくてはならない。それ以外にもインターフェースが違いすぎて、私はプログラミングなら、決してA級ライセンスのレベルではないにしても、都内どこでも普通にドライブできるベテラン・ドライバに相当すると思うが、Visual Studio Community 2017だと全く操作の分からない車をいきなり運転したようなもので、これでは到底公道を運転できない。さて、何か月かかけて慣らし運転するとしよう。それにしても私はなんで、若い時からプログラムのインターフェースに慣れるのが苦手なのだろうか。プログラミングそのものはまあ苦手感はないとおもうのだが。

森肇教授の「カオス 流転する自然」岩波科学ライブラリーに、レスラー系のアトラクタというものが示されていたので、これも描いてみることにした。なお、レスラー系といっても、アントニオ猪木であるわけがなく、創案者はドイツ人のレスラー(Roessler)という人らしい。"oe"は本当はウムラウトなのだが、それは往々にして文字化けするのでここでは使わない。その方程式は、Wikiも併せて参照すると、以下のとおりである。

ここで定数a = 0.2、b = 0.2、c = 5.7でルンゲ・クッタで描いた軌跡が以下のとおりである。

なんでも、cが特定の値より大きいときカオスになるらしくて、そこでどのあたりが臨界値か調べてみた。

定数a = 0.2、b = 0.2、c = 4.7でも以下のとおり定常軌道である。

定数a = 0.2、b = 0.2、c = 2.0だと以下のとおり定常軌道である。

定数a = 0.2、b = 0.2、c = 3.0だと以下のとおり定常軌道である。。

定数a = 0.2、b = 0.2、c = 4.1だと以下のとおり定常軌道である。。

定数a = 0.2、b = 0.2、c = 4.2だと以下のとおり定常軌道かもしれないが微妙である。よって、4.2より少し大きいところに臨界値があると思われる。

プログラムのソースコードは、以下のとおりである。

//---------------------------------------------------------------------------

#include <vcl.h>
#include <math.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double fnx(double y,double z); 
double fny(double a,double x,double y);
double fnz(double b,double c,double x,double z);
double ROTx(double TH1, double TH2, double x, double y, double z);
double ROTy(double TH1, double TH2, double x, double y, double z);
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
const double H = 0.00001;
const pi = 3.14156;
double x = 2.0,y = 2.0,z = 13.0;
double a = 0.2, b = 0.2, c = 4.25;
double a1,a2,a3,a4;
double b1,b2,b3,b4;
double c1,c2,c3,c4;
double px, py;
int N = 0;

 while (N++ < 50000000 )  {
  a1 = fnx(y,z); 
  b1 = fny(a,x,y);
  c1 = fnz(b,c,x,z);
  a2 = fnx(y+H*b1/2.0,z+H*c1/2.0);
  b2 = fny(a,x+H*a1/2.0,y+H*b1/2.0);
  c2 = fnz(b,c,x+H*a1/2.0,z+H*c1/2.0);
  a3 = fnx(y+H*b2/2.0,z+H*c2/2.0);
  b3 = fny(a,x+H*a2/2.0,y+H*b2/2.0);
  c3 = fnz(b,c,x+H*a2/2.0,z+H*c2/2.0);
  a4 = fnx(y+H*b3,z+H*c3);
  b4 = fny(a,x+H*a3,y+H*b3);
  c4 = fnz(b,c,x+H*a3,z+H*c3);
  x += H*(a1+2.0*a2+2.0*a3+a4)/6.0;
  y += H*(b1+2.0*b2+2.0*b3+b4)/6.0;
  z += H*(c1+2.0*c2+2.0*c3+c4)/6.0;
  if ( N > 100 )                {
   px = 400.0 + 15.0*ROTx(pi/4.0, pi/4.0,x,y,z);
   py = 400.0 + 15.0*ROTy(pi/4.0, pi/4.0,x,y,z);
   Canvas->Pixels[px][py] = clBlack;
  }
 }
}
//---------------------------------------------------------------------------
double fnx(double y,double z)  {
 return -y-z;
}

double fny(double a,double x,double y)  {
 return x + a*y;
}

double fnz(double b,double c,double x,double z)  {
 return b + x*z - c*z;
}

double ROTx(double TH1, double TH2, double x, double y, double z)  {
 return cos(TH1)*cos(TH2)*x + sin(TH1)*y + cos(TH1)*sin(TH2)*z;
}

double ROTy(double TH1, double TH2, double x, double y, double z)  {
 return -sin(TH1)*cos(TH2)*x + cos(TH1)*y - sin(TH1)*sin(TH2)*z;
}

エノンのアトラクタというかエノン写像の続きだが、Wikiによると、ロジスティック写像のように分岐した後カオスになる図が存在するようなので、描いてみることにした。下記のとおりである。

この図は、下記のようなエノン写像において、b = 0.3と固定し、aを1.0から1.42まで次第に増やしていく。
xi+1 = 1 - axi2 + yi
yi+1 = byi

すなわち、図において横軸はaであり、縦軸はxの値である。aが1.42で打ち切られているのは、それよりaを増やすと式が発散してしまうからである。

この図からは、ロジスティック写像よりも複雑な分岐の様子がわかる。なお図において、縦軸と横軸にそれぞれ0.1刻みで目盛りをつけてみた。

なお参考までに、縦軸をyにしたときの図は、こちらである。

見て取れるように、yの方も分岐はあるのだが、値の幅が相対的に小さくなっている。

最後に例によって、プログラムのソースコードを付ける。

//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
 double x,y,a,b,tx,ty,px,py;
 int N = 0;
 int i;
 b = 0.3;
 for ( i = 0; i < 422; i++ )      {
   a = 1.0 + (double) 0.5*i/500.0; // パラメータaの値のセット
   x = y = 0.1;  // 初期値
   N = 0;
   while ( N++ < 800 )     {
     tx = 1.0 - a*x*x + y;
     ty = b*x;
	 x = tx; y = ty;
	 if ( N > 100 )          {
	   px = 100 + i;
	   py = 250.0 - 150.0 * x;
	   Canvas->Pixels[px][py] = clBlack;
	 }
  }
 }
 // ここで座標線を引く
 Canvas->Pen->Color = clRed;
 Canvas->MoveTo(1,250); Canvas->LineTo(550,250);
 Canvas->MoveTo(100,2); Canvas->LineTo(100,600);
 for ( i = 200; i <= 500; i += 100 )    {
	Canvas->MoveTo(i,25); Canvas->LineTo(i,475);
 }
 for ( i = 1; i <= 15; i++ )             {
   Canvas->MoveTo(80,250+i*15); Canvas->LineTo(120,250+i*15);
 }
 for ( i = 1; i <= 15; i++ )             {
   Canvas->MoveTo(80,250-i*15); Canvas->LineTo(120,250-i*15);
 }
}
//---------------------------------------------------------------------------

カオスを勉強しているうちに、大学時代に授業に出たことがある森肇教授が書かれた啓蒙書があることを知った。それは、「カオス 流転する自然」岩波科学ライブラリーである。

大学時代、森教授は統計力学の授業を担当していて、大変偉い方だという評判だった。なんでも非平衡の統計力学が専門だと伺ったのだが、当時不真面目な学生だった私はほとんど授業に出なかった。勿体ないというか申し訳ないので、せめてもの償いにこの本を買って勉強させていただくことにした。

その本を読んでいると、エノン写像というのが出てきた。なんでも、以下のような式らしい。
xi+1 = 1 - 1.4xi2 + 0.3yi
yi+1 = yi

ただ、ここにはあまり詳しいことが載ってないので、Wikiで確認して、微妙に違う以下の式でアトラクタを描いてみた。
xi+1 = 1 - 1.4xi2 + yi
yi+1 = 0.3yi

すると、下記のようになかなか格好いい、ナイキの商標のような軌跡になるのだった。

なお、例によってプログラムのソースコードも示す。

//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
 double x,y,a,b,tx,ty,px,py;
 int N = 0;
 int i,j;
 b = 0.3;
 a = 1.4;
   x = y = 0.1;  // 初期値
   while ( N++ < 10000 )     {
     tx = 1.0 - a*x*x + y;
     ty = b*x;
     x = tx; y = ty;
     if ( N > 100 )         {
       px = 300.0 + 150*x;
       py = 300.0 - 300.0*y;
       Canvas->Pixels[px][py] = clBlack;
    }
  }

}
//---------------------------------------------------------------------------

3体問題で重心座標のことを学んだので、一応2体問題にも重心座標を適用してみることにした。

下記の例をとる。

m1 = 1.1;
m2 = 1.0;

// 位置と速度の初期値
x1 = 3.2;  y1 = 0.0; u1 = 0.01; v1 = 0.2;
x2 = -3.0;  y2 = 0.0; u2 = 0.0; v2 = -0.21;

これの普通の座標での軌道は、以下のとおりであった。

それに重心座標を適用すると以下のようになる。

見て取れるように、まったくなんの変哲もない二つの楕円軌道になるだけなのだが、当たり前のようでもきちんと確認しておきなさいというのが私が幾多の警察小説を読んで学んだことなのであった。

最後に一言。本当はカープ優勝の記事を書こうと思っていたのに、なんで巨人が阪神に勝ってくれんのんじゃ、と広島弁で愚痴をいうのだった。

ステートフル・ブートで述べたLME288というリンク・エラーであるが、ウェブで調べていて意外な解決策があることが分かった。

すなわち、要するに問題を起こしているのは、C:\Program Files(x86)\Embarcadero\RAD Studio\11.0\bin\ilink32.exeというプログラムなので、エクスプローラでとりあえず、そのプログラムがあるフォルダを開く。

そして、ilink32.exeをマウスの右ボタンでクリックして、「プロパティ」のセクションを選び、「互換性」のタブで、Windows 7を選ぶと、なんとLME288というリンク・エラーが解消された。

このような設定は、パソコンを立ち上げる度に実行しなくてはならないのか、それともパソコンを立ち上げる度にilink32.exeをマウスの右ボタンでクリックする必要があるのかどうか分からないが、ともかく問題はilink32.exeのWindows 10との互換性にあったようで、それが分かればまあいいとしたものである。

で、立ち上げなおしても依然としてLME288というリンク・エラーは起きなかったので、互換モードの設定は一応、persistentのようである。とりあえずよかった。

それにしても、今までどうしてこんな解決策がウェブで見つけられなかったのか、不思議でならない。こうしてまた私のウェブ不信症が嵩じるのだった。

これで解決と打ち切ろうとしたら、まだ綾があった。すなわち、ちゃんとコンパイル・リンクできる状態でパソコンを立ち上げなおし、コマンド・ラインでコンパイル・リンクしてみると、問題なく動作できる。ところが本体のGUIのコンパイラの方ではLME288のエラーが出た。これはおかしいと、ilink32.exeの互換性設定メニューにいって、よくみると互換性判定機能というのがあるのでその機能を走らせてみると、なんと正しい互換性はWindows 8なのだという。それはそうだよね、C++ Builder XE4はたしか、2013リリースだものねと思ったのだが、ではどうして今までWindows 7互換モードで動いていたのか、という疑問も出てくる。

これでは突然文句をいってくるクレーマーだよね、と思ってしまった。

それだけではない。システムを立ち上げなおすと、まるでサイコロを振るようにLME288が起きてしまう。それを解決するには、暗中模索のように何度もシャットダウンや再起動しなくてはならず、困惑すること頻りなのだった。

棋譜解析

user-pic
0

棋友が、先日対戦した将棋の棋譜を送ってきた。いままでは私がざっと並べて適当にコメントを述べていたのだが、折角私のパソコンに激指14が導入されているので、それを使って解析することにした。

実は激指14の棋譜解析機能は、王位戦第3局などのような記事でも使っている。それでこの際、その手順のようなものをマニュアル風に述べてみる。

プロの棋譜は基本的にウェブからダウンロードできる。それには大抵、プロ棋士の解説やコメントがあり、それも大いに参考させて頂く。ところがメールで送られてきた棋譜はMS-WORDに打ち込まれたもので、そのままでは激指14が処理できないので、棋譜をいったん印刷してそれをみながら、激指14上で、人間対人間の対戦モードで対戦というかたちで棋譜を辿った。

こうして辿ってみるだけでも、なるほどそう指すものなのか、とか、これはどうみても悪手ではないかと感じたりするので、それはノートにメモっておく。また棋譜を辿る途中で、ときどきありえない、あるいはそこに駒が移動できないだろうという明白な誤記を発見することがあるが、私の将棋の常識に照らして、ここはこの手の間違いだよね、と分かる場合はそのまま進めて、あえて棋友に聞かない。こうして最後まで棋譜を辿ったら一応、棋譜を名前をつけて保存する。

そして、激指14のメニューの「解析」→「棋譜解析」を選ぶと、「解析の棋力を指定して下さい。」とプロンプトされるので、迷わず最強である「Pro+」を選ぶ。すると、初手から最終手までの解析が始まるのだが、「Pro+」だと数十分くらい時間がかかったりするので、本を読んだり、コーヒーを飲んだりして待つ。

解析が終わると、個々の手に疑問手、悪手、敗着などのラベルが付けられることがある。その各々の手が何手目かと、どのラベルか、激指14が代わりにそこからどんな手を読んでいるかかが別ウインドウで表示されるので、それをメモする。

そして、ラベルがついた手と、それに対する激指14の読みを手数昇順で、テキストファイルなどで書き留める。そこに私なりのコメントをつけるとともに、棋友から、ここはどうだったかと聞かれている箇所があったなら、そこにラベルがついていない場合、その局面まで局面を進めて、「解析」→「次の一手」を選ぶ。すると激指14は次の一手とそこからの読み、そしてそこでの形勢を示してくれる。なお、「解析」→「次の一手」の場合は、強さを聞いてこないので、あらかじめメニューの「ファイル」→「設定」で「次の一手」の棋力をPro+に設定しておく。なお、面白いのは同じ局面で「次の一手」をさせるとその都度微妙に手の読みと形勢判断が揺れる。このあたり、確率を入れているからだろう。

ともかくこうやってコンテンツ集めて、なけなしの作文力でまとめて棋友に送ったら、そこそこ喜んでもらえたので、とりあえずよかった。自分の棋力があてにならないのでパソコン将棋に頼っているのが情けない、と思わないのが情けないかも。

なんとか3体問題の面白そうな軌道をみてみたいとウェブで検索していて、中村 健という方のこちらのページをみつけて、あまりにもいろいろなほしい情報が開示されているので感激した。そこで、中村さんに感謝しつつ、この情報を使わせていただくことにした。

まず、ここで教わったことは重心座標を使う、ということである。勘の悪い私は、重心座標といわれても、は?としばらく意味が分からなかったが、ああ、3体の重心を座標の原点とみなして描画するのだな、とやっとわかった。確かに、私が試みたある3体の軌道では、3体が一緒になってx座標が増える方向に移動して、とうとう視界からみえなくなってしまった。こういうのも重心座標だと、きちんと追える。

すなわち、下記のように、まず重心(xc,yc)を求め、x1 - xcのような計算で、重心からの相対的な距離をプロットする。ここでLは表示のスケーリングであり、この値を適当に調整して、見やすいように、且つ画面からはみださないように表示される軌道のサイズを決める。

 xc = (m1*x1 + m2*x2 + m3*x3)/(m1+m2+m3);
 yc = (m1*y1 + m2*y2 + m3*y3)/(m1+m2+m3);
 px = 320.0 + L*(x1 - xc);
 py = 290.0 - L*(y1 - yc);
 Canvas->Pixels[px][py] = clRed;
 px = 320.0 + L*(x2 - xc);
 py = 290.0 - L*(y2 - yc);
 Canvas->Pixels[px][py] = clBlue;
 px = 320.0 + L*(x3 - xc);
 py = 290.0 - L*(y3 - yc);
 Canvas->Pixels[px][py] = clPurple;

そこで、中村さんのホームページから気まぐれで、いくつかピックアップして描画してみた。

まず、面白そうだなと思ったのが、周期解の近くとあった、この軌道。単なる気分で、もとの質量がm1 = m2 = m3 = 3だったのを2.5にしてみたら、変な二層構造になってしまった。この場合は、L = 30とした。

質点の番号mxyuv
1(赤)2.5-300-0.5
2(青)2.53000.5
3(紫)2.505.190.50

次は単純にみえる8の字解。この場合、L = 150とした。単純なのは結果であって、これは今世紀に入って、C. Moore, A. Chenciner & R. Montgomery, C. Sim´o らによって発見されたもので、そのパラメータを決定するために相当な努力があった。

質点の番号mxyuv
1(赤)1/3-1000
2(青)1/31000.5
3(紫)1/300-0.60.923

次は、中村さんがBroucke解と呼ぶもの。私が示した軌道は、早めに打ち切ったものである。理由は、途中から計算誤差のせいだと思うが、突然軌道が発散してしまったからである。それに対処するためにHの刻み値を小さくしたら、軌道を描くのが低速になりすぎて、うんざりした。なお、L = 130とした。

質点の番号mxyuv
1(赤)1/3-1000.3407
2(青)1/31000.3407
3(紫)1/3000.1-0.6814

次は、基本に戻って、新周期解の一種。L = 60とした。

質点の番号mxyuv
1(赤)3-300-0.5
2(青)33000.5
3(紫)305.180.4460

ところで、3体問題の軌道を描いてみての感想は、一見安定しているようでいて、わずかな計算誤差などが、軌道の崩壊につながる、ということである。そのため、私が使っている4次のルンゲクッタ法で、刻みの値Hとして極めて小さいものを使ったりしているのだが、するとどんどん軌道の進み速度が低下する。これはやはり、8次のルンゲクッタ法を使う必要があるのか、と思いつつ、何を調べればそれが載っているのかわからず、困っている。

別のこととして、パラメータの決定がある。たとえば、レヴィの図形ジュリア集合IFSのアトラクタのどれにおいても、適切なパラメータを使用しないと、描画がほとんど意味をなさなくなるので、私はウェブや文献などから学んだ値をほぼそのまま使っているのだが、こういうのをオリジナルに決定するのって、きっとすごく大変なんだよね、と自分に言い聞かせるだけで、どうするのかまったく想像もつかない。なんか情けないが能力のなさはどうしようもない。ともかく、重複をおそれず、ソースコードもリストしておく。

//---------------------------------------------------------------------------

#include <vcl.h>
#include <math.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double fx(double x1, double y1, double x2, double y2);
double fy(double x1, double y1, double x2, double y2);
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{

 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Run1Click(TObject *Sender)
{
const double H = 0.000001;  // 刻み
double x1,y1,u1,v1;
double x2,y2,u2,v2;
double x3,y3,u3,v3;
double a1,a2,a3,a4;  // dx1
double b1,b2,b3,b4;  // dy1
double c1,c2,c3,c4;  // du1
double d1,d2,d3,d4;  // dv1
double e1,e2,e3,e4;  // dx2
double f1,f2,f3,f4;  // dy2
double g1,g2,g3,g4;  // du2
double h1,h2,h3,h4;  // dv2
double p1,p2,p3,p4;  // dx3
double q1,q2,q3,q4;  // dy3
double r1,r2,r3,r4;  // du3
double s1,s2,s3,s4;  // dv3
double m1,m2,m3;       // 質量
double px, py;
double xc, yc;  //重心
double pi = 3.14159;
double TH= 2.0*pi/3.0;
double L = 60.0;
int N = 0;           // 繰り返し回数
m1 = 3.0;
m2 = 3.0;
m3 = 3.0;


// 位置と速度の初期値
x1 = -3.0;  y1 = 0.0; u1 = 0.0; v1 = -0.5;
x2 = 3.0; y2 = 0.0; u2 = 0.0; v2 = 0.5;
x3 = 0.0; y3 = 5.18; u3 = 0.446; v3 = 0.0;

while ( N++ < 1000000000 )        {

 a1 = u1;
 b1 = v1;
 c1 = -m2*fx(x1,y1,x2,y2) - m3*fx(x1,y1,x3,y3);
 d1 = -m2*fy(x1,y1,x2,y2) - m3*fy(x1,y1,x3,y3);
 e1 = u2;
 f1 = v2;
 g1 = -m1*fx(x2,y2,x1,y1) - m3*fx(x2,y2,x3,y3);
 h1 = -m1*fy(x2,y2,x1,y1) - m3*fy(x2,y2,x3,y3);
 p1 = u3;
 q1 = v3;
 r1 = -m1*fx(x3,y3,x1,y1) - m2*fx(x3,y3,x2,y2);
 s1 = -m1*fy(x3,y3,x1,y1) - m2*fy(x3,y3,x2,y2);
 a2 = u1+H*c1/2.0;
 b2 = v1+H*d1/2.0;
 c2 = -m2*fx(x1+H*a1/2.0,y1+H*b1/2.0,x2+H*e1/2.0,y2+H*f1/2.0) 
     - m3*fx(x1+H*a1/2.0,y1+H*b1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 d2 = -m2*fy(x1+H*a1/2.0,y1+H*b1/2.0,x2+H*e1/2.0,y2+H*f1/2.0) 
      - m3*fy(x1+H*a1/2.0,y1+H*b1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 e2 = u2+H*f1/2.0;
 f2 = v2+H*h1/2.0;
 g2 = -m1*fx(x2+H*e1/2.0,y2+H*f1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m3*fx(x2+H*e1/2.0,y2+H*f1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 h2 = -m1*fy(x2+H*e1/2.0,y2+H*f1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m3*fy(x2+H*e1/2.0,y2+H*f1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 p2 = u3+H*r1/2.0;
 q2 = v3+H*s1/2.0;
 r2 = -m1*fx(x3+H*p1/2.0,y3+H*q1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m2*fx(x3+H*p1/2.0,y3+H*q1/2.0,x2+H*e1/2.0,y2+H*f1/2.0);
 s2 = -m1*fy(x3+H*p1/2.0,y3+H*q1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m2*fy(x3+H*p1/2.0,y3+H*q1/2.0,x2+H*e1/2.0,y2+H*f1/2.0);
 a3 = u1+H*c2/2.0;
 b3 = v1+H*d2/2.0;
 c3 = -m2*fx(x1+H*a2/2.0,y1+H*b2/2.0,x2+H*e2/2.0,y2+H*f2/2.0) 
      - m3*fx(x1+H*a2/2.0,y1+H*b2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 d3 = -m2*fy(x1+H*a2/2.0,y1+H*b2/2.0,x2+H*e2/2.0,y2+H*f2/2.0) 
      - m3*fy(x1+H*a2/2.0,y1+H*b2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 e3 = u2+H*f2/2.0;
 f3 = v2+H*h2/2.0;
 g3 = -m1*fx(x2+H*e2/2.0,y2+H*f2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m3*fx(x2+H*e2/2.0,y2+H*f2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 h3 = -m1*fy(x2+H*e2/2.0,y2+H*f2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m3*fy(x2+H*e2/2.0,y2+H*f2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 p3 = u3+H*r2/2.0;
 q3 = v3+H*s2/2.0;
 r3 = -m1*fx(x3+H*p2/2.0,y3+H*q2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m2*fx(x3+H*p2/2.0,y3+H*q2/2.0,x2+H*e2/2.0,y2+H*f2/2.0);
 s3 = -m1*fy(x3+H*p2/2.0,y3+H*q2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m2*fy(x3+H*p2/2.0,y3+H*q2/2.0,x2+H*e2/2.0,y2+H*f2/2.0);
 a4 = u1+H*c3;
 b4 = v1+H*d3;
 c4 = -m2*fx(x1+H*a3,y1+H*b3,x2+H*e3,y2+H*f3) 
      - m3*fx(x1+H*a3,y1+H*b3,x3+H*p3,y3+H*q3);
 d4 = -m2*fy(x1+H*a3,y1+H*b3,x2+H*e3,y2+H*f3) 
      - m3*fy(x1+H*a3,y1+H*b3,x3+H*p3,y3+H*q3);
 e4 = u2+H*f3;
 f4 = v2+H*h3;
 g4 = -m1*fx(x2+H*e3,y2+H*f3,x1+H*a3,y1+H*b3) 
      - m3*fx(x2+H*e3,y2+H*f3,x3+H*p3,y3+H*q3);
 h4 = -m1*fy(x2+H*e3,y2+H*f3,x1+H*a3,y1+H*b3) 
      - m3*fy(x2+H*e3,y2+H*f3,x3+H*p3,y3+H*q3);
 p4 = u3+H*r3;
 q4 = v3+H*s3;
 r4 = -m1*fx(x3+H*p3,y3+H*q3,x1+H*a3,y1+H*b3) 
      - m2*fx(x3+H*p3,y3+H*q3,x2+H*e3,y2+H*f3);
 s4 = -m1*fy(x3+H*p3,y3+H*q3,x1+H*a3,y1+H*b3) 
      - m2*fy(x3+H*p3,y3+H*q3,x2+H*e3,y2+H*f3);


 x1 += H*(a1 + 2.0*a2 + 2.0*a3 + a4)/6.0;
 y1 += H*(b1 + 2.0*b2 + 2.0*b3 + b4)/6.0;
 u1 += H*(c1 + 2.0*c2 + 2.0*c3 + c4)/6.0;
 v1 += H*(d1 + 2.0*d2 + 2.0*d3 + d4)/6.0;
 x2 += H*(e1 + 2.0*e2 + 2.0*e3 + e4)/6.0;
 y2 += H*(f1 + 2.0*f2 + 2.0*f3 + f4)/6.0;
 u2 += H*(g1 + 2.0*g2 + 2.0*g3 + g4)/6.0;
 v2 += H*(h1 + 2.0*h2 + 2.0*h3 + h4)/6.0;
 x3 += H*(p1 + 2.0*p2 + 2.0*p3 + p4)/6.0;
 y3 += H*(q1 + 2.0*q2 + 2.0*q3 + q4)/6.0;
 u3 += H*(r1 + 2.0*r2 + 2.0*r3 + r4)/6.0;
 v3 += H*(s1 + 2.0*s2 + 2.0*s3 + s4)/6.0;
 
 xc = (m1*x1 + m2*x2 + m3*x3)/(m1+m2+m3);
 yc = (m1*y1 + m2*y2 + m3*y3)/(m1+m2+m3);
 px = 320.0 + L*(x1 - xc);
 py = 290.0 - L*(y1 - yc);
 Canvas->Pixels[px][py] = clRed;
 px = 320.0 + L*(x2 - xc);
 py = 290.0 - L*(y2 - yc);
 Canvas->Pixels[px][py] = clBlue;
 px = 320.0 + L*(x3 - xc);
 py = 290.0 - L*(y3 - yc);
 Canvas->Pixels[px][py] = clPurple;
}
}
//---------------------------------------------------------------------------
double fx(double x1, double y1, double x2, double y2) {
 double r,R;
 r = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2);
 R = r * sqrt(r);
 return (x1 - x2)/R;
}

double fy(double x1, double y1, double x2, double y2) {
 double r,R;
 r = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2);
 R = r * sqrt(r);
 return (y1 - y2)/R;
}

3体問題

user-pic
0

さて、2体問題に続き、いよいよ懸案の3体問題の記述に着手することにした。私が若い頃、ルンゲクッタ法で2体問題までは記述できていたのに、3体問題で挫折したのは、3体問題だと、記述が面倒になるからなのだった。それは難しいというよりも、たんに面倒というだけなのだが、当時使用していたコンピュータ言語はBASICで、関数を定義して呼び出すということさえできない。それで細かい記入ミスをしてしまい、私の散漫な注意力では書き間違いをみつけきれなかったようである。

2つの質点の間の引力を記述する関数は、fx(x1,y1,x2,y2)とfy(x1,y1,x2,y2)として定義したので、それをそのまま使えばよい。また、面倒な変数の置き換えは、Perlを使うことにした。

そこで例えば、下記のような内容のテキストファイルをRK.txtとして用意し、これに下記ソースコード1のPerlプログラムをかけることで、うまく変換でき、その結果のRKX.txtを利用することで、ソースコード2のC++のプログラムを得た。これは難なくコンパイルできた。ソースコード2をちらっとでも見ていただくと、練達のプログラマは除き、これを手書きだけで正しく書くことが若干煩わしいことは分かっていただけるだろう。こうして3体問題がなんとか記述できたのだが、適当に位置と速度を与えても、理解困難な不規則な軌道をとって、あっという間に発散してしまう。そこで、とりあえず原点を中心とした円に三角形を形成するように同じ質量の質点を配置し、円に対する接線の方向の初期速度を与えることで、下記のような安定解を得た。

これは、当然に昔から知られていて、ラグランジュの正三角形解と呼ばれているそうである。これは安定解といっても、実は針の先端で踊る天使のように?不安定で、一方の質点の質量を1%増やしただけでも、下記のように発散してしまう。

既知の3体問題の解はいろいろ知られているようであるが、ソースコードのせいで長くなったので、それは別稿とする。

RK.txt

 a1 = u1;
 b1 = v1;
 c1 = -m2*fx(x1,y1,x2,y2) - m3*fx(x1,y1,x3,y3);
 d1 = -m2*fy(x1,y1,x2,y2) - m3*fy(x1,y1,x3,y3);
 e1 = u2;
 f1 = v2;
 g1 = -m1*fx(x2,y2,x1,y1) - m3*fx(x2,y2,x3,y3);
 h1 = -m1*fy(x2,y2,x1,y1) - m3*fy(x2,y2,x3,y3);
 p1 = u3;
 q1 = v3;
 r1 = -m1*fx(x3,y3,x1,y1) - m2*fx(x3,y3,x2,y2);
 s1 = -m1*fy(x3,y3,x1,y1) - m2*fy(x3,y3,x2,y2);

ソースコード1

open FI, "<RK.txt" or die "zzz";
open FO, ">RKX.txt";

$S1{x1} = "x1+H*a1/2.0";
$S1{y1} = "y1+H*b1/2.0";
$S1{u1} = "u1+H*c1/2.0";
$S1{v1} = "v1+H*d1/2.0";
$S1{x2} = "x2+H*e1/2.0";
$S1{y2} = "y2+H*f1/2.0";
$S1{u2} = "u2+H*f1/2.0";
$S1{v2} = "v2+H*h1/2.0";
$S1{x3} = "x3+H*p1/2.0";
$S1{y3} = "y3+H*q1/2.0";
$S1{u3} = "u3+H*r1/2.0";
$S1{v3} = "v3+H*s1/2.0";

$S2{x1} = "x1+H*a2/2.0";
$S2{y1} = "y1+H*b2/2.0";
$S2{u1} = "u1+H*c2/2.0";
$S2{v1} = "v1+H*d2/2.0";
$S2{x2} = "x2+H*e2/2.0";
$S2{y2} = "y2+H*f2/2.0";
$S2{u2} = "u2+H*f2/2.0";
$S2{v2} = "v2+H*h2/2.0";
$S2{x3} = "x3+H*p2/2.0";
$S2{y3} = "y3+H*q2/2.0";
$S2{u3} = "u3+H*r2/2.0";
$S2{v3} = "v3+H*s2/2.0";

$S3{x1} = "x1+H*a3";
$S3{y1} = "y1+H*b3";
$S3{u1} = "u1+H*c3";
$S3{v1} = "v1+H*d3";
$S3{x2} = "x2+H*e3";
$S3{y2} = "y2+H*f3";
$S3{u2} = "u2+H*f3";
$S3{v2} = "v2+H*h3";
$S3{x3} = "x3+H*p3";
$S3{y3} = "y3+H*q3";
$S3{u3} = "u3+H*r3";
$S3{v3} = "v3+H*s3";

while (<FI>)     {

 printf FO "%s",$_;
}

seek FI,0,0;
$N = 2;
while (<FI>)    {
 chomp;
   foreach $x (sort keys %S1)  {
     $y = $S1{$x};
     $_ =~ s/$x/$y/g;
   }
 printf FO "%s%d%s\n",substr($_,0,2),$N,substr($_,3,120);
}

seek FI,0,0;
$N++;
while (<FI>)    {
 chomp;
   foreach $x (sort keys %S2)  {
     $y = $S2{$x};
     $_ =~ s/$x/$y/g;
   }
 printf FO "%s%d%s\n",substr($_,0,2),$N,substr($_,3,120);

}
seek FI,0,0;
$N++;
while (<FI>)    {
 chomp;
   foreach $x (sort keys %S3)  {
     $y = $S3{$x};
     $_ =~ s/$x/$y/g;
   }
 printf FO "%s%d%s\n",substr($_,0,2),$N,substr($_,3,120);

}

ソースコード2

//---------------------------------------------------------------------------

#include <vcl.h>
#include <math.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double fx(double x1, double y1, double x2, double y2);
double fy(double x1, double y1, double x2, double y2);
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{

 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Run1Click(TObject *Sender)
{
const double H = 0.000001;  // 刻み
double x1,y1,u1,v1;
double x2,y2,u2,v2;
double x3,y3,u3,v3;
double a1,a2,a3,a4;  // dx1
double b1,b2,b3,b4;  // dy1
double c1,c2,c3,c4;  // du1
double d1,d2,d3,d4;  // dv1
double e1,e2,e3,e4;  // dx2
double f1,f2,f3,f4;  // dy2
double g1,g2,g3,g4;  // du2
double h1,h2,h3,h4;  // dv2
double p1,p2,p3,p4;  // dx3
double q1,q2,q3,q4;  // dy3
double r1,r2,r3,r4;  // du3
double s1,s2,s3,s4;  // dv3
double m1,m2,m3;       // 質量
double px, py;
double pi = 3.14159;
double TH= 2.0*pi/3.0;
int N = 0;           // 繰り返し回数
m1 = 1.0;
m2 = 1.0;
m3 = 1.0;


// 位置と速度の初期値
x1 = 3.0*cos(TH);  y1 = 3.0*sin(TH); u1 = -0.1*sin(TH); v1 = 0.1*cos(TH);
x2 = 3.0*cos(2.0*TH); y2 = 3.0*sin(2.0*TH); u2 = -0.1*sin(2.0*TH); v2 = 0.1*cos(2.0*TH);
x3 = 3.0; y3 = 0.0; u3 = 0.0; v3 = 0.1;

while ( N++ < 100000000 )        {

 a1 = u1;
 b1 = v1;
 c1 = -m2*fx(x1,y1,x2,y2) - m3*fx(x1,y1,x3,y3);
 d1 = -m2*fy(x1,y1,x2,y2) - m3*fy(x1,y1,x3,y3);
 e1 = u2;
 f1 = v2;
 g1 = -m1*fx(x2,y2,x1,y1) - m3*fx(x2,y2,x3,y3);
 h1 = -m1*fy(x2,y2,x1,y1) - m3*fy(x2,y2,x3,y3);
 p1 = u3;
 q1 = v3;
 r1 = -m1*fx(x3,y3,x1,y1) - m2*fx(x3,y3,x2,y2);
 s1 = -m1*fy(x3,y3,x1,y1) - m2*fy(x3,y3,x2,y2);
 a2 = u1+H*c1/2.0;
 b2 = v1+H*d1/2.0;
 c2 = -m2*fx(x1+H*a1/2.0,y1+H*b1/2.0,x2+H*e1/2.0,y2+H*f1/2.0) 
     - m3*fx(x1+H*a1/2.0,y1+H*b1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 d2 = -m2*fy(x1+H*a1/2.0,y1+H*b1/2.0,x2+H*e1/2.0,y2+H*f1/2.0) 
      - m3*fy(x1+H*a1/2.0,y1+H*b1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 e2 = u2+H*f1/2.0;
 f2 = v2+H*h1/2.0;
 g2 = -m1*fx(x2+H*e1/2.0,y2+H*f1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m3*fx(x2+H*e1/2.0,y2+H*f1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 h2 = -m1*fy(x2+H*e1/2.0,y2+H*f1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m3*fy(x2+H*e1/2.0,y2+H*f1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 p2 = u3+H*r1/2.0;
 q2 = v3+H*s1/2.0;
 r2 = -m1*fx(x3+H*p1/2.0,y3+H*q1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m2*fx(x3+H*p1/2.0,y3+H*q1/2.0,x2+H*e1/2.0,y2+H*f1/2.0);
 s2 = -m1*fy(x3+H*p1/2.0,y3+H*q1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m2*fy(x3+H*p1/2.0,y3+H*q1/2.0,x2+H*e1/2.0,y2+H*f1/2.0);
 a3 = u1+H*c2/2.0;
 b3 = v1+H*d2/2.0;
 c3 = -m2*fx(x1+H*a2/2.0,y1+H*b2/2.0,x2+H*e2/2.0,y2+H*f2/2.0) 
      - m3*fx(x1+H*a2/2.0,y1+H*b2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 d3 = -m2*fy(x1+H*a2/2.0,y1+H*b2/2.0,x2+H*e2/2.0,y2+H*f2/2.0) 
      - m3*fy(x1+H*a2/2.0,y1+H*b2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 e3 = u2+H*f2/2.0;
 f3 = v2+H*h2/2.0;
 g3 = -m1*fx(x2+H*e2/2.0,y2+H*f2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m3*fx(x2+H*e2/2.0,y2+H*f2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 h3 = -m1*fy(x2+H*e2/2.0,y2+H*f2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m3*fy(x2+H*e2/2.0,y2+H*f2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 p3 = u3+H*r2/2.0;
 q3 = v3+H*s2/2.0;
 r3 = -m1*fx(x3+H*p2/2.0,y3+H*q2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m2*fx(x3+H*p2/2.0,y3+H*q2/2.0,x2+H*e2/2.0,y2+H*f2/2.0);
 s3 = -m1*fy(x3+H*p2/2.0,y3+H*q2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m2*fy(x3+H*p2/2.0,y3+H*q2/2.0,x2+H*e2/2.0,y2+H*f2/2.0);
 a4 = u1+H*c3;
 b4 = v1+H*d3;
 c4 = -m2*fx(x1+H*a3,y1+H*b3,x2+H*e3,y2+H*f3) 
      - m3*fx(x1+H*a3,y1+H*b3,x3+H*p3,y3+H*q3);
 d4 = -m2*fy(x1+H*a3,y1+H*b3,x2+H*e3,y2+H*f3) 
      - m3*fy(x1+H*a3,y1+H*b3,x3+H*p3,y3+H*q3);
 e4 = u2+H*f3;
 f4 = v2+H*h3;
 g4 = -m1*fx(x2+H*e3,y2+H*f3,x1+H*a3,y1+H*b3) 
      - m3*fx(x2+H*e3,y2+H*f3,x3+H*p3,y3+H*q3);
 h4 = -m1*fy(x2+H*e3,y2+H*f3,x1+H*a3,y1+H*b3) 
      - m3*fy(x2+H*e3,y2+H*f3,x3+H*p3,y3+H*q3);
 p4 = u3+H*r3;
 q4 = v3+H*s3;
 r4 = -m1*fx(x3+H*p3,y3+H*q3,x1+H*a3,y1+H*b3) 
      - m2*fx(x3+H*p3,y3+H*q3,x2+H*e3,y2+H*f3);
 s4 = -m1*fy(x3+H*p3,y3+H*q3,x1+H*a3,y1+H*b3) 
      - m2*fy(x3+H*p3,y3+H*q3,x2+H*e3,y2+H*f3);


 x1 += H*(a1 + 2.0*a2 + 2.0*a3 + a4)/6.0;
 y1 += H*(b1 + 2.0*b2 + 2.0*b3 + b4)/6.0;
 u1 += H*(c1 + 2.0*c2 + 2.0*c3 + c4)/6.0;
 v1 += H*(d1 + 2.0*d2 + 2.0*d3 + d4)/6.0;
 x2 += H*(e1 + 2.0*e2 + 2.0*e3 + e4)/6.0;
 y2 += H*(f1 + 2.0*f2 + 2.0*f3 + f4)/6.0;
 u2 += H*(g1 + 2.0*g2 + 2.0*g3 + g4)/6.0;
 v2 += H*(h1 + 2.0*h2 + 2.0*h3 + h4)/6.0;
 x3 += H*(p1 + 2.0*p2 + 2.0*p3 + p4)/6.0;
 y3 += H*(q1 + 2.0*q2 + 2.0*q3 + q4)/6.0;
 u3 += H*(r1 + 2.0*r2 + 2.0*r3 + r4)/6.0;
 v3 += H*(s1 + 2.0*s2 + 2.0*s3 + s4)/6.0;
 px = 290.0 + 70.0*x1;
 py = 290.0 - 70.0*y1;
 Canvas->Pixels[px][py] = clRed;
 px = 290.0 + 70.0*x2;
 py = 290.0 - 70.0*y2;
 Canvas->Pixels[px][py] = clBlue;
 px = 290.0 + 70.0*x3;
 py = 290.0 - 70.0*y3;
 Canvas->Pixels[px][py] = clPurple;
}
}
//---------------------------------------------------------------------------
double fx(double x1, double y1, double x2, double y2) {
 double r,R;
 r = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2);
 R = r * sqrt(r);
 return (x1 - x2)/R;
}

double fy(double x1, double y1, double x2, double y2) {
 double r,R;
 r = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2);
 R = r * sqrt(r);
 return (y1 - y2)/R;
}

このアーカイブについて

このページには、過去に書かれたブログ記事のうちコンピュータカテゴリに属しているものが含まれています。

次のカテゴリはスピリチュアルです。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。