〔プログラムの説明〕

図に示すように，2 次元平面の中に n 個（ 3 ≦ n ≦ 20 ）の点が与えられたときに， これらの点すべてを含む円のうち，半径が最小である円の中心座標と半径を求めるプログラムである。

(1)　各点の位置は，x 座標と y 座標で与えられる。 各座標値はともに 10⁴ を超えない正の値である。

(2)　ここでは，次のアルゴリズムによって，求める円の中心を探し出している。

　�@　仮中心点を原点（0.0，0.0）に，移動係数を 0.5 に設定する。

　�A　n 個の点の中で，仮中心点から最も遠い位置にある点（最遠点）の方向に向かって， 一定の距離（以下，移動距離と呼ぶ）だけ仮中心点を移動する。 このときの移動距離は，移動前の仮中心点から最遠点までの距離に移動係数を乗じた値とする。 これを一定回数（ここでは 100 回）繰り返した後に �B へ移る。
　なお，繰返しの途中で，仮中心点の移動距離が一定の値（10⁻⁹ ） 以下になった場合には，移動前の仮中心点を円の中心座標とみなし，処理を終了する。

　�B　移動係数を 倍にして，�A へ戻る。

(3)　関数 solvcircle の仕様は，次のとおりである。

　　void solvcircle( double pt_x[], double pt_y[],
　　　　　　　　　　double *cx, double *cy, double *r );

　・機能：引数で受け取った点の座標値を基に，これらの点を内包する最小円の 中心座標（x 座標と y 座標）と半径を求める。

　・引数：pt_x　対象となる各点の x 座標を格納した配列変数

pt_y　対象となる各点の y 座標を格納した配列変数

cx　求まった円の中心点の x 座標を格納する変数へのポインタ

cy　求まった円の中心点の y 座標を格納する変数へのポインタ

r　求まった円の半径を格納する変数へのポインタ

(4)　関数 distance の仕様は，次のとおりである。

　　double distance( double x1, double y1,

　　　　　　　　　　double x2, double y2 );

　・機能：引数で指定された２点間の距離を求める。

　・引数：x1，y1　一方の点の x 座標（x1）と y 座標（y1）

　　x2，y2　他方の点の x 座標（x2）と y 座標（y2）

　・返却値：2 点（x1，y1）と（x2，y2）の距離

(5)　関数 main は，4 個の点座標が与えられたときの関数 solvcircle の使用例を示している。

(6)　プログラムでは，解を 10⁻⁸ までの精度で求める。

#include <stdio.h>
#include <math.h>
#define CONVRGV   1.0e-9
#define EOD       -1.0
#define LOOPMAX   100
#define ON        1
#define OFF       0
void solvcircle( double[], double[], double*, double*, double* );
double distance( double, double, double, double );
main()
{
    double pos_x[] = { 10.0,  50.0,  50.0, 100.0, EOD };
    double pos_y[] = { 50.0,  10.0, 100.0,  50.0, EOD };
    double x, y, r;
    solvcircle( pos_x, pos_y, &x, &y, &r );
    printf( "    円の中心座標 ( %12.8f , %12.8f )  半径 = %12.8f\n",
            x, y, r );
}
void solvcircle( double pt_x[], double pt_y[],
                    double *cx, double *cy, double *r )
{
    double mvrate=0.5, maxdist, dist;
    double xsv=0.0, ysv=0.0;
    int k, i, t, cvgflg;
    *cx = *cy = 0.0;  /* 仮中心点の初期座標 */
    cvgflg = OFF;
    while ( cvgflg == OFF ) {
        for ( t=0; t < LOOPMAX; t++ ) {
            maxdist = 0.0;
            for ( i=0; pt_x[i] != EOD; i++ ) {
                dist = distance( *cx, *cy, pt_x[i], pt_y[i] );
                if ( dist > maxdist ) {
                    maxdist = dist;
                    k = i;
                }
            }
            *cx += ( pt_x[k] -  ) * mvrate;
            *cy += ( pt_y[k] -  ) * mvrate;
            if (  <= CONVRGV ) {
                *cx = xsv;
                *cy = ysv;
                cvgflg = ON;
                break;
            }
            xsv = *cx;
            ysv = *cy;
        }
        mvrate /= 2.0;
    } 
    *r = ;
}
double distance( double x1, double y1, double x2, double y2 )
{
    return sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );
}

ア　*cx

イ　*cy

ウ　dist

エ　distance( pt_x[0], pt_y[0], *cx, *cy )

オ　distance( pt_x[0], pt_y[0], xsv, ysv )

カ　maxdist

キ　pt_x[0]

ク　pt_y[0]

ア　distance ( pt_x[k], pt_y[k], *cx, *cy )

イ　distance ( pt_x[k], pt_y[k], xsv, ysv )

ウ　distance ( xsv, ysv, *cx, *cy )

エ　maxdist

オ　maxdist - dist