#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#include "nrutil.h"

double func(double x[])
{
	double a0, a1, a2;

	a0 = (x[0] - 1.0); // (x-1)
	a1 = (x[1] - 2.0); // (y-2)
	a2 = (x[2] - 3.0); // (z-3) のこと．
	return fabs(a0 * a0 * a0) + a1 * a1 + fabs(a2);
}

int main(void)
{
	double **p, *y, ftol; // 二次元配列用ポインタ（入力用変数），
	int i, j, ndim, mpts, nfunk;

	ndim = 3; // 変数の次数：ここでは，x,y,z
	mpts = ndim + 1;
	ftol = 1.0e-12; // 許容誤差

	p = matrix(mpts, ndim); // 2次元のメモリ確保,,これが，p[i]を可能にしているのだろう！．
	y = vector(mpts);		// 1次元のメモリ確保
	for(i = 0; i < mpts; i++){
		//printf("first\n"); // 4回動いてる
		for(j = 0; j < ndim; j++){ // 3回動いてる
			//printf("second\n");
			p[i][j] = (double)rand() / (double)RAND_MAX; 
			// ポインタpに値を代入
			// p[i][j]の値を最大1に正規化している．
			//printf("%lf\n",p[i][j]);
		}
	}

	/* *********************************************** */
	/*   [input]に関する処理を表示するために存在する, 極論すれば無くても良い*/
	/* *********************************************** */
	printf("f(x,y,z) = |(x-1)^3| + (y-2)^2 + |z-3|\n");
	printf("極小点：x=1, y=2, z=3\n\n");
	printf("[ input ]\n");
	printf("\t\t x        y        z\n");
	for(i = 0; i < mpts; i++)
	{
		printf("function: ");
		for(j = 0; j < ndim; j++) printf("%lf ",p[i][j]);
		//printf("%lf\n",p[i]);
		printf(" = %lf\n", func(p[i]));// 関数の極値を表示,p[i]=0
	}
	printf("tolerance = %e\n\n", ftol);
	/* *********************************************** */

	
	amoeba(p, y, ndim, ftol, func, &nfunk);

	printf("[ output ]\n");
	printf("count=%d\n", nfunk);
	printf("\t\t x        y        z\n");
	for(i = 0; i < mpts; i++)
	{
		printf("function: ");
		for(j = 0; j < ndim; j++) printf("%lf ",p[i][j]);  // p[i][j]は収束後の変数値
		printf(" = %lf\n", y[i]); // y[i]は収束後の関数極小値
	}

	free_matrix(p);
	free_vector(y);

	return 0;
}