// Очень ВАЖНЫЙ параметр остановки счёта !!!
#define EPS 1.e-6

// общий итерационный процесс поиска решения
// модификация Исаева-Сонина
int resh(int n1,int n2,ext_value1 p[])
{
	int i,j,l=0;
	ext_value1 q[n2],pn[n1]; // массивы невязок и пробных параметров пристрелки
	double df0,dfn,w,dp;
	double A[n2*n2]; // матрица производных
	double b[n2];   // столбец правой части
	double gamma=1; // коэффициент Исаева-Сонина

// цикл требуется для правильной инициализации размерности и сохранении значения alpha:
for(i=0;i<n1;i++) pn[i]=p[i];

// printf("resh: f (iter=0) start\n"); 
	f(n1,n2,p,q); // вычисление функции невязок и производных
// printf("resh: f (iter=0) end\n");
	
do{
	printf("---> l=%d : gamma=%lf\n",l,gamma);
	gamma=1.0;

    printf("Базовая точка (p)\n");
	for(i=0;i<n2;i++) {
//		printf("(%d) ",i );
		printf("%le  ",p[i].get() );
//		for(j=0;j<n;j++) printf("%lf ",p[i].get(j) );
//		printf("\n");
	}
	printf("\n");
    // вычисляем декартову норму невязок
    // и сохраняем градиент в свободном столбце СЛАУ
    df0=0.0;
    for(i=0;i<n2;i++) {
		b[i]=q[i].get();
		df0+=b[i]*b[i];
	}
	df0=sqrt(df0);
	
//	printf("декартова норма градиента (b)\n");
//	printf("df0=%le\n",df0);
// Внимание!!! 
// При слишком маленьком EPS хороший ответ считается плохим
// Нужно обдумать ещё одну остановку счёта !!!
// Внимание : условие остановки счёта дублируется ниже
	if( df0<EPS ) return l; // это условие выхода - возвращаем число итераций 
	
	// матрица Якоби (вторых производных) сохраняется как матрица СЛАУ:
	for(i=0;i<n2;i++) for(j=0;j<n2;j++) A[i*n2+j]=q[i].get(j);
	// решение СЛАУ (матрица испорчена, в b оказывается ответ)
	gauss(A,b,n2);
	// вычисляем пробное смещение:
	l++; // увеличили счётчик числа итераций
	
	printf("Смещение (dp):\n");
	for(i=0;i<n2;i++) printf("%lf  ",b[i]);
	printf("\n");
	    while(1)
	    {
			for(i=0;i<n2;i++) pn[i]=p[i]-gamma*b[i];
			f(n1,n2,pn,q);
			dfn=0.0;
			for(i=0;i<n2;i++) {
				w=q[i].get();
				dfn+=w*w;
			}
			dfn=sqrt(dfn);
			dp=0.0; // вычисляем "абсолютно-относительную норму" смещения
			for(i=0;i<n2;i++) {
				w=fabs( p[i].get() );
				if( w>1 ) w=fabs(b[i])/w;
				else w=fabs(b[i]);
				dp+=w*w;
			}
			dp=sqrt(dp);
// Внимание !!! 
// Здесь написано условие выхода по малости градиента или малости смещения			
			if( dfn<EPS || (gamma>0.8 && dp<EPS) ) {
				for(i=0;i<n2;i++) p[i]=pn[i];
				return l; // это условие выхода - возвращаем число итераций
				// ответ сохранён в массиве p[n]
			}
			if( dfn<df0 ) {
				for(i=0;i<n2;i++) p[i]=pn[i];
				break; // новая точка лучше норме вектор-функции невязок,
				// но недостаточно хороша - отправились на новую итерацию
				// можно оптимизировать, так как условие выхода проверяется дважды
				// пока оптимизацию отложил
			}
			if( gamma > EPS ) gamma/=2; // дробление шага Исаева-Сонина со спорным условием малости
			else return -29; // "плохая" остановка счёта (плохая функция)
		}
// конец цикла выбора нового смещения;
// начинается новая итерация
} while(l<15);	
// это плохая остановка счёта: сделали 15 итераций и ответ найти не удалось
	return -l;
}