//擬牛頓(QuasiNewton)法求解非線性方程組一組實根
#include "Comm.h"		//公共部分頭文件
#include "Matrix.h"		//矩陣類及相關函數等的定義
#include "LinearEquation.h"	//線性方程(組)求解頭文件
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#include <iostream>		//輸入輸出流頭文件

const int  N=2;         //方程組變量的個數
const int  M=20;        //迭代的最大次數
const double eps=1e-5;  //eps為精度

template <class _Ty> 
void FunctionValueRSN(valarray<_Ty>& x, valarray<_Ty>& y)
{//輸入方程組
/*
	x^2+y^2-5=0
	(x+1)*y-(3*x+1)=0
*/
	y[0]=x[0]*x[0]+x[1]*x[1]-5;
	y[1]=(x[0]+1)*x[1]-(3*x[0]+1);
}


template <class _Ty>
inline int 
RootQuasiNewton(_Ty eps, _Ty t, _Ty h, valarray<_Ty>& x, int k)
{
	 int i, j, l;
    _Ty am, z, beta, d;
	
	int n = x.size();		//方程組中各方程個數，也是未知數個數

	matrix<_Ty> a(n, n);	//定義二維數組對象a
	valarray<_Ty> b(n);		//定義一維數組對象b
	valarray<_Ty> y(n);		//定義一維數組對象y
	    
	l = k;
	am = 1.0 + eps;
    while(am>=eps)
    {
		FunctionValueRSN(x, b);
        am = 0.0;
        for(i=0; i<n; i++)
        {
			z = Abs(b[i]);
            if(z>am) am = z;
        }
        if(am>=eps)
        {
			l--;
            if(l==0)
			{
				cout << "Fail(0)!" << endl;
				return(0);
			}

            for(j=0; j<n; j++)
            {
				z = x[j];
				x[j] += h;
                FunctionValueRSN(x, y);
                for(i=0; i<n; i++) 
					a(i,j) = y[i];
                x[j] = z;
            }

			if(LE_TotalChoiceGauss(a, b) == 0)
			{
				cout << "Fail(-1)!" << endl;
				return(-1);					//矩陣a奇異
			}

            beta = 1.0;

            for(i=0; i<n; i++) beta -= b[i];
            
			if(FloatEqual(beta,0))
			{
				cout << "Fail(-2)!" << endl;
				return(-2);					//beta=0
			}

            d = h / beta;
            
			for(i=0; i<n; i++) x[i]=x[i]-d*b[i];

            h = t * h;
        }
    }
    return(k-l);					//正常結束，返回迭代次數
}

void main()
{
	double a[N] = {1, 1};              //在(1，1)附近求方程的解,維數“2”是方程組變量的個數
	valarray<double> va(a, N);
	double t=0.01,h=0.001;  
	RootQuasiNewton(eps, t, h, va,  M);
	cout<< "方程組的解為 :  ";
	ValarrayPrint(va);
}