/*
 * 各種反復法による平方根の計算
 * 基本プログラム
 * コンパイル：cc sqrt.c -o sqrt -lm
 */

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

// 表示オプション
#define TRACE	// 途中経過を表示する

// 精度（打ち切り誤差の最大値）
#define EPS	1.0e-5

// 反復回数の限度
#define MAX	1000

// 反復回数
int count = 0;

/*
 * 問題：x = sqrt(A) を計算せよ．ただし，A > 1 とする．
 * → x^2 = A
 * → f(x) = x^2 - A = 0 の解 x を求めよ．
 * （f(x) = A - x^2 = 0 としても同じハズ．）
 */
#define A	2.0

// 関数 f(x)
double f(double x)
{
	return (x*x - A);	// x>0に対して単調増加となる定義方法
//	return (A - x*x);	//    　〃　　単調減少　　　〃
	// 注意：二分法では単調増加／減少によって処理方法が変化
}

// 導関数 f'(x) = df/dx （ニュートン法で使用）
double df(double x) {
	return (2.0*x);		// x>0に対して単調増加となる定義方法
//	return (-2.0*x);	//    　〃　　単調減少　　　〃
}

// 写像の初期値
#define X0	A

// 二分法の初期値
#define XL	1.0	// 下限値
#define XU	A	// 上限値

// 写像 x = g(x)
#define PROBLEM	"実験(1) ニュートン法"
// ニュートン法の漸化式：x(i+1)=x(i)-f(x(i))/f'(x(i))
double g(double x)
{
	return (x - f(x)/df(x));
}
/*
*/

/*
#define PROBLEM	"実験(2) 反復写像の不適切な使用例"
// 問題の方程式を変形：x=√A → x^2=A → x=A/x
double g(double x)
{
	return (A/x);
}
// これだと収束しません!! Lipschitz条件が成立してますか？
*/

/*
#define PROBLEM	"実験(3) 反復写像の適切な使用例"
// 実験(2)の式を更に変形：x=A/x → x+x=x+A/x → x=(x+A/x)/2
double g(double x)
{
	return ((x + A/x)/2.0);
}
// ちなみにこれは，ニュートン法と同じ計算となっている．
// 他にも複数の変形版が考えられる．
*/


// エラー出力
void fatal(char *msg)
{
	fprintf(stderr, "!!! 解析失敗 !!! %s\n", msg);
	exit(1);
}

// 反復法（反復写像，ニュートン法）
// x0: 初期値
// 注意：反復過程でとりうるxの範囲内で，
//	Lipschitz条件 |g'(x)|<1 が成立すれば収束する．
//	i.e. 不成立だと収束の保証なし．
double mapping(double x0)
{
	double	x;	// x_(i+1)

	count = 0;
#ifdef TRACE
	printf("%3d : %f\n", count, x0);
#endif
	while (1) {
		x = g(x0);
		count++;
#ifdef TRACE
		printf("%3d : %f\n", count, x);
#endif
		if (fabs(x - x0) < EPS) break;
		if (count >= MAX) fatal("mapping(): not converged.");

		x0 = x;
	}
	return (x);
}

// 課題：二分法
// x0, x1: 初期値（下限，上限）
double bisect(double x0, double x1)
{
	double	x;	// 中点

	count = 0;
#ifdef TRACE
	printf("%3d : %f %f\n", count, x0, x1);
#endif
	while (1) {
		x = (x0 + x1)/2.0;
		count++;
#ifdef TRACE
		printf("%3d : %f\n", count, x);
#endif
		if (fabs(x - x0) < EPS) break;
		if (count >= MAX) fatal("bisect(): not converged.");

		if (f(x) < 0.0) {	// f(x) が負ということは...
			x0 = x;	// 中点を下限にする（f(x)が単調増加の場合）
		} else {
			x1 = x;	// 中点を上限にする（　　　　〃　　　　　）
		}
	}
	return (x);
}

int main()
{
	double	a;
	printf("平方根の計算：\n");
	printf("精度（打ち切り誤差の最大値）：%e\n\n", EPS);

	printf("数学ライブラリ関数：\n");
	printf("  sqrt(%f) = %f\n\n", A, sqrt(A));

	printf("%s：\n", PROBLEM);
	printf("  sqrt(%f) = %f\n", A, mapping(X0));
	printf("  count = %d\n\n", count);

	printf("二分法：\n");
	printf("  sqrt(%f) = %f\n", A, bisect(XL,XU));
	printf("  count = %d\n\n", count);
	return (0);
}